synthetik-technologies / blastfoam

A CFD solver for multi-component compressible flow with application to high-explosive detonation, explosive safety and air blast
Other
215 stars 96 forks source link

Instability in 3D run after using rotateFields #43

Closed johng12 closed 2 years ago

johng12 commented 2 years ago

Dear blastFoam Team,

An instability seems to develop in my simulation when I use rotateFields to map from a 2D cylindrical geometry (axisymmetric wedge) to the full 3D domain. I am guessing one or more fields are not getting correctly initialized, but I'm not sure which.

Some additional details: 2D Case (detonation+shock propagation): Serial run with 2 active phases - water and tnt w/linear activation 3D Case: Parallel run with 3 active phases - water, tnt (basic JWL EOS only), air (free-surface)

At t=0 in the 3D run, everything looks correct. After the first time step, I see a large pressure spike develop at the original detonation location (i.e. charge center). See attached figures.

The command I am using with rotateFields is given here: mpirun -np 16 rotateFields -sourceTime 0.00027 -centre '(-0.588 0 0)' -fields '(p rho U e T rhoU rhoE lambda.tnt alpha.air alpha.tnt alpha.water rho.water rho.air rho.tnt)' ../2D_shock/ -parallel

From what I can tell, only a few of the fields listed are actually getting mapped to the 3D domain, but that may not be a problem. Any help or guidance you might be able to provide would be extremely helpful.

Thanks, John Gilbert rotateFieldExample_t0 rotateFieldExample_t1 rotateFieldExample_t2

jheylmun commented 2 years ago

Hi John, There were several bugs in the rotateFields utility, however they should be fixed now. Please let me know if you are running in to any other problems.

Thanks, Jeff

johng12 commented 2 years ago

Hi Jeff,

Thanks for the update. I've installed v5.0, but I'm still seeing some strange behavior after using rotateFields. I need to package up a good test case for you.

In the meantime, I tried running the mappedBuidling3D/wedge tutorial and I noticed the velocity field seems to be asymmetrical after the mapping from 2D -> 3D. I think it's most easily seen in the z-component of velocity. I've included two screenshots of Uz below, the first is at the end of the 2D run, and the 2nd is at the start of the 3D run. I would expect Uz to be symmetric about the axis of rotation. Can you tell me if this is expected behavior?

Thanks, John

axisymmetricCharge_wedge_Uz_orig

axisymmetricCharge_3D_Uz_orig

jheylmun commented 2 years ago

Hi John, I thought I had fixed this, but I guess it was just for the 1D case. I am assuming I just need a special condition for rotation from 2D to 3D. I will look in to it.

Thanks, Jeff

jheylmun commented 2 years ago

Alright, I limited the axis which the rotation is about, so it should only be rotating around the z axis. Based on the tutorial case it looks like it should be correct now. Does the fix help you cases as well?

johng12 commented 2 years ago

Hi Jeff,

Sorry it's taken me so long to get back to you on this. I am testing out the updated 'rotateFields' utility, but I'm still running into some issues. The figures below show pressure contours (1) immediately after mapping from a 1D-axisymmetric case, but prior to taking any time steps, and (2) a few time steps after the start of the simulation. Note the odd behavior at the center of the charge.

pContour_mapped_0 22ms pContour_mapped_0 23ms

I saw in the User Guide, that you recommend mapping the e field for multiphase simulations, so I thought that might be my issue. However, adding the e field to the -additionalFields list resulted in the following error message:


Calulating map and rotation tensors [100] [100] [100] --> FOAM FATAL ERROR: [100] Unknown patch type for patch outlet for e Perhaps the library is missing? include the necessary library with 'libs ("lib*.so")' in the controlDict


If I go into the e file, the patch type is listed as gradientEnergy:


outlet
{
    type            gradientEnergy;
    gradient        uniform 0;
    value           uniform 206730.612245;
}

Any thoughts?

Thanks, John

jheylmun commented 2 years ago

Hey John, I haven't updated the user guide since the improvements to rotate fields were made, mostly because I wanted to make sure it was all working first, I will be updating it once all of this has been straightened out though. I have made some significant improvements to how the energy is initialize though so you shouldn't need to map the "e" field any more. If you do want to map "e" though, in the controlDict, add the line

libs ("libblastThermodynamics.so");

it will include the additional boundary condition types that are used for "e".

I would recommend using the mappedBuilding3D case as a template. The option "-uniform" may be what is missing since that is responsible for handling the detonation points (aka not re-activating the original detonation point).

Thanks, Jeff