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

1D Stefan problem validation of new solver

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   September 24, 2020, 11:13
Default 1D Stefan problem validation of new solver
  #1
New Member
 
Join Date: Mar 2020
Posts: 12
Rep Power: 4
Huihui Xiao is on a distinguished road
Hello everyone,

I developed a new evaporation solver based on some excellent solvers (interFoam, interPhaseChangeFoam and evapVOFHardt), and when I validated the 1D Stefan problem, the velocity field is not correct (see attached), and some key code can be found as follows:
1. Source terms in alphaEqn (the same as the source terms in interPhaseChangeFoam)

// explicit source term
volScalarField Su
(
IOobject
(
"Su",
runTime.timeName(),
mesh
),
- dotM/rho1
);


//zeroField Sp;

// implicit source term
volScalarField Sp
(
IOobject
(
"Sp",
runTime.timeName(),
mesh
),
fvc::div(phi) - dotM*(1.0/rho2-1.0/rho1)
);

2. UEqn with new source terms
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
- fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U) // new
+ MRF.DDt(rho, U)
+ turbulence->divDevRhoReff(rho, U)
==
fvOptions(rho, U)
);

3. TEqn

Info << "Solve TEqn" << endl;
{
volScalarField k = alpha1*k1 + alpha2*k2; // one-field formulation
volScalarField rhoCp = alpha1*rho1*cp1 + alpha2*rho2*cp2;
// one-field formulation; rho:density; cp: specific heat capacity

surfaceScalarField alphaPhi = (rhoPhi-phi*rho2)/(rho1-rho2);
// hint: rhoPhi = alpha1*rho1*phi + (1-alpha1)*rho2*phi
// rho1 must not be equal to rho2
surfaceScalarField rhoCpPhi = alphaPhi*(rho1*cp1-rho2*cp2)+phi*rho2*cp2;
// rhoCpPhi = alpha1*rho1*cp1*phi + (1-alpha1)*rho2*cp2*phi

fvScalarMatrix TEqn
(
fvm::ddt(rhoCp,T)
+ fvm::div(rhoCpPhi,T)
- fvm::laplacian(k,T)
==
fvm::Sp(fvc::ddt(rhoCp) + fvc::div(rhoCpPhi),T)
// new source term 01 (Coupled to mass source term)
+ fvm::Sp(TEqnSou,T)
- TEqnSou*Tsat
// new source term 02 (cooling due to evaporation)
);

TEqn.relax();
TEqn.solve();

Info<< "min/max(T) = " << min(T).value() << "; "
<< max(T).value() <<endl;
}

4. pEqn
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAUf, p_rgh) == fvc::div(phi) - dotM*(scalar(1.0)/rho2 - scalar(1.0)/rho1) //new
);

5. source terms (modified from evapVOFHardt)

volScalarField jLoc = (T-Tsat)/(Rint*hEvap);

dotM = jLoc*magGradAlpha;

forAll(alpha1.ref(),iCell)
{
if (magGradAlpha.ref()[iCell] > 0.0)
{
dotM[iCell] *= 1.0;
// only set source terms in interface cells
}

else
{
dotM[iCell] *= 0.0;
}
}

// the second term of the TEqn
TEqnSou = -magGradAlpha/Rint;

6. main program
#include "sourceTerms.H"
evapRate = fvc::domainIntegrate(dotM);

#include "TEqn.H"

Info<< "Evaporation rate is = " << evapRate.value() << "kg/s at t = " << runTime.timeName()
<< "s" << endl;

If anyone know how to fix the problem, it will be best, and if I have something wrong in my code, please let me know.

Many thanks,
Huihui
Attached Images
File Type: jpg test.jpg (32.8 KB, 33 views)
Huihui Xiao is offline   Reply With Quote

Old   February 8, 2021, 10:10
Default
  #2
New Member
 
hamdani hakim
Join Date: Mar 2019
Posts: 16
Rep Power: 5
kimou is on a distinguished road
Hi,

You work on 1D twho-phase-flow probleme, , Me too, i work on 1D two-phase-flow but compressible configuration . I hope that you resolve your problem



Please, I would like to know how you chose your discretizations (fvSchemes and fvSolution) for your problem ??

best regret,
Attached Images
File Type: png Capture du 2021-02-08 15-51-54.png (72.3 KB, 22 views)
kimou is offline   Reply With Quote

Reply

Tags
evaporation, vof

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
Problem with validation of regionCouple solver tladd OpenFOAM Verification & Validation 0 May 22, 2014 18:09
Is there a problem in the Euler solver? Combas SU2 4 March 28, 2014 15:48
[ANSYS Meshing] Problem with axisymmetric solver Rasel ANSYS Meshing & Geometry 0 June 20, 2012 04:35
mesh.update problem in a new FSI solver ICL OpenFOAM 0 October 8, 2011 14:16
patching problem unsteady solver yellow-stuff Main CFD Forum 0 September 25, 2009 01:26


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