Hello,
I am working on mod
Hello,
I am working on modeling of mass-transfer between two immiscible liquids. Recently I started using OpenFoam. I have tried some simple simulations by manipulating the dambreak problem,(for example a denser liquid drop falling in a less denser liquid and also the rayleigh-taylor instability problem). I intend to use VOF for my project and it would be very helpful to know more about the method implemented in the present version of OpenFoam, i.e. is it CICSAM or a different one ? Also I would like to know if the option in interFoam - movingMesh is the same as Adaptive grid refinement around the interFace ? Thanks in advance. |
> i.e. is it CICSAM or a diff
> i.e. is it CICSAM or a different one
It's a new technique I developed a few years ago to resolve some of the fundamental problems of CICSAM and other traditional VOF interface compression methods. The differences have already been debated at length, have a look through previous threads on the subject. > Also I would like to know if the option in interFoam - movingMesh is the same as Adaptive grid refinement around the interFace ? No. |
I was wondering what's the ac
I was wondering what's the actual and preferred expression for the compressive term in gamma equation and a very brief explanation of how it is derived. The formula for 'phir' in the code differs from what Henrik has mentioned in (Eq. (4.15)). In the code:
phir=cGamma()*mag(phi/mesh.magSf())*interface.nHatf() while Henrik's suggests (if I have formulated it correctly): phir=cGamma*nhatf*max(nhat * phi / mesh.magSf()**2 ) Are they really different (am I missing something) and if yes, which one is better? 2) where is the smoothing function for 'gamma' in interFoam? Does it have a significant effect on VOF performance or just affects surface tension prediction? thanks |
I have tried various options
I have tried various options for the compression term in the gamma equation and the one I recommend is the one that is currently in interFoam. It is not derived, it is selected and there are other choices. If you would like to find out which option is best for you try them out on your case.
There is no "smoothing" function for gamma in interFoam, it was found to be detrimental to the overall performance of the interface capturing. |
Hi Henry,
Can I get any ref
Hi Henry,
Can I get any references for the VOF method you have implemented ? Thanks |
Sorry, I haven't written many
Sorry, I haven't written many papers; it's too time consuming and stops me doing the interesting work.
|
There is a decent description
There is a decent description of the implemented Interface-Capturing Methodology (probably a bit out of date now) in a PhD Thesis by dr. Henrik Rusche:
Computational Fluid Dynamics of Dispersed Two-Phase Flows at High Phase Fractions, Henrik Rusche, Imperial College of Science, Technology & Medicine, December 2002. Hrv |
That describes an old version
That describes an old version of my method, not the one currently in interFoam.
|
I'm currently having a look in
I'm currently having a look into the gammaEqn.H from OF 1.4.1
Well, some parts look like the equation (4.15) from Rusches PhD thesis. But still I'm not sure about some issues: - Where is the time derivative of gamma? - What is interface.nHatf()? Is it the normal vector to the cell face? - What doese MULES::explicitSolve01(gamma, phi, phiGamma) mean? Thanks so far! |
Dear Sebastian,
the time de
Dear Sebastian,
the time derivative is not there because of the usage of MULES. Look at OpenFOAM-1.5.x\src\finiteVolume\fvMatrices\solvers\MULES\MU LES.H for details: " MULES: Multidimensional universal limiter with explicit solution. Solve a convective-only transport equation using an explicit universal multi-dimensional limiter. Parameters are the variable to solve, the normal convective flux and the actual explicit flux of the variable which is also used to return limited flux used in the bounded-solution. " Best regards Christian |
The time derivative is account
The time derivative is accounted for within the MULES solver. In the gammaEqn.H only explicit fluxes are calculated, which are needed in MULES.
The interface.nHatf() represents a cell face unit interface normal flux. It is evaluated from the dot product of the cell face surface vector and the interface unit normal calculated at the cell face: nHatf_ = nHatfv & Sf, where nHatfv = gradGammaf/(mag(gradGammaf) + deltaN_). and deltaN is a stabilization factor for the case of gradGammaf = 0. For the full implementation look in src/transportModels/interfaceProperties/interfaceProperties.C Regards. |
Hi,
I have two questions rela
Hi,
I have two questions related to the interFoam solver: 1) I used the utility "barycenter" posted by Sebastian. It works fine with axisymmetric cases. In 2D or 3D cases and using the utility, it doesn't work (the position of the bubble is in decrease when time increase). Is there a signification for that? I find also that the position of the bubble increase when the value of Min(gamma) is negative. Is there a relation with Min(gamma) and the barycenter of the bubble?! MULES: Solving for gamma Liquid phase volume fraction = 0.99921 Min(gamma) = -2.20718e-11 Max(gamma) = 1 2) I want to switch to OF-1.5. In the file interFoam.C there is an additional line comparing to OF-1.4.1: p = pd + rho*gh; Is it a construction of the pressure p field or what? if yes, how can I get values of this field for different times. Best regards |
Hi,
I have two questions rela
Hi,
I have two questions related to the interFoam solver: 1) I used the utility "barycenter" posted by Sebastian. It works fine with axisymmetric cases. In 2D or 3D cases and using the utility, it doesn't work (the position of the bubble is in decrease when time increase). Is there a signification for that? I find also that the position of the bubble increase when the value of Min(gamma) is negative. Is there a relation with Min(gamma) and the barycenter of the bubble?! MULES: Solving for gamma Liquid phase volume fraction = 0.99921 Min(gamma) = -2.20718e-11 Max(gamma) = 1 2) I want to switch to OF-1.5. In the file interFoam.C there is an additional line comparing to OF-1.4.1: p = pd + rho*gh; Is it a construction of the pressure p field or what? if yes, how can I get values of this field for different times. Best regards |
Thanks so far for your respons
Thanks so far for your responses.
I will try to get them step by step (and respond if there are any questions). First of all: Is there a reference for MULES in some kind of citable form? |
Dear Edin.
Thanks for your
Dear Edin.
Thanks for your response. I have read your answer carefully and looked the code up. But I still have some questions related to the code. Is this the correct analytical representation of the code? http://www.cfd-online.com/OpenFOAM_D...es/1/10363.png I'm not sure what to do with the deltaN. In the code it's calculated like this: 1e-8/average(gamma.mesh().V()), 1/3) What does gamma.mesh().V() mean? See you (maybe at SLA ...) |
Sebastian,
I will come up t
Sebastian,
I will come up to you. The coefficient cGamma is not in the analytical transport equation, but it is used in modeling the relative velocity. gamma.mesh().V(), for a non-moving mesh simply takes the volumes of control cells throughout the domain, so the average of it to the power of 1/3 gives a representative cell dimension (lenght). Regards. |
Dear Edin.
Thanks for your
Dear Edin.
Thanks for your answere. So, the equation will take this form? http://www.cfd-online.com/OpenFOAM_D...es/1/10372.png |
On a side comment, I plan on u
On a side comment, I plan on using interFoam to do 2-D flows. Does anyone know how to implement proper BC?
|
Hello All,
One new question
Hello All,
One new question related to VOF solvers (interFoam). I modelled a ship sailing in a basin with water and air. The watersurface at the Inlet and internal has been set at Z=4 m. After a few seconds of physical time (so not CPU time) the solver crashes on very low deltaT. Looking the results in paraview one can see what happens: http://www.cfd-online.com/OpenFOAM_D...es/1/10436.jpg The picture shows the outlet of the domain. It looks like the watersurface is attracted towards Z=0 at the outlet patch. Maybe it might eventually go down further, provided that the solver should not crash. However, after that I translated the geometry such that the initial waterlevel is at Z=0 and kept the rest unchanged. Than the solution runs smooth as expected. BC's on the outlet: gamma and U: zeroGradient pd: fixedValue uniform 0 My questions: - is there any mechanism that forces a surface towards Z=0? - Has the problem described here something to do with continuity and/or mass conservation? - I set BC's on the top of the domain at fixedValue (vesselspeed) for U and zeroGradient for pd, this to reduce strong air vortices which influence solution stability. Should I better set the BC's on the top at atmosphere? Hope someone can shine his/hers light, not because I got stuck (I found a work around) but out of curiousity, Brgds, Mark Mark |
If this is an outlet than mayb
If this is an outlet than maybe the water is simply "leaking" out of it.
Maybe you can try fixedValue 0 for U at the out- and inlet? |
Hello Sebastian,
The inlet
Hello Sebastian,
The inlet has fixedValue U uniform (5 0 0) (approx). I can try to set this BC on outlet as well but I think this BC is too strong, it may influence the results too much. Isn't it very common to set fixedValue U on inlet and zeroGradient on outlet and for p the reverse? Brgds and nice days, |
Has anyone any comment on the
Has anyone any comment on the "bug" side of things, as opposed to the work arounds?
I have been working on similar flow setups, with the added complication that the air/water interface does not line up with the principle coordinate system. So the air/water interface is not at z=constant. So far I have not been able to achieve a stable simulation, which blows up, about after 1.5s of physically simulated time. Of course, I may have other issues with the case setup, but elimination of other possibilities is always a help. So if the interfoam or rasinterfoam solver buggy, if the air/water interface is not at z=0??? Thanks. albcem |
Albcem,
One thing which com
Albcem,
One thing which comes to my mind: is your gravity vector set correctly? Brgds, Mark |
Hello Mark,
Thanks for the
Hello Mark,
Thanks for the tip. However I double checked the consistency of the gravity direction. It checks out correctly. In my latest boundary condition settings, here is what I impose: U: atmospheric_wind_profile-boat_velocity at the air side, -boat_velocity at the water side. These inlet and outlet conditions are imposed at all vertical walls. Zero-gradient conditions at the top and bottom faces of the rectangular domain should allow for the flow to preserve conservation laws. (I get 1e-14 continuity error in each time-step) pd: Impose the pressure only at the top of the atmosphere following the Wigley hull example. k-epsilon: Impose inlet & inletOutlet types with non-zero k & epsilon values calculated according to Fluent manual recommendations using reference lengths typical of tested boats. These are imposed at vertical walls. At top and bottom patches, I have zero-gradient. ---- For all variables, I have also tested zeroGradient, inletOutlet types at several of the vertical and horizontal walls as sense allowed so. All have been fruitless. Possibilities on where things may be going wrong: The wind profile I impose, obeying a basic 1/7 power law, does not go to zero exactly at the air/water interface which may be causing high gradients there. However, I observe descent convergence on the simpleFOAM solution. So I have doubts on this argument. I use tetrahedral mesh, with seemingly fine enough resolution at the air/water interface (40mm->320mm variation). Can this be an issue? I have read papers in which tetrahedral mesh with such resolution seems to do fine... I may be introducing too much or incorrect turbulence into the flow messing things up. I will try a laminar solution to see if the blow up is delayed or eliminated - As the simpleFLOW solution converges I doubt this one as well... Out of ideas, anything else I can check? Would it help to upload the case somewhere? Albcem |
Hello Albcem,
I do not know
Hello Albcem,
I do not know your BC's. However the air side of boat-like case setups often show high air velocities, causing instabilities. In my cases I am not interested in the air side, I just set the same uniform fixedValue on U on the air side as on the water side. I only distinghuish these two by gamma. Attached a 0 directory containing my BC's. Only the gamma file is to be filled in for the actual case setup. I wrote a small utility for this (see elsewhere on this forum). It sets gamma 1 on internal field and Inlet for Z<=0 and gamma=0 elsewhere. I also experienced instable solutions using tet meshes. Since OF-1.5 I use snappyHexMesh, which I strongly suggest. However, we started with a question about Z=0 reference in VOF calculations. What is actually your problem and case setup? Brgds, Mark |
Hello Albcem,
I do not know
Hello Albcem,
I do not know your BC's. However the air side of boat-like case setups often show high air velocities, causing instabilities. In my cases I am not interested in the air side, I just set the same uniform fixedValue on U on the air side as on the water side. I only distinghuish these two by gamma. http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif pd http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif U http://www.cfd-online.com/OpenFOAM_D...hment_icon.gif gamma Attached a 0 directory containing my BC's. Only the gamma file is to be filled in for the actual case setup. I wrote a small utility for this (see elsewhere on this forum). It sets gamma 1 on internal field and Inlet for Z<=0 and gamma=0 elsewhere. I also experienced instable solutions using tet meshes. Since OF-1.5 I use snappyHexMesh, which I strongly suggest. However, we started with a question about Z=0 reference in VOF calculations. What is actually your problem and case setup? Brgds, Mark |
Hello.
Happy New Year to al
Hello.
Happy New Year to all of you! I have some further questions regarding the VOF methode in interFoam. Concerning the gamma-equation (gammaEqn.H): phir is limited to max(phic)! Is this the maximum phic in the cell or in the whole computational domain? Rusche wrote about limiting the phir in the transition area. This looks more like limiting it tho the whole domain ... Concerning the U-equation (UEqn.H): Am I right that fvVectorMatrix UEqn ( fvm::ddt(rho, U) + fvm::div(rhoPhi, U) - fvm::laplacian(muf, U) - (fvc::grad(U) & fvc::grad(muf)) ) is an equation for face values (before interpolation to the cell centers)? So ( fvc::interpolate(interface.sigmaK())*fvc::snGrad(g amma) - ghf*fvc::snGrad(rho) - fvc::snGrad(pd) ) * mesh.magSf() has to be of this kind as well to satisfy the momentum predictor? Since the snGrad-Schemes are face-interpolation-Schemes I think fvc::interpolate is one as well? So what is fvc::reconstruct used for? The Programmers and User Guide keep quiet about these two fvc-Schemes. Greetings. S. |
Let me put my question in anot
Let me put my question in another way:
How is the interface.sigmaK() discretisized? Is this a source term? Then I suppose it is calculated explicitly, since it is called within fvc::interpolate? But then it is a volume field? Shouldn't it be a surface field? |
Hello Sebastian,
I am also
Hello Sebastian,
I am also studying the code of interfoam. As I know of the interfoam solver. The sigmaK() is the product of the curvature and the surface tension coefficient. Since the surface tension effect is treated by CSF model, interfoam discretisized it as body force term. You can read the PhD Thesis by dr. Henrik Rusche P.154~P.155 for further infomation. Jim |
Thank you Jim.
Then, how is
Thank you Jim.
Then, how is this body force discretisized? As far as I got it, sigmaK() is explicitly calculated (fvc) at the cell faces. Does it need further treatment? Greetings. S. |
Hallo Sebastian,
As far as
Hallo Sebastian,
As far as I see interface.sigmaK() (product of surface tension and curvature) returns a volScalarField. These values are then interpolated to the cell faces with fvc::interpolate(). Different schemes are here possible for example linear (= central differencing). Michael |
Thank you Michael.
Ok, sigm
Thank you Michael.
Ok, sigmaK() returns a volScalarField following interfaceProperties.H: tmp<volscalarfield> sigmaK() const { return sigma_*K_; } And K_ is called with calculateK() in interfaceProperties.C: void interfaceProperties::calculateK() { const fvMesh& mesh = gamma_.mesh(); const surfaceVectorField& Sf = mesh.Sf(); // Cell gradient of gamma volVectorField gradGamma = fvc::grad(gamma_); // Interpolated face-gradient of gamma surfaceVectorField gradGammaf = fvc::interpolate(gradGamma); //gradGammaf -= // (mesh.Sf()/mesh.magSf()) // *(fvc::snGrad(gamma_) - (mesh.Sf() & gradGammaf)/mesh.magSf()); // Face unit interface normal surfaceVectorField nHatfv = gradGammaf/(mag(gradGammaf) + deltaN_); correctContactAngle(nHatfv.boundaryField()); // Face unit interface normal flux nHatf_ = nHatfv & Sf; // Simple expression for curvature K_ = -fvc::div(nHatf_); } As far is I see the gradGamma at the cell centers is interpolated to the cellfaces and named gradGammaf. This gradGammaf is normalized to unit length and named nHatfv. nHatfv is "dot producted" with the cell face vector Sf and named nHatf_. The curvature K_ is than declared as -div(nHatf_). If nHatf_ is calculated at the cell face how can sigma_*K_ be a volScalarField? I'm still not getting it. Where am I thinking wrong? |
Hallo Sebastian,
there are di
Hallo Sebastian,
there are different versions of fvc::div(). It is possible to have surfaceFields or volFields inside the brackets (). But as far as I see all have volFields as return value. Here we put in the surfaceScalarField nHatf_ and get back the volScalarField K_. It seems right that K_ is a volScalarField because most other variables (like pressure) are also definded as volFields. The mathematics behind the fvc::div() operation can be seen in the Programmer's Guide 2.4.5 / 2.4.10 (look at surfaceIntegrate). Maybe a better reference for what is done here is the PhD thesis of Onno Ubbink (equation 3.47). Greetings Michael |
Ok, Michael, Thank you.
It
Ok, Michael, Thank you.
It dawns me slowly - Ubbink (3.47) makes it some sense to me. Since all the gradients in in natfv are explicitly interpolated the fvc::div(nHatf_) is just a tensor operation and does not need any approximation? And the "outer" fvc::interpolation() is meant for the curvation at the cell center. But why interpolate it from the neighbours, if can be calculated from -sigma*div(nHatf_) ? |
Ok, Michael, Thank you.
It
Ok, Michael, Thank you.
It dawns me slowly - Ubbink (3.47) makes it some sense to me. Since all the gradients in in natfv are explicitly interpolated the fvc::div(nHatf_) is just a tensor operation and does not need any approximation? And the "outer" fvc::interpolation() is meant for the curvation at the cell center. But why interpolate it from the neighbours, if can be calculated from -sigma*div(nHatf_) ? What happens to the cell volume if the force due to surface tension is treated like a source term? |
Ok, Michael, Thank you.
It
Ok, Michael, Thank you.
It dawns me slowly - Ubbink (3.47) makes it some sense to me. Since all the gradients in in natfv are explicitly interpolated the fvc::div(nHatf_) is just a tensor operation and does not need any approximation? And the "outer" fvc::interpolation() is meant for the curvation at the cell center. But why interpolate it from the neighbours, if can be calculated from -sigma*div(nHatf_) ? It the force due to surface tension is treated like a source term, what happens to the cell volume? |
Hello again.
Thinking about
Hello again.
Thinking about the interface compression with phir I wondered about something. If the artificial velocity is directed into the direction of the interface normal, would this lead to a velocity directed to one side of the interface and compress it only on this side? Is there an illustration of the behaviour of such a scheme to the interface? S. |
Well, the answers are getting
Well, the answers are getting thinner.
Lets try an easier one: Using VOF means using blended density and viscosity in the Navier-Stokes-Equations: rho = gamma*rho1+(1-gamma)*rho2 Is that applicable to the velocity, too? u = gamma*u1+(1-gamma)*u2 |
Is that applicable to the velo
Is that applicable to the velocity, too?
u = gamma*u1+(1-gamma)*u2 That's a good question. But one should be clear about the two-fluid model concept. In two-fluid model (VOF here), only one velocity for both phases. NO u1 and u2. |
Then, why is Jasak stating thi
Then, why is Jasak stating this in his slide
http://powerlab.fsb.hr/ped/kturbo/Op..._12Jan2005.pdf on page 4? |
All times are GMT -4. The time now is 06:36. |