Author |
Message |
Jared Barber (jared_barber)
Member Username: jared_barber
Post Number: 30 Registered: 01-2007
| Posted on Wednesday, April 16, 2008 - 07:02 pm: | |
Hello, I have moving particles hitting walls. In real incompressible flow, lubrication theory tells us (for these particular particles and walls) that the particles should never hit the walls. To prove that this should be the case, I have taken the "penalty" used in "viscous.pde" (grad^2(p) = penalty*div(u)) and ramped it up to 10000 (we set the penalty to be a constant in our equations since the length scale is O(1) always). In this case, the particles don't hit the wall suggesting the reason the particles are hitting the wall in penalty = 100 case is because the incompressibility condition in these gaps (where the particles hit a wall) is not being enforced as well as it could be. The ultimate goal, then, is to enforce more incompressibility in regions where particles are in danger of hitting walls. One way to do this is to assign penalty = 10000 everywhere for all of our runs. The problem is that this increases computational time a lot. Our simulations go from 1 min per run to 10 min per run. That is what happens if you take our code and replace penalty = 100 with penalty = 10000. Another way to do find out what happens when penalty = 10000, however, is to take and make this into a staged problem by setting penalty = staged(100,10000). One would hope this would speed up the computation since a really good guess for the 10000 case would be obtained from the 100 case. Unfortunately, this does not seem to work as well as it could, 20 min down to 8 min (e.g.). So, the other idea I came up with is to take the penalty and make it a function of time. That is penalty = (e.g.) 100*10^(2*t). Then you solve for time from t = 0 to 1. Surpisingly enough this works well for our problem (which involves 42 scalar equations, a lot) taking approx 15-20 min of computation down to 4 min. I have attached a similar smaller problem (with only 1 scalar equation and a much simpler geometry...because of this just setting penalty = 10000 works faster than trying the attached time approach...this is not the case for our other problem). By looking at this problem, you can see how it works. One problem is I don't know the accuracy of the solution returned. For instance, if I insert the correct initial values for u, v, and p (obtained from running Lube_Flow with k = .01) the values from each method (setting penalty = 10000 or using the time approach) actually agree relatively well whereas if I don't insert the correct initial values, the values can differ a bit (compare u0 for instance). I'm not sure why this is the case since there is no explicit time dependence in this problem (give a penalty, u, v, and p depend only on the value of that penalty and not on any values of any parameters from previous or subsequent time steps). So the question is, how accurate will this time-dependent penalty approach be in the end? As a side note, is there anyway to obtain the distance from (x,y) in the domain to the nearest boundary? Thanks, Jared
|
Robert G. Nelson (rgnelson)
Moderator Username: rgnelson
Post Number: 1101 Registered: 06-2003
| Posted on Friday, April 18, 2008 - 02:53 pm: | |
Rather than address the time dependent problem, I went back to the beginning and tried to speed up the staged version. If you look at the U plots, you will see there are some twists at the exit planes. This indicates a conflict between the boundary condition and the needs of the PDE. Such a conflict can cause the solution to struggle, because there is really no solution that satisfies both. I replaced your outflow condition with the "tautological" ones, Natural(U)=sigma12 Natural(V)=sigma22 at the top Natural(U)=-sigma12 Natural(V)=-sigma22 at the bottom. This improves the smoothness of the solution at the exit, but there is still a whip in the contours, indicating that P=0 is not really a consistent BC either. I made a four stage run, starting with high viscosity and low penalty, then lowering the viscosity and raising K, so that K*penalty goes from 100 to 100000. I added histories to track U0 and dx(U) as the penalty increases. These plots show the consistency of the result at the higher penalties. This script completes in 6 seconds on my 1.5GHz Athlon.
|
Jared Barber (jared_barber)
Member Username: jared_barber
Post Number: 31 Registered: 01-2007
| Posted on Sunday, May 11, 2008 - 12:00 pm: | |
Hello, I have an even more challenging problem in terms of incompressibility. I found out that using a space-dependent penalty actually seems to the job ok, for arcs approaching lines. By the way, once again, is there any quick way to find out the d(x,y) where d(x,y) is the minimum distance from the point (x,y) to all segments of the boundary? Now, the new problem is that I want to keep corners from hitting lines. I have attached one time step from my full-fledged code in a zipped file so that what I am talking about can be seen. In the simulation, a "corner" of one object is very close to a "wall" of a different object. In my code it is realistic to assume that the "corner" will oppose approach to the "wall" (experiencing forces on the order of 1/(distance to wall)^3 pushing it away from the wall). As such there are a select few equations with very large coefficients (that are trying to enforce slow approach) while all of the rest have O(1) coefficients. I want to solve this system but it seems FlexPDE is having difficulties doing this for me. Way #1: Try to solve the system straightforward. (To do this, just open "Norm_Flow.pde" and run it.) Problem: this takes in excess of 20 minutes (I have yet to have had the patience to let it finish.) In addition, if I coax it by using some nice transfer data (including initial mesh and good initial guess for u, v, and p), I once got it to solve in 16 minutes. Even then it seems that while all the scalar equations should equal 0, they actually equal O(1000). In addition, the corner is still approaching the wall much faster than it should be. Way #2: Turn on the corner forces using a staged variable ("lub_on" in my code, just change lub_on in c_vs_defs.txt file to be staged(0,1) or staged(0,0.5,1) or anything like that, same results). Problem, the nodal equations corresponding to the corner node and the colliding wall nodes are not 0 but rather O(1000000). So the question is, how do I force the equations governing the movement of the corner and the wall near the corner to be 0? It seems if I use staging, those equations are allowed to be O(1000000). If I don't use staging, it seems it tries not to let it be O(1000000) but it takes a long long time for it to get an answer which is still not that good in the end. Is there some sort of request I can make to enforce those special equations more than the rest? (And why does staging allow the residual on such scalar equations to be so big when unstaged runs don't?) Please let me know. Also, I am fairly certain that I have no small errors lurking beneath the mounds of include files I have included but I can always try to make a simpler example to look at if anyone would think that would be more helpful. Hoping for a couple suggestions, Jared |
Robert G. Nelson (rgnelson)
Moderator Username: rgnelson
Post Number: 1119 Registered: 06-2003
| Posted on Sunday, May 11, 2008 - 04:52 pm: | |
As I believe I mentioned before, I think the process of scaling equations is not working well when there are many GLOBAL equations. Frankly, we never expected anyone to have so Many.... We have discovered that the problem behaves much better if you turn matrix scaling off and use a much smaller overshoot: SELECT SCALING=OFF OVERSHOOT=1e-6 The conjugate-gradient linear system solver iterates until the solution updates are resolved to ERRLIM*OVERSHOOT. If some components of the solution are much larger or smaller than the others, requiring a tighter convergence helps resolve the smaller components. With these changes, the problem you sent completes in 18 minutes on our 3GHz Pentium, with global Ux_x and Vx_x values like 1e-9. |
Jared Barber (jared_barber)
Member Username: jared_barber
Post Number: 32 Registered: 01-2007
| Posted on Sunday, May 11, 2008 - 09:52 pm: | |
Hello, I am very sorry but I completely forgot to turn the "problem" on. I'm not sure what possessed me. Anyway, I have attached a new version with the particular "lub_on" parameter set to 1 rather than 0. I think if you run this you may discover different results. (Hopefully along the lines of what I spoke of before.) Jared
|
Marek Nelson (mgnelson)
Moderator Username: mgnelson
Post Number: 34 Registered: 07-2007
| Posted on Monday, May 12, 2008 - 09:00 pm: | |
I can't explain why, but it appears that by introducing the 'lub_on' components the matrix is ill-conditioned enough that FlexPDE does not converge. By gradually introducing the 'lub_on' factor in stages, as one might do in a non-linear problem, FlexPDE is able to converge. I staged the problem with lub_on = STAGED(0, 0.01, 0.1, 0.25, 0.5, 1) and the problem ran in about 30 min. (I also used the previously mentioned SELECT SCALING=off) I didn't do this, but you can use HISTORY plots to display the changes to your global variables as the stages are applied. This would help you to pick more meaningful staged values of 'lub_on'. I just picked a few arbitrarily. Also, now that the problem converges, you may want to decrease your error limit. You are using errlim=0.01 which is fairly loose. |
|