
[Sponsors] 
June 18, 2012, 04:58 
Solver for hypersonic flow in OpenFoam

#1 
Member
Charles
Join Date: May 2011
Location: France
Posts: 77
Rep Power: 15 
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 steadystate solver which is fine for hypersonic flows ? Thanks, Charles 

June 19, 2012, 02:56 

#2 
Member
Join Date: May 2009
Posts: 32
Rep Power: 17 
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 multicore 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 

June 19, 2012, 03:26 

#3 
Member
Charles
Join Date: May 2011
Location: France
Posts: 77
Rep Power: 15 
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 

June 19, 2012, 04:34 

#4 
Senior Member
Olivier
Join Date: Jun 2009
Location: France, grenoble
Posts: 272
Rep Power: 17 
hello,
you may try to use coEuler / SLTS time scheme, then your time step will be higer (at least ~ x5 time faster ). regards, olivier 

June 20, 2012, 02:54 

#5 
Member
Join Date: May 2009
Posts: 32
Rep Power: 17 
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.


July 23, 2012, 07:21 

#6 
Member
Alexander
Join Date: Mar 2009
Posts: 49
Rep Power: 17 
If you are interesting in local time stepping (use "pseudotime") 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 :)


July 24, 2012, 02:52 

#7 
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 368
Rep Power: 20 
Hi Alexander,
I'm interested in rhoCentralFoam+pseudotime stepping for accelerating convergence to steady state. Can you please post your code? Thanks V. 

July 25, 2012, 07:16 

#8 
Member
Alexander
Join Date: Mar 2009
Posts: 49
Rep Power: 17 
There are three files in attachment (modified rhoCentralFoam.C, createFields.H, new setrDeltaT.H).
Below I post diff's rhoCentralFoam.C Code:
 ../../OpenFOAM2.0.0/applications/solvers/compressible/rhoCentralFoam/rhoCentralFoam.C<>20110616 16:48:41.000000000 +0400 +++ rhoCentralFoam.C<>20120725 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++; Code:
 ../../OpenFOAM2.0.0/applications/solvers/compressible/rhoCentralFoam/createFields.H<>20110616 16:48:41.000000000 +0400 +++ createFields.H<>20120725 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 +); + Code:
 ../../OpenFOAM2.0.0/applications/solvers/compressible/rhoCentralFoam/setrDeltaT.H<>19700101 03:00:00.000000000 +0300 +++ setrDeltaT.H<>20120725 18:50:06.000000000 +0400 @@ 0,0 +1,28 @@ +{ + + // Set the reciprocal timestep 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; + +} Code:
ddtSchemes { default localEuler rDeltaT; } Code:
 ../OpenFOAM2.0.0/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C 20110616 16:48:41.000000000 +0400 +++ /opt/openfoam200/src/finiteVolume/fields/fvPatchFields/derived/advective/advectiveFvPatchField.C 20120408 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() = 

July 25, 2012, 08:26 

#9 
Senior Member
Vesselin Krastev
Join Date: Jan 2010
Location: University of Tor Vergata, Rome
Posts: 368
Rep Power: 20 
Thanks a lot, I'll give a try to the modified code as soon as possible.
Regards V. 

August 27, 2012, 09:05 

#10 
New Member
Join Date: Aug 2009
Location: Stuttgart, Germany
Posts: 20
Rep Power: 16 
Hi Alexander,
I originally had a slightly different solution for the local timestepping 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 timestep 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 Last edited by ndr; August 28, 2012 at 02:20. 

August 28, 2012, 05:25 

#11 
Member
Alexander
Join Date: Mar 2009
Posts: 49
Rep Power: 17 
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. 

August 28, 2012, 08:24 

#12 
New Member
Join Date: Aug 2009
Location: Stuttgart, Germany
Posts: 20
Rep Power: 16 
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. 

August 28, 2012, 09:21 

#13 
Member
Alexander
Join Date: Mar 2009
Posts: 49
Rep Power: 17 
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?.. 

August 28, 2012, 10:43 

#14 
New Member
Join Date: Aug 2009
Location: Stuttgart, Germany
Posts: 20
Rep Power: 16 
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. 

August 29, 2012, 08:22 

#15 
New Member
Join Date: Aug 2009
Location: Stuttgart, Germany
Posts: 20
Rep Power: 16 
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. 

August 29, 2012, 11:06 

#16 
Member
Alexander
Join Date: Mar 2009
Posts: 49
Rep Power: 17 
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 

August 29, 2012, 14:13 

#17 
New Member
Join Date: Aug 2009
Location: Stuttgart, Germany
Posts: 20
Rep Power: 16 
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 speedup is still considerably high though and for quasi steadystate 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.... 

August 30, 2012, 03:20 

#18 
Member
Alexander
Join Date: Mar 2009
Posts: 49
Rep Power: 17 
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?


August 30, 2012, 09:31 

#19 
New Member
S. Mitra
Join Date: Oct 2011
Location: US
Posts: 1
Rep Power: 0 
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. 

October 21, 2012, 01:42 

#20 
Super Moderator

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 ?


Tags 
hypersonic wake 
Thread Tools  Search this Thread 
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Solver for an incompressible, turbulent flow with heat transfer  tH3f0rC3  OpenFOAM Running, Solving & CFD  9  June 17, 2019 06:12 
nano fluif flow in openfoam  anijdon  OpenFOAM  0  February 18, 2011 16:22 
What solver to use for inviscid flow simulation over missile  mecbe2002  OpenFOAM  0  April 27, 2010 11:10 
What solver to use for simulating flow behavior in inlet manifold  mecbe2002  OpenFOAM  4  April 9, 2010 08:01 
OpenFOAM Training in Europe and USA  hjasak  OpenFOAM  0  August 8, 2008 05:33 