I have a unconfined circular impinging jet problem that has variable fluid properties and buoyancy effects present in the flow field. The jet is facing vertically upward along the axial direction. I have much trouble with this simulation when I add in the effects of buoyancy and variable fluid properties. The b.c. on the vertical stationary free surface is found from the continuity eqn. to be du/dr = -u/r in axisymmetric form. The hortzontal boundary is also a stationary free surface whose b.c. becomes dv/dz = 0, again from continuity eqn. The code is a collocated form of the SIMPLE method using the Power Law method for the convection terms. The flow is steady and laminar with a Re number through the jet of about 100. When I have no buoyancy and variable fluid properties the solution looks very good. When the other effects are added this problem turns into a nightmare! If the solution converges it does not look physically realistic. The program has already solved this problem with a coarser grid so I don't think there's a bug in the code, but now that I refine the mesh I have this problem. Even before I had much trouble solving this with the coarser grid but I managed to get to the solution. If anyone has any insight into this problem could you please contact me or response in this discussion group.

(1). I don't think it is a troubled CFD simulation. It is fairly typical. Recently, I have been doing a code check out for 3-D flows. Each time I change the mesh, I get a different result. And I also have the slow convergence problem. (and it is a commercial code.) So, welcome to the amazing world of CFD. (2). After having said that, I think we need to talk about the logic first. (3). It is fine to say that the code is bug free. (it is better to have a friend say that, after he has read through the listing.) So, this is definitely something you can do. (4). For your problem, the mesh arrangement should be fairly straightforward. Since you didn't mention anything about the mesh, it is a good place to look into. I mean, the mesh can be Cartesian (even though it is axisymmetrical.) with mesh lines parallel to z-axis and mesh lines normal to z-axis. With this system, you need to define the boundary locations, which you didn't describe in your message. (5). In other words, after you have defined your computational mesh, you need to define the computational domain and the boundary locations. (5). Cartesian mesh is not the only one for this type of problem. You can also use a streamline oriented mesh which closely follow the streamlines from the jet, through the stagnation region to the exit plane. (6). In other words, the mesh system selected will also change the boundary location and boundary conditions. (7). Once the mesh system is properly defined, the fundamental task is to solve the governing equations with the proper boundary conditions. In this area, it is important to select a method which you are sure that it can give you a converged answer, unless you are interested in making it a post-doctoral research topics. This is essential. The word "collocated" and "SIMPLE" are not brands in fashion. If it can not solve your problem, it is no good. (8). To find out whether the code is good or not ( or the method is suitable or not), the only way to do is to run some test cases with known solutions. (9). At this point, I would strongly suggest that you run some simple test cases first to isolate the problem. The following cases may be useful to locate the problem...a fully developed laminar flow at Re=100, a free jet flow with a fully-developed profile (parabolic) at Re=100, a confined jet problem with a fully-developed jet profile at Re=100, then the free jet impingment problem with a fully-developed jet profile. In each case, make sure that the mesh independent solution is obtained. (10). Once you have done these test cases, I think, you should be able to tell whether the code is good for your problem or not. The body force and the variable density are two extra factors which could change your solution and convergence. At that point, you can examine its impact on the solution through the coefficients and the source term of the governing equations. (11). As I always said, CFD is not just running a CFD code. Most of the time, you are getting the wrong answer. But if you follow the logic mentioned above, it is possible to make the road to the converged answer less trouble, even if it is not always possible to make it trouble-free. (12).The basic requirement to get involved in CFD is one must have the courage, it is applicable to students, engineers, managers, and company CEOs.

John: From your experience what are the limitations of this SIMPLE algorithm with the collocated grid arrangement? I have tested this code with many problems (i.e., driven cavity with and without buoyancy, 2D and axisymmetric channel flows, even this impinging jet problem gives good results when compared to the original paper without buoyancy and variable fluid properties). I have done another impinging jet problem where there were no stationary free surfaces except down stream of the nozzle where the exit was located. That problem also had buoyancy w/ variable fluid properties and gave good results. Is it possible that this collocated grid may have trouble with the stationary free surface b.c.'s? Please let me know what you think.

WIth the addition of buoyancy forces the coupling between the momentum and energy equations gets stronger. A relaxation of the equations, source terms or coefficients might be necessary.

(1). I have not figured out the exact definition of this "stationary free surface b.c.'s" yet. I have not seen it before. (2). First of all, you need to find out whether such boundary condition is consistent with the governing equations or not. (3). Then the next step is to see whether it requires the update from iteration to iteration. And ask a question about whether it will lead to convergence or not. Is it possible to get oscillating solutions? If the boundary condition is too weak, it will float between iterations. In this case, replace it by a stronger boundary conditions. (4). Or you could relax the updating frequency, that is to update the b.c. once every few iterations. In this way, you can keep a stronger control over the change at the boundary condition. (5). I don't know how you handle the variable properties or the body force term. But since you were able to get good solutions in other cases, perhaps the b.c. is the problem area.

John - Thanks for carrying on a discussion of my problem with me. To clear up the definition of what I mean by a 'stationary free surface' I hope I have not misled anyone since this is what I call the boundary, not the author of the original paper it was published in. In the original paper this boundary was referred to as a 'free boundary'. It is actually an imaginary bounday where the velocity is not initially known. The continuity eqn. is used to determine what the velocity is at every iteration at this imaginary boundary. The author does not say how he handles this numerically. I have to use my experience and judgement on its implementation into the code. Below is a sketch of the computational domain:

^ Z,V-vel

|

| flat plate surface

|--------------------------------------|

| |

|<----axis of symmetry |

| ^ | |

| | |<--nozzle wall |

| | | |

| | | free surface----->|

| | |----------------------------------|----->R,U-vel

| ^

inlet | velocity free surface

I hope this sketck helps visualize the problem. I use the b.c. du/dr = -u/r on the vertical free boundary and dv/dz = 0 on the horizontal free boundary. Currently, I use a first order backward difference for the vertical and first order forward difference for the hortzontal boundary in order to compute velocity at every iteration. Does this sound reasonable to you? What other way could I impose a 'stronger' b.c. at these free boundaries? Would fixing the boundaries for every 5 or 10 iterations before updating them be a better idea? Please let me know what you think.

Sincerely, Chuck

Sorry about the picture, Ithought it would look better than it does! I can't seem to sketch it with this browser.

(1). It's good enough for me. It is not very complicated. (2). I think the computational domain can be redefined and the boundary condition improved. (3). With the current configuration, the problem has two inlets to the computational domain, one is the small jet, where you have specified the velocity profile. The other is a very large inlet (from the jet wall to the outer vertical boundary), where the condition is totally unknown (not specified). So, the problem becomes a co-axial jet problem, with the outer jet very large and the velocity un-specified. It is hard to mix a lot of unknown with a small jet, and expect to get some converged solution. (4). So, you must change this condition from dv/dz=0 to v=V0, where V0 is a given velocity for the outer jet flow (a very, very weak jet) and it can be V0 >= 0.0. Actually, it is easier to set the V0 and U0 to zero or any fixed small number. When V0 and U0 are zero, the outer jet is now replaced by a wall. When V0 is a small given value, you have a low speed outer jet. In either case, the boundary is properly specified, and it is a strong boundary condition. (5). You can study the consequence of this boundary condition (different V0 value) and its location on the final solution. (you can move this boundary further down to reduce the influence on the top wall. (6). The outer boundary (the vertical cylindrical surface, also a vertical line in this case) can be handled in two ways. (7). The first approach is to increase the diameter of the outer boundary, r= R,outer where R,outer is much larger than the jet diameter. The advantage of this is when R,outer is large, du/dr will approach zero. This will explicitly remove the u from the equation, because u will be small, and r will be large, thus u/r will be even smaller. Now you have a zero gradient outer boundary condition. (8). The other approach is to put the whole problem inside a large cylinder with the top of the cylinder extended over the disk wall. This open top can be used as exit boundary. The boundary condition will be simply dv/dz=0.0 . Make sure that the current top disk is smaller than the size of the outer cylinder, so that the fluid can flow through the gap. Something like a co-axial jet problem inside a pipe with a disk wall at the center. (9). In this way, you can replace the very weak boundary into a well-defined problem. (10). The guideline is that always try to define the computational domain from the actual testing setup point of view. And try to keep the boundary further away from the point of interest.