# conditional solving of transport equation

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

 January 20, 2012, 08:41 conditional solving of transport equation #1 Member   Join Date: Nov 2011 Location: Berlin Posts: 31 Rep Power: 14 hi people, I would like to set up a solver for a solidification process with solving a transport equation for a species only in the liquid phase of the sytem. There is a volScalarField alpha which defines the state of phase (0 complete solid. How can I define a solver which works only in the non solid part of the domain like: // definition of eq. {liqEqn = ....} // conditional solving of liqEqn if (alpha != 0) liqEqn.solve(); Is there somewhere a similar case/tutorial/documentation reference? thank you for advice, dzi (I use OF 2.01)

 January 22, 2012, 13:10 #2 Senior Member   David Gaden Join Date: Apr 2009 Location: Winnipeg, Canada Posts: 437 Rep Power: 21 It is very difficult to only use portions of the mesh for the matrix solution in OpenFOAM. You'd probably have to create a new temporary mesh, and create new variables on it - then you'd have to create new boundaries where it is cut-off... Rather than that, you probably want to work with the full mesh, and modify the matrix so that the portion from the solid cells reduce to a trivial equation. I'm thinking you could create a custom preconditioner for your matrix. I don't know exactly what you want to do to the matrix to achieve this, though. __________________ ~~~ Follow me on twitter @DavidGaden

 January 23, 2012, 03:35 #3 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 29 I think the easiest solution is to multiply all terms in the equation with alpha. __________________ *On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer.

January 24, 2012, 10:51
#4
Senior Member

Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
Quote:
 Originally Posted by dzi hi people, I would like to set up a solver for a solidification process with solving a transport equation for a species only in the liquid phase of the sytem. There is a volScalarField alpha which defines the state of phase (0 complete solid. How can I define a solver which works only in the non solid part of the domain like: // definition of eq. {liqEqn = ....} // conditional solving of liqEqn if (alpha != 0) liqEqn.solve(); Is there somewhere a similar case/tutorial/documentation reference? thank you for advice, dzi (I use OF 2.01)
You can simply define a flag in this way:

volScalarField solveEq(pos(alpha-alphaCutOff));
volScalarField doNotEq(1-solveEq);

which is 1 when alpha > alphaCutOff, and 0 elsewhere. Then there are two possible solutions, depending on what you are trying to do (momentum equation or scalar equation?)

- Momentum equation:

Code:
```fvMatrix UEqn
(
solveAlpha*
(
)+
doNotSolve*
(
fvm::Sp(coeff,yourVariable)   // Set the variable to zero or to a value here
)
);```
- Scalar equations: you can use what done above, or drop elements from the matrix and fix the value of the solution directly using
Code:
`myEqn.setValues(...)`
See for example OpenFOAM/OpenFOAM-2.1.x/applications/solvers/multiphase/bubbleFoam/wallDissipation.H for an example.

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.

 January 25, 2012, 03:47 #5 Member   Join Date: Nov 2011 Location: Berlin Posts: 31 Rep Power: 14 thank you for the replies, for me the easiest solution is the suggestion from akidess to multiply all, or parts of the equation with alpha. Looks like it can be solved and I get something out which makes sense. The other suggestions also sound interesting. I will try if I come to a limit with the first solution, but on the first glance they seem to be more sophisticated. Thank you again for helping on this topic! dzi

January 25, 2012, 10:19
#6
Senior Member

Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,912
Rep Power: 36
Quote:
 Originally Posted by dzi thank you for the replies, for me the easiest solution is the suggestion from akidess to multiply all, or parts of the equation with alpha. Looks like it can be solved and I get something out which makes sense.
This depends on what you are solving for, but in general, multiplying by zero both sides of the equation without doing anything else will introduce a (0 = 0) equation and make the linear system singular.

If you are solving the momentum equation, you will also have to deal with the problem of the central coefficient going to zero, which will lead to a segmentation fault when you calculate H/A. You can check how this problem is addressed in compressibleTwoPhaseEulerFoam.

Best,
__________________
Alberto Passalacqua

GeekoCFD - A free distribution based on openSUSE 64 bit with CFD tools, including OpenFOAM. Available as in both physical and virtual formats (current status: http://albertopassalacqua.com/?p=1541)
OpenQBMM - An open-source implementation of quadrature-based moment methods.

To obtain more accurate answers, please specify the version of OpenFOAM you are using.

 January 30, 2013, 01:57 Solidification in OpenFOAM #7 Member   ,... Join Date: Apr 2011 Posts: 92 Rep Power: 14 Hi I am solving a binary alloy solidification problem with OpenFOAM. For the species equation, which seems to be the most critical equation especially at the region where channels form, I have used zero grad boundary conditions; which causes a flow into the wall. Have you guys also used grad(CL) = 0, and grad(C) = 0 at the boundaries?

 December 3, 2015, 08:04 conditional solving of transport equation #8 Member   Alexander Nekris Join Date: Feb 2015 Location: France Posts: 32 Rep Power: 11 Hello dear foamers, I want to revive this old thread with a new question. I wonder if it is possible in OpenFOAM to solve a transport equation only for cells whose variable values are larger than a certain value. Like: solve Eqn. only for values >0… Let me explain: I simulate a fluid consisting of two gas fractions gas1 and gas2. Gas2 does not cover the whole simulation domain but is present in places. So, basically, there are regions where the concentration of gas2 is zero. Consequently, the partial density rho of gas2 in those regions is zero too. In order to simulate a thermal non-equilibrium I need to solve the energy equation for both gases separately. I want to solve the energy equation for gas2 only for regions, where the density of gas2 is larger than a certain value (e.g. 10e-9 or something…). My energy equation is taken from rhoCentralFoam, where the internal energy is solved in two steps: Step1: the predictor step without diffusion term: solve ( fvm::ddt(rhoE) + fvc::div(phiEp) - fvc::div(sigmaDotU) ); Than calculating internal energy and correcting boundary conditions: e = rhoE/rho - 0.5*magSqr(U); e.correctBoundaryConditions(); thermo.correct(); rhoE.boundaryField() = rho.boundaryField()* ( e.boundaryField() + 0.5*magSqr(U.boundaryField()) ); Step2: the diffusion correction equation. if (!inviscid) { solve ( fvm::ddt(rho, e) - fvc::ddt(rho, e) - fvm::laplacian(turbulence->alphaEff(), e) ); thermo.correct(); rhoE = rho*(e + 0.5*magSqr(U)); } The predictor step is working just fine, but since I have to divide rhoE by rho of gas2 I understandably get the error, because there are cells in the simulation domain with rho = 0. What I did is: forAll(e, celli) { if (rho [celli] > 0.000001) { e [celli] = rhoE [celli]/rho [celli] - 0.5*magSqr(U[celli]); } else { e [celli] = e [celli]*0.0; } } But I have my suspicion, that this technique is highly inefficient concerning the computational speed. The next problems occur by executing: e.correctBoundaryConditions(); and thermo.correct(); When executing the command “correctBoundaryConditions” the solver crashes with the following report: #0 Foam::error:rintStack(Foam::Ostream&) at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? ... … ... Floating point exception (core dumped) I took a look at “src/OpenFOAM/fields/GeometricFields/SlicedGeometricField/SlicedGeometricField.C”, where I suppose the function “correctBoundaryConditions” is defined, but my OpenFOAM knowledge is too little to understand the file. My questions are: - Is there a better way than using the “forAll” loop, like I did? - How could I handle the printStack -error coming up by calling the correctBoundaryConditions-function? Many thanks in advance!

 December 3, 2015, 09:43 #9 Senior Member   Olivier Join Date: Jun 2009 Location: France, grenoble Posts: 272 Rep Power: 17 Hello, I guess you should also loop over boundaryMesh, not only (internal) cells. regards, olivier

 December 3, 2015, 10:31 #10 Member   Alexander Nekris Join Date: Feb 2015 Location: France Posts: 32 Rep Power: 11 Hello Olivier, hello all. Olivier, thank you for the quick reply. I just did it like this: Code: ``` forAll(e.boundaryField(), patchi) { fvPatchScalarField& ePatch = e.boundaryField()[patchi]; forAll(ePatch, facei) { if (rho[facei] > 0.000001) { e.correctBoundaryConditions(); } else { e[facei] = e[facei]*0.0; } } }``` It compiles just fine. The same problem appears with the command thermo.correct(). What I did there is: Code: ``` forAll(e, celli) { if (e[celli] > 0.000001) { thermo.correct(); } }``` It compiles also fine. I guess, that the function e.correctBoundaryConditions() is only executed for the particular cell “I” at each forAll-loop and not over all boundary cells. Is it right? Because otherwise I would get the same error as without using forAll. The next challenge is the corrector step: Code: ``` solve ( fvm::ddt(rho, e) - fvc::ddt(rho, e) - fvm::laplacian(turbulence->alphaEff(), e) );``` Here I have the same error as with the correctBoundaryConditions()-function, namely: #0 Foam::error:rintStack(Foam::Ostream&) at ??:? #1 Foam::sigFpe::sigHandler(int) at ??:? ... … ... Floating point exception (core dumped) Is there some technique to solve this equation conditionally (for e > 0.00...)? Thanks in advance. Regards Alex

February 22, 2017, 10:26
#11
Member

Kalpana Hanthanan Arachchilage
Join Date: May 2015
Location: Orlando, Florida, USA
Posts: 30
Rep Power: 10
Hi Alberto,

I know it's an old post. However, I have tried the method you mentioned below. My problem is a simulation of a nano-droplet. the particle movement is simulated as a concentration field which is only defined in liquid phase. So the transport equation has to be solved only in liquid phase. I have tried several methods, but particles tend to diffuse into vapor side. I tried the method you mentioned, by multiplying solveEq by the transport equation. As expected particles doesn't diffuse in to vapor side. but, it's not mass conserved.

Do you have any idea, how to make it mass conserved?

Kalpana

Quote:
 Originally Posted by alberto You can simply define a flag in this way: volScalarField solveEq(pos(alpha-alphaCutOff)); volScalarField doNotEq(1-solveEq); which is 1 when alpha > alphaCutOff, and 0 elsewhere. Then there are two possible solutions, depending on what you are trying to do (momentum equation or scalar equation?) - Momentum equation: Code: ```fvMatrix UEqn ( solveAlpha* ( //Put your equation here )+ doNotSolve* ( fvm::Sp(coeff,yourVariable) // Set the variable to zero or to a value here ) );``` - Scalar equations: you can use what done above, or drop elements from the matrix and fix the value of the solution directly using Code: `myEqn.setValues(...)` See for example OpenFOAM/OpenFOAM-2.1.x/applications/solvers/multiphase/bubbleFoam/wallDissipation.H for an example. Best,

 March 4, 2017, 17:09 #12 Member   Paolo Capobianchi Join Date: Sep 2014 Posts: 35 Rep Power: 11 Hi kal1943335, I am trying to tackle a similar problem, but different in physics. I have used a similar approach as the one suggested by Alberto and I have also noticed loss of mass. I think that is a matter of boundary conditions at the the liqui-gas interface. If you try to solve a PDE in a defined region, e.g. in the liquid phase, you should specify an interfacial condition for the variable you are solving. I was thinking about somethig like a zero flux. In other words, if C is the variable of interest, I would impose: n . grad(C) = 0 where n is the local unit vector normal to the interface. However, in order to do so, you should know the positon of the interface at each instant. Maybe a possible strategy would be to onsider an iso-volume fraction surface and apply the BC's there. I am not totally sure it would work. I am still trying to sort this problem out. Best, Paolo

March 4, 2017, 19:43
#13
Senior Member

Arjun
Join Date: Mar 2009
Location: Nurenberg, Germany
Posts: 1,273
Rep Power: 34
Quote:
 Originally Posted by pablitobass Hi kal1943335, I am trying to tackle a similar problem, but different in physics. I have used a similar approach as the one suggested by Alberto and I have also noticed loss of mass. I think that is a matter of boundary conditions at the the liqui-gas interface. If you try to solve a PDE in a defined region, e.g. in the liquid phase, you should specify an interfacial condition for the variable you are solving. I was thinking about somethig like a zero flux. In other words, if C is the variable of interest, I would impose: n . grad(C) = 0 where n is the local unit vector normal to the interface. However, in order to do so, you should know the positon of the interface at each instant. Maybe a possible strategy would be to onsider an iso-volume fraction surface and apply the BC's there. I am not totally sure it would work. I am still trying to sort this problem out. Best, Paolo
I think the best way to solve this thing is to set up everything for whole mesh and not on parts. This way the issue of mass loss does not come.

Once you set up this way what you end up is number of points where Ap=0, and Src = 0.
This Ap = 0 is a problem for linear solver, but since Src = 0 visit the Ap of the matrix and set Ap = 1.0E-20 for all the Ap which are less than 1E-20.

Then next step is to solve this linear system with BiCG type linear solver and be done with it.

March 5, 2017, 16:16
#14
Member

Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 11
Quote:
 Originally Posted by arjun I think the best way to solve this thing is to set up everything for whole mesh and not on parts. This way the issue of mass loss does not come. Once you set up this way what you end up is number of points where Ap=0, and Src = 0. This Ap = 0 is a problem for linear solver, but since Src = 0 visit the Ap of the matrix and set Ap = 1.0E-20 for all the Ap which are less than 1E-20. Then next step is to solve this linear system with BiCG type linear solver and be done with it.
Hi Arjun ,

Thanks for your reply. I will try the solution you proposed which look appealing. However, I am afraid that in my case I should solve the transport equation in a sub-region of the domain. I have to deal with a PDE defined at the interface between two different fluids. My idea was to define the transported quantity in a small but finite region across the interface, detect at each time step the cells at the boundary and then apply the BC's there. At the moment I have only found a way to detect those cell. Now I am stuck on how to impose the conditions I want. I am not sure it will be an easy task to accomplish. I will keep you posted if I will find a solution to the problem.

Best,
Paolo

March 6, 2017, 12:19
#15
Member

Kalpana Hanthanan Arachchilage
Join Date: May 2015
Location: Orlando, Florida, USA
Posts: 30
Rep Power: 10
Thank you Arjun for your reply. I will try that method as well. I have tried this problem in several ways.

Here the diffusion coefficient, D is zero in vapor side and it's updated with alpha in every iteration.

1. I implemented my transport equation as follows.

fvm::ddt(C)
+fvm::div(phi,C)
-fvm::laplacian(D,C)
== 0

this method actually diffuse the field C into vapor side.

2.
volScalarField liquidField = max(min(alpha1,1),1e-80);

liquidField*fvm::ddt(C)
+liquidField*fvm::div(phi,C)
-liquidField*fvm::laplacian(D,C)
== 0

Still it diffused into vapor side.

3.
Here, I used an approach similar to alphaEqn. I used MULES to solve for C.

so the equation is like,

surfaceScalarField phiC
(
fvc::flux
(
phiAlpha,
C,
CScheme
)
+ fvc::flux
(
C,
CrScheme
)
);
MULES::explicitSolve(geometricOneField(),C, phi, phiC, SpC ,SuC, 1, 0);

laplacian term is explicitly calculated and included in SuC.

magGradAlphaFCutNorm is used instead of (1-alpha1)*alpha1 in alphaEqn as (1-alpha1)*alpha1 term didn't do the interface compression properly. I must note my C field is in the order of 1e-3. so the compression term has to be changed. As alpha varies from 0 to 1, the diffusion is not much dominant. but, as my variable is from 0 to 1e-3, the diffusion into vapor side is dominant. It provide good results when i change the order of magGradAlphaFcutNorm, but with time, it generate more mass in the liquid field. It's because, the compression term isn't physical.

And Paolo, please check the interTrackFoam implementation in openfoam extended versions. they solved similar problem. but, they used a moving mesh to simulate the interface and they discretized the interface surface, for surfactant transport equation.
I'm not sure whether it'll help you or not. but, have a look. And you have mentioned that applying zeroflux boundary condition at the isovolume surface. Do you have any idea how to implement such boundary condition. if so, that would be awesome.

Thank you very much for your inputs.

Kalpana.

Quote:
 Originally Posted by arjun I think the best way to solve this thing is to set up everything for whole mesh and not on parts. This way the issue of mass loss does not come. Once you set up this way what you end up is number of points where Ap=0, and Src = 0. This Ap = 0 is a problem for linear solver, but since Src = 0 visit the Ap of the matrix and set Ap = 1.0E-20 for all the Ap which are less than 1E-20. Then next step is to solve this linear system with BiCG type linear solver and be done with it.

March 7, 2017, 17:28
#16
Member

Paolo Capobianchi
Join Date: Sep 2014
Posts: 35
Rep Power: 11
Quote:
 Originally Posted by kal1943335 And Paolo, please check the interTrackFoam implementation in openfoam extended versions. they solved similar problem. but, they used a moving mesh to simulate the interface and they discretized the interface surface, for surfactant transport equation. I'm not sure whether it'll help you or not. but, have a look. And you have mentioned that applying zeroflux boundary condition at the isovolume surface. Do you have any idea how to implement such boundary condition. if so, that would be awesome. Thank you very much for your inputs. Kalpana.
Hi Kalpana,

I wasn't aware about that solver. I will give a look on it.

At the moment I still have no idea on how to implement the bc's. Hopefully I'm wrong, but a first glance seems a really though task... In case I will find a way, I'll let you know.

Thanks for the suggestion.

Best,
Paolo

May 8, 2018, 18:57
#17
New Member

Pei Li
Join Date: Nov 2015
Posts: 10
Rep Power: 10
Quote:
 Originally Posted by kal1943335 Thank you Arjun for your reply. I will try that method as well. I have tried this problem in several ways. Here the diffusion coefficient, D is zero in vapor side and it's updated with alpha in every iteration. 1. I implemented my transport equation as follows. fvm::ddt(C) +fvm::div(phi,C) -fvm::laplacian(D,C) == 0 this method actually diffuse the field C into vapor side. 2. volScalarField liquidField = max(min(alpha1,1),1e-80); liquidField*fvm::ddt(C) +liquidField*fvm::div(phi,C) -liquidField*fvm::laplacian(D,C) == 0 Still it diffused into vapor side. 3. Here, I used an approach similar to alphaEqn. I used MULES to solve for C. so the equation is like, surfaceScalarField phiC ( fvc::flux ( phiAlpha, C, CScheme ) + fvc::flux ( -fvc::flux(-phir, magGradAlphaFcutNorm, CrScheme), C, CrScheme ) ); MULES::explicitSolve(geometricOneField(),C, phi, phiC, SpC ,SuC, 1, 0); laplacian term is explicitly calculated and included in SuC. magGradAlphaFCutNorm is used instead of (1-alpha1)*alpha1 in alphaEqn as (1-alpha1)*alpha1 term didn't do the interface compression properly. I must note my C field is in the order of 1e-3. so the compression term has to be changed. As alpha varies from 0 to 1, the diffusion is not much dominant. but, as my variable is from 0 to 1e-3, the diffusion into vapor side is dominant. It provide good results when i change the order of magGradAlphaFcutNorm, but with time, it generate more mass in the liquid field. It's because, the compression term isn't physical. And Paolo, please check the interTrackFoam implementation in openfoam extended versions. they solved similar problem. but, they used a moving mesh to simulate the interface and they discretized the interface surface, for surfactant transport equation. I'm not sure whether it'll help you or not. but, have a look. And you have mentioned that applying zeroflux boundary condition at the isovolume surface. Do you have any idea how to implement such boundary condition. if so, that would be awesome. Thank you very much for your inputs. Kalpana.
Hi Kalpana. Could you explain more about how to make 'CScheme' work? I tried to use it in fvm::div(), but error: ‘CScheme’ was not declared in this scope always occurs. Thanks a lot.

May 8, 2018, 19:26
#18
Member

Kalpana Hanthanan Arachchilage
Join Date: May 2015
Location: Orlando, Florida, USA
Posts: 30
Rep Power: 10
Quote:
 Originally Posted by HectorLee Hi Kalpana. Could you explain more about how to make 'CScheme' work? I tried to use it in fvm::div(), but error: ‘CScheme’ was not declared in this scope always occurs. Thanks a lot.
Hi Hector,
CScheme and CrScheme are defined similar to interFoam. in alphaEqn.H, they defined words alphaScheme and alpharScheme as follows,

word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");

So I defined CScheme and CrScheme and included these two entries in fvSchemes.

As far as I understood, these words refer to the corresponding name in fvSchemes under div terms and flux is calculated using those discretization schemes. More information on this method is explained in Jasak's dissertation.

Kalpana

 May 15, 2018, 00:40 #19 New Member   Pei Li Join Date: Nov 2015 Posts: 10 Rep Power: 10 Thanks a lot, Kalpana. Now it is more clear with your explanation and the errors has disappeared.

 January 24, 2019, 03:54 #20 Member   AJAY BHANDARI Join Date: Jul 2015 Location: INDIA Posts: 57 Rep Power: 10 Hi all, I think my post best fits here. I want to solve the continuity equation on the mesh elements of my computational domain. But on one mesh element i want to solve the same continuity equation with some additional terms. Can anybody please suggest something in this regard. How should i start with. Any help will be appreciated. Best Ajay

 Tags conditional, solidification, solving