CFD Online Discussion Forums

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   InterFoam: Add an equation to only one phase (https://www.cfd-online.com/Forums/openfoam-programming-development/110512-interfoam-add-equation-only-one-phase.html)

voingiappone December 13, 2012 01:08

InterFoam: Add an equation to only one phase
 
1 Attachment(s)
Hello again everybody,

as I already explained here, I want to insert a tracer in my simulation to calculate the mixing time of the reactor.
I want to consider only the liquid phase as I am going to insert a certain amount of NaCl as tracer and follow the time evolution of the concentration with a probe.

Now I was just able to add the tracer using as a base the wiki page on how to add temperature to icoFoam. However, if I run the simulation I get my tracer propagating out of the liquid phase through the interface into the gas.

Watching this thread (http://www.cfd-online.com/Forums/ope...interfoam.html) I made an idea of how to read physical properties for different phases but I cannot figure out how to make my equation to be solved only for the liquid phase.
To put things simple, now it looks like it is just overlapping the liquid-gas system (i have a global equation, a global diffusion coefficient, etc). See attached screenshot (tracer concentration profile on the right).

Is there anybody who can help me with this task?

Thanks!

Bernhard December 13, 2012 03:04

Did you try to determine the diffusion coefficient based on the liquid fraction alpha1? There is no way to solve an equation for only one of the phases using interFoam.

voingiappone December 13, 2012 05:13

Quote:

Originally Posted by Bernhard (Post 397285)
Did you try to determine the diffusion coefficient based on the liquid fraction alpha1? There is no way to solve an equation for only one of the phases using interFoam.

Bernhard,

thank you for the precious hint! I thought that the equations were solved for both the phases... but then, where the alpha1 parameter meaning goes??
So, you are right and I should take that into account, which brings us (me) to a big problem... I don't have any idea on how to modify the C++ code.

I tried to add alpha1 in the ::div and ::ddt functions but they seem to be programmed to accept only 2 arguments, hence I need to define a new variable (or a function) which already incorporates the calculations.

It is such a simple task (a matter of multiplying, summing and subtracting) but I don't know how to implement it in the code. I also tried copy/pasting other functions (also from the thread mentioned above) but to no avail. I'm stuck (again).

Thanks!

Bernhard December 13, 2012 05:42

For the implementation you have linked some threads for the temperature implementation, which will be similar. However, you will probably not be able to keep all the concentration in the liquid phase only by numerical diffusion. This is a complicated issue, which is probably described in literature a lot.

nimasam December 13, 2012 05:49

bernhard is right, there is no way :), to solve your equation just for one phase
but to reduce numerical diffusion, in each time step you should update your concentration with alpha1

voingiappone December 13, 2012 07:02

@Bernhard

I was already aware of that.... or at least I was imagining that such a behaviour could cause problems.

@nimasam

I would like to do that but I lack the ability to create functions, define new variables or manipulate existing ones.
To put it simple, I don't know C++...:confused: Copy/pasting a tutorial and adapt it to the new solver is on thing; to write your fuctions is another.


Regarding the diffusion problems... well, my aim in this study is to model a phtobioreactor to assess whether a swirling flow arises and if it does, to calculate the frequency. This is not something like, you create your model, you implement it and then validate it through experiments. It is the reverse: I would like to explore different reactor configurations and when I find the good one, I will build it and see if the effect is the one that I expect. For this task the standard interFoam suffices and performs perfectly (I am already writing a paper with the validated results). The mixing time issue is something I would like to add in the workflow as I suspect that it can be used as a significant parameter... not sure though.
However, I expect a Tmix of more than 1 minute so, the small effect of the diffusion between phases can become relevant. Just need to test it...

If somebody knows how I can multiply my concentration for alpha1 before executing the divs and ddts I would be very grateful. for those who know C++ should be a trivial task...

Tnks!

Bernhard December 13, 2012 07:33

Quote:

Originally Posted by voingiappone (Post 397318)
If somebody knows how I can multiply my concentration for alpha1 before executing the divs and ddts I would be very grateful. for those who know C++ should be a trivial task...

It is easier for people to reply on this request if you already show a code snippet of what you are trying now, and explain why it doesn't work (show compiling errors).

voingiappone December 13, 2012 08:00

Quote:

Originally Posted by Bernhard (Post 397324)
It is easier for people to reply on this request if you already show a code snippet of what you are trying now, and explain why it doesn't work (show compiling errors).

I suppose you're right but I actually don't have many things to show.

The equation I would like to implement is really basic: if I neglect the diffusion (as the coefficient is 10^-9 order):

d(rho*c)+div(rho*c*u)=0

which in OF should translate to:

fvm::ddt(rho, c)
+ fvm::div(rho, c, U)

If I get it right. Now I guess that if you don't take into account alpha1, interFoam just solves the eq for the continuum....
The first problem here is that "div(rho, c, U)" is not recognized and the compialtion fails (I don't have the trace as I am already at home). If I remove c or U (ie I leave only two terms) the file compiles.

So, here I would like to add the alpha1 dependence in these derivatives but, as they only accept a limited number of parameters I suppose I should assign the calculations to another variable and than pass it to the main equation. For example:

CAlpha == c*alpha1
CAlphaRho == CAlpha*rho

and then:

fvm::ddt(rho, CAlpha)
+ fvm::div(CAlphaRho, U)

This should work (should get c=~0 for the gas and solve normally when alpha1=~1).

Of course this code does not compile, complaining that neither CAlpha nor CAlphaRho have been declared for this scope. As I don't know how to define variables and subsequently how to use them in calculations, I don't know how and where to insert which code.

Sorry for even not being able to better explain the situation...

Bernhard December 13, 2012 08:08

If you want to have convection of terms, you need the velocity at the phases, so you could use div(phi,c) or div(phiRho,c) (see http://openfoamwiki.net/index.php/Ma...ver_is_writing )

Some strange things are happening if you are defining CAlpha == c*alpha1. You should try to get your solver correct without these steps, and thus solving for c like you solve for T. As you only have convection, physically, as there is no mass flowing across the interface, you will also have no transport of c. Then you could try to multiply c after the solution step by alpha1 ( c *= alpha1 ), but I don't know how this will work out, and at least some concentration will be lost.

voingiappone December 13, 2012 08:32

Quote:

Originally Posted by Bernhard (Post 397330)
If you want to have convection of terms, you need the velocity at the phases, so you could use div(phi,c) or div(phiRho,c) (see http://openfoamwiki.net/index.php/Ma...ver_is_writing )

This makes sense... it is the same in the Temperature equation. Sorry for pointing to a bad formula.

Quote:

Some strange things are happening if you are defining CAlpha == c*alpha1. You should try to get your solver correct without these steps, and thus solving for c like you solve for T. As you only have convection, physically, as there is no mass flowing across the interface, you will also have no transport of c. Then you could try to multiply c after the solution step by alpha1 ( c *= alpha1 ), but I don't know how this will work out, and at least some concentration will be lost.
So, if I got the point right I should solve the equation "as is" (using phi) and then multiply the obtained c for alpha1.... :

Code:


fvScalarMatrix CEqn
(
fvm::ddt(rho, c)
+ fvm::div(phi, c)
);         

CEqn.solve();

c *= alpha1;

The multiplication step will fit just like that? Need to add something to the code to define it? alpha1 and c are existing variables so I suppose I don't need anything more than what I have.

Loosing some of the concentration in the multiplication step may be a minor drawback (I can look for a constant asymptotic value for the concentration and then normalize to the max). If the solver "eats" some of the concentration for each step, I will never get an asymptotic value though....

I will test this solution tomorrow. I hope this fix will work!

Bernhard, thank you so much for your time!

michielm December 13, 2012 12:47

A possible solution, though it is a bit of a dodgy method, would be to add in the lost concentration in each timestep.

You could do this by checking in each timestep what the difference between the initial concentration and the current concentration is. The difference is the amount lost to numerical diffusion and you could add this back into your domain.

For example by elevating all cell concentrations by an equal amount, but this might give you a small artificial increase of mixing speed (this is actually the same as your normalization). The other, probably more physical, option is to add the concentration at the interface (so in cells where 0<alpha1<1) because that is the location where the concentration is lost as well.

Again: this is not at all a physical way of dealing with it, but since the numerical diffusion is unphysical as well you might get away with a solution of this kind.

voingiappone December 14, 2012 02:39

1 Attachment(s)
Quote:

Originally Posted by Bernhard (Post 397330)
Then you could try to multiply c after the solution step by alpha1 ( c *= alpha1 ), but I don't know how this will work out, and at least some concentration will be lost.

Thank you very much for the formula! The solver compiles flawlessly and I was able to run a test simulation with no problems.
Regarding the lost concentration, I am afraid that this may be a major problem. I was quite optimistic but looking at the results, I changed my mind.
I saw 10 ml of 1 g/l solution dilute to 1 mg/l when the tracer did't occupy more than one third of the total volume (1 liter). This means it should dilute up to 1/30th, more or less but it almost vanished (1 mg/l).
This comes from the fact that the bubbles moving through the tracer plume tend to "eat" it like a worm while the liquid should open and close leaving the bubble pass over instead. You can see three bubbles I evidenced in a screenshot I attach.
This brings us to what michielm said:

Quote:

A possible solution, though it is a bit of a dodgy method, would be to add in the lost concentration in each timestep
This is probably what I need to do to reduce the effect so I started analyzing the problem. A little bit deeper.
Your first suggestion to blindly add back the lost concentration does not fits well my way of seeing the problem: if I add back the deltaC in each time step, the actual concentration drop due to the mixing gets lost in the process and I will end up with a constant profile. Or maybe I didn't get your point right. However this is a really sexy statement:

Quote:

The other, probably more physical, option is to add the concentration at the interface (so in cells where 0<alpha1<1) because that is the location where the concentration is lost as well.
I totally agree with this: being able to do that could lead to a (somewhat) acceptable solution. How do we do this???

I spent some time in thinking about this and I came to an idea which I will try to explain to see what you think about that. Hope not to sound naive ;)

First of all we have the same equation I compiled (with the alpha1 inclusion) solved to obtain the concentration and then we multiply for alpha1:
Code:

solve(CEqn)  ----> C'(t)  ----> C'(t)*alpha1 = C(t)
If everything works as expected, in the first step, no problem occurs during the solution (alpha1 is still out of play). So in theory if I make a SUM of all the tracer quantities in the domain, it should give me the same amount of tracer I injected in the first time step. This for each time step, due to the mass conservation + lack of outlets in the reactor:


The tracer amount that vanishes in the time step (compared with the overall mass balance) should be then added back in the domain as michielm noted. At this point, adding it blindly in the domain (in every cell the same amount) will only screw things up... so I thought: what if I was able to sum the tracer back only in those cells where alpha1>0.5? And, wouldn't be nice to be able not to add anything back where alpha1=1 (i.e. do not artificially increase the max concentration)? So I schematized the interface situation like this:


That's the profile of the interface and thus the profile of the concentration C after the multiplication. The part that's removed by alpha1<0.5 is the one that vanishes. So, I was wondering if it was possible to add it back in the blue dashed zone. Again, I don't want to add it back blindly, otherwise my C value in the alpha1=1 zone would get bigger than before the correction! So I am planning to use a sort of weighting system based on alpha1, using something like a gaussian profile like the one shown here:


First you check whether the alpha1 value is grater than 0.5 and then you do your final correction basing on the profile of that curve:


Of course you loose more then what you put back in the domain but at least you add it back to the right place.

This is more or less what I was thinking to do... what do you think about it? And of course, how do I implement it without screwing up the code?

Thank you and sorry for the long post!

Bernhard December 14, 2012 02:55

This is certainly an approach to keep the total mass of concentration at the end of each time step constant. However, there is no physical argument at all why you should do this. Your actually correcting the concentration with alpha purely. By the way, is the mass balanced without this correction? Now you want to partially undo this with a complicated approach. You can also think of multiplying with a function looking like Heaviside(alpha1-0.5), which will be doing something better already. Maybe this is what you are planning to do already? Also, you can never assume alpha1 across the interface anti-symmetric, so it will never be an optimal solution.

But please, do not try to solve this problem independently. I am quite convinced that people have been looking for solutions to this issue already quite a few times. Search in literature to find out what has been done, and compare the results of different approaches. This problem will take you a long time to solve, and then to only find out that it has been done in the 1990s will be frustrating:

Have a look into how the interface is kept sharp in interFoam. Some artificial interface compression scheme is often used. Maybe it is possible to use this artificial convection in your equations as well, as it will keep the concentration away from the interface. Good luck, it is going to be tough :)

voingiappone December 14, 2012 03:50

Quote:

Originally Posted by Bernhard (Post 397476)
This is certainly an approach to keep the total mass of concentration at the end of each time step constant. However, there is no physical argument at all why you should do this. Your actually correcting the concentration with alpha purely. By the way, is the mass balanced without this correction?

You're right, no physical meaning in all this. I came to this conclusion just following what michielm suggested which is "non physical".
By the way, if I do directly multiplying by alpha1, the result is not mass balanced. And I dare to say it is badly unbalanced.

Quote:

Now you want to partially undo this with a complicated approach. You can also think of multiplying with a function looking like Heaviside(alpha1-0.5), which will be doing something better already. Maybe this is what you are planning to do already? Also, you can never assume alpha1 across the interface anti-symmetric, so it will never be an optimal solution.
Actually I was planning to define the gaussian function (which I generically indicated with FUNC in the previous post) to look like the heaviside but without the step shape by defining an half gaussian for the + side and the reverse for the - side (to give that fringe look to the interface).


Quote:

This problem will take you a long time to solve, and then to only find out that it has been done in the 1990s will be frustrating:
Have a look into how the interface is kept sharp in interFoam. Some artificial interface compression scheme is often used. Maybe it is possible to use this artificial convection in your equations as well, as it will keep the concentration away from the interface. Good luck, it is going to be tough :)
Will definitely need to go through some literature and the suggestion of the interface compression can be useful too. I will see if I can get to some point.
Meanwhile thanks for your helpful support!

olivierG December 14, 2012 07:22

hello,

You may take a look at the two Phase Euler solver ... if your case is applicable for Euler-Euler, then it will be easier to solve your scalar trough one phase only.
If you stay with VOF, you may try to add 2 scalar instead of one, and balance from a <-> b depending of alpha: you will be conservative for (a,b), then it should be easier to rebalance your scalar at interface, and stay conservative.

regards,
olivier

voingiappone December 15, 2012 02:01

Quote:

Originally Posted by olivierG (Post 397527)
You may take a look at the two Phase Euler solver ...

Hello Olivier, thanks for your suggestion.
Unfortunately I have tried the Euler-Euler solver many times and no matter what I do and how I initialize the domain, I never end up with a realistic fluid patter. Not to mention the gas hold-up. So I was forced to leave it for the more complicated (and time consuming) VOF. At least this one reproduces the behavior of my bubble column / air-lift system flawlessly.

Quote:

If you stay with VOF, you may try to add 2 scalar instead of one, and balance from a <-> b depending of alpha: you will be conservative for (a,b), then it should be easier to rebalance your scalar at interface, and stay conservative.
Sounds intriguing but I cannot figure out the exact mechanism.... I define two parameters (let's say C1 and C2) then I solve for both the same equation but in the first I multiply for alpha and in the latter I multiply for (1-alpha).
Then? How do I balance from a to b using alpha?

Thanks a lot.

Luca

Cyp December 15, 2012 04:51

hello !

What about the thermodynamic equilibrium at the interface ? According to your figure you have a partitioning at the gaz/liquid interface.

Have you read this topic http://www.cfd-online.com/Forums/ope...interface.html ?

You can also adapt the example to an equilibrium concentration at the interface and solve only the diffusion in one phase. (it is the single way to uncouple the equations in both phases)

Regards,
Cyp

voingiappone December 16, 2012 02:54

Quote:

Originally Posted by Cyp (Post 397654)
hello !

What about the thermodynamic equilibrium at the interface ? According to your figure you have a partitioning at the gaz/liquid interface.

Have you read this topic http://www.cfd-online.com/Forums/ope...interface.html ?

You can also adapt the example to an equilibrium concentration at the interface and solve only the diffusion in one phase. (it is the single way to uncouple the equations in both phases)

Regards,
Cyp

Hello Cyp!

Thanks a lot for the post you linked! That looks like the perfect solution! Defining C on the whole domain and then apply the derivatives was actually something I thought about but definitely not in the way you wrote there. That actually makes sense and moreover should be what I am looking for!
I just have to find a proper model to fit the interface... using H does not fit as I want C in the gas strictly =0 (H= infinite, right?).

Will think about that but meanwhile let me thank you for putting me on the (probably) right track!

Luca

Cyp December 16, 2012 05:21

If we get back to the diffusion example I mentioned in the other topic, we can merely consider the equation in the beta phase :

\nabla \cdot  D_{\beta}\nabla C_{\beta} = 0,

Both phases are connected through a thermodynamic equilibrium condition at the interface:
C_{\beta}= C_{eq} \texttt{  at  } \mathcal{A}_{\beta \gamma}


What you look for is an partial differential equation that govern
C = \alpha C_{\beta}
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}

Following the same line of development of the other topic, you can express the derivative of C as:
\nabla C = \alpha \nabla C_{\beta} + C_{\beta}\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}
+ D_{\beta}\nabla C_{\beta}\nabla \alpha + \nabla \cdot D C_{\beta}\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}\nabla \alpha + D_{\beta}\nabla C_{\beta}\cdot\nabla \alpha

This additional term represents the interfacial jump condition and the flux continuity. From the thermodynamic equilibrium relationship you can write:
\nabla\cdot D \nabla C = \nabla \cdot D C_{eq}\nabla \alpha + D_{\beta}\nabla C_{\beta}\cdot\nabla \alpha

I think it is something like that you are looking for. May be you should impose a zero flux at the interface (you get rid of D_{\beta}\nabla C_{\beta}\cdot\nabla \alpha in the previous equation)...

Best regards,
Cyp

voingiappone December 16, 2012 23:36

Hello Cyp, and thanks for the development of the equations!

Taking another way I got to the same result but in my case it was inserted in the context of the equation I shown some posts above. I left it untouched (by actually removing diffusion) as I was not able to compile the file with my poor C++ knowledge.

But let's say we can do like you kindly calculated and we end up with the equation (no flux at the interface):
Quote:

Originally Posted by Cyp (Post 397752)

\nabla\cdot D \nabla C = \nabla \cdot D C_{eq}\nabla \alpha

This represents the interfacial jump condition however I can leave it untouched as a function of C_{\beta} right?
Now I should add it to the main equation to take into account the time dependence:

\frac{\partial C}{\partial t}+\nabla \phi C + \nabla\cdot D \nabla C = 0

which, applying your calculations, should become:

\frac{\partial C}{\partial t}+\nabla \phi C + \nabla \cdot D C\nabla \alpha= 0

But the \nabla \phi C part should be treated the same way to include the alpha dependence right? Or am I missing something?

BTW, also thanks for showing me the [math] tag! Really useful!

Cyp December 17, 2012 03:46

Quote:

Originally Posted by voingiappone (Post 397830)
Hello Cyp, and thanks for the development of the equations!

Taking another way I got to the same result but in my case it was inserted in the context of the equation I shown some posts above. I left it untouched (by actually removing diffusion) as I was not able to compile the file with my poor C++ knowledge.

But let's say we can do like you kindly calculated and we end up with the equation (no flux at the interface):


This represents the interfacial jump condition however I can leave it untouched as a function of C_{\beta} right?
Now I should add it to the main equation to take into account the time dependence:

\frac{\partial C}{\partial t}+\nabla \phi C + \nabla\cdot D \nabla C = 0

which, applying your calculations, should become:

\frac{\partial C}{\partial t}+\nabla \phi C + \nabla \cdot D C\nabla \alpha= 0

But the \nabla \phi C part should be treated the same way to include the alpha dependence right? Or am I missing something?

BTW, also thanks for showing me the [math] tag! Really useful!


Hello,

Actually, \nabla \cdot D C_{eq}\nabla \alpha is an additional source term. If you re-do the equations development starting from
\frac{\partial C}{\partial t}+\nabla \cdot \phi C + \nabla\cdot D \nabla C

You will obtain (neglecting the flux term), something like :
\frac{\partial C}{\partial t}+\nabla \cdot \phi C + \nabla\cdot D \nabla C = \nabla \cdot D C_{eq}\nabla \alpha

The rhs must be treated as an explicit term in your C equation. C_{eq} is evaluated from the thermodynamics. It depends on C_{\gamma} or be constant, it's up to you. As a starting point, I suggest you to consider it as a constant value.

Regards,
Cyp

voingiappone December 18, 2012 05:33

Quote:

Originally Posted by Cyp (Post 397850)
Hello,

Actually, \nabla \cdot D C_{eq}\nabla \alpha is an additional source term. If you re-do the equations development starting from
\frac{\partial C}{\partial t}+\nabla \cdot \phi C + \nabla\cdot D \nabla C

You will obtain (neglecting the flux term), something like :
\frac{\partial C}{\partial t}+\nabla \cdot \phi C + \nabla\cdot D \nabla C = \nabla \cdot D C_{eq}\nabla \alpha

The rhs must be treated as an explicit term in your C equation. C_{eq} is evaluated from the thermodynamics. It depends on C_{\gamma} or be constant, it's up to you. As a starting point, I suggest you to consider it as a constant value.

Regards,
Cyp

Howdy Cip,

got the point! I'll see if I can find some time to implement the equation in the code and run a test simulation so I can get back here and report the result.

Thanks!

Luca

voingiappone December 18, 2012 23:17

1 Attachment(s)
Cyp,

you're the man! The problem was (almost) completely solved using your approach! I had some problems after inserting the alpha laplacian in the .C file as I had to specify to solve it as an explicit term but after figuring out how to fix that, the simulation runs smoothly.

The tracer is still coming out from the liquid phase but only to a minor extent, due to the "froth" like interface smearing you have when a lot of bubbles hit the surface but that's something that really happens in the reactor so I think that won't bother me that much.

I attach a screenshot of the simulation where you can see that.

Man, thank you so much for your help! I can now move on to my next problem: coupling my algae growth curve in the system, but that will be done afterwards and not using openFOAM.

Luca

vishal_s October 5, 2013 02:51

Hello luca,I am trying to simulate jet(liquid) in a different density medium (air) . The jet gets vaporized once it reaches it's specific vaporization temperature/pressure. I read ur thread but I am bit confused, how to add a transport scalar to the 2nd phase (air) to take into account vapor in interPhasechange Foam. What changes I have to make in the solvers??

Please help

voingiappone October 8, 2013 22:45

Quote:

Originally Posted by vishal_s (Post 455146)
take into account vapor in interPhasechange Foam

Hello vishal, first of all let me tell you that I don't use that particular software so I may not be the best person to ask to....
For what I was able to read on the forum, that solver is used when a liquid phase vaporizes into its vapor, that is you are trying to make water to vaporize in air (three phases). I don't know if you can do that with that specific solver.

It may look like you can use my approach in your situation but I think it is very different: in my case I have a low concentration parameter that changed the phase, not the liquid itself so its effect on liquid/gas properties is negligible. For this reason you have to keep in mind something:

- adding a new equation (at least in interFoam) creates a sort of new layer over with the new scalar is solved. What I mean is that your main equations get solved beforehand and then, the new one is solved using (where needed) p and U values from the main ones. This makes the new equation "phase insensitive" that's in the end the reason why I created this thread.

- when you add a new parameter with the alpha1 dependence you get two of them because:

C = C_{liq} \cdot alpha1 + C_{gas} \cdot (1 - alpha1)

for this reason, if you add a new scalar it won't get solved together with the main equations and it will have to be adjusted to consist of just one phase (because you actually only need three). So, I think this is not the right approach. Moreover in your case, the liquid phase vanishes as it vaporizes so major changes that cannot be neglected happen in the momentum equation!

The only thing I can think of is that you solve the problem using the 2 phases available and make adjustments for each time step to the gas phase. What i mean is that once solved the time step, your newly added equation can average the gas properties with those of air and then save the calculated data in a dedicated file. You may not want to average using alpha1 (or you get the liquid in the way again) but maybe with vapor pressure may work...?? Something like:

C_{air} = \frac{P^{*}(T)-P}{P^{*}(T)} \cdot (1 - alpha1)

This can be the air concentration in the vapor cloud (not in the liquid) where at the interface you have a pressuse equal to the vapor pressure for the given temperature and as you move away from the interface you have a lower pressure, up to zero where air concentration in air becomes, obviously 1. Doing this you can see the air concentration distribution. Then you can use this parameter to average gas properties between the air and the vapor values and write a dimensioned scalar file in each time step folder if you wish.
However this approach will solve the jet flow phase change using ONLY the vapour transport properties stored in the transportProperties file. If they are similar to air properties (and they may be) the solution will not be that different from the right one.

I am sorry if this is not what you are looking for but if you need help let's discuss it here. If you like it I can tell you how to do it.

Luca

vishal_s October 18, 2013 01:58

Thanks Luca. I will give a try and let you know.

Jonas Ansoni April 9, 2014 17:45

1 Attachment(s)
Quote:

Originally Posted by voingiappone (Post 398036)
Howdy Cip,

got the point! I'll see if I can find some time to implement the equation in the code and run a test simulation so I can get back here and report the result.

Thanks!

Luca

Hello Luca!

I'm having the same problem you got, as shown in attached Fig. I'm implementing the equation, but I'm having trouble to insert the term \nabla \cdot D C_{eq}\nabla\alpha.

How have you implemented this term in the .C file?

Thank you in advance

voingiappone April 10, 2014 00:27

Hi Jonas,

you did not explain which problem you are facing so I just try to guess. The only problem I had to implement the code in the C file was that you use an already calculated field (alpha1) in a laplacian so you have to calculate it as an explicit:

Code:

fvc::laplacian(DC, alpha1)
This was the only problem I had that did not let me compile the file.

Let me know if you solved the problem.

cheers

voingiappone April 10, 2014 00:28

Sorry, I forgot to specify that DC is a new variable defined as the product between the diffusion coefficient and the equilibrium concentration.

lzhou November 25, 2014 11:38

Hello Luca,

I am doing a simulation with interFoam and I would also like to add a concentration equation (the concentration of sediment actually) into the solver. The problem is similar as yours which is that I want the concentration in the air to be zero and the equation only be solved in the water zone. Reading this thread gives me a lot of help. I just want to ask how did you define the equilibrium concentration at the interface ? As the concentration in air is zero and the concentration in water is unknown and to be calculated ... Maybe I don't understand correctly...Can you please help me with this problem ?

Thanks !

voingiappone November 25, 2014 14:19

Hi Izhou!

First of all, if you're dealing with suspended solids, I seem to recall that a multiphase interFoam is available where you can model them. Did you give it a shot?

Anyway, if you want to proceed with this approach, I will tell you what I did. I went through different approaches, the most correct of which is the one where you set the eq. concentration equal to the bulk concentration, to be solved for. As you correctly noticed. I was not able at first to compile or run the code, so I changed the approach and set it as a constant. I eventually succeeded but used that with the phase change.

I was going to use this to calculate the mixing time and validating with an experiment so I made some assumptions:

- whatever the initial pulse concentration, the interface concentration would be 0.
- whatever the initial concentration, the interface concentration at t=t_mix would be exactly equal to the bulk concentration

As you can see in the various screenshots, the liquid gas interface in an air lift reactor is on the top part (bubbles excluded) so I decided to approximate the interface concentration with the perfectly mixed concentration, that is the concentration attained after t_mix by using the initial NaCl in the pulse.

This is an approximation but the error appears to be minimal when compared with a mixing time experiment in the actual reactor. If your concentrations are low enough and the instruments are sensitive enough to measure small pulses, you can actually verify this as I did and the result will not differ from the simulation. However, I don't know your situation and your simulation conditions so, the only thing is trying to set it up and verify if it works for you. If it doesn't, you should add C_eq=C as a variable and solve for it as usual. This is the case when you do have the phase change and the interface concentration is solved. I ended up with the equation (in the .cpp file):

Code:

fvm :: laplacian ( alpha1 , C )
Of course, this is only the part referring to the interface equilibrium. Either of the two approaches should do the trick.
Let me know if you made it run successfully.

Cheers

Luca

lzhou December 2, 2014 10:02

5 Attachment(s)
Hi Luca,

Thanks very much for your kindly help and I am sorry for the delay of my response.

I have tried both of the ways you remanded and I have also looked into the "multiPhaseInterFoam". I'm not sure this solver can be used to solve my problem as it seems to me that the mixture of water and sediment is not the same as the mixture of water and oil.. Meanwhile, there is another solver called "settlingFoam" which seems to be interesting, but it does not include the free surface effect so I think I still have to use your method to solve this problem.:)

My case is the water jet under a sluice gate as showed in fig.case. And I added the concentration equation into interFoam solver to calculate the suspended load concentration of sediment.

I first set the Ceq as the concentration of sediment at steady state. During the development (fig.C1_development) it seems to cause a little error at the interface but after reaching steady state the result is fine(fig.C1_steady). So I thought the problem during the development is because the Ceq should be the same as the concentration at the interface. Then I set the value of Ceq as the C value at the interface as follow:
Code:

surfaceScalarField Cf = fvc::interpolate(C);

fvScalarMatrix CEqn
    (
        fvm::ddt(C)
        + fvm::div(phiC, C)
        - fvm::laplacian(nuEff, C)
        ==
        - fvc::laplacian(nuEff*(Cf), alpha1)
    );
     
  CEqn.solve();

With this equation, the problem during the development is solved but the final result show more convection into the air phase. I didn't find any other way to have a better result...

Later in my case I have to add the falling velocity of the sediment into the concentration equation and the concentration at steady state will then no longer be uniform, so for now I will use the second way to implement the concentration equation...But I'll continue to look at this problem.

And there is another question. I don't know what is the value of DC,which is the mass diffusivity, so I just use the effective viscosity...Is this reasonable ? Do you have any suggestions ?

Thanks again for your reply and have a nice day !

Lu ZHOU

AbbasRahimi July 28, 2016 10:14

Quote:

Originally Posted by voingiappone (Post 521069)
Hi Izhou!

First of all, if you're dealing with suspended solids, I seem to recall that a multiphase interFoam is available where you can model them. Did you give it a shot?

Anyway, if you want to proceed with this approach, I will tell you what I did. I went through different approaches, the most correct of which is the one where you set the eq. concentration equal to the bulk concentration, to be solved for. As you correctly noticed. I was not able at first to compile or run the code, so I changed the approach and set it as a constant. I eventually succeeded but used that with the phase change.

I was going to use this to calculate the mixing time and validating with an experiment so I made some assumptions:

- whatever the initial pulse concentration, the interface concentration would be 0.
- whatever the initial concentration, the interface concentration at t=t_mix would be exactly equal to the bulk concentration

As you can see in the various screenshots, the liquid gas interface in an air lift reactor is on the top part (bubbles excluded) so I decided to approximate the interface concentration with the perfectly mixed concentration, that is the concentration attained after t_mix by using the initial NaCl in the pulse.

This is an approximation but the error appears to be minimal when compared with a mixing time experiment in the actual reactor. If your concentrations are low enough and the instruments are sensitive enough to measure small pulses, you can actually verify this as I did and the result will not differ from the simulation. However, I don't know your situation and your simulation conditions so, the only thing is trying to set it up and verify if it works for you. If it doesn't, you should add C_eq=C as a variable and solve for it as usual. This is the case when you do have the phase change and the interface concentration is solved. I ended up with the equation (in the .cpp file):

Code:

fvm :: laplacian ( alpha1 , C )
Of course, this is only the part referring to the interface equilibrium. Either of the two approaches should do the trick.
Let me know if you made it run successfully.

Cheers

Luca

Hi Luca,

I'm struggling with the same problem you did. I tried to used the steps you have provided to develop a new scalarTransportFoam to avoid diffusion of tracer into other phase. No success so far! Is there any chance you share the equations you put in scalarTransportFoam solver?

Thanks and regards,
Abbas

nimasam July 29, 2016 04:15

Dear Abbas

look following papers:
1- A unified single-field model framework for Volume-Of-Fluid
simulations of interfacial species transfer applied to bubbly flows

2-Numerical simulation of species transfer across fluid interfaces in
free-surface flows using OpenFOAM

AbbasRahimi July 29, 2016 09:59

Thank you Nima for response. Have you tried the work in any of these papers?

nimasam July 30, 2016 05:03

yes, i did, it works :)

AbbasRahimi July 30, 2016 21:39

Quote:

Originally Posted by nimasam (Post 612039)
yes, i did, it works :)

Is there any chance you share your code with me? My email is rahimi.abas@gmail.com

kuria May 1, 2018 08:06

Hey voingiappone

It would be great if you could upload the code you used to fix this problem :)


Thank you!

roshanak.rabiee878 November 26, 2020 16:05

mixing time tracer
 
hi
I was wondering if you solved your problem or not.
I had the same case, I want to calculate the mixing time. I created a new solver based on interFOAM and added a new transport equation. my problem is two-point?
I need to inject this tracer inside of domain, not at the inlet
time of injection is not at the beginning, it applied at a specific time and specific duration
How Can I do that? can you help me?
Quote:

Originally Posted by voingiappone (Post 397276)
Hello again everybody,

as I already explained here, I want to insert a tracer in my simulation to calculate the mixing time of the reactor.
I want to consider only the liquid phase as I am going to insert a certain amount of NaCl as tracer and follow the time evolution of the concentration with a probe.

Now I was just able to add the tracer using as a base the wiki page on how to add temperature to icoFoam. However, if I run the simulation I get my tracer propagating out of the liquid phase through the interface into the gas.

Watching this thread (http://www.cfd-online.com/Forums/ope...interfoam.html) I made an idea of how to read physical properties for different phases but I cannot figure out how to make my equation to be solved only for the liquid phase.
To put things simple, now it looks like it is just overlapping the liquid-gas system (i have a global equation, a global diffusion coefficient, etc). See attached screenshot (tracer concentration profile on the right).

Is there anybody who can help me with this task?

Thanks!


sadra2003 June 7, 2022 08:33

Quote:

Originally Posted by Jonas Ansoni (Post 485038)
Hello Luca!

I'm having the same problem you got, as shown in attached Fig. I'm implementing the equation, but I'm having trouble to insert the term \nabla \cdot D C_{eq}\nabla\alpha.

How have you implemented this term in the .C file?

Thank you in advance


Hello Jonas,

I am trying to simulate a single bubble movement in a solution of water and sugar with interIsoFoam solver, OF2112. I modified the solver and coupled the density, surface tension and viscosity of solution to the concentration of sugar which is different in various parts of the domain. In order to solve the distribution of sugar (a passive scalar) in the geometry, I added a new equation to the solver as below:

fvScalarMatrix CEqn
(
fvm::ddt(C)
+ fvm::div(phi, C)
- fvm::laplacian(dc,C)
==
fvOptions(C)
);


CEqn.relax();
fvOptions.constrain(CEqn);
CEqn.solve();
fvOptions.correct(C);

Now, the problem is, sugar concentration penetrates inside the bubble which is not correct. I would like to know how can I prevent sugar entering the bubble? I would be more than happy if you share your opinion with me.


All times are GMT -4. The time now is 18:02.