Unphysical temperatures in chtMultiRegionFoam
5 Attachment(s)
Dear everyone, I am trying to solve a somewhat straightforward case where a fluid is flowing around a solid zone, where the solid has a quite complex geometry (it's a 3-d printed structure). I attach pictures of the geometry with a short explanation: the solid region intersects the external walls of the fluid zone, which is cylindrical. In the picture where you only see the foam, you can see the external walls of the foam colored white. In the picture with both foam and fluid, you can see the fluid walls colored white and the foam external walls colored green. Thus the external cylinder contains two separate external walls, i.e. those of the fluid region and those of the solid foam. You will see them mentioned as 'walls' in B.C of both regions. The rest is quite straightforward. I am using standard B.C for pressure and velocity, namely fixed pressure inlet and zerogradient on the outlet and a mass flow rate for inlet velocity with a zero gradient on the outlet. I am applying arbitrary B.Cs to both the foam and fluid, just to test whether i can achieve a stable solution with chtMultiRegionFoam. I am however obtaining unphysical temperatures in the fluid zone after short simulated times (e.g. 0.01s if I use backward time discretization schemes in the solid zone, or 0.001s if i use steadyState time schemes in the solid zone).
The problem seems to be in the energy equation and velocity solution coupling. Given my B.C., there is no possible way in which the fluid should reach temperatures above 500K, since the external walls and inlet are at 500K themselves and that should be the maximum, yet it does. My crash is always due to Equation of state exceeding the range of 300-2000K at some point. When this happens, the Courant number exceeds 1 for obvious reasons, but nothing strange occurs beforehand. After some work, it seems that the issue is for some reason close to the inlet, and I attach a picture of a Threshold function of ParaFoam where I have isolated all cells above 500K after the instabilities incurred in one of my simulations. Checking the velocity field, I also get velocities that are way too high and unphysical solutions (e.g. 12m/s in the zone of unphysical temperature, whereas my inlet is at uniform 0.2m/s), also visible in one of the pictures attached. I also attach an example log of what the solution looks like when instabilities have already arised (although not making the equation of state crash because T is below 2000K). No solver is maxing out its iterations thus I have not yet understood where the issue lies. As you can see fluid zone max T is already above 500K, which is already unphysical with my given B.Cs. Attachment 93342 Attachment 93343 Attachment 93344 Attachment 93345 Attachment 93346 As you can see, they are all at the inlet of the fluid zone, which does not even touch the solid foam (as there is a small buffer zone between both fluid inlet and outlet and the foam itself). I have already tried different time and space discretization schemes, and the last I tried should be one of the most stable, also attached below. Ignore the extra entries that you see for some variables, those are for a custom solver I am not currently using and they do not create any issue in the solution of the solver. I have also obtained stable solutions by isolating the fluid zone mesh and running it on rhoSimpleFoam with a fixed foam_to_fluid wall temperature. It would seem that in cht my velocity-energy coupling is messed up for some reason. B.C for the fluid Code:
Code:
internalField uniform 400; Code:
ddtSchemes Code:
ddtSchemes Code:
solvers Code:
|
I have not read everything but you are using linear schemes. This is highly unstable and can lead to such temperature swings way above what is physically possible. linear schemes are second order but not limited to anything. Hence they explode very easily. Switch those to upwind for first order (recommended until you get a first stable solution) or some second order accurate scheme with some limiter. Like linearUpwind or limitedLinear. Or even vanLeer etc etc. And if you are solving a steady state try the bounded argument as well "bounded Gauss upwind" for example.
Code:
div(phi,U) Gauss linear; |
Dear Bloerb, thank you for your insights. I have tried upwind schemes and they haven't solved my issue. You are right about the steadyState schemes but I tried with and without and there is no big difference. However, I seem to have identified a further key element of what causes my problems: velocity solution in the PIMPLE algorithm. Even if I run my solver for a single time step, I already get a velocity field which is messed up (as the one in picture above, with high velocities near the inlet). I have encountered this issue only with PIMPLE, as if I run a single time step with SIMPLE with exactly the same settings (fVsolution and fvSchemes), and even if I let it run until steady state, I get very reasonable and physical results both for velocity and temperature. For the comparison between PIMPLE and SIMPLE I isolated my fluid mesh region, imposed fixed temperature B.Cs on the wall and run rhoPimpleFoam and rhoSimpleFoam. I have tried a variety of schemes (although not all yet), and rhoSimpleFoam works perfectly well, while rhoPimpleFoam gives unreasonable velocity fields after already a few time steps. If someone has any experience of these two giving large differences on the same case setup, let me know please
|
Are you solving for a steady state or not? That is the question you need to answer. Because you are mixing a lot of different concepts. If you are, you need to add relaxation factors. You can't just switch between steadyState and backward in one domain. That is absolutely nonphysical.
You are currently not using any kind of relaxation. The difference between SIMPLE and PIMPLE however is that SIMPLE needs relaxation and you typically do not converge each time step in OpenFOAM solvers. It is hence no wonder that simpleFoam does not blow up but PimpleFoam does. Solving these problems is always the same approach. First switch to first order schemes. Trying schemes is useless. upwind is the most stable scheme there is. If it blows up it is not you schemes. So upwind for everything. Next run checkMesh and make sure your mesh is not beyond good and evil. I am going to assume you are indeed interested in the unsteady behaviour. Try these settings. If they do not do the trick, you could post your U and p_rgh boundary conditions and your thermoPhysicalProperties. solid Code:
ddtSchemes Code:
ddtSchemes Code:
PIMPLE |
1 Attachment(s)
Quote:
Two last questions I have now are 1) how can SimpleFoam avoid convergence, keeping track of iterations and residuals at each step? Since it provides final residuals and iterations, I would normally think it works exactly like PimpleFoam in terms of equations convergence (i.e. stopping when a certain residual is reached). Is this really not the case, and the solution of equations just stops too early to capture my 'unstable' velocity behavior? 2) can the relaxation of equations affect my solution meaningfully? In the example you gave, there is zero effect of relaxation given that the relaxation factors are 1, but if I were to change that, or even relax the fields instead of the equations as you proposed, could it possibly lead me to physical solutions of the velocity field? I don't think so (but might be wrong), because my very issue is that I reach a converged velocity field, albeit a non-physical one. Attachment 93394 p field Code:
dimensions [1 -1 -2 0 0 0 0]; initial U field to potentialFoam Code:
dimensions [0 1 -1 0 0 0 0]; |
To solve for the steady state behaviour use chtMultiRegionSimpleFoam. If you are using the foundation version this solver does not exist. The behaviour was added to chtMultiRegionPimpleFoam. But you need to be careful with that. Otherwise it'll blow up there are a few topics out there about that. You basically need to set the correct relaxation factors.
PIMPLE is SIMPLE+PISO. So the first nOuterCorrectors are SIMPLE and the last is PISO. You have not specified any nOuterCorrectors, so they are 1. Hence you are solving in pure PISO mode. There are two kinds of relaxation factors U and p_rgh for the first loops and UFinal p_rghFinal for the last loop. SIMPLE and PISO are basically the same. The main difference is that simple solves the system just once and piso solves it nCorrector times. SIMPLE uses relaxation and is hence not bound to a time step or CFL number. But solving a system once and using relaxation to make sure it does not blow up looses you the time accuracy. Hence you can only use it to get to a steady state. Or you do dozens of simple loops per time step instead of one. That is the whole point. Relaxation will smooth out the solution. So if a steady state is possible you'll reach it quicker. But if there oscilations for example, the relaxation will smooth those out to reach that steady state. So the answer is simple. Switch to chtMultiRegionSimpleFoam or use the steady state settings for the foundation version. Which would mean adding relaxation of 0.7 and 0.3 for UFinal and p_rghFinal. |
Quote:
|
My typical advice is to not use any relaxation for h or e at all or 0.99 if needed. Because it is rarely if ever necessary. At least assuming everything else was set up properly. And for U, p and turbulence properties use their respective values from the tutorials. The h default relaxation in many tutorials is 0.7 and this will take ages to converge in cht settings. I personally never had to use relaxation ever and have solved multi million cell cases with dozens of fluid and solid regions coupled together with complex flow on awful tet meshes. I found it only necessary when starting a transient run without any steady state initialization of the flow field.
If you can subcycle the fluid/solid regions depends on your OpenFOAM version. The easiest would be to look at the source code of the chtMultiRegionFoam solver. Typically it is the fluid region that is slowing everything down. And there have been some implementations for extra loops added in there. But since the current foundation version replaced the solver construct entirely with something more flexible I can't say. What's possible in your version. |
Quote:
-operate in simpleFoam mode of the solver by leaving the PIMPLE {} subdictionary in system/fvSolution empty. -check relaxation factors of your U,p_rgh fields. They gave me some trouble in obtaining a stable solution. -Use upwind schemes initially before switching to second-order accuracy. -use steadyState time discretization schemes. -Initialize your flow field with a reasonable solution. This was the most crucial factor for me, as using a U field solution from potentialFoam was not able to yield a stable solution in my case. I thus isolated my fluid region and applied a fixed wall temperature B.C. (instead of coupling fluids and solids) and ran rhoSimpleFoam( simpleFoam works equally well) until convergence (you still need relaxation with simpleFoam). By starting chtMultiRegionFoam with the latter converged velocity field, the simulation ran very smoothly. For the moment, my issues are solved. I hope i don't have to come back to ask again for help too soon. |
All times are GMT -4. The time now is 18:27. |