Weird result from 'phase field solver'. Need help.
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:
/*---------------------------------------------------------------------------*\ Code:
//Cahn-Hilliard Eq'n solving Code:
// properties update Code:
volScalarField rAU(1.0/UEqn.A()); transportProperties Code:
/*--------------------------------*- C++ -*----------------------------------*\ 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 |
I found where the problem occurs during calculation but not exactly.
4 Attachment(s)
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 |
Is the problem due to the part of solving chemical potential and Cahn-Hilliard eq'n?
3 Attachment(s)
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()); 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 |
Possible reason is found: numerical instability in Cahn-Hilliard equation
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 |
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? |
phase field out of bounds
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 |
All times are GMT -4. The time now is 22:07. |