CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Equilibrium thru Interface

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree15Likes

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 23, 2018, 09:40
Default
  #21
Member
 
Peng Liang
Join Date: Mar 2014
Posts: 59
Rep Power: 10
tjliang is on a distinguished road
Quote:
Originally Posted by Cyp View Post
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
Thans a lot Cyp, then what does n mean in your equation? Do you think that I should use interfoam to start for my plasma-liquid interaction? I see a lot of other standard multiphase solvers in Openfoam and I am quite confused as to which one to choose.

Bests,

Peng
tjliang is offline   Reply With Quote

Old   April 23, 2018, 09:49
Default
  #22
Cyp
Senior Member
 
Cyprien
Join Date: Feb 2010
Location: Stanford University
Posts: 299
Rep Power: 16
Cyp is on a distinguished road
Quote:
Originally Posted by tjliang View Post
Thans a lot Cyp, then what does n mean in your equation? Do you think that I should use interfoam to start for my plasma-liquid interaction? I see a lot of other standard multiphase solvers in Openfoam and I am quite confused as to which one to choose.

Bests,

Peng

\mathbf{n}_{\beta \gamma} is the normal vector to the gas/liquid interface. In your problem, I will say it depends. If you have an immiscible interface between your two fluids. Yes, you can go for a Volume of Fluid solver as a starting point.
Cyp is offline   Reply With Quote

Old   July 2, 2018, 15:29
Default
  #23
New Member
 
Rajesh Singh
Join Date: Jun 2010
Posts: 15
Rep Power: 14
rsingh7 is on a distinguished road
Dear FOAMer,

I have imported StarCCM+ mesh to openfoam 5.0 and quality of mesh was also good. The CCM mesh is trimmed mesh with polyhedral cells. The flow simulation for interFoam with species transport equation (Haroun formulation)diverges after some iterations. Below is report of the checkMesh. Help for resolving this issue would be appreciated.


Thanks in advance

*******************************
Create time

Create polyMesh for time = 0

Time = 0

Mesh stats
points: 2709670
faces: 7273254
internal faces: 6550582
cells: 2217107
faces per cell: 6.23508
boundary patches: 10
point zones: 0
face zones: 1
cell zones: 1

Overall number of cells of each type:
hexahedra: 1558763
prisms: 17795
wedges: 1051
pyramids: 73
tet wedges: 51
tetrahedra: 184093
polyhedra: 455281
Breakdown of polyhedra by number of faces:
faces number of cells
4 4
5 43
6 76272
7 247451
8 4264
9 70800
10 3318
11 356
12 16481
13 700
14 204
15 35328
16 60

Checking topology...
Boundary definition OK.
Cell to face addressing OK.
Point usage OK.
Upper triangular ordering OK.
Face vertices OK.
Number of regions: 1 (OK).

Checking patch topology for multiply connected surfaces...
Patch Faces Points Surface topology
floweriod1 110302 113851 ok (non-closed singly connected)
flow:sheet1 229147 166554 ok (non-closed singly connected)
flow:sheet2 228457 166118 ok (non-closed singly connected)
flowutg 9532 9665 ok (non-closed singly connected)
flow:wallt 9576 9998 ok (non-closed singly connected)
floweriod2 110292 113851 ok (non-closed singly connected)
flow:inl 5071 5594 ok (non-closed singly connected)
flow:wallb 10551 11107 ok (non-closed singly connected)
flowutl 5792 6040 ok (non-closed singly connected)
flow:ing 3952 4069 ok (non-closed singly connected)

Checking geometry...
Overall domain bounding box (-0.0133625 -0.013 -0.0851) (0.0133628 0.013 0.007)
Mesh (non-empty, non-wedge) directions (1 1 1)
Mesh (non-empty) directions (1 1 1)
Boundary openness (4.18881e-14 -2.88815e-15 1.82866e-17) OK.
Max cell openness = 2.78469e-16 OK.
Max aspect ratio = 9.20971 OK.
Minimum face area = 1.10391e-11. Maximum face area = 6.4e-07. Face area magnitudes OK.
Min volume = 5.01959e-16. Max volume = 5.12e-10. Total volume = 3.4427e-05. Cell volumes OK.
Mesh non-orthogonality Max: 57.8894 average: 11.7876
Non-orthogonality check OK.
Face pyramids OK.
Max skewness = 3.2775 OK.
Coupled point location match (average 0) OK.

Mesh OK.

End
rsingh7 is offline   Reply With Quote

Old   June 7, 2022, 08:55
Default
  #24
Member
 
sadra mahmoudi
Join Date: Feb 2021
Location: Austria
Posts: 34
Rep Power: 3
sadra2003 is on a distinguished road
Quote:
Originally Posted by Cyp View Post
@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,

Hello Cyp,

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.
sadra2003 is offline   Reply With Quote

Old   June 7, 2022, 08:57
Default a question
  #25
Member
 
sadra mahmoudi
Join Date: Feb 2021
Location: Austria
Posts: 34
Rep Power: 3
sadra2003 is on a distinguished road
Quote:
Originally Posted by Astrodan View Post
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

Mara 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

Hello Astrodan,

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.
sadra2003 is offline   Reply With Quote

Reply

Tags
drop, interface position, interfoam, openfoam

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Wind turbine simulation Saturn CFX 58 July 3, 2020 02:13
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


All times are GMT -4. The time now is 07:24.