[Santos-Dumont]

Hi,

I'm trying to simulate the wake behind an object at hypersonic speeds (around M=15) with OpenFoam.
The solver SonicFoam seems inappropriate for hypersonic speeds because the shock is too close to the object. So I used rhoCentralFoam which gives nice results.

The problem is that I have a very huge domain as I'm interested in the far wake behind an object and rhoCentralFoam is a transient solver. My first case (half of a sphere with a domain lenght of 2000 times the sphere diameter) needed 10 days of calculations on 7 CPU cores. Even if the results are nice the calculation time is too important.

Which solver would you recommend for such a case ? Is there a steady-state solver which is fine for hypersonic flows ?

Thanks,
Charles

[KrisT]

I don't think there are too many options here:

In principle you could save time by using adaptive mesh (maybe you already do?), but when I last tested (on rhoCentralFoam), the capability in OpenFOAM it couldn't handle multi-core simulations in an efficient manner, i.e., there was no redistribution of cells after refinement. Which resulted in that - when using several cores - one didn't gain much (on single core it was working well though).

I would be very interested to hear if there has been improvements in OpenFOAM in this respect (did my tests on 1.5 or 1.6).

Apart from this, try to find more cores ;-)

/jon

[Santos-Dumont]

Hi Kris,

I think I could try coarsening the mesh or adapt it but I might reduce the precision of the results as I am still trying to get a mesh convergence on that case (using another CFD software).

I had no particular trouble running it in parallel. I didn't make a measurement on the actual speed gained by the parallel run but I tried to run it on a single core and it was way, way much slower.

Charles

[olivierG]

hello,

you may try to use coEuler / SLTS time scheme, then your time step will be higer (at least ~ x5 time faster ).

regards,
olivier

[KrisT]

Hi! Sorry, but I meant you should try the _adaptive_ mesh refinement, so that mesh is updated during the calculation (according to some measure - involving gradient of density or so). That is, mesh is refined around shocks, and more coarse where not much happens in the field. But as I said, not sure if this is efficient when running in parallel.

[sahas]

If you are interesting in local time stepping (use "pseudo-time") there is some information in the forum about its implementation in rhoCentralFoam. I have done it for rhoCentraFoam2.0, and I can post here the code if demanded :-)

[vkrastev]

Hi Alexander,
I'm interested in rhoCentralFoam+pseudo-time stepping for accelerating convergence to steady state. Can you please post your code?

Thanks

V.

[sahas]

There are three files in attachment (modified rhoCentralFoam.C, createFields.H, new setrDeltaT.H).
Below I post diff's

rhoCentralFoam.C

--- ../../OpenFOAM-2.0.0/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C<---->2011-06-16 16:48:41.000000000 +0400
+++ rhoCentralFoam.C<-->2012-07-25 18:46:33.000000000 +0400
@@ -35,6 +35,7 @@
 #include "turbulenceModel.H"
 #include "zeroGradientFvPatchFields.H"
 #include "fixedRhoFvPatchScalarField.H"
+#include "fvcSmooth.H"
.
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
.
@@ -149,8 +150,9 @@
         // estimated by the central scheme
         amaxSf = max(mag(aphiv_pos), mag(aphiv_neg));
.
-        #include "compressibleCourantNo.H"
         #include "readTimeControls.H"
+        #include "setrDeltaT.H"
+        #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
.
         runTime++;

createFields.H

--- ../../OpenFOAM-2.0.0/applications/solvers/compressible/rhoCentralFoam/createFields.H<------>2011-06-16 16:48:41.000000000 +0400
+++ createFields.H<---->2012-07-25 18:48:04.000000000 +0400
@@ -111,3 +111,19 @@
         thermo
     )
 );
+
+volScalarField rDeltaT
+(
+    IOobject
+    (
+        "rDeltaT",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh,
+    1/dimensionedScalar("maxDeltaT", dimTime, GREAT),
+    zeroGradientFvPatchScalarField::typeName
+);
+

setrDeltaT.H

--- ../../OpenFOAM-2.0.0/applications/solvers/compressible/rhoCentralFoam/setrDeltaT.H<>1970-01-01 03:00:00.000000000 +0300
+++ setrDeltaT.H<------>2012-07-25 18:50:06.000000000 +0400
@@ -0,0 +1,28 @@
+{
+
+    // Set the reciprocal time-step from the local Courant number
+    rDeltaT.dimensionedInternalField() = max
+    (
+      1/dimensionedScalar("maxDeltaT", dimTime, maxDeltaT),
+      fvc::surfaceSum(amaxSf)().dimensionedInternalField()
+              /(maxCo*mesh.V())
+    );
+
+    Info<< "    Flow        = "
+            << gMin(1/rDeltaT.internalField()) << ", "
+            << gMax(1/rDeltaT.internalField()) << endl;
+
+    // Limit the largest time scale
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    rDeltaT.max(1/maxDeltaT);
+
+
+    // Spatially smooth the time scale field
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+    scalar rDeltaTSmoothingCoeff = 0.1;
+    fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff);
+
+    Info<< "    Overall     = " << min(1/rDeltaT).value()
+        << ", " << max(1/rDeltaT).value() << nl << endl;
+
+}

To activate Local Time Stepping (LTS) you should specify the next option in system/fvSchemes:

ddtSchemes
{
    default         localEuler rDeltaT;
}

Apropos if you need the "waveTransmissive"boundary condition with LTS, you should correct the file OpenFOAM-2.0.0/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C in the next way

--- ../OpenFOAM-2.0.0/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C    2011-06-16 16:48:41.000000000 +0400
+++ /opt/openfoam200/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C    2012-04-08 21:32:23.000000000 +0400
@@ -30,6 +30,7 @@
 #include "EulerDdtScheme.H"
 #include "CrankNicholsonDdtScheme.H"
 #include "backwardDdtScheme.H"
+#include "localEulerDdtScheme.H"
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -226,6 +227,7 @@
         (
             ddtScheme == fv::EulerDdtScheme<scalar>::typeName
          || ddtScheme == fv::CrankNicholsonDdtScheme<scalar>::typeName
+         || ddtScheme == fv::localEulerDdtScheme<scalar>::typeName
         )
         {
             this->refValue() =

Also in the version 2.0 there is a bug with Prandl number, which is alway 1 (see bug http://www.openfoam.org/mantisbt/view.php?id=475)

[vkrastev]

Thanks a lot, I'll give a try to the modified code as soon as possible.

Regards

V.

[ndr]

Hi Alexander,

I originally had a slightly different solution for the local time-stepping implementation in rhoCentralFoam for OF 2.1.x, only my rDeltaT was calculated a bit different.

However, when I tried to validate the solver for a simple flat plate case with a laminar Mach 2 flow, the results were not completely satisfying. So I gave your solver a try and interestingly nearly the same behaviour occurred: Especially the pressure field is fluctuating wildly downstream of the shock wave which originates at the plate's leading edge. I attached the contour plots of the pressure field both for the case using LTS and for exactly the same setup using the normal timestepping. Without LTS only one shock wave is captured as it should be. For the boundary conditions: The flow is coming from the left with uniform velocity, pressure and temperature. I used wave transmissive for the right outlet and the top, the plate at the bottom starts at 50 cm (before that a symmetry plane is used). For the reconstruction schemes I used the Gamma scheme.

As I had used the original setup for validation of the boundary layer I also had a look into that. And as you can see in the velocity and temperature profiles across the BL at a length of 1 m (red is without LTS, blue is with LTS) those profiles are not really identical... The red profiles are validated against literature and an analytical solution, so they are definitely correct.
I already tried to change the time-step smoothing, but that had no effect at all.

Do you have an idea where this fluctuation may come from? Or did you maybe observe similar problems when testing your solver?

Regards

Nils

[sahas]

Hi Nils,

what Courant number "maxCo" do you use? I found that sometimes there are fluctuations in the fields when Courant number is greater than approximately 0.1.

[ndr]

Hi Alexander,

ok, I didn't realize I had to take a Courant number that low. I ran my calculations as usual with a maxCo of 0.4, but I also tried 0.2 without any improvement in the results.
But as you suggested I will also try 0.1 (and 0.05) and will let you know if that helps with the fluctuations.

[sahas]

Also, the result may strongly depend on interpolation scheme used. For example, vanLeer sometimes gives oscillations, whereas upwind does not. As well you may try SFCD scheme.
P.S. Of course, interpolation scheme is weakly related to ddt scheme, but who knows exactly?..

[ndr]

Yes, I already noticed that vanLeer seems to be unstable for certain problems, especially when dealing with hypersonic boundary layers.
But until now the Gamma scheme always worked fine for me, so I think I will try to keep as many parameters constant as possible (when comparing localEuler and Euler) so that differences in the results can easily be traced back to the ddt scheme.

[ndr]

Hi Alexander,

thanks for the hint regarding the maximum Courant number, that did the trick. If I use maxCo 0.1 or lower it works perfectly well.
It's a bit of a drawback that the number of iterations now increases considerably compared to the original Euler case, which was calculated at maxCo = 0.4 and without any problems. But the LTS still seems to converge faster to a steady state by a factor of about 2.

[sahas]

Glad to hear good news, Nils

As I understand, the "drawback" is a kind of "phantom". Look: when you use Euler scheme you limit the time step corresponding to Courant number 0.4 in some cell (may be or greatest, or smallest, or anything else). But in another cell Courant number corresponding to chosen time step may be much less than 0.4 (for example, 0.04) - and this cell may be source of oscillations. Whenever you use localEuler, in all cells (almost) you have Courant number 0.4, also in "problem" cell, and the critical Courant number for it cannot be greater than 0.1. This is just my theory, I did not have "hard look" on it

[ndr]

Yes, I was aware of that. However, due to my grid structure the size of the cells in flow direction (which should be mainly responsible for Courant number calculation as there is hardly any crossflow) was approximately the same all over the grid, only the wall normal distance was refined further.

So when I ran the Euler calculation I indeed did not get a Courant of 0.4 all over the grid, but the discrepancy between the cells was not that big. If I remember correctly the minimum Courant number was between 0.05 and 0.1, but that occured in a region far away from the wall and any shocks, so those cells should not be critical for any oscillations.

But it seems that somewhere around Co = 0.1 there is the upper limit for LTS, as I tested my old version of the solver again and it showed the same behaviour. The speed-up is still considerably high though and for quasi steady-state calculations it seems to be the only reasonable choice for rhoCentralFoam, all other simulations may take ages to compute depending on your geometry and flow type....

[sahas]

By the way, inside boundary layer the diffusive term has significant value so I think one should consider also the diffusion Courant number (mu*dt/dx^2) - may be, it is responsible for arising oscillations?

[SMitra]

Hi,

You may try to use normalized properties for this compressible problem. That way, you may be able to increase the time step that you need to solve the problem. This will give you quicker qualitative results for your problem.
For correct quantitative solution, you will have to use the proper properties.

Hope this helps.

PS: For normalized prop: look at "thermophysicalproperties" in the ForwardStep problem in sonicFoam tutorial.

[praveen]

Originally Posted by sahas

By the way, inside boundary layer the diffusive term has significant value so I think one should consider also the diffusion Courant number (mu*dt/dx^2) - may be, it is responsible for arising oscillations?

The diffusive terms in momentum and energy equation are treated implicitly. So diffusive time step should not be a problem. Is the mesh made of quadrilaterals or triangles ?