[hadisafaei]

Hi All.
I'm trying to simulate impact of a droplet to a solid surface with the velocity of 100 m/s. I want to investigate effect of air compressibility, So i use compressibleInterFoam solver.
after some timesteps I get strange results. About time t=6e-8s the velocity suddenly raises up to 1000m/s !!!!
I guess this is because of continuty equation. so i check total continuty of phases by:
mdot = fvc::ddt(rho)+fvc::div(rhophi )
and summation in whole domain. It is Ok.
then i check continuty of each phases by:
mdot_phase1 = fvc::ddt(alpha1 , rho1)+fvc::div(alpha1rho1 , phi)
mdot_phase2 = fvc::ddt(alpha2 , rho2)+fvc::div(alpha2rho2 , phi)
and summation in whole domain. It is not Ok.
what could be the reasons ?
Can someone give me a hint on how can I get better results?

[hadisafaei]

this is my case
setup.zip

[richpaj]

Hello there,

I took a look at your case and thought some of the following points might
be of use:

1) the impact speed is very high (100m/s) in relation to the speed of sound for air at your chosen temperature and pressure, this immediately suggests that you should be using the transonic (i.e. implicit) formulation for the pEqn, you can switch this on in fvSolution with

PIMPLE
{
momentumPredictor no;// <-----possibly also
nCorrectors 3;
nNonOrthogonalCorrectors 0;
transonic yes;//<------
}

You will also need the following (or similar) in fvSchemes

div(phid1,p_rgh) Gauss upwind;
div(phid2,p_rgh) Gauss upwind;


and the specification of something like a GAMG solver for the new pEqn in fvSolution.

2) I'm not sure what "air compressibility" effects you are hoping to see in this simulation incidentally, especially since this flow is dominated by the inertia
of the falling droplet. Perhaps you want to try and track shock fronts
in the droplet itself? In which case you'd need maxDeltaT << droplet-diameter/sonic-speed.

You could monitor the impact force with a suitable function object in controlDict if only to confirm that the total momentum
delivered to the target wall roughly agrees with that of the incoming droplet.

3) Since this is a compressible flow, you might like to consider using a compressible thermo for both phases so that for phase1 you should probably use a "perfectFluid" EOS instead of the current "rhoConst".

4) The singularity in (p,U and the T) that develops around t=6e-08s most likely originates from the nature of the wedge geometry itself.
Instead of using a sharp wedge you might like to try a "blunt wedge", and if the singularity persists to simply increase the mesh refinement in both directions (so ensuring cells of unit aspect ratio).

You could also increase the initial distance of the droplet from the target to allow for more
"solution development" (i.e. ambient p, U adjustment) before the point of impact.

Alternatively, you might try different matrix-solver schemes in fvSolution for p
(ditto for the convection schemes used in the pEqn).


5)
At t=0

Liquid phase volume fraction = 0.00240267 Min(alpha.phase1) = 0 Min(alpha.phase2) = -1.15853e-05

and at t=1e-06s I see the following:

Liquid phase volume fraction = 0.00233543 Min(alpha.phase1) = 0 Min(alpha.phase2) = -1.69835e-05

so that the mass continuity error is not so bad (one part in 10^4) and this
despite the period of (p,U) overshoot.

A straighforward integration of the "explicit" form

fvc::ddt(alpha2 , rho2)+fvc::div(alpha2rho2 , phi)

is unlikely to provide an accurate assessment for mass continuity.


6) Another issue is possibly your choice of boundary conditions. If the endTime is fairly
short (~ domain-length/sonic-speed) then what you currently have is probably fine,
but for longer durations you might like to consider the effect of any reflected pressure waves.
(cf. waveTransmissive for zero reflection).

I hope some of the above might be of help.

Good luck!

[hadisafaei]

Thank you for your hints
I applied some changes to the case and here is the results.


1)I used transonic formulation for the pEqn. Unfortunately the results were worse than the previous state and velocity grown up to 1700m/s!!!


2)I have two purposes:
First I want to see effect of air compressibility to the air traping ( in droplet) in point of impact.
Second I want to simulate the impact of a droplet with an air cavity inside it(.i.e a hollow droplet) to a solid surface Due to high velocity of impact, effect of compressiblity of air is important.
I checked (maxDeltaT << droplet-diameter/sonic-speed) It was true.


3) The results were slightly better by using perfectFluid EOS.


4)I don’t know what you mean “ bulent wedge”. Instead I used a full 3D box domain but singularity point persists yet. Mesh refinment get me worse results!!!
Other changes ( Increasing initial distance of the droplet from the substrate and using different matrix-solver schemes) didn't Make significant changes in results. Finally after some time steps, this singular point spread in whole domain and then and the solution were diverge. I really do not know what is the origin of this singular point.
My guess this singular point occurs because of pressure boundary condition (fixedfluxpressure) ) in lower wall and also relation between pressure and density(EOS). To prove it I used "rhoConst" for both phases and saw the singular point was eliminated.




5)In compressibleInterFoam formulation, total mass continuity is obtained from mass continuity of each phase .pressure equation is obtained from total mass continuity equation which has been discussed in old posts. This formulation seen in this paper too:
http://www.sciencedirect.com/science...45793013001229

But in compressibleInterFoam Unlike its formulations, first pressure equation is solved then density of each phases is obtained by EOS of each phase and finally mixture density is calculated ( line 127 in pEqn.H (OpenFoam2.3))
This mixture density didn’t satisfy mass continuty and so compressibleInterFoam forced mass continuty manually by solve(fvm::ddt(rho) + fvc::div(rhoPhi)); ( line 91 in compressibleInterFoam.C )
Now my question is:
Is it necessary mass continuity for each phase to be established ? and how?
how can I eliminate singular point?
6) before shock wave has arrived the geometry boundary the solution were diverge.


thanks for any help

[richpaj]

Hello there,

Originally Posted by hadisafaei

1)I used transonic formulation for the pEqn. Unfortunately the results were worse than the previous state and velocity grown up to 1700m/s!!!

that's odd, when I ran some tests some time ago I thought the velocity
decreased to something like 800m/s. It's possible I also made some other changes but frustratingly I can't seem to locate the modified case I made.

Originally Posted by hadisafaei

I don’t know what you mean “ bulent wedge”. Instead I used a full 3D box domain but singularity point persists yet. Mesh refinment get me worse results!!!

I've attached an image of a blunt wedge (it eliminates the singular point along the axis). In order to compute a 3D geometry the computing resources required would be significant, whilst smaller & coarser meshes are liable to "evaporate" (owing to progressively large mass continuity errors).

Originally Posted by hadisafaei

My guess this singular point occurs because of pressure boundary condition (fixedfluxpressure) ) in lower wall and also relation between pressure and density(EOS). To prove it I used "rhoConst" for both phases and saw the singular point was eliminated.

The singularity is a combination of jumps in p-rho-U-rhoPhi
which might become manifest in one of the explicit terms in the pEqn.
I think it's possible to track through the code and find out where exactly
by creating relevant fields in createFields.H.
(I don't believe the singularity derives from the pressure BCs, fixedFluxPressure simply adjusts
the wall flux to account for gravity if it is present).

Originally Posted by hadisafaei

5)..
This mixture density didn’t satisfy mass continuty and so compressibleInterFoam forced mass continuty manually by solve(fvm::ddt(rho) + fvc::div(rhoPhi)); ( line 91 in compressibleInterFoam.C )
Now my question is:
Is it necessary mass continuity for each phase to be established ? and how?
how can I eliminate singular point?

I'm sure there are many approaches to this but in compressibleInterFoam
the following seems to be occurring:

i) "solve(fvm::ddt(rho) + fvc::div(rhoPhi))" smooths the estimate for rho
and correctly updates the BCs for rho.

ii)"alphaEqns" acts to enforce correct phase continuity.

iii) Phase flux rhoPhi is updated in "alphaEqnsSubCycle.H"

iv) In TEqn.H, mixture.correct() updates all the various phase properties


One approach to resolving your singularity might be to monitor some relevant
diagnostic fields esp. rhoPhi (use fvc::reconstruct(..) ) and some
of the source terms in the pEqn. (Incidentally, I presumable the case runs
fine at lower speeds?).

Sorry I can't do more than offer a few comments, a few too many deadlines
to deal with currently.


Good luck.