[Victor_R]

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.

U_field_inlet.jpg

Foam_geometry.png

T_instability_inlet.png

Fluid+foam.png

log_solution_example.jpg

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

  
 inlet
    {
        type           fixedValue;
        value          uniform 500;
    }

    outlet
    {
        type            zeroGradient;
    }

    fluid_to_foam
    {
        type            compressible::turbulentTemperatureCoupledBaffleMixed;
        value           $internalField;
        Tnbr            T;
        kappaMethod     fluidThermo;
    }

    walls
    {
        type            fixedValue;
        value           uniform 500;
    }

B.C for the solid:

internalField   uniform 400;

boundaryField
{
    foam_to_fluid
    {
        type            compressible::turbulentTemperatureCoupledBaffleMixed;
        value           $internalField;
        Tnbr            T;
        kappaMethod     solidThermo;
    }

    walls
    {
        type            fixedValue;
        value           uniform 500;
    }
}

fvSchemes for fluid ( i have tried upwind already):

ddtSchemes
{
    default         backward;
}

gradSchemes
{
    default         cellLimited Gauss linear 0.5;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss linear;
    div(phi,K)      Gauss linear;
    div(phi,h)      Gauss linear;
    div(phi,k)      Gauss linear;
    div(phi,omega)  Gauss linear;
    div(phi,Yi_h)   Gauss limitedLinear 0.1;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear limited corrected 0.5;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         limited corrected 0.5;
}

fvSChemes

ddtSchemes
{
    default       steadyState;// backward;
}

gradSchemes
{
    default         leastSquares;
}

divSchemes
{
    default         none;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fvSolution for the fluid

solvers
{
    rho
    {
        solver          diagonal;
    }

    rhoFinal
    {
        $rho;
    }

    p_rgh
    {
        solver           PCG;
        preconditioner   FDIC;
        tolerance        1e-7;
        relTol           0.01;
    }

    p_rghFinal
    {
        $p_rgh;
        relTol           0;
    }

    "(U|ki|h|epsilon|Yi|omega)"
    {
        solver           smoothSolver;
        smoother         GaussSeidel;
        tolerance        1e-7;
        relTol           0.1;
    }

    "(U|k|epsilon|Yi|omega)Final"
    {
        $U;
        relTol           0;
    }

    h
    {
        solver           PBiCG;
        preconditioner   DILU;
        tolerance        1e-7;
        relTol           0.1;
    }

fvSolution for the foam:

solvers
{
    "(h|e|rho|Yi|Ysi|p)"
    {
        solver           PCG;
        preconditioner   DIC;
        tolerance        1e-7;
        relTol           0.1;
    }

    "(h|e|rho|Yi|Ysi|p)Final"
    {
        $h;
        relTol           0;
    }
}

PIMPLE
{
    nNonOrthogonalCorrectors 0;

    correctMeshPhi no;
}

relaxationFactors
{
    equations
    {
        ".*"            1;
    }
}

[Bloerb]

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.

    div(phi,U)      Gauss linear;
    div(phi,K)      Gauss linear;
    div(phi,h)      Gauss linear;
    div(phi,k)      Gauss linear;
    div(phi,omega)  Gauss linear;

If your mesh is big try GAMG for p (and h) instead of PCG for a significant speed up (likely 50 to 80% faster solution time). And change nNonOrthogonalCorrector to 1. Why are your ddtSchemes steadyState and backward in fluid and solid? You are either using a steady state solver or you aren't.

[Victor_R]

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

[Bloerb]

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

ddtSchemes
{
    default       Euler;
}

gradSchemes
{
    default          Gauss linear;
}

divSchemes
{
    default         none;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluid

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss upwind;
    div(phi,K)      Gauss upwind;
    div(phi,h)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,omega)  Gauss upwind;
    div(phi,Yi_h)   Gauss upwind;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

and

PIMPLE
{
    momentumPredictor   yes;
    nCorrectors         3;
    nNonOrthogonalCorrectors 1;
}

relaxationFactors
{
    equations
    {
        "h.*"           1;
        "U.*"           1;
    }
}

[Victor_R]

Originally Posted by Bloerb

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

ddtSchemes
{
    default       Euler;
}

gradSchemes
{
    default          Gauss linear;
}

divSchemes
{
    default         none;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluid

ddtSchemes
{
    default         Euler;
}

gradSchemes
{
    default         Gauss linear;
    grad(U)         cellLimited Gauss linear 1;
}

divSchemes
{
    default         none;

    div(phi,U)      Gauss upwind;
    div(phi,K)      Gauss upwind;
    div(phi,h)      Gauss upwind;
    div(phi,k)      Gauss upwind;
    div(phi,omega)  Gauss upwind;
    div(phi,Yi_h)   Gauss upwind;
    div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

and

PIMPLE
{
    momentumPredictor   yes;
    nCorrectors         3;
    nNonOrthogonalCorrectors 1;
}

relaxationFactors
{
    equations
    {
        "h.*"           1;
        "U.*"           1;
    }
}

Dear Bloerb, thanks again. Indeed I am mixing things and my lack of expertise plays into this. I am interested in the steady-state behavior at the moment, but might be interested in the transient evolution of the system in the future. You now made it clear that using steadyState schemes needs special care to actually reach a converged solution at each time step (i.e. relaxation), and it makes perfect sense. You also said SIMPLE needs relaxation, but I assume you were talking about PIMPLE, which is inherently transient. I did not know that SIMPLE is much more permissive under the convergence point of view at each iteration, and thus I now understand why it doesn't blow up. It also makes a lot of sense since we only care about its final solution. Now that I implemented your suggestions (upwind, schemes and relaxation) and see that they don't work, I start to think that the solution obtained with SIMPLE is incorrect/ partial. The mesh is fine from checkMesh. I still have the same issues, obtaining non physical velocity fields after a few iterations. Find my U and p boundary conditions below, together with a log of my first 3 iterations. I see that i max out the iterations in the first two time steps (indicating poor convergence) but after that it never maxed out again. Regarding my U initial field, I am using a solution from potentialFoam which I have checked and seems to make sense. Since my mesh is complex i need to start with potentialFoam. The initial U condition given to potentialFoam is reported below and agrees with what i see in post processing results from potentialFoam with paraFoam. The issues start with rhoPimpleFoam.

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.


image_2023-02-17_193603363.png
p field

dimensions      [1 -1 -2 0 0 0 0];

internalField   uniform 101325;

boundaryField
{
    inlet
    {
        type            zeroGradient;
    }

    outlet
    {
        type            fixedValue;
        value           $internalField;
    }

    fluid_to_foam
    {
        type            zeroGradient;
    }

    walls
    {
        type            zeroGradient;
    }

initial U field to potentialFoam

dimensions      [0 1 -1 0 0 0 0];

internalField   uniform (0 0 0);

boundaryField
{
    outlet
    {
        type            zeroGradient;
    }
    fluid_to_foam
    {
        type            noSlip;
    }
    inlet
    {
        type            flowRateInletVelocity;
        volumetricFlowRate      2.1e-5;
    }
    walls
    {
        type            noSlip;
    }

[Bloerb]

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.

[Victor_R]

Originally Posted by Bloerb

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.

Dear Bloerb, thank you very much for your contribution. I implemented your changes and tested it by isolating the fluid region, using rhoSimpleFoam vs. rhoPimpleFoam. Now with under relaxation and outercorrectors, finally PIMPLE (from rhoPimpleFoam) seems to be working fine. My final issue is using (as you said) the foundation version, and the chtMultiRegionFoam solver to couple solid and fluid. Sadly, in this solver one cannot manipulate iterations on the fluid region independently of the solid (and maybe it's meant to ensure faster convergence of the coupled fluid-solid solution). Still, I would prefer to perform several outer corrections exclusively on my fluid zone before my solid zone is solved again (at the cost of heavy computational effort) but this doesn't seem to be possible. Since I had so much trouble stabilizing my fluid region with a fixed B.C. in rhoPimpleFoam, I imagine fluid-solid B.C. that change after each iteration are not going to help in getting a stable solution. Update: they don't , even operating in SIMPLE mode my simulation blows up. I am thus trying to understand which kind of relaxation I should use on my solid zone to reach a stable coupled solution, do you have any tip to offer?. Also, are you aware if chtMultiRegion only reads PIMPLE controls from the system/fvSolution (as it seems to ignore any nOuterCorrectors or other PIMPLE control I put in in the system/region/fvSolution file) or it also exclusively reads relaxation factors from system/fvSolution ignoring entries in the system/region/fvSolution file?

[Bloerb]

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.

[Victor_R]

Originally Posted by Bloerb

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.

Dear Bloerb, I've finally managed to stabilize my simulations. Thank you for your key contributions. The take away lessons you gave and that worked to use chtMultiRegionFoam (OF9 foundation version), although probably well-known good practices, are:

-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.