[Gfs-devel] Update on Mass Transfer Across the Interface

[Gael Guedon]

Hi all,

it has been a while we did not discuss about the mass transfer issue. I finally come to some improvement recently regarding the volume source using a 1D version of the single phase Euler solver of Gerris.

I succeeded in reducing significantly the spurious oscillations in velocity and pressure following the observations and analyses presented in the following references:

[1] M. F. Lai, A projection method for reacting flow in the zero Mach number limit. PhD thesis, University of California, Berkeley, 1993.

[2] W. J. Rider. Approximate projection methods for incompressible flows: Implementation, variants and robustness. Technical Report LA-UR-2000, Los Alamos National Laboratory, 1995.

[3] W. J. Rider, Filtering non-solenoidal modes in numerical solutions of incompressible flows. Int. J. Numer. Meth. Fluids, 28:789-814, 1998.

Three main modifications are applied:

1) an additional projection is performed to compute the scalar gradient field of the divergent part of velocity;

2) the projection is modified so that pressure corrects for the solenoidal part of velocity only;

3) a vertex-divergence projection filter is applied to improve coupling between the divergent part of velocity and the volume source (as suggested in [1]).

I hope finalizing a document which summarizes the approach and the results obtained soon. I am still not sure if this may work at the end because I am working with a 1D version... however it seems to be promising.

Modifying the projection should not be a problem and I am more concerned about the vertex-centred divergence on adaptive grid which may need appropriate interpolations (especially if we would like to include solid).

As a first step I will try to implement a 2D version on uniform grid to see if the method is consistent.

José, as I remembered you were working on the correction of the VOF transport, I am curious to know at which point you are now.

I will keep you informed.

Cheers,

Gael

[Jose M. López-Herrera Sánchez]

Hi all,

2012/7/26 Gael Guedon <gaelg...@gmail.com>

I hope finalizing a document which summarizes the approach and the results obtained soon.

I am willing to see your memo.
 

José, as I remembered you were working on the correction of the VOF transport, I am curious to know at which point you are now.

I apologize by not answering you before but I have been on holidays the last two weeks of July. I have implemented a new source type (GfsSourceInterfaceFlux) what serves to impose a flux of a variable perpendicular to an interface defined by a TracerVOF. It usage would be:

SourceInterfaceFlux Variable intensity_function TracerVOF

what serves to compute a source of the type:

int_{interf} intensity_function \delta_{kronecker} \mathbf{n} ds

It uses for the calculation of the normal of the interface the available interface reconstruction data (m \cdot x = \alpha) so not smoothing is performed. It can be used to either correct the VOF transport equation  what writes

F_t + v \cdot \nabla F = int_{interf} (\dot{m}/\rho) \delta_{kronecker} \mathbf{n} ds

Therefore to correct the VOF transport equation it should be added:

SourceInterfaceFlux T (\dot{m}/\rho)   T

and, in principle,  for adding a not-smoothing volume source by means of: 

SourceInterfaceFlux P function_2  T

I have focused on its usage as a corrector of the VOF transport equation so I have compared to the gael implementation. In reality both treatment are equivalent since we both solve the same equation,

F_t + U_I \cdot \nabla F = 0

U_I is the velocity of the interface being U_I = v + v_{relative}. v is the fluid velocity and v_{relative} is the relative velocity of the interface. v clearly depends on the volume source imposed (or its absence see D. Lakehal et al. / Int. J. Heat and Fluid Flow 23 (2002) 242–257) .  What is not clear (unless for me) is what is really v (or what it should be). Clearly is an average velocity but I can not figure out a clear physical interpretation of it. I suspect that this unclear aspect contribute to the difficulties of interpreting the pressure values discussed in the thread:

  https://mail.google.com/mail/u/0/?shva=1#search/volume+source/1372ed2329af2eff

To check the performance of both implementation I have use an ellipse that grow at a constant rate being the fluid velocity zero, so the form has to be maintained.

Some comments:

• The Gael implementation works fairly well but I see difficulties if TracerVOFConcentration are wanted to be used.
  
• My implementation has the strength that does not imply large changes of the existing code but in the present form (because it is naively implemented) it is prone to suffer  under/over shoots. The splitting character of the VOF advection scheme does not give any problem.
  

 I am dubious about what line I should follow, i.e. to abandon my line centering on the gael's line working into make it compatible with the TracerVOFConcentration or to keep working in mine trying to skip over/under shoots problems (I have in mind to shift the plane using the value of alpha to see if a neighbouring full cell has to be changed to be an interface cell). Certainly, the implementation of gael is more mature than mine. Any feedback is welcomed!

Attached you find:

• A patch with the GfsSourceInterfaceFlux
• A patch which contains the implementations of Gael
• simulations files of the ellipse test.

Finally a question:
It is about the available T_x, T_y and T_z there is any reason to be different of zero in full cells? This is why you, Gael, prefers define the yours T_nx, T_ny, ...? 
 
Cheers

Jose

[Gael Guedon]

Hi José,

> I am willing to see your memo.

Just finished the document! You can download it at the following link:

https://dl.dropbox.com/u/18801312/GuedonMassTransfer_2012_08_02.pdf

> I apologize by not answering you before but I have been on holidays the last two weeks of July.

> I have implemented a new source type (GfsSourceInterfaceFlux) what serves to impose a flux

> of a variable perpendicular to an interface defined by a TracerVOF.

No problem. I will try to have a look at your implementation.. too many things to do before going

in vacation..

> I have focused on its usage as a corrector of the VOF transport equation so I have compared

> to the gael implementation. In reality both treatment are equivalent since we both solve the

> same equation,
>
> F_t + U_I \cdot \nabla F = 0
>
> U_I is the velocity of the interface being U_I = v + v_{relative}. v is the fluid velocity and

> v_{relative} is the relative velocity of the interface. v clearly depends on the volume source

> imposed (or its absence see D. Lakehal et al. / Int. J. Heat and Fluid Flow 23 (2002) 242–257).

> What is not clear (unless for me) is what is really v (or what it should be). Clearly is an average

> velocity but I can not figure out a clear physical interpretation of it. I suspect that this unclear

> aspect contribute to the difficulties of interpreting the pressure values discussed in the thread:
>
> https://mail.google.com/mail/u/0/?shva=1#search/volume+source/1372ed2329af2eff

I am not sure how to interpret v in case of volume source too. However I do not think it is linked

with the pressure distribution. If there only a volume source, the values of pressure correspond

actually to the jumps in dynamic pressure and should have nothing to do with the volume source.

Effectively, in the ideal case where a constant volume source is applied, from one time step to

another, the pressure should correct only for the advection term which then corresponds to the

jump in dynamic pressure (i.e. the divergence of velocity at the previous time should cancel

out with the divergence of velocity at the new time). In practice, the projection is approximated,

thus there is a remaining imbalance in velocity divergence which generates oscillations near

the boundaries of the volume source where the imbalance is highest. The exact solution of

pressure should therefore be proportional to rho*u^2 which is verified.

> To check the performance of both implementation I have use an ellipse that grow at a constant

> rate being the fluid velocity zero, so the form has to be maintained.
>
> Some comments:
>
> The Gael implementation works fairly well but I see difficulties if TracerVOFConcentration are

> wanted to be used.
> My implementation has the strength that does not imply large changes of the existing code but

> in the present form (because it is naively implemented) it is prone to suffer  under/over shoots.

> The splitting character of the VOF advection scheme does not give any problem.

I am not sure about this but I have the feeling that the under/over shoots are due to a term of the

form:

F \nabla \cdot \mathbf{u}

that is remove during the Lagrangian advection of VOF. In my implementation, the smoothing

should compensate the variation of this term by varying the fluxes at faces, which in your

implementation are constant since the velocity jump is not smoothed over 2-3 cells, i.e. it does

not follow the position of the interface.

> I am dubious about what line I should follow, i.e. to abandon my line centering on the gael's

> line working into make it compatible with the TracerVOFConcentration or to keep working in

> mine trying to skip over/under shoots problems (I have in mind to shift the plane using the

> value of alpha to see if a neighbouring full cell has to be changed to be an interface cell).

> Certainly, the implementation of gael is more mature than mine. Any feedback is welcomed!

I think you could try to study the term I presented above if not done yet. However I am not sure

if I am interpreting well your comments about under/over shoots. In any case, as suggested

by Stephane, we should adopt an implementation based on the geometrical VOF.

> Attached you find:
>
> A patch with the GfsSourceInterfaceFlux
> A patch which contains the implementations of Gael
> simulations files of the ellipse test.

I will try to have a look.. it will be hard before september..

> Finally a question:
> It is about the available T_x, T_y and T_z there is any reason to be different of zero in full

> cells? This is why you, Gael, prefers define the yours T_nx, T_ny, ...? 

I think you could have some issues when interpolating to faces the normal to the interface,

however it could be handled quite easily by simply using the value at the interfacial cell.

I will keep you updated about the "possibly improved" implementation of the volume source.

Cheers,

Gael

[Gael Guedon]

Hi all,

my previous mail is a little bit messy so do not forget to expand trimmed content

to see all the comments if your e-mail application is hiding quotes.

Sorry

Gael

[Jose M. López-Herrera Sánchez]

Hi Gael,

Just finished the document! You can download it at the following link:

https://dl.dropbox.com/u/18801312/GuedonMassTransfer_2012_08_02.pdf

I have just downloaded. Good work! I suspect I will need some time to understand it.

I am not sure about this but I have the feeling that the under/over shoots are due to a term of the

form:

F \nabla \cdot \mathbf{u}

One of the reasons to use for testing an ellipse in a fluid at rest was to keep only the source additional test. So for the test case I am solving

F_t= \int_{Interface} function \delta \mathbf{n} ds

and the term you mention is null in this case. Observe that I do not modify the velocity field but the VOFtracer field. That is because I put the blame in my naive implementation not in any other cause. 

I think you could try to study the term I presented above if not done yet. However I am not sure

if I am interpreting well your comments about under/over shoots. In any case, as suggested

by Stephane, we should adopt an implementation based on the geometrical VOF.

I going to follow the  geometrical approach that Stephane suggested.


Cheers (and good holidays!)

José

PS, Your email was perfectly readable so do not worry!

[Gael Guedon]

Hi José,

> One of the reasons to use for testing an ellipse in a fluid at rest was to keep only the source additional test. So for the test case I am solving

>
> F_t= \int_{Interface} function \delta \mathbf{n} ds
>
> and the term you mention is null in this case. Observe that I do not modify the velocity field but the VOFtracer field. That is because I put the blame in my naive implementation not in any other cause.

sorry I misunderstood the test case. I just had a quick look at your patch and if I understood well, you are adding the flux due to mass transfer just after the flux update of the direction-split advection scheme.

Did you tried to add the flux due to mass transfer outside the advection scheme so that the normal to the interface and area of the interface is not updated meanwhile? I have no idea but maybe this is the update that generate over/under-shots...

Another thing that come to my mind is that the fluxes are not computed at face-centre so some imbalance from one cell to another could occur.

> I going to follow the  geometrical approach that Stephane suggested.

Maybe you could try to modify my implementation so that it is using the geometrical VOF. This should requires some checks on neighbors, i.e. if the neighbor is full cell, face-value of "mass transfer" is cell-value, otherwise interpolate (not sure if this has to be done this way though..).

> Cheers (and good holidays!)

Thank you!

Cheers,

Gael

[Jose M. López-Herrera Sánchez]

Hi Gael,

sorry I misunderstood the test case. I just had a quick look at your patch and if I understood well, you are adding the flux due to mass transfer just after the flux update of the direction-split advection scheme.

In effect, On each interface cell, I am adding/removing an additional flux of tracer beside to that already existing across the cell faces.
 

Did you tried to add the flux due to mass transfer outside the advection scheme so that the normal to the interface and area of the interface is not updated meanwhile? I have no idea but maybe this is the update that generate over/under-shots...

I think it not make any difference, If I understand well the Stephane's scheme, the geometrical stuff for each interface cell (vector m and alpha) is updated AFTER the tracer distribution is corrected. In other words, first tracer distribution at instant n+1 is calculated for all cells then vector m and alpha are updated to instant n+1 from the T^(n+1) distribution.
 

Another thing that come to my mind is that the fluxes are not computed at face-centre so some imbalance from one cell to another could occur.

 

 I do not think so.  The flux I am imposing is across the interface NOT across the cell faces.

Cheers

Jose

[Jose M. López-Herrera Sánchez]

Hi all,

I have been working a little bit on the movement of the interface due to mass transfer. I have followed two lines more or less in parallel: the one of gael which consist into add an additional face velocity with origin in the mass transfer and the mine in which such mass transfer is set as an additional flux in each interface. I have used to test the codes the following problems: 

• A inclined line which grows/shrink at a constant rate due to a mass transfer in a fluid at rest. 
• A corner   which grows/shrink at a constant rate in a fluid at rest. 
• An ellipse  which grows/shrink at a constant rate in a fluid at rest. 

All the stuff (patches, simulation files etc) is available at:

https://www.dropbox.com/sh/vopko5nrdd8lrex/KBcBM6u9QU

Gael approach:

Following the suggestion of Stephane I have modified the implementation of Gael to avoid smoothing using reconstructed VOF normal. I find very little difference between the Gael's smooth approach and the mine using the reconstructed VOF normal. I have left for a future the other suggestion of Stephane (to create a new Tracer class).

This approach gives very good results except close to the boundaries. I am pretty sure that this is due to the present implementation. The boundaries (and solids if present) have to be transparent to the additional velocities with origin in the mass transfer. I have tried to solve this by adding the flag  FTT_TRAVERSE_BOUNDARY_FACES to (masstransfer.c line  295)

  gfs_domain_face_traverse (domain, dimension == 2 ? FTT_XY : FTT_XYZ,

   FTT_PRE_ORDER, FTT_TRAVERSE_LEAFS, -1,

   (FttFaceTraverseFunc) correct_normal_velocity, &par)

but I have not succeed (I obtained a warning)  because I using it in a bad manner (I have not been able to find an example of using of this flag).  So any tip on its usage is welcome !

In general the results with this approach are excellent. Curiously if the interface is growing, the corners are a little bit smoothed.  This creates a much more better-looking reconstructed interface (see the ellipse case). On the contrary the error on the mass volume is a little bit higher than mine. I put the blame (maybe naively) on the average performed in faces shared by interfacial cells in this method.

This approach has, in my opinion, some weaknesses. First, I found difficult to treat well an eventual collapsing. Second,  if a TracerVOFConcentration is present I found difficult how to  evolve it properly (conservatively).

Jose approach:      

The implementation is far to be finished since it is prepared for uniform grid and it not takes into account the metrics (maybe for this I would need some help). This method has the advantage that cover the collapse and the extension to deal with a TracerVOFConcentration would be relatively easily.

As conclusion, I still doubt what approach would be better although I feel inclined to my approach. So If any of you can find some time for send me a feedback...

Cheers

Jose
 

[Gael Guedon]

Hi José,

thank you for improving my approach and compare it with yours. Sorry for the very late feedback, I have been busy last month and could not work on mass transfer in Gerris... so unfortunately I don't really have updates...

I gave a look at the results and it seems that your approach is more robust at boundaries and during collapse and nucleation. If you think it is also compatible with TracerVOFConcentration, maybe it is more suited than my approach. For the ellipse case, it seems to me that the error observed in your approach may come from the accuracy of the fluxes approximation, though I don't know yet how to improve it..

In the next days I should try to implement the velocity filter I mentioned last time and see if it works. Then I should be able to have a look at the interface transport.

I will keep you informed!

Cheers,

Gael

------------------------------------------------------------------------------
Live Security Virtual Conference
Exclusive live event will cover all the ways today's security and
threat landscape has changed and how IT managers can respond. Discussions
will include endpoint security, mobile security and the latest in malware
threats. http://www.accelacomm.com/jaw/sfrnl04242012/114/50122263/
_______________________________________________
Gfs-devel mailing list
Gfs-...@lists.sourceforge.net
https://lists.sourceforge.net/lists/listinfo/gfs-devel

[Jose M. López-Herrera Sánchez]

Hi Gael,

. I hope I will have some time to work a little bit on the mass transfer.

2012/10/11 Gael Guedon <gaelg...@gmail.com>

Hi José,

thank you for improving my approach and compare it with yours. Sorry for the very late feedback, I have been busy last month and could not work on mass transfer in Gerris... so unfortunately I don't really have updates...

Do not worry!, I am late answering you also. As you I have been a little busy this last month. Me neither have updates (yet).
 

I gave a look at the results and it seems that your approach is more robust at boundaries and during collapse and nucleation. If you think it is also compatible with TracerVOFConcentration, maybe it is more suited than my approach. For the ellipse case, it seems to me that the error observed in your approach may come from the accuracy of the fluxes approximation, though I don't know yet how to improve it..

Yes, I agree with you that maybe the fluxes can be better accounted for. Although,  I have not a clear idea on how to do it.
 

In the next days I should try to implement the velocity filter I mentioned last time and see if it works. Then I should be able to have a look at the interface transport.

Ok perfect!

Cheers
José

[Gaurav Tomar]

Hi

Perhaps time to revisit my suggestion on shifting the mass flux to the
neighbor full-cell cell-face even before the code solves the pressure
Poisson eqn. IMHO, it is not just about the correct motion of the
tracer of each kind but also the velocity field.

I would suggest one should evaluate this possibility of accurately
putting the volume source AND sink term (source in one and sink in the
other: not equal obv.).

Another issue:
In VOF, in an attempt to conserve volume we loose area accuracy.
Adaptive mesh refinement does come to rescue but for a uniform grid
one should be able to estimate the error in the computation of the
mass transfer.

[PS: Congrats to SP on creating an open forum through gerris where
people can discuss fluid mechanics science and code development more
openly than ever. ]

Cheers.

> ------------------------------------------------------------------------------
> Everyone hates slow websites. So do we.
> Make your web apps faster with AppDynamics
> Download AppDynamics Lite for free today:
> http://p.sf.net/sfu/appdyn_sfd2d_oct

> _______________________________________________
> Gfs-devel mailing list
> Gfs-...@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/gfs-devel
>

------------------------------------------------------------------------------
WINDOWS 8 is here.
Millions of people. Your app in 30 days.
Visit The Windows 8 Center at Sourceforge for all your go to resources.
http://windows8center.sourceforge.net/
join-generation-app-and-make-money-coding-fast/

[Jose M. López-Herrera Sánchez]

Hi Gaurav,

Another issue:
In VOF, in an attempt to conserve volume we loose area accuracy.
Adaptive mesh refinement does come to rescue but for a uniform grid
one should be able to estimate the error in the computation of the
mass transfer.

True, indeed, it is what I observe. About the estimation you suggest. How do you think it can be done?

 

Cheers

José

[Stephane Popinet]

Hi Gael,

> Just finished the document! You can download it at the following link:
>
> https://dl.dropbox.com/u/18801312/GuedonMassTransfer_2012_08_02.pdf

Looks like the link is dead. Could you please update it?

cheers

Stephane

------------------------------------------------------------------------------
Monitor your physical, virtual and cloud infrastructure from a single
web console. Get in-depth insight into apps, servers, databases, vmware,
SAP, cloud infrastructure, etc. Download 30-day Free Trial.
Pricing starts from $795 for 25 servers or applications!
http://p.sf.net/sfu/zoho_dev2dev_nov

[Gael Guedon]

Hi Stephane,

sorry I removed the file by mistake some time ago.. The link should work now

https://dl.dropbox.com/u/18801312/GuedonMassTransfer_2012_08_02.pdf

I would be very thankful if you could give any feedback on this.

Cheers,

Gael

[Gael Guedon]

Hi all,

I have some news on this issue :-)

I finished the preliminary implementation of the vertex divergence projection filter as well as the modified projection, as proposed in the document I provide you.

By the way, there were two small errors in the document so here is an updated version

https://dl.dropbox.com/u/18801312/GuedonMassTransfer_2012_11_16.pdf

I provide a patch for the implementation, which may require reviews because maybe some steps can be simplified/optimized..

Here is a brief description of what I did:

Vertex divergence filter (GfsEventVertexDivFilter):

- Implemented as an event

- Executed at each event_do

- Uses a regular 3x3x3 stencil

- Store gradient of velocity in order to compute values at virtual finer cells

- Filtering deactivated at corner boundaries

- If volume source present (as a GfsVariablePoissonDivergent, see below), remove the divergent part of the velocity before filtering

- Does not consider metrics (should not be important) and "alpha" (did not check influence yet)

Modified projection and volume source implementation (GfsVariablePoissonDivergent):

- Implemented using a new variable derived from GfsVariablePoisson

- GfsSimulation modified so that the Poisson equation for the potential is solved at the right moment

- GfsSimulation modified in order to correct the velocity field before the approximate projection

- The mac projection is not modified (uses gfs_divergent_hook to add the volume source)

- Need at least one boundary with Dirichlet condition for the potential otherwise does not converge

I made some tests using a sharp volume source (band of 1 cell thickness) and the results are identical to the results obtained from the 1D code in MATLAB (see document).

Which is very promising for MassTransfer simulation (at least I hope so..).

I still need to do some further tests on sharp sources to see if this is working properly and stable (2D, 3D and unsteady sources...)

Regarding 2D test, I re-used the already existing one (cylinder source in the test-case list). I changed the time-step to 0.0001 so that the instabilities observed by Daniel appear.

I computed analytically the theoretical pressure considering the volume source and added the related ErrorNorm and pressure distribution plots to the case.

I made another parameter file in which the filtering is used and put it in a sub-folder "./filtered" in a similar manner as for the Reynolds test-case.

Briefly, the pressure distribution in the filtered case is very clean, without oscillations.

Here is a series of commands you can use to run the test case

tar -zxvf source.tar.gz

cd source

sh source.sh source.gfs

cd filtered

sh ../source.sh filtered.gfs

I also provide a preview of the pressure distributions obtained with and without filtering so that you can see the improvements.

Any feedback on this is welcome! thanks

I will now have a look to the transport of the interface.

Cheers,

Gael

[Gael Guedon]

Sorry there was a mistake in the script to launch the test case... please find attached the correct version.

I also added two other tests where adaption is not uniform within the cylindrical volume source.

Cheers,

Gael

[Daniel Fuster]

Hi Gael

I have not looked into too much detail, but I found a different
analytical solution for the pressure (for the limiting case of Rout ->
infty). I was happy with it because it was the solution for the Pmac
field and also the solution I recover when I solve the Poisson
equation once more to get the correct pressure from the velocity
field. See

https://groups.google.com/forum/?fromgroups=#!topic/gerris-devel/2vHNuGtY2DM

Indeed, from your equation in the gfs file I find the same external
pressure evolution than you. However, the internal equation you have
is

P = P_{out} - { 3 s^2 r^2 \over 8 R_{out}^2 } + { s^2 R_c^4 \over 8
R_{out}^2 } ( 1 + { 2 R_{out}^2 \over R_c^2 } ) \; \; if \; \; r \le
R_c

which for Rout -> infty is constant (I guess it is a typo)

P = P_{out} + { s^2 R_c^2 \over 4 }

I find

P = P_{out} - { s^2 r^2 \over 8 }

Maybe I did something wrong (quite possible indeed) but it would be
nice to check because

best
Daniel

2012/11/22 Gael Guedon <gaelg...@gmail.com>:

> ------------------------------------------------------------------------------
> Monitor your physical, virtual and cloud infrastructure from a single
> web console. Get in-depth insight into apps, servers, databases, vmware,
> SAP, cloud infrastructure, etc. Download 30-day Free Trial.
> Pricing starts from $795 for 25 servers or applications!
> http://p.sf.net/sfu/zoho_dev2dev_nov
> _______________________________________________
> Gfs-devel mailing list
> Gfs-...@lists.sourceforge.net
> https://lists.sourceforge.net/lists/listinfo/gfs-devel
>

--
www.danielfuster.com

[Gael Guedon]

Hi Daniel,

thank you for the quick feedback!

which for Rout -> infty is constant (I guess it is a typo)

yes you are right, I made a mistake.. You should read in the description

P = P_{out} - { 3 s^2 r^2 \over 8 } + { s^2 R_c^4 \over 8 R_{out}^2 } ( 1 + { 2 R_{out}^2 \over R_c^2 } ) \; \; if \; \;  r \le R_c

and

return r > R0 ? C1*(1.-Rout*Rout/(r*r)) : 

              -3.*sour*sour*r*r/8.+C1*(1.+2.*Rout*Rout/(R0*R0));

in the analytic solution.

I am aware of your different solution and I think I have an explanation.

In Gerris, the following equation is solved

rho ( ddt (U) + div (UU) ) = - grad (P) + SourceU

which is equivalent in solving

ddt (rhoU) + div (rhoUU) = - grad (P) + SourceU

together with

ddt (rho) + div (rhoU) = SourceM

for the case of spatially and temporally constant density.

So for the analytic solution, I have the Bernoulli one outside the cylinder.

Inside the cylinder I have integrated the following equation

grad (P) = - U . grad (U) - U div (U) = - U . grad (U) - U s

which recovers what Gerris is solving for.

In your case, I guess you are integrating this equation inside the cylinder

(i.e. Bernoulli solution)

grad (P) = - U . grad (U)

that you recover with the final pressure projection present in your patch

in which the p->sources are omitted and the advection term is estimated

as U . grad (U).

The problem is that I am still not sure which solution represent the "real"

pressure... at least I have been able to recover what Gerris should solves

for...

My guess is that Gerris is solving for a simplification of the "complete" Navier-Stokes

equations even in presence of mass/volume source so the pressure obtained

from Gerris should be the "real" one.

I hope someone is able to make things clear.

Thank you for your help!

Best

Gael

[Jose M. López-Herrera Sánchez]

Hi Gael, 

First, let me congratulate you by the great work your are doing (and sharing with us). Your results look very promising!  I am a remorseful by leaving a little bit aside the mass transfer problem we were discussing meanwhile you were working hard on it;  My alibi is that until Christmas I have a heavy period of lectures in my university. I hope to be able to contribute more intensively again to the issue once the courses finish. 

Cheers

José

 

2012/11/22 Gael Guedon <gaelg...@gmail.com>

[Daniel Fuster]

Hi Gael

Yes, I fully agree with your explanation. The set of mathematical
equations I was expecting (and I want) to solve is

rho ( ddt (U) + U \cdot (grad U) ) = - grad (P) [1]
div (U) = SourceP [2]

when sourceP = 0, those are the equations written in the JCP Popinet 2009.

I think that, as you say, what it indeed coded for the cell centered values is
ddt (rhoU) + div (rhoUU) = - grad (P) [3]
div (U) = SourceP [4]

so that whenever I specify a Source P, we find that for rho=cte

div (rhoUU) = rho( U . grad (U) + U div (U)) = rho(U . grad (U) + U SourceP )

which implies that what gerris is solving is

rho ( ddt (U) + U \cdot (grad U) ) = - grad (P) - rho U SourceP
div (U) = SourceP

If we agree that sourceP is the term one should expect from the
solution 3-4, then I am fine with the current implementation (and your
correction to obtain a smooth pressure:). In the end, the solution of
the equations is more physical. In any case, for my problem I want to
solve for 1,2. Now that I know what happens, I can try to figure out a
solution. Thanks a lot

Best

Daniel

2012/11/22 Gael Guedon <gaelg...@gmail.com>:

[Gael Guedon]

Hi José and Daniel,

thank you for the support and discussions, this is always rewarding and helpful!

José I understand your alibi, I am quite busy with lectures as well (surely less than you).

Daniel, I agree that it depends on what equations are wanted to be solved.

I just would like to add a comment on the equations you wrote

I think that, as you say, what it indeed coded for the cell centered values is
ddt (rhoU) + div (rhoUU) = - grad (P)                  [3]
div (U) = SourceP                                            [4]

Actually, what is coded is

rho (ddt (U) + div (UU)) = - grad (P) + SourceU

div (U) = SourceP

because if you have a spatial or temporal variation of density, it is not taken into account.

The fact that you have div (UU) instead of U . grad (U) is due to the implementation of the Godunov scheme (see JCP Popinet 2003).

So as a result, this implies that Gerris is solving for the following equivalent equations (considering rho=cte)

ddt (rhoU) + div (rhoUU) = - grad (P) + SourceU

div (U) = SourceP

I just reversed your reasoning so that one should not confuse the fact that even if you assign a variable density, some terms will not be taken into account.

Cheers,

Gael 

[Stephane Popinet]

Hi Gael, Daniel,

> Actually, what is coded is
>
> rho (ddt (U) + div (UU)) = - grad (P) + SourceU
> div (U) = SourceP
>
> because if you have a spatial or temporal variation of density, it is not
> taken into account.
> The fact that you have div (UU) instead of U . grad (U) is due to the
> implementation of the Godunov scheme (see JCP Popinet 2003).
>
> So as a result, this implies that Gerris is solving for the following
> equivalent equations (considering rho=cte)
>
> ddt (rhoU) + div (rhoUU) = - grad (P) + SourceU
> div (U) = SourceP

I agree, although in the case where SourceP = 0, which is what is done
normally, what is actually solved can although be expressed as (for
rho non constant):

ddt (rhoU) + div (rhoUU) = - grad (P) + SourceU

ddt (rho) + U div (rho) = 0
div(U) = 0

that is

U ddt (rho) + UU div (rho) + rho ddt(U) + rho div(UU) = - grad (P) + SourceU
ddt (rho) + U div (rho) = 0
div(U) = 0

or again

rho (ddt(U) + div(UU)) = - grad (P) + SourceU
ddt (rho) + div (rho U) = 0
div(U) = 0

which is what is implemented in Gerris i.e. no terms are missing even
in the case of variable density.

But yes, if SourceP != 0, then one has to be careful.

cheers

Stephane

PS: Gael, shall I wait until this is fully clarified before reviewing
your patches?

[Gael Guedon]

Hi Stéphane,

thanks a lot for clarifying the equations solved by Gerris!

I think you can start reviewing my patch if you have time. The most delicate part is the modification of the approximate projection.

I used a new object for this (GfsVariablePoissonDivergent) and this may not be the most elegant implementation...

Let me know if some points are obscure.

Thank you for your precious time and help!!

Best,

Gael

[Gael Guedon]

Hi all,

I had some time to analyze the equations to solve in case of mass source or mass transfer, and I would like to have your feedback

For the moment in Gerris we have

rho [ddt (U) + U \cdot grad (U)] + u [ddt (rho) + rho div (U) + U \cdot grad (rho)] = - grad (P) + SourceU

ddt (rho) + U grad (rho) = 0

div (U) = SourceP

which is equivalent to

rho [ddt (U) + div (UU)] = - grad (P) + SourceU
ddt (rho) + div (rhoU) = rho SourceP
div (U) = SourceP

So when a SourceP is given, a mass source/sink is solved for. The patch I provided is suitable for such situation.

If a mass transfer is desired (eventual volume source/sink but no mass source/sink), then what should be solved is

rho [ddt (U) + U \cdot grad (U)] + u [ddt (rho) + rho div (U) + U \cdot grad (rho)] = - grad (P) + SourceU

ddt (rho) + U grad (rho) = - rho SourceP

div (U) = SourceP

which is equivalent to

rho [ddt (U) + U \cdot grad (U)] = - grad (P) + SourceU

ddt (rho) + div (rhoU) = 0

div (U) = SourceP

So depending on the desired problem, one should have the opportunity to correct the equations solved.

I think this can be done by modifying slightly the advection source term before the approximate projection, for example by removing an estimation of U div (U) at time n+1/2.

The patch I provided should also be useful in the case of mass transfer because it is dealing with another issue.

When I will have some time available, I will have a look at a possible correction of the advection term.

Cheers,

Gael