|
[Sponsors] |
Weird result from 'phase field solver'. Need help. |
|
LinkBack | Thread Tools | Search this Thread | Display Modes |
May 7, 2021, 02:55 |
Weird result from 'phase field solver'. Need help.
|
#1 |
Member
Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 10 |
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 <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; } // ************************************************************************* // 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(); 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); } 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); 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; 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 |
|
May 12, 2021, 00:12 |
Additional information
|
#2 |
Member
Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 10 |
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(); 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); } 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); } } } 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; } ); 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; } 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; } // ************************************************************************* // 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 |
|
May 12, 2021, 09:26 |
I found where the problem occurs during calculation but not exactly.
|
#3 |
Member
Jun
Join Date: Nov 2015
Posts: 57
Rep Power: 10 |
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. Please, help me. Any comment would be appreciated. Best regards, Jun |
|
May 21, 2021, 02: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: 10 |
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(); 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 |
|
June 10, 2021, 02:28 |
Possible reason is found: numerical instability in Cahn-Hilliard equation
|
#5 |
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 |
|
March 29, 2022, 23:36 |
|
#6 |
New Member
Xiao Mingkun
Join Date: Jul 2021
Posts: 3
Rep Power: 5 |
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? |
|
December 7, 2022, 13:09 |
phase field out of bounds
|
#7 |
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 |
|
Tags |
phase field method, pimplefoam, surface tension, two-phase |
Thread Tools | Search this Thread |
Display Modes | |
|
|
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 |