| Author | Message | |||
     Dmitriy Reznikov (barabus) New member Username: barabus Post Number: 1 Registered: 04-2007 |
In all fluid flow examples/code samples, as well as in Backstrom book, the continuity equation is replaced with the following derived equation : div(grad(p)) = 2*ro*[dx(u)*dy(v) - dy(u)*dx(v)] - C*mu/L0^2*(dx(u)+dy(v)) where L0 is a characteristic dimension for a particular problem, and C is a "constant, a suffitiently large number chosen to enforce continuity". This equation gives me a lot of headache, because the solution progress, as well as wether the problem will converge at all, is very much dependent on the value of this constant C. Backstrom writes "FlexPDE wouldn't accept PDEs of o(<2)" and this is the reason that this equation is derived. Is there any other way, a better way, to enforce continuity. What if to enable [ select FIRSTPARTS on ], can I have continuity equation in its native form : dx(u)+dy(v)=0 ( for incompressible 2D ) ? Will this work ? How do I know that this would be a right way of doing this ? | |||
| Dmitriy Reznikov (barabus) New member Username: barabus Post Number: 2 Registered: 04-2007 |
| |||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 814 Registered: 06-2003 |
The embarassment of the Navier-Stokes equations is that you have in principle three equations in two unknowns (u,v). The system is overdetermined in (u,v) and underdetermined in pressure. One way out of this dilemma is to replace the continuity equation with an equation relating pressure to velocity, resulting in a system of three equations in three unknowns (u,v,p). My preference in this matter is to use the "Penalty Pressure" approach described in the comments to the example problem "Samples | Steady_State | Fluids | Viscous.pde". In essence, this defines a fictitious equation of state that is "slightly compressible" rather than strictly incompressible. This approach corresponds to the second term of Professor Backstrom's equation, and yes, there is an arbitrariness in the assignment of the compressibility C. However, the strictly incompressible equation is a mathematical abstraction and is a very harsh equation to solve numerically. The artificial compressibility coefficient allows some control of the tradeoff between computational efficiency and density variations. The first term in Professor Backstrom's equation is derived from a combination of the two momentum equations, and I believe it to be a tautology that contributes nothing to the solution. Professor Backstrom does not agree with me, and I am not an expert in NS equations. So I leave it up to you. In selecting a value for C (in truth a "penalty" coefficent), there are two extremes to be avoided: 1) if C is too small, mass will not be conserved; 2) if C is too large, the pressure equation will dominate the numerics and the computation will be expensive and slow. Incidentally, the factor C*mu/L0^2 us really the "penalty" coefficient. Professor Backstrom uses the mu/L0^2 factor because it is dimensionally correct and imparts some degree of adaptability to the term. It has no other fundamental justification and is in that sense arbitrary. A more serious issue is the boundary condition on the pressure equation. The form used by Prof Backstrom is justifiable from arguments of momentum conservation when a stream hits a wall. But it is numerically troublesome, and contributes to instabilities in some cases. I always use NORMAL(P)=0, and will continue to do so until someone proves to me that it is significantly in error. I suggest that you try this BC and see if some of your troubles don't go away. FlexPDE accepts PDEs of o(<2). The problem is that the equation dx(u)+dy(v)=0, when used as a defining equation for P, blows up because it has no dependence on P. It's like saying "given a+b=0 find c". | |||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 815 Registered: 06-2003 |
PS A typo: The boundary condition in paragraph 7 should read NATURAL(P)=0, not NORMAL(P)=0. | |||
| Dmitriy Reznikov (barabus) Junior Member Username: barabus Post Number: 3 Registered: 04-2007 |
Thank you very much, Dr.Nelson, I am trying to implement your comments. Specifically, taking out the (mu/L0^2) multiplier and disabling the (2*ro*[dx(u)*dy(v) - dy(u)*dx(v)]) in the "continuity" equation. I will post later on the results of these. | |||
| Dmitriy Reznikov (barabus) Member Username: barabus Post Number: 4 Registered: 04-2007 |
Ok, I tried, but I actually find the natural(p)=mu*(normal(vector(1,0))*del2(u)+normal(vector(0,1))*del2(v)) boundary condition at the walls to work better for me than natural(p)=0. Another thing is no matter how I'm trying to get a right value for the C constant, I cannot get the solution to converge. Dear Robert, could you please take a quick look at my problem if you will find some obvious mistakes ? The problem is very simple: the air enters on the left and is supposed to exit on the right + the internal piece of metal is supposed to be cooled.
| |||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 818 Registered: 06-2003 |
As I presume you are aware, Newton's method is only safely convergent when the starting point is sufficiently close to the solution that there are no local minima between it and the solution point. The solver must be able to "slide down" the functional until it hits the bottom. Unfortunately, this condition is not met for high-Reynolds number flow unless the starting velocity distribution is very close to the correct solution. The equations allow a multitude of eddy patterns, all of which introduce local minima in the path of the Newton search. One way out of this dilemma is to use the FlexPDE facility of STAGING to work toward the high-Re solution in a sequence of low-Re preliminary solutions. The easiest way to do this is to stage the density, starting at zero so the nonlinear terms disappear entirely, and staging up to the final value. The final value may not be reachable. Partly this is because the finite element equations become ill-conditioned at high-Re. The attached script is a modification of yours in which I have staged the density through a series of arbitrary jumps starting at zero. This script seems to run almost to completion, but appears to be failing to converge at the last stage. FlexPDE applies a Petrov-Galerkin weighting method to produce upwind differencing at high velocities. You could try SELECT UPFACTOR=2 (or perhaps more) to increase the weighting on the upwind correction terms. Don;t let this number get too large, or the solution could be distorted.
| |||
| Robert G. Nelson (rgnelson) Moderator Username: rgnelson Post Number: 819 Registered: 06-2003 |
PS. I didn't mean to imply in an earlier post that Backstrom's factor (mu/L0^2) should be removed. It is still useful in scaling the term and making the "C" factor more consistent across a range of problems. |
my_current_problem