Equilibrium thru Interface

hfsf · November 19, 2012, 03:23

Hi Foamers!
I'm having a little bit of trouble trying to remake the work done at the work of Waheed et al (2002) (http://www.sciencedirect.com/science...17931002001242), where some axissymetric wedge section of a drop inside a cilinder is taken, this drop containing solute, and the outer phase is solute-free. So, essencially I want to computate the solute transfer.

The main problem is that i need to establish a equilibrium of the solute concentration inside and outside the drop, what I believe I could do by phisically stablishing a interface boundary in my domain and using groovyBC to ser the dragient expression. So, essencially, my domain would consist of 2 regions: the drop and the outer phase. By now, the mesh construction is OK, and I could stablish this phisical boundary and the paraFoam recognizes it.

But running my case in my own implementation of interfoam (interfoam + solute transport), i noticed that the interface bondary (so called defaultFaces) is not needed. In fact, I can remove it from the boundary listed in the 0 directory with no problem! That's not right, my solver just ignores my equilibrium condition. What should I do? My master depends on this series of simulation

I'm attaching my case files. Thanks in advance, I hope someone out there could help me.

PS: In order to save space, after download and extraction please run blockMesh.

akidess · November 19, 2012, 04:28

I think you need two domains with non-identical vertices, which you then couple using some sort of a baffle boundary condition. Have a look at the tutorials for chtMultiRegionFoam, where this is done for the temperature equation.

- Anton

Cyp · November 19, 2012, 06:51

Do you have a partitioning relationship between the concentration in both side of the interface ?

Something like :

$C_{g} = H C_{l} \texttt{ at } \mathcal{A}_{gl}$

If it is the case, you can add a transport equation in interFoam solver using the theory of distributions (see Haroun (2012) for exemple)). The concentration jump will be considered as an additional term in your transport equation.

Best,
Cyp

hfsf · November 19, 2012, 07:07

Dear Cyp and Anton,

Thanks for the fast reply. I will check your suggetions, thanks in advance.

And yes Cyp, i have a partitioning condition like the one you showed. But how could I implement this on the solver to consider this jump only thru the interface boundary? I will check the reference you suggested.

Cyp · November 19, 2012, 08:27

I can explain it to you through a simple example. Consider only the diffusion between two phases (beta and gamma for instance) :

In the beta-phase you have
$\nabla \cdot D_{\beta}\nabla C_{\beta} = 0,$

and in the gamma-phase
$\nabla \cdot D_{\gamma} \nabla C_{\gamma} = 0.$

Both phases are connected through a flux continuity at the interface
$\textbf{n}_{\beta\gamma} \cdot D_{\beta} \nabla C_{\beta} = \textbf{n}_{\beta\gamma} \cdot D_{\gamma} \nabla C_{\gamma} \texttt{ at } \mathcal{A}_{\beta \gamma},$

and the thermodynamic equilibrium condition reads:
$C_{\beta}=H C_{\gamma} \texttt{ at } \mathcal{A}_{\beta \gamma}$

What you look for is an partial differential equation that govern
$C = \alpha C_{\beta} + (1-\alpha) C_{\gamma}$
where $\alpha$ is the phase indicator provided from the VOF solution. With such a formulation, C is defined on the whole domain. In the same manner, you can defined a diffusion field as
$D = \alpha D_{\beta} + (1-\alpha) D_{\gamma}$

Now you express the derivative of C :
$\nabla C = \alpha \nabla C_{\beta} + (1-\alpha) \nabla C_{\gamma} + (C_{\beta} - C_{\gamma})\nabla \alpha$

multiplying this relation by D and applying the divergence operator, you get :

$\nabla \cdot D \nabla C = \alpha \nabla \cdot D_{\beta}\nabla C_{\beta} + (1-\alpha) \nabla \cdot D_{\gamma} \nabla C_{\gamma}$
$+ (D_{\beta}\nabla C_{\beta}-D_{\gamma} \nabla C_{\gamma})\nabla \alpha + \nabla \cdot D (C_{\beta} - C_{\gamma})\nabla \alpha$

Just keep in mind that according to the distribution theory you have : $\textbf{n}_{\beta\gamma} = -\nabla \alpha$ . Consequently, the previous equation reduces to:

$\nabla\cdot D \nabla C = \nabla \cdot D (C_{\beta} - C_{\gamma})\nabla \alpha,$

This additional term represents the interfacial jump condition. If there is a continuity, you can get rid of it. However, if you have a partitioning relation, you have to consider it.

At the interface, we have $C_{\beta}=H C_{\gamma}$ .

Consequently, $C_{\beta} - C_{\gamma}= (1-H)C_{\gamma}$

more over, $C = \alpha C_{\beta} + (1-\alpha) C_{\gamma} = (\alpha H + (1-\alpha))C_{\gamma}$

So $C_{\beta} - C_{\gamma}= \frac{(1-H)}{\alpha H + (1-\alpha)} C$

So your diffusion equation becomes :
$\nabla\cdot D \nabla C = \nabla \cdot D \frac{(1-H)}{\alpha H + (1-\alpha)} C \nabla \alpha,$

With such a formulation, you will automaticly have a jump condition at the interface between beta and gamma.

You can also optimised the solution with $D = \frac{D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$

I let you adapt this exemple to the advection-diffusion equation.

Best regards,
Cyp

voingiappone · June 28, 2013, 03:55

Quote:

Originally Posted by Cyp

So your diffusion equation becomes :
$\nabla\cdot D \nabla C = \nabla \cdot D \frac{(1-H)}{\alpha H + (1-\alpha)} C \nabla \alpha,$

With such a formulation, you will automaticly have a jump condition at the interface between beta and gamma.

You can also optimised the solution with $D = \frac{D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$

I let you adapt this exemple to the advection-diffusion equation.

Best regards,
Cyp

Cyp,

I don't actually know if you are still following this tread but I have a couple of questions.
It came for me too the time to implement the phase jump condition so I came back to this useful thread. In my previous case (a non-volatile tracer) I just had the Laplacian of a bunch of constants and alpha1: no problem in the solution if you insert the explicit laplacian (alpha1 already calculated). But in this case it is different because we have C an alpha1 simultaneously in the laplacian.

For what I can see in the other terms OpenFOAM always expects in the laplacian a dimensionedScalar and a volScalarField. So, I gathered all the constant terms in the fraction and calculated them before the C equation:

$H = D \frac{(1-H)}{\alpha H + (1-\alpha)}$

Being this a function of alpha1 I had to define it in the Createfields.H as another volScalarField. Again, no problem. What I actually cannot understand is: how do I formulate this in C++??

$\nabla\cdot D \nabla C = \nabla \cdot H \cdot C \nabla \alpha,$

should become

Code:

fvm::laplacian(DC, C)

But if I code it this way, my solver won't compile. For what I see, the compiled does not like the fact that C is in the Laplacian.
Moreover my low C++ knowledge prevents me for finding an alternative formulation.
Do you have any hint?

Thanks!

block · March 11, 2016, 10:53

Quote:

Originally Posted by voingiappone

Cyp,

I don't actually know if you are still following this tread but I have a couple of questions.
It came for me too the time to implement the phase jump condition so I came back to this useful thread. In my previous case (a non-volatile tracer) I just had the Laplacian of a bunch of constants and alpha1: no problem in the solution if you insert the explicit laplacian (alpha1 already calculated). But in this case it is different because we have C an alpha1 simultaneously in the laplacian.

For what I can see in the other terms OpenFOAM always expects in the laplacian a dimensionedScalar and a volScalarField. So, I gathered all the constant terms in the fraction and calculated them before the C equation:

$H = D \frac{(1-H)}{\alpha H + (1-\alpha)}$

Being this a function of alpha1 I had to define it in the Createfields.H as another volScalarField. Again, no problem. What I actually cannot understand is: how do I formulate this in C++??

$\nabla\cdot D \nabla C = \nabla \cdot H \cdot C \nabla \alpha,$

should become

Code:

fvm::laplacian(DC, C)

But if I code it this way, my solver won't compile. For what I see, the compiled does not like the fact that C is in the Laplacian.
Moreover my low C++ knowledge prevents me for finding an alternative formulation.
Do you have any hint?

Thanks!

Hi CyprienThank you very much for the explanation, can i please have a reference for these equation? I would really like to get in depth to understand the equations better.

Kind Regards

Rimsha

Cyp · March 11, 2016, 13:57

@Article{Haroun2010,
Title = {Volume of fluid method for interfacial reactive mass transfer: Application to stable liquid film },
Author = {Y. Haroun and D. Legendre and L. Raynal},
Journal = {Chemical Engineering Science },
Year = {2010},
Number = {10},
Pages = {2896 - 2909},
Volume = {65},
Abstract = {A volume of fluid method is developed in order to simulate reactive mass transfer in two-phase flows and is applied to study reactive laminar liquid film. The thermodynamic equilibrium of chemical species at the interface is considered using Henry's law. The chemical species concentration equation is solved using primitive variables and local fluxes are locally directly calculated at the interface. The present treatment of jump discontinuity of chemical concentration is consistent with a volume of fluid approach and the difficulty to calculate accurate local mass flux across interface is overcome. For plane interface, the precision of the numerical simulation is found to be very satisfactory while for curved interface a special procedure has been developed to reduce the development of spurious fluxes at the interface. The algorithm is validated for different cases by comparison with available solutions. The method is then applied to study non-reactive and reactive mass transfer in a falling liquid film. The results show that the liquid side mass transfer is well predicted by the Higbie (1935) theory when the transfer is controlled by the film advection provided that adequate parameters are considered, i.e. the actual velocity at interface and not the average liquid film velocity. For situations controlled by diffusion, the Sherwood number tends to a constant value characteristic of purely diffusive situations. For the reactive mass transfer, first and second order irreversible chemical reactions in the liquid phase are considered. The numerical results are compared respectively, with Danckwerts (1970) and Brian et al. (1961) solutions and good agreement is observed. The proposed Volume of Fluid method is shown to be well adapted to deal with interfacial reactive mass transfer problems. },
DOI = {http://dx.doi.org/10.1016/j.ces.2010.01.012},
File = {Haroun2010.pdf:ARTICLES/Haroun2010.pdf:PDF},
ISSN = {0009-2509},
Keywords = {\{CFD\}},
URL = {http://www.sciencedirect.com/science/article/pii/S0009250910000291}
}

Cheers,

block · March 11, 2016, 14:01

Quote:

Originally Posted by Cyp

@Article{Haroun2010,
Title = {Volume of fluid method for interfacial reactive mass transfer: Application to stable liquid film },
Author = {Y. Haroun and D. Legendre and L. Raynal},
Journal = {Chemical Engineering Science },
Year = {2010},
Number = {10},
Pages = {2896 - 2909},
Volume = {65},
Abstract = {A volume of fluid method is developed in order to simulate reactive mass transfer in two-phase flows and is applied to study reactive laminar liquid film. The thermodynamic equilibrium of chemical species at the interface is considered using Henry's law. The chemical species concentration equation is solved using primitive variables and local fluxes are locally directly calculated at the interface. The present treatment of jump discontinuity of chemical concentration is consistent with a volume of fluid approach and the difficulty to calculate accurate local mass flux across interface is overcome. For plane interface, the precision of the numerical simulation is found to be very satisfactory while for curved interface a special procedure has been developed to reduce the development of spurious fluxes at the interface. The algorithm is validated for different cases by comparison with available solutions. The method is then applied to study non-reactive and reactive mass transfer in a falling liquid film. The results show that the liquid side mass transfer is well predicted by the Higbie (1935) theory when the transfer is controlled by the film advection provided that adequate parameters are considered, i.e. the actual velocity at interface and not the average liquid film velocity. For situations controlled by diffusion, the Sherwood number tends to a constant value characteristic of purely diffusive situations. For the reactive mass transfer, first and second order irreversible chemical reactions in the liquid phase are considered. The numerical results are compared respectively, with Danckwerts (1970) and Brian et al. (1961) solutions and good agreement is observed. The proposed Volume of Fluid method is shown to be well adapted to deal with interfacial reactive mass transfer problems. },
DOI = {http://dx.doi.org/10.1016/j.ces.2010.01.012},
File = {Haroun2010.pdf:ARTICLES/Haroun2010.pdf:PDF},
ISSN = {0009-2509},
Keywords = {\{CFD\}},
URL = {http://www.sciencedirect.com/science/article/pii/S0009250910000291}
}

Cheers,

Brilliant, thank you very Cyprien, highly appreciated

Rimsha

Astrodan · May 24, 2016, 05:54

Note that although it was ages ago there is a mistake in the derivation done by Cyp:

Quote:

Originally Posted by Cyp

At the interface, we have $C_{\beta}=H C_{\gamma}$ .

Consequently, $C_{\beta} - C_{\gamma}= (1-H)C_{\gamma}$

is wrong, it should be $C_{\beta} - C_{\gamma}= (H-1)C_{\gamma}$

leading to the result of

$\nabla\cdot D \nabla C = \nabla \cdot D \frac{(H-1)}{\alpha H + (1-\alpha)} C \nabla \alpha,$

which, by the way, describes the view from "the other phase" compared to Haroun, who ends up at
$\nabla\cdot D \nabla C = \nabla \cdot D \frac{(1-H)}{\alpha + H (1-\alpha)} C \nabla \alpha,$
which would be equal to assuming $C_\gamma = H C_\beta$ (which in fact is just $H_1 = \frac{1}{H_2}$ ).
Note that Cyps answer was a mixup of both solutions of nominator and denominator.

Edit: Thanks anyway Cyp, this was generally a great help.
And while at it, has anyone here managed to reproduce the calculations in the paper with OpenFOAM? I kind of struggle with different Henry coefficients (some work, some produce an exception).

Astrodan · June 14, 2016, 11:14

Finally my solver works now the way it was supposed to be. For anyone who cares I will attach my modified solver as well as two testcases. All modifications were made by made, but based on the publication by Haroun et al (2010), cited above, as well as Nieves-Remacha et al (2015):

OpenFOAM Computational Fluid Dynamic Simulations of Two-Phase Flow and Mass Transfer in an Advanced-Flow Reactor

María José Nieves-Remacha, Lu Yang, and Klavs F. Jensen
Industrial & Engineering Chemistry Research 2015 54 (26), 6649-6659
DOI: 10.1021/acs.iecr.5b00480

The Attachements:

interHarounFoam:

This is the solver. it is a modified version of interFoam, implementing both Equations found in Nieves-Remacha (2015), Supplemental. This means it implements as well the Haroun approach, as an approach by Marshall et al. (2012). It is merely for testing purpose, the diffusion coefficients are fixed to $1e^{-5} m^2 s^{-1}$ in the createFields.H, and the fields used are named H (for Haroun) and M (for Marshall).

NievesRemacha2015:

This is a sample case for the first example in the supplemental by Nieves-Remacha et al. It contains all relevant data to run interHarounFoam and generate one set of the results.

Haroun2010:

A sample case to reproduce the example shown in Figure 1 a/b of Haroun et al. (2010). This case can be run as it is, too.

A couple of further notes:

Both sample cases contain a visual.pvsm, which can be loaded in paraview (either using File > Load State) or starting paraview with the parameter --state=visual.pvsm.
You will have to fix the path to your case, though.
Since diffusion coefficients and Herny coefficient are fixed in the solver, you will need to modify the createFields.H and re-compile the solver every time you want to try other values. Sorry about that. I am working on a library, which does the same thing, only better, which I intend to "publish" here soon.
In the fvSchemes make sure to use a Gauss linear scheme (or Gauss cubic, for what it matters) to solve div(phi,C). Most other schemes will fail, somewhere between extremely wrong results to crashing the solver.
Nieves-Remacha implement both equations as balances using the cellSurface (surfaceScalarField phiC), calculating first phiC ~ snGrad(alpha) and then div(phiC). Mathematically, the same could be achieved by having a laplacian(c, alpha) (something like that). Here OF allows to use a volScalarField c, which can also be created based on Harouns equation. However, this doesn't work either. I couldn't be bothered with the numerics or maybe even a simpler, mathematical reason for this, but you might care about this if you ever want to implement it on your own. Took me some days..

I hope this helps someone, I worked on this for quite a while.

May the foam be with you,
Timm

katete · May 11, 2017, 10:03

Hello,

I don't know if anyone is still following this conversation but I have a question concerning the harmonic averaging of the diffusivity as it is mentioned in Haroun et al. (2010):

As I know, the harmonic averaging is usually performed with a factor 2.
So does anyone know why in this case it is
$D = \frac{D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$

instead of
$D = \frac{2 D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$

I would be very glad if you could help me out with this.

Thanks and best regards!
Katharina

Astrodan · May 11, 2017, 10:17

I'm not sure how the harmonic mean is derived exactly, but if you look at wikipedia, there is a definition for the weighted harmonic mean:

$H = \left( \frac{\sum w_i x_i^{-1}}{w_i} \right)^{-1}$

which you can rewrite for the harmonic function given.

Apart from that, simple calculation would show you that a factor of two leads to invalid results, as $\alpha = 1$ or $\alpha = 0$ will lead the the corresponding diffusion coefficient twice as big as it should be in the pure phase.

katete · May 11, 2017, 10:40

Thank you so much for the quick reply, Astrodan. Especially your last sentence makes it very clear why there is no need for an additional factor.

katete · May 11, 2017, 10:55

Sorry, I have to ask another thing:

So, the equation

$D = \frac{D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$

is a substitute for the usual linear average of the diffusion:

$D = \alpha D_{\gamma} + (1-\alpha)D_{\beta}$

If we have $\alpha = 0$ for gas and $\alpha = 1$ for liquid, the usual linear averaging leads to
$D = D_{\gamma}$ for $\alpha = 1$ and $D = D_{\beta}$ for $\alpha = 0$ .

However, if I use the harmonic averaging, I get:
$D = D_{\gamma}$ for $\alpha = 0$ :
$D = \frac{D_{\beta} D_{\gamma}}{0* D_{\gamma} + 1*D_{\beta}}= D_{\gamma}$

and $D = D_{\beta}$ for $\alpha = 1$ .

Am I getting this wrong and making a mistake somewhere? Shouldn't the result be the same for the case of pure phases?

phani45 · September 9, 2017, 18:59

Hi Adrean,

I have implemented your equation and was able to get the concentration distribution for 1D diffusional problem which was matching exactly with Analytical solution. But when I try to calculate the fluxes. I'm getting very high values as a result my Mass transfer coefficient valuations are wrong.

can you please brief me, how did you calculate the fluxes across interface?

My code is as attached in the attachment.......also you can have a look at the C profile from the picture attached

Eagerly waiting for your reply

jawfi · September 16, 2017, 04:54

Quote:

Originally Posted by Astrodan

Finally my solver works now the way it was supposed to be. For anyone who cares I will attach my modified solver as well as two testcases. All modifications were made by made, but based on the publication by Haroun et al (2010), cited above, as well as Nieves-Remacha et al (2015):

OpenFOAM Computational Fluid Dynamic Simulations of Two-Phase Flow and Mass Transfer in an Advanced-Flow Reactor

María José Nieves-Remacha, Lu Yang, and Klavs F. Jensen
Industrial & Engineering Chemistry Research 2015 54 (26), 6649-6659
DOI: 10.1021/acs.iecr.5b00480

The Attachements:

interHarounFoam:

This is the solver. it is a modified version of interFoam, implementing both Equations found in Nieves-Remacha (2015), Supplemental. This means it implements as well the Haroun approach, as an approach by Marshall et al. (2012). It is merely for testing purpose, the diffusion coefficients are fixed to $1e^{-5} m^2 s^{-1}$ in the createFields.H, and the fields used are named H (for Haroun) and M (for Marshall).

NievesRemacha2015:

This is a sample case for the first example in the supplemental by Nieves-Remacha et al. It contains all relevant data to run interHarounFoam and generate one set of the results.

Haroun2010:

A sample case to reproduce the example shown in Figure 1 a/b of Haroun et al. (2010). This case can be run as it is, too.

A couple of further notes:

Both sample cases contain a visual.pvsm, which can be loaded in paraview (either using File > Load State) or starting paraview with the parameter --state=visual.pvsm.
You will have to fix the path to your case, though.
Since diffusion coefficients and Herny coefficient are fixed in the solver, you will need to modify the createFields.H and re-compile the solver every time you want to try other values. Sorry about that. I am working on a library, which does the same thing, only better, which I intend to "publish" here soon.
In the fvSchemes make sure to use a Gauss linear scheme (or Gauss cubic, for what it matters) to solve div(phi,C). Most other schemes will fail, somewhere between extremely wrong results to crashing the solver.
Nieves-Remacha implement both equations as balances using the cellSurface (surfaceScalarField phiC), calculating first phiC ~ snGrad(alpha) and then div(phiC). Mathematically, the same could be achieved by having a laplacian(c, alpha) (something like that). Here OF allows to use a volScalarField c, which can also be created based on Harouns equation. However, this doesn't work either. I couldn't be bothered with the numerics or maybe even a simpler, mathematical reason for this, but you might care about this if you ever want to implement it on your own. Took me some days..

I hope this helps someone, I worked on this for quite a while.

May the foam be with you,
Timm

Hi Timm,

Sorry but how can I open the modified solvers you attached ? I am interested to use these equations and modify them but I wasn't able to load cases with your modified solvers ? your help is appreciated !

Thanks,

bbita · October 20, 2017, 20:18

Dear Foamer,

I am going to use the same approach for the multiphase case. Any idea about the mass flux?

Thanks,

Astrodan · November 8, 2017, 11:26

Quote:

Originally Posted by jawfi

Sorry but how can I open the modified solvers you attached ? I am interested to use these equations and modify them but I wasn't able to load cases with your modified solvers ? your help is appreciated !

As far as I see there is the interHarounFoam.zip file attachted, which you will have to download and buid. The other are the sample cases in which the solver runs.

Quote:

Originally Posted by phani45

I have implemented your equation and was able to get the concentration distribution for 1D diffusional problem which was matching exactly with Analytical solution. But when I try to calculate the fluxes. I'm getting very high values as a result my Mass transfer coefficient valuations are wrong.

can you please brief me, how did you calculate the fluxes across interface?

My code is as attached in the attachment.......also you can have a look at the C profile from the picture attached

I don't fully understand your code, since it kind of misses context/comments. However, Haroun simply used Ficks Law (J = - D dC/dX) to calculate the flux, while I have no clue so far what your equations are about. Try that first, and then see if the results are reasonable. phiCf, by the way, is only a corrector to achieve a jump condition at the interface, and no actual occuring flux.

I think that for some meshes a high flux might occur simply due to numerical interpolation errors at the mesh points in immediate vicinity of the interface. Check that, too.

Quote:

Originally Posted by bbita

I am going to use the same approach for the multiphase case. Any idea about the mass flux?

Well, tell us what you got so far, what fails, and where your problem is exactly. No one can answer based on a question like that. Apart from that, check the source papers and above posts..

bbita · January 21, 2018, 13:03

Hi Foamers,

I have the transport equation for concentration which is validated with the analytical solution. Now I am trying to add phase change through Sp and Su. I used CEqn.flux() and add that term as Su and Sp to alphaEqn.H but it didn't work. Any idea how I should solve this problem with phase change.

Thanks,

Bita

November 19, 2012, 04:28		#2
akidess Senior Member Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 30	I think you need two domains with non-identical vertices, which you then couple using some sort of a baffle boundary condition. Have a look at the tutorials for chtMultiRegionFoam, where this is done for the temperature equation. - Anton hfsf likes this. __________________ On twitter @akidTwit Spend as much time formulating your questions as you expect people to spend on their answer.

November 19, 2012, 06:51		#3
Cyp Senior Member Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 299 Rep Power: 18	Do you have a partitioning relationship between the concentration in both side of the interface ? Something like : $C_{g} = H C_{l} \texttt{ at } \mathcal{A}_{gl}$ If it is the case, you can add a transport equation in interFoam solver using the theory of distributions (see Haroun (2012) for exemple)). The concentration jump will be considered as an additional term in your transport equation. Best, Cyp hfsf likes this. Last edited by Cyp; November 19, 2012 at 09:27.

November 19, 2012, 08:27		#5
Cyp Senior Member Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 299 Rep Power: 18	I can explain it to you through a simple example. Consider only the diffusion between two phases (beta and gamma for instance) : In the beta-phase you have $\nabla \cdot D_{\beta}\nabla C_{\beta} = 0,$ and in the gamma-phase $\nabla \cdot D_{\gamma} \nabla C_{\gamma} = 0.$ Both phases are connected through a flux continuity at the interface $\textbf{n}_{\beta\gamma} \cdot D_{\beta} \nabla C_{\beta} = \textbf{n}_{\beta\gamma} \cdot D_{\gamma} \nabla C_{\gamma} \texttt{ at } \mathcal{A}_{\beta \gamma},$ and the thermodynamic equilibrium condition reads: $C_{\beta}=H C_{\gamma} \texttt{ at } \mathcal{A}_{\beta \gamma}$ What you look for is an partial differential equation that govern $C = \alpha C_{\beta} + (1-\alpha) C_{\gamma}$ where $\alpha$ is the phase indicator provided from the VOF solution. With such a formulation, C is defined on the whole domain. In the same manner, you can defined a diffusion field as $D = \alpha D_{\beta} + (1-\alpha) D_{\gamma}$ Now you express the derivative of C : $\nabla C = \alpha \nabla C_{\beta} + (1-\alpha) \nabla C_{\gamma} + (C_{\beta} - C_{\gamma})\nabla \alpha$ multiplying this relation by D and applying the divergence operator, you get : $\nabla \cdot D \nabla C = \alpha \nabla \cdot D_{\beta}\nabla C_{\beta} + (1-\alpha) \nabla \cdot D_{\gamma} \nabla C_{\gamma}$ $+ (D_{\beta}\nabla C_{\beta}-D_{\gamma} \nabla C_{\gamma})\nabla \alpha + \nabla \cdot D (C_{\beta} - C_{\gamma})\nabla \alpha$ Just keep in mind that according to the distribution theory you have : $\textbf{n}_{\beta\gamma} = -\nabla \alpha$ . Consequently, the previous equation reduces to: $\nabla\cdot D \nabla C = \nabla \cdot D (C_{\beta} - C_{\gamma})\nabla \alpha,$ This additional term represents the interfacial jump condition. If there is a continuity, you can get rid of it. However, if you have a partitioning relation, you have to consider it. At the interface, we have $C_{\beta}=H C_{\gamma}$ . Consequently, $C_{\beta} - C_{\gamma}= (1-H)C_{\gamma}$ more over, $C = \alpha C_{\beta} + (1-\alpha) C_{\gamma} = (\alpha H + (1-\alpha))C_{\gamma}$ So $C_{\beta} - C_{\gamma}= \frac{(1-H)}{\alpha H + (1-\alpha)} C$ So your diffusion equation becomes : $\nabla\cdot D \nabla C = \nabla \cdot D \frac{(1-H)}{\alpha H + (1-\alpha)} C \nabla \alpha,$ With such a formulation, you will automaticly have a jump condition at the interface between beta and gamma. You can also optimised the solution with $D = \frac{D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$ I let you adapt this exemple to the advection-diffusion equation. Best regards, Cyp nimasam, voingiappone, hfsf and 4 others like this. Last edited by Cyp; November 19, 2012 at 09:26.

March 11, 2016, 13:57		#8
Cyp Senior Member Cyprien Join Date: Feb 2010 Location: Stanford University Posts: 299 Rep Power: 18	@Article{Haroun2010, Title = {Volume of fluid method for interfacial reactive mass transfer: Application to stable liquid film }, Author = {Y. Haroun and D. Legendre and L. Raynal}, Journal = {Chemical Engineering Science }, Year = {2010}, Number = {10}, Pages = {2896 - 2909}, Volume = {65}, Abstract = {A volume of fluid method is developed in order to simulate reactive mass transfer in two-phase flows and is applied to study reactive laminar liquid film. The thermodynamic equilibrium of chemical species at the interface is considered using Henry's law. The chemical species concentration equation is solved using primitive variables and local fluxes are locally directly calculated at the interface. The present treatment of jump discontinuity of chemical concentration is consistent with a volume of fluid approach and the difficulty to calculate accurate local mass flux across interface is overcome. For plane interface, the precision of the numerical simulation is found to be very satisfactory while for curved interface a special procedure has been developed to reduce the development of spurious fluxes at the interface. The algorithm is validated for different cases by comparison with available solutions. The method is then applied to study non-reactive and reactive mass transfer in a falling liquid film. The results show that the liquid side mass transfer is well predicted by the Higbie (1935) theory when the transfer is controlled by the film advection provided that adequate parameters are considered, i.e. the actual velocity at interface and not the average liquid film velocity. For situations controlled by diffusion, the Sherwood number tends to a constant value characteristic of purely diffusive situations. For the reactive mass transfer, first and second order irreversible chemical reactions in the liquid phase are considered. The numerical results are compared respectively, with Danckwerts (1970) and Brian et al. (1961) solutions and good agreement is observed. The proposed Volume of Fluid method is shown to be well adapted to deal with interfacial reactive mass transfer problems. }, DOI = {http://dx.doi.org/10.1016/j.ces.2010.01.012}, File = {Haroun2010.pdf:ARTICLES/Haroun2010.pdf:PDF}, ISSN = {0009-2509}, Keywords = {\{CFD\}}, URL = {http://www.sciencedirect.com/science/article/pii/S0009250910000291} } Cheers, the_ichthyologist likes this.

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
Wind turbine simulation	Saturn	CFX	60	July 17, 2024 06:45
An error has occurred in cfx5solve:	volo87	CFX	5	June 14, 2013 18:44
RPM in Wind Turbine	Pankaj	CFX	9	November 23, 2009 05:05
Convective Heat Transfer - Heat Exchanger	Mark	CFX	6	November 15, 2004 16:55
Replace periodic by inlet-outlet pair	lego	CFX	3	November 5, 2002 21:09

November 19, 2012, 07:07		#4
hfsf Member Hanniel Freitas Join Date: Jan 2012 Location: Natal, Brazil Posts: 34 Rep Power: 14	Dear Cyp and Anton, Thanks for the fast reply. I will check your suggetions, thanks in advance. And yes Cyp, i have a partitioning condition like the one you showed. But how could I implement this on the solver to consider this jump only thru the interface boundary? I will check the reference you suggested.

May 11, 2017, 10:03		#12
katete New Member Join Date: Jan 2016 Posts: 10 Rep Power: 10	Hello, I don't know if anyone is still following this conversation but I have a question concerning the harmonic averaging of the diffusivity as it is mentioned in Haroun et al. (2010): As I know, the harmonic averaging is usually performed with a factor 2. So does anyone know why in this case it is $D = \frac{D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$ instead of $D = \frac{2 D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$ I would be very glad if you could help me out with this. Thanks and best regards! Katharina

May 11, 2017, 10:17		#13
Astrodan Member Timm Severin Join Date: Mar 2014 Location: Munich Posts: 63 Rep Power: 12	I'm not sure how the harmonic mean is derived exactly, but if you look at wikipedia, there is a definition for the weighted harmonic mean: $H = \left( \frac{\sum w_i x_i^{-1}}{w_i} \right)^{-1}$ which you can rewrite for the harmonic function given. Apart from that, simple calculation would show you that a factor of two leads to invalid results, as $\alpha = 1$ or $\alpha = 0$ will lead the the corresponding diffusion coefficient twice as big as it should be in the pure phase.

May 11, 2017, 10:40		#14
katete New Member Join Date: Jan 2016 Posts: 10 Rep Power: 10	Thank you so much for the quick reply, Astrodan. Especially your last sentence makes it very clear why there is no need for an additional factor.

May 11, 2017, 10:55		#15
katete New Member Join Date: Jan 2016 Posts: 10 Rep Power: 10	Sorry, I have to ask another thing: So, the equation $D = \frac{D_{\beta} D_{\gamma}}{\alpha D_{\gamma} + (1-\alpha)D_{\beta}}$ is a substitute for the usual linear average of the diffusion: $D = \alpha D_{\gamma} + (1-\alpha)D_{\beta}$ If we have $\alpha = 0$ for gas and $\alpha = 1$ for liquid, the usual linear averaging leads to $D = D_{\gamma}$ for $\alpha = 1$ and $D = D_{\beta}$ for $\alpha = 0$ . However, if I use the harmonic averaging, I get: $D = D_{\gamma}$ for $\alpha = 0$ : $D = \frac{D_{\beta} D_{\gamma}}{0* D_{\gamma} + 1*D_{\beta}}= D_{\gamma}$ and $D = D_{\beta}$ for $\alpha = 1$ . Am I getting this wrong and making a mistake somewhere? Shouldn't the result be the same for the case of pure phases?

October 20, 2017, 20:18		#18
bbita Member Join Date: May 2017 Posts: 44 Rep Power: 9	Dear Foamer, I am going to use the same approach for the multiphase case. Any idea about the mass flux? Thanks,

January 21, 2018, 13:03		#20
bbita Member Join Date: May 2017 Posts: 44 Rep Power: 9	Hi Foamers, I have the transport equation for concentration which is validated with the analytical solution. Now I am trying to add phase change through Sp and Su. I used CEqn.flux() and add that term as Su and Sp to alphaEqn.H but it didn't work. Any idea how I should solve this problem with phase change. Thanks, Bita