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

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

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

Like Tree8Likes
  • 1 Post By mykkujinu2201
  • 1 Post By mykkujinu2201
  • 1 Post By mykkujinu2201
  • 2 Post By mykkujinu2201
  • 3 Post By mykkujinu2201

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   May 7, 2021, 02:55
Default Weird result from 'phase field solver'. Need help.
  #1
Member
 
Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 10
mykkujinu2201 is on a distinguished road
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, C are chemical potential, mixing energy density scale, interface thickness and order parameter, respectively.
Order parameter indicates each phase and interface. C=1 \rightarrow phase 1, C=-1 \rightarrow phase 2, -1<C<1 \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},M 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, p, \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
arescao likes this.
mykkujinu2201 is offline   Reply With Quote

Old   May 12, 2021, 00:12
Default Additional information
  #2
Member
 
Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 10
mykkujinu2201 is on a distinguished road
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
Attached Images
File Type: png PressureAlongCenterLine.png (25.8 KB, 29 views)
File Type: png Umagnitude.png (37.3 KB, 34 views)
File Type: png Residual_continuity.png (93.5 KB, 21 views)
arescao likes this.
mykkujinu2201 is offline   Reply With Quote

Old   May 12, 2021, 09:26
Default I found where the problem occurs during calculation but not exactly.
  #3
Member
 
Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 10
mykkujinu2201 is on a distinguished road
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
Attached Images
File Type: png PressureAlongCenterLine_100.png (38.3 KB, 8 views)
File Type: png PressureAlongCenterLine_100_blow.png (45.0 KB, 7 views)
File Type: png Umagnitude_100.png (62.6 KB, 27 views)
File Type: png Pressure_100_blow.png (30.0 KB, 14 views)
arescao likes this.
mykkujinu2201 is offline   Reply With Quote

Old   May 21, 2021, 02:05
Default Is the problem due to the part of solving chemical potential and Cahn-Hilliard eq'n?
  #4
Member
 
Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 10
mykkujinu2201 is on a distinguished road
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 M can be either constant or variable (1-C^2).

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 n 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 n+1 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
Attached Images
File Type: png residuals.png (122.6 KB, 11 views)
File Type: png velocityField.png (32.1 KB, 23 views)
File Type: png pressure.png (32.5 KB, 11 views)
acgnipper and arescao like this.
mykkujinu2201 is offline   Reply With Quote

Old   June 10, 2021, 02:28
Default Possible reason is found: numerical instability in Cahn-Hilliard equation
  #5
Member
 
Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 10
mykkujinu2201 is on a distinguished road
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.
mykkujinu2201 is offline   Reply With Quote

Old   March 29, 2022, 23:36
Default
  #6
New Member
 
Xiao Mingkun
Join Date: Jul 2021
Posts: 3
Rep Power: 5
xiaomk is on a distinguished road
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?
xiaomk is offline   Reply With Quote

Old   December 7, 2022, 13:09
Default phase field out of bounds
  #7
Member
 
Thomas Flint
Join Date: Jan 2016
Posts: 60
Rep Power: 10
tom_flint2012 is on a distinguished road
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
tom_flint2012 is offline   Reply With Quote

Reply

Tags
phase field method, pimplefoam, surface tension, two-phase

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
PEMFC model with FLUENT brahimchoice FLUENT 22 April 19, 2020 16:44
Access phase fraction field in a phase change model hend OpenFOAM Programming & Development 0 October 17, 2019 06:03
Sharing links for two phase solver packages developed by openfoam community swap_9068 OpenFOAM Programming & Development 1 April 2, 2017 06:43
fluent divergence for no reason sufjanst FLUENT 2 March 23, 2016 17:08
''unknown radialModelType type Gidaspow'' PROBLEM WITH THE BED TUTORIAL AndoniBM OpenFOAM Running, Solving & CFD 2 March 25, 2015 19:44


All times are GMT -4. The time now is 06:05.