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

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

 May 7, 2021, 01:55 Weird result from 'phase field solver'. Need help. #1 Member   Jun Join Date: Nov 2015 Posts: 57 Rep Power: 8 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 , where , , , are chemical potential, mixing energy density scale, interface thickness and order parameter, respectively. Order parameter indicates each phase and interface. phase 1, phase 2, interface. 2) Cahn-Hilliard equation , where , are velocity and mobility, respectively. 3) Continuity equation 4) Momentum equation , where , , , are density, pressure, viscosity, gravitational acceleration and surface tension force, respectively. Viscosity and density are calculated using the equations below. , where each subscript indicates each phase. Mixing energy density scale, is related to surface tension as , where is surface tension. Surface tension force, is calculated as . 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 . 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 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.

May 11, 2021, 23:12
#2
Member

Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 8
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,
can be rewritten as
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);

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

fvVectorMatrix UEqn
(
fvm::ddt(rho,U) + fvm::div(rhoPhi,U)
+ MRF.DDt(rho,U)
==
fvOptions(rho, U)
);

UEqn.relax();

fvOptions.constrain(UEqn);
if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(

(
+ 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*f_stf*mesh.magSf());

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

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)

2. Gas phase (phase 2)

3. Mobility

4. Infterface thickness

5. Surface tension

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;
}

{
//	default Gauss linear;
default         Gauss linear skewCorrected 0.5;//pointCellsLeastSquares;
}

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;

}

laplacianSchemes
{
default         Gauss linear corrected;
}

interpolationSchemes
{
default         linear;
}

{
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 [m] and [m]. Actually, his result shows much closer to [m].
However, my result is [m], which is nearly [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
 PressureAlongCenterLine.png (25.8 KB, 17 views) Umagnitude.png (37.3 KB, 18 views) Residual_continuity.png (93.5 KB, 15 views)

May 12, 2021, 08:26
I found where the problem occurs during calculation but not exactly.
#3
Member

Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 8
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 .
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.

Any comment would be appreciated.

Best regards,
Jun
Attached Images
 PressureAlongCenterLine_100.png (38.3 KB, 5 views) PressureAlongCenterLine_100_blow.png (45.0 KB, 5 views) Umagnitude_100.png (62.6 KB, 12 views) Pressure_100_blow.png (30.0 KB, 8 views)

May 21, 2021, 01:05
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: 8
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

2. Cahn-Hilliard equation is solved as below.

3. Chemical potential of time is calculated as
.

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.

Best regards,
Jun
Attached Images
 residuals.png (122.6 KB, 8 views) velocityField.png (32.1 KB, 10 views) pressure.png (32.5 KB, 7 views)

 June 10, 2021, 01:28 Possible reason is found: numerical instability in Cahn-Hilliard equation #5 Member   Jun Join Date: Nov 2015 Posts: 57 Rep Power: 8 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.

 March 29, 2022, 22:36 #6 New Member   Xiao Mingkun Join Date: Jul 2021 Posts: 3 Rep Power: 3 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?

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