Weird result from 'phase field solver'. Need help.

mykkujinu2201 · May 7, 2021, 01:55

Dear Forum,

Phase field solver i made does not work well.

Currently, I am interested in bubble simulation in micro scale using OpenFOAM.

I first considered to use interFoam (VOF) but there was spurious velocity issue.
Thus, I changed to Coupled Level-Set Volume-of Fluid method (CLSVOF).
However, even for CLSVOF, spurious velocity issue is not solved.

Finally, I decide to use phase Field method, in which Cahn-Hilliard equation is additionally solved and surface tension is calculated not using curvature of interface but using chemical potential.

I read the references as below.

1.Jamshidi, F. CFD Simulation of an Air Bubble in Microfluidics using Volume of Fluid and Phase Field Methods in OpenFOAM. https://publikationen.bibliothek.kit.edu/1000078984 (2016) doi:10.5445/IR/1000078984.
2.Cai, X., Wörner, M. & Deutschmann, O. Implementation of a Phase Field Method in OpenFOAM® for Simulation of Spreading Droplets and Verification by Test Problems. 16.
3.Jamshidi, F. et al. On suitability of phase-field and algebraic volume-of-fluid OpenFOAM® solvers for gas–liquid microfluidic applications. Computer Physics Communications 236, 72–85 (2019).

Their method consists of equations below.

1) Chemical potential equation
$\phi=\lambda\left(\frac{1}{\epsilon^2}C(C^2-1)-\nabla^2C\right)$ , where $\phi$ , $\lambda$ , $\epsilon$ ,

are chemical potential, mixing energy density scale, interface thickness and order parameter, respectively.
Order parameter indicates each phase and interface.

$\rightarrow$ phase 1,

$\rightarrow$ phase 2,

$\rightarrow$ interface.

2) Cahn-Hilliard equation
$\frac{\partial C}{\partial t}+\left(\underline{u}\cdot\nabla\right)C=M\nabla^2 \phi$ , where $\underline{u}$ ,

are velocity and mobility, respectively.

3) Continuity equation
$\nabla \cdot \underline{u} = 0$

4) Momentum equation
$\frac{\partial (\rho \underline{u})}{\partial t}+\left(\underline{u}\cdot\nabla\right)\underline{u}=-\nabla p + \nabla \cdot \left[\mu\left(\nabla \underline{u} + \nabla^T \underline{u}\right)\right] +\rho \underline{g} +\underline{f}_{st}$ , where $\rho$ ,

, $\mu$ $\underline{g}$ , $\underline{f}_{st}$ are density, pressure, viscosity, gravitational acceleration and surface tension force, respectively.

Viscosity and density are calculated using the equations below.
$\mu=\frac{1}{2}\left[(C+1)\mu_1+(C-1)\mu_2\right]$
$\rho=\frac{1}{2}\left[(C+1)\rho_1+(C-1)\rho_2\right]$ ,
where each subscript indicates each phase.

Mixing energy density scale, $\lambda$ is related to surface tension as $\sigma\epsilon=\frac{2\sqrt{2}}{3}\lambda$ , where $\sigma$ is surface tension.

Surface tension force, $\underline{f}_{st}$ is calculated as $\underline{f}_{st}=-C\nabla \phi$ .

I followed the instructions in the references to make solver with OpenFOAM.
Based on pimpleFoam, I made phase field solver.

1. main file

Code:

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    pimpleFoam

Description
    Large time-step transient solver for incompressible, turbulent flow, using
    the PIMPLE (merged PISO-SIMPLE) algorithm.

    Sub-models include:
    - turbulence modelling, i.e. laminar, RAS or LES
    - run-time selectable MRF and finite volume options, e.g. explicit porosity

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulentTransportModel.H"
#include "pimpleControl.H"
#include "fvOptions.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "postProcess.H"

    #include "setRootCase.H"
    #include "createTime.H"
    #include "createMesh.H"
    #include "createControl.H"
    #include "createTimeControls.H"
    #include "createFields.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"

    turbulence->validate();

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

    Info<< "\nStarting time loop\n" << endl;

    while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"

        runTime++;

        Info<< "Time = " << runTime.timeName() << nl << endl;

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "CahnHilliardEq.H"
            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }
/*
            if (pimple.turbCorr())
            {
                laminarTransport.correct();
                turbulence->correct();
            }
*/
        }

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //

CahnHilliardEq.H

Code:

//Cahn-Hilliard Eq'n solving
//mixing energy
dimensionedScalar lambda(scalar(3.0)*sigma*eps/(Foam::sqrt(scalar(2.0))*scalar(2.0)));

//Chemical potential
pot=lambda/(eps*eps)*C*(C-scalar(1.0))*(C+scalar(1.0))-fvc::laplacian(lambda,C);
pot.correctBoundaryConditions();

//order parameter transport
fvScalarMatrix CEqn
(
    fvm::ddt(C) + fvm::div(phi,C) - fvc::laplacian(M,pot) 
);

CEqn.relax();
CEqn.solve();
C.correctBoundaryConditions();

UEqn.H

Code:

// properties update
mu==(scalar(0.5)*((scalar(1.0)+C)*mu1+(scalar(1.0)-C)*mu2));
rho==(scalar(0.5)*((scalar(1.0)+C)*rho1+(scalar(1.0)-C)*rho2));
rho.correctBoundaryConditions();
mu.correctBoundaryConditions();
surfaceScalarField rhoPhi=(fvc::interpolate(rho)*phi);

pot==lambda/(eps*eps)*C*(C*C-scalar(1.0))-lambda*fvc::laplacian(C);

f_st==-C*fvc::grad(pot);
f_st.correctBoundaryConditions();

// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);

fvVectorMatrix UEqn
(
    fvm::ddt(rho,U) + fvm::div(rhoPhi,U)
  + MRF.DDt(rho,U)
  - fvm::laplacian(mu,U)-(fvc::grad(U)&fvc::grad(mu))
== 
   f_st
   +rho*g
   +fvOptions(rho, U)
);

UEqn.relax();

fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
    solve(UEqn == -fvc::grad(p));

    fvOptions.correct(U);
}

pEqn.H

Code:

volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
surfaceScalarField phiHbyA
(
    "phiHbyA",
    fvc::flux(HbyA)
  + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi)
);

MRF.makeRelative(phiHbyA);

adjustPhi(phiHbyA, U, p);

tmp<volScalarField> rAtU(rAU);

if (pimple.consistent())
{
    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
    phiHbyA +=
        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
    HbyA -= (rAU - rAtU())*fvc::grad(p);
}

// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAtU(), MRF);

// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
    // Pressure corrector
    fvScalarMatrix pEqn
    (
        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
    );

    pEqn.setReference(pRefCell, pRefValue);

    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

    if (pimple.finalNonOrthogonalIter())
    {
        phi = phiHbyA - pEqn.flux();
    }
}

#include "continuityErrs.H"

// Explicitly relax pressure for momentum corrector
p.relax();

U = HbyA - rAtU()*fvc::grad(p);
U.correctBoundaryConditions();
fvOptions.correct(U);

However, when I tested the solver for static bubble case, wrong result came out.

transportProperties

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "constant";
    object      transportProperties;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
phase1
{
//	transportModel  Newtonian;
	rho 	rho [1 -3 0 0 0 0 0] 997;
	mu	mu  [1 -1 -1 0 0 0 0] 997.e-6;
}

phase2
{
//	transportModel  Newtonian;
	rho 	rho [1 -3 0 0 0 0 0] 1.2e-0;
	mu	mu  [1 -1 -1 0 0 0 0] 1.8372e-5;
}

parameter
{
	eps	eps [0 1 0 0 0 0 0] 1.e-5;
	M	M [-1 3 1 0 0 0 0] 1.e-11;//1.e-11;	
}

properties
{
	g	g	[0 1 -2 0 0 0 0] (0.0 0.0 0.0);
	sigma	sigma [1 0 -2 0 0 0 0] 7.283e-2;
}
	
//dummy	
transportModel  Newtonian;
nu	1.e-4;

All the boundary conditions are cyclic boundary conditions.

I doubt about surface tension force term but I do not know the exact reason.

If you have any comments about it, it would be really helpful.

Best regards,
Jun

mykkujinu2201 · May 11, 2021, 23:12

Dear forums,

I think I need to update the progress to solve this issue.

With the code above, I got diverged solution right after few iteration.

I tried to modify the code more and I found that using force flux is more stable than directly adding force term in UEqn.
Thus, I modified UEqn.H and pEqn.H comparing with buoyantBoussinesqPimpleFoam.
I also found that by incompressibility,
$\vec{f}_{st}=-C\nabla \phi$ can be rewritten as $\vec{f}_{st}=\phi\nabla C$
Also, I modified the Cahn-Hilliard Eq'n part.

CahnHilliardEq.H

Code:

//Cahn-Hilliard Eq'n solving
//mixing energy
dimensionedScalar lambda(scalar(3.0)*sigma*eps/(Foam::sqrt(scalar(2.0))*scalar(2.0)));


//Chemical potential
pot==lambda/(eps*eps)*C.oldTime()*(C.oldTime()-scalar(1.0))*(C.oldTime()+scalar(1.0))-lambda*fvc::laplacian(C.oldTime());
pot.correctBoundaryConditions();

//order parameter transport
fvScalarMatrix CEqn
(
    fvm::ddt(C) + fvm::div(phi,C) - M*fvc::laplacian(pot) 
);

CEqn.relax();
CEqn.solve();
C.correctBoundaryConditions();

pot==lambda/(eps*eps)*C*(C-scalar(1.0))*(C+scalar(1.0))-lambda*fvc::laplacian(C);
pot.correctBoundaryConditions();

UEqn.H

Code:

// properties update
dimensionedScalar mu1(nu1*rho1);
dimensionedScalar mu2(nu2*rho2);
mu==(scalar(0.5)*((scalar(1.0)+C)*mu1+(scalar(1.0)-C)*mu2));
mu.correctBoundaryConditions();
volScalarField nu(mu/rho);
nu.correctBoundaryConditions();
rho==(scalar(0.5)*((scalar(1.0)+C)*rho1+(scalar(1.0)-C)*rho2));
rho.correctBoundaryConditions();
surfaceScalarField muf(fvc::interpolate(mu));
surfaceScalarField rhoPhi(fvc::interpolate(rho)*phi);

surfaceScalarField f_stf=fvc::interpolate(pot)*fvc::snGrad(C);
// Solve the Momentum equation
MRF.correctBoundaryVelocity(U);

fvVectorMatrix UEqn
(
    fvm::ddt(rho,U) + fvm::div(rhoPhi,U)
  + MRF.DDt(rho,U)
  - fvm::laplacian(muf,U)-(fvc::grad(U)&fvc::grad(muf))//fvc::div((mu*dev(T(fvc::grad(U)))))//(fvc::grad(U)&fvc::grad(mu))
== 
fvOptions(rho, U)
);

UEqn.relax();

fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
//    solve(UEqn == -fvc::grad(p)+fvc::reconstruct(f_stf));
	solve
	(
	    UEqn
	    == 
	    fvc::reconstruct
	    (
	   
	      (
		-fvc::snGrad(p)
		+ f_stf
	      )*mesh.magSf()
	    )
	);

    fvOptions.correct(U);
}

pEqn.H

Code:

{
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));

U == rAU*UEqn.H();
U.correctBoundaryConditions();

//surfaceScalarField phiSurf(-rAUf*fvc::interpolate(C)*fvc::snGrad(pot)*mesh.magSf());
surfaceScalarField phiSurf(rAUf*f_stf*mesh.magSf());

surfaceScalarField phiHbyA
(
    "phiHbyA",
    (
        fvc::flux(HbyA)
      + rAUf*fvc::interpolate(rho)*fvc::ddtCorr(U, phi)
    )
    + phiSurf
);

adjustPhi(phiHbyA, U, p);
MRF.makeRelative(phiHbyA);

constrainPressure(p, rho, U, phiHbyA, rAUf, MRF);

// Non-orthogonal pressure corrector loop
while (pimple.correctNonOrthogonal())
{
    // Pressure corrector
    fvScalarMatrix pEqn
    (
////        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
    );
    pEqn.setReference(pRefCell, pRefValue);
    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));

    if (pimple.finalNonOrthogonalIter())
    {
        phi = phiHbyA - pEqn.flux();
        p.relax();
        U = HbyA + rAU*fvc::reconstruct((phiSurf - pEqn.flux())/rAUf);
        U.correctBoundaryConditions();
        fvOptions.correct(U);
    }
}
}

I tested the 2D static bubble case following the validation case of Jamshidi.
The domain is a square and the domain size is (0.001m x 0.001m).
The bubble has a diameter of 0.0005m, and it locates at the center of the domain.
Mesh is generated as (50 x 50), so total 2500 cell points are allocated.
I forgot to mention at the previous post but I initialized order parameter as below.

funkySetFieldsDict

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  1.7.1                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      funkySetFieldsDict;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

expressions
(
	initialC
	{
		field C;
		expression "tanh((sqrt(sqr(pos().x-5.e-4)+sqr(pos().y-5.e-4))-2.5e-4 )/(sqrt(2)*1.e-5))";
		keepPatches 1;
	}

);

Unlike VOF, it has finite interface thickness as initial condition.

All of the boundary conditions are periodic (cyclic).
The other properties are as below.
1. Liquid phase (phase 1)
$\rho_1: 997.0 [kg/m^3]$
$\nu_1: 1.0\times10^{-1} [m^2/s]$
2. Gas phase (phase 2)
$\rho_2: 1.2 [kg/m^3]$
$\nu_2: 1.531\times10^{-5} [m^2/s]$
3. Mobility
$M: 1.0\times10^{-11} [m^3s/kg]$
4. Infterface thickness
$\epsilon: 1.0\times10^{-5} [m]$
5. Surface tension
$\sigma: 7.283\times10^{-2} [kg/s^{-2}]$

Also, I present fvScheme and fvSolution which I used.

fvSchemes

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  2.3.0                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSchemes;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

ddtSchemes
{
    default         backward;//CrankNicolson 1.0;//backward;//CrankNicolson 1.0;
    ddt(C)	bounded Euler;
} 

gradSchemes
{
//	default Gauss linear;
    default         Gauss linear skewCorrected 0.5;//pointCellsLeastSquares;
//    grad(p)	Gauss linear; 
//    grad(C)     fourth;
//    grad(pot)	Gauss linear;
//    grad(U) Gauss cubic;
//    grad(magSqr(f_st)) Gauss linear;
}

divSchemes
{
    default        Gauss linear;
    div(phi,C)   Gauss Gamma 0.5;
    div(rhoPhi,U) Gauss Gamma 0.5;//Gauss Gamma 0.5;
    div(phi,U) Gauss Gamma 0.5;//Gauss Gamma 0.5;
//    div(phi,grad(pot)) Gauss linear;
//    div((mu*dev(T(grad(U))))) Gauss linear;
  
}

laplacianSchemes
{
    default         Gauss linear corrected;
}

interpolationSchemes
{
    default         linear;
}

snGradSchemes
{
    default         corrected;
}

fluxRequired
{
    default         no;
    U;
    p;
    C;
    pcorr;
}


// ************************************************************************* //

fvSolution

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v2012                                 |
|   \\  /    A nd           | Website:  www.openfoam.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     2.0;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvSolution;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

solvers
{
    C
    {

        solver          GAMG;
	nVcycles	10;
	tolerance	1e-12;
	relTol 		0;
	smoother		DILUGaussSeidel;
	nSmoothingSteps	4;
	nPreSweeps	2;
	nPostSweeps	1;
	nFinestSweeps	1;
	cacheAgglomeration 	true;
	nCellsInCoarestLevel	10;
	agglomerator	faceAreaPair;
	mergeLevels	1;

    }
    CFinal
    {
        $C;
        relTol          0;
    }

    p
    {

 	solver	PCG;
	preconditioner
	{
		preconditioner          GAMG;
		nVcycles	10;
		tolerance	1e-10;
		relTol 		0;
		smoother		DICGaussSeidel;
		nSmoothingSteps	4;
		nPreSweeps	0;
		nPostSweeps	2;
		nFinestSweeps	1;
		cacheAgglomeration	true;
		nCellsInCoarestLevel	10;
		agglomerator	faceAreaPair;
		mergeLevels	1;
	}
	tolerance	1.e-12;
	relTol		0;
	minIter		2;
//	maxIter		50;

    }

    pFinal
    {
 	solver	PCG;
	preconditioner
	{
		preconditioner          GAMG;
		nVcycles	10;
		tolerance	1e-10;
		relTol 		0;
		smoother		DICGaussSeidel;
		nSmoothingSteps	4;
		nPreSweeps	0;
		nPostSweeps	2;
		nFinestSweeps	1;
		cacheAgglomeration	true;
		nCellsInCoarestLevel	10;
		agglomerator	faceAreaPair;
		mergeLevels	1;
	}
	tolerance	1.e-12;
	relTol		0;
	minIter		2;
//	maxIter		50;
  
    }

    U
    {
//	type		coupled;
        solver          GAMG;
	nVcycles	10;
	tolerance	1e-12;
	relTol 		0;
	smoother		DILUGaussSeidel;
	nSmoothingSteps	4;
	nPreSweeps	2;
	nPostSweeps	1;
	nFinestSweeps	1;
	cacheAgglomeration	false;
	nCellsInCoarestLevel	100;
	agglomerator	algebraicPair;
	mergeLevels	1;
    }

    UFinal
    {
	$U;
        relTol          0;
    }

}


PIMPLE
{
    correctPhi true; 
    checkMeshCourantNo true;
    momentumPredictor yes;
    nCorrectors     1;
    nOuterCorrectors 1000;
    nNonOrthogonalCorrectors 0;
    pRefCell        0;
    pRefValue       0;
    residualControl
    {
	U
	{
	    tolerance 1e-10;
	    relTol 0;
	}
	p
	{
	    tolerance 1e-12;
	    relTol 0;
	}
	C
	{
	    tolerance 1e-12;
	    relTol 0;
	}

    }
}


relaxationFactors
{
	fields
	{
		p 0.3;
		pFinal 1;
	}
	equations
	{
		"(U|C)"	0.7;
		"(U|C)Final"	1;

	}
}

// ************************************************************************* //

Fortunately, the code did not blow at all and I got the result as attached.
Theoretically, the pressure difference between the outside and the inside of the bubble should be 280 [Pa] so thi is wrong. However, this is due to the number of the grid according to Jamshidi's thesis.
Though it is wrong, the pressure difference is similar to that of Jamshidi. (nearly 50 [Pa])
According to his result about spurious current, the maximum spurious current should be between $1.0\times10^{-6}$ [m] and $1.0\times10^{-5}$ [m]. Actually, his result shows much closer to $1.0\times10^{-6}$ [m].
However, my result is $9.2\times10^{-05}$ [m], which is nearly $1.0\times10^{-04}$ [m].
Also, I found strange behavior of residuals and time step continuity.
After run, I did not modified fvSchemes or fvSolution but some peaky trends occur at residuals during the calculation.
Similarly, global continuity error also showed abrupt increases or decreases during the simulation and, consequently, cumulative continuity error increases up to nearly 0.013.

I think this is not the problem of phase field method itself but the problem of some settings in fvScheme, fvSolution or even at pEqn.H.

If you have any comment, it would be really appreciated.

Best regards,
Jun

mykkujinu2201 · May 12, 2021, 08:26

Dear forum,

I kept monitoring the simulation.
I tested for 100 by 100 case.
The code and the set-up are exactly as same as above except for the grid.

I found that, as time goes, the pressure difference between the inside and the outside of the bubble gets closer to the theoretical value (280 [Pa]) and spurious velocity also gets lower.
Nearly after [MATH] t=0.01, the theoretical Laplace pressure is achieved and the spurious velocity is nearly $2\times 10^{-4} [m/s]$ .
However, after few iterations, pressure equation is not conserved at all, and consequently, the code gave me weird velocity field and pressure field.

I tried several times with different schemes but I got the same result.
It seems that something happens after pressure gets saturated.
I tried to find similar cases in other solvers but I could not find the reason or even a similar case.
If you know any similar error, possible reason, any other opinion or comments.
Please, help me.

Any comment would be appreciated.

Best regards,
Jun

mykkujinu2201 · May 21, 2021, 01:05

Dear Forum,

I keep trying to resolve the issue but I could not resolve it at all.
I first doubted pEqn.H and UEqn.H but now I suspect the part of solving Cahn-Hilliard equation and chemical potential.

Recently, I read Kim's paper (Kim, J. Phase-Field Models for Multi-Component Fluid Flows. Communications in Computational Physics 12, 613–661 (2012)).
Here, it is written that the mobility

can be either constant or variable

.

I tested with using the variable mobility with the definition above.
As a result, the code ran without any problem.
The pressure along the centerline shows pressure jump as similar as theoretical result.
However, the maximum spurious velocity converged to 3.2e-03 [m/s], which I expected as the order of 1.e-05 [m/s].
See the attached pictures.

The mobility affects order parameter and chemical potential, and those two affect capillary force term.
Therefore, now I suspect Cahn-Hilliard equation and chemical potential.
Cahn-Hilliard equation and chemical potential are solved as below.
1. Chemical potential of time

is calculated as
$\phi_n=\lambda\left(\frac{1}{\epsilon^2}C_n(1-C_n^2)-\nabla^2{C_n}\right)$

2. Cahn-Hilliard equation is solved as below.
$\frac{\partial C}{\partial t}+\left(\vec{u}\cdot C\right)=\nabla \cdot (M(1-C_{n+1}^2) \nabla \phi_n)$

$\rightarrow$ $C_{n+1}$

3. Chemical potential of time

is calculated as
$\phi_{n+1}=\lambda\left(\frac{1}{\epsilon^2}C_{n+1}(1-C_{n+1}^2)-\nabla^2{C_{n+1}}\right)$ .

In the code, the process is written as

Code:

pot==lambda/(eps*eps)*C.oldTime()*(C.oldTime()-scalar(1.0))*(C.oldTime()+scalar(1.0))-fvc::laplacian(lambda,C.oldTime());
pot.correctBoundaryConditions();

fvScalarMatrix CEqn
(
    fvm::ddt(C) + fvm::div(phi,C) == fvc::laplacian(M*(1.0-C*C),pot)
);

CEqn.relax();
CEqn.solve();
C.correctBoundaryConditions();

pot==lambda/(eps*eps)*C*(C-scalar(1.0))*(C+scalar(1.0))-fvc::laplacian(lambda,C);
pot.correctBoundaryConditions();

Because the original 4th order is not directly solved in OpenFOAM,
I chose to solve in a segregated manner but now I think I organized the equation wrong.

Please, give me any advice or comment.

Best regards,
Jun

mykkujinu2201 · June 10, 2021, 01:28

Dear forum,

For several weeks, I have tested the case with varying grid, time step size and mobility.

What I found is that grid and time step size have no effect on the result.

At a similar moment, the code blows.

However, from the mobility test, I found that reduced mobility relaxed the sudden increase of spurious velocity and the code did not blow.

I used 0.01 times the previous mobility.

In addition, I have found that the laplace term of chemical potential in Cahn-Hilliard equation is vulnerable to numerical instability.

Is there any method to cure the numerical instability of Cahn-Hilliard equations that is applicable in OpenFOAM?

Please, let me know.

It would be appreciated.

Best regards,
Jun

xiaomk · March 29, 2022, 22:36

Hi, Jun

I am interested in your achievement of phase field method in OpenFOAM, and I have two questions,

1, how is the adherent contact angle achieved in your codes?

For the adherent contact angle, refer to Eq. (2.7) in Cai X. (2016) Interface-Resolving Simulations of Gas-Liquid Two-Phase Flows in Solid Structures of Different Wettability, dissertation.

2, what is the boundary condition of the order parameter C?

tom_flint2012 · December 7, 2022, 12:09

Hi Jun,

I'm fairly experienced with phase field methods for microstructure evolution (Allen cahn),. and I'm just starting to play around with can-hilliard type models like this.

I implemented something that looks almost identical to what you have here. However my phase field goes out of bounds, from -1.1 to + 1.1 either side of the interface, did you ever come across anything like this, and how did you fix it?

I could just clip the field to -1, +1 but I don't think that's the correct thing to do but I'm happy to be told otherwise.

Thanks,

Tom

June 10, 2021, 01:28	Possible reason is found: numerical instability in Cahn-Hilliard equation	#5
mykkujinu2201 Member Jun Join Date: Nov 2015 Posts: 57 Rep Power: 10	Dear forum, For several weeks, I have tested the case with varying grid, time step size and mobility. What I found is that grid and time step size have no effect on the result. At a similar moment, the code blows. However, from the mobility test, I found that reduced mobility relaxed the sudden increase of spurious velocity and the code did not blow. I used 0.01 times the previous mobility. In addition, I have found that the laplace term of chemical potential in Cahn-Hilliard equation is vulnerable to numerical instability. Is there any method to cure the numerical instability of Cahn-Hilliard equations that is applicable in OpenFOAM? Please, let me know. It would be appreciated. Best regards, Jun tom_flint2012, XinXin and arescao like this.

December 7, 2022, 12:09	phase field out of bounds	#7
tom_flint2012 Member Thomas Flint Join Date: Jan 2016 Posts: 60 Rep Power: 10	Hi Jun, I'm fairly experienced with phase field methods for microstructure evolution (Allen cahn),. and I'm just starting to play around with can-hilliard type models like this. I implemented something that looks almost identical to what you have here. However my phase field goes out of bounds, from -1.1 to + 1.1 either side of the interface, did you ever come across anything like this, and how did you fix it? I could just clip the field to -1, +1 but I don't think that's the correct thing to do but I'm happy to be told otherwise. Thanks, Tom

Similar Threads
Thread	Thread Starter	Forum	Replies	Last Post
PEMFC model with FLUENT	brahimchoice	FLUENT	22	April 19, 2020 15:44
Access phase fraction field in a phase change model	hend	OpenFOAM Programming & Development	0	October 17, 2019 05:03
Sharing links for two phase solver packages developed by openfoam community	swap_9068	OpenFOAM Programming & Development	1	April 2, 2017 05:43
fluent divergence for no reason	sufjanst	FLUENT	2	March 23, 2016 16:08
''unknown radialModelType type Gidaspow'' PROBLEM WITH THE BED TUTORIAL	AndoniBM	OpenFOAM Running, Solving & CFD	2	March 25, 2015 18:44

March 29, 2022, 22:36		#6
xiaomk New Member Xiao Mingkun Join Date: Jul 2021 Posts: 3 Rep Power: 4	Hi, Jun I am interested in your achievement of phase field method in OpenFOAM, and I have two questions, 1, how is the adherent contact angle achieved in your codes? For the adherent contact angle, refer to Eq. (2.7) in Cai X. (2016) Interface-Resolving Simulations of Gas-Liquid Two-Phase Flows in Solid Structures of Different Wettability, dissertation. 2, what is the boundary condition of the order parameter C?