
[Sponsors] 
January 20, 2012, 09:41 
conditional solving of transport equation

#1 
Member
Join Date: Nov 2011
Location: Berlin
Posts: 31
Rep Power: 7 
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<alpha<1). alpha is the liquid fraction, alpha = 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, 14:10 

#2 
Senior Member
David Gaden
Join Date: Apr 2009
Location: Winnipeg, Canada
Posts: 436
Rep Power: 14 
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 cutoff... 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, 04:35 

#3 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 1,182
Rep Power: 21 
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, 11:51 

#4  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,910
Rep Power: 27 
Quote:
volScalarField solveEq(pos(alphaalphaCutOff)); volScalarField doNotEq(1solveEq); 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 ) ); Code:
myEqn.setValues(...) 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. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

January 25, 2012, 04:47 

#5 
Member
Join Date: Nov 2011
Location: Berlin
Posts: 31
Rep Power: 7 
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, 11:19 

#6  
Senior Member
Alberto Passalacqua
Join Date: Mar 2009
Location: Ames, Iowa, United States
Posts: 1,910
Rep Power: 27 
Quote:
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. OpenQBMM  An opensource implementation of quadraturebased moment methods. To obtain more accurate answers, please specify the version of OpenFOAM you are using. 

January 30, 2013, 02:57 
Solidification in OpenFOAM

#7 
Member
,...
Join Date: Apr 2011
Posts: 92
Rep Power: 6 
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, 09:04 
conditional solving of transport equation

#8 
New Member
Alexander Nekris
Join Date: Feb 2015
Location: France
Posts: 12
Rep Power: 4 
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 nonequilibrium 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. 10e9 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/</SPAN>fields/</SPAN>GeometricFields/</SPAN>SlicedGeometricField/</SPAN>SlicedGeometricField.C</SPAN>”, 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 correctBoundaryConditionsfunction? Many thanks in advance! </SPAN> 

December 3, 2015, 10:43 

#9 
Senior Member
Olivier
Join Date: Jun 2009
Location: France, grenoble
Posts: 266
Rep Power: 10 
Hello,
I guess you should also loop over boundaryMesh, not only (internal) cells. regards, olivier 

December 3, 2015, 11:31 

#10 
New Member
Alexander Nekris
Join Date: Feb 2015
Location: France
Posts: 12
Rep Power: 4 
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 forAllloop 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 

Tags 
conditional, solidification, solving 
Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
High Courant Number @ icoFoam  Artex85  OpenFOAM Running, Solving & CFD  11  February 16, 2017 14:40 
SimpleFoam k and epsilon bounded  nedved  OpenFOAM Running, Solving & CFD  15  December 2, 2016 03:40 
transsonic nozzle with rhoSimpleFoam  Unseen  OpenFOAM Running, Solving & CFD  7  April 16, 2014 03:38 
Low Mach number Compressible jet flow using LES  ankgupta8um  OpenFOAM Running, Solving & CFD  7  January 15, 2011 14:38 
Could anybody help me see this error and give help  liugx212  OpenFOAM Running, Solving & CFD  3  January 4, 2006 19:07 