sandbox/zaleski/vdw.h
Equations and code
This solver implements a scheme for a diffuse interface model. The intent is to use it for the meso-scale modelling of various physical problems involving interfaces, such as contact lines, thin sheet break up etc.
Basics
The Van der Waals phase field equation can be written as in Nadiga, B. T., & Zaleski, S. (1995). Investigations of a two-phase fluid model. arXiv preprint comp-gas/9511003. hereafter NZ95. The model reads \displaystyle \partial_t \rho + \partial_i(\rho u_i) = 0 and \displaystyle \partial_t (\rho u_i ) + \partial_j ( \rho u_i u_j ) = \partial_j \sigma_{ij} with \rho the density, \mathbf{u} = (u_i)_i the velocity. We will note \mathbf{q} = \rho \mathbf{u} the momentum. The solver will start with the basic basilisk framework
#include "run.h"
#include "timestep.h"The variables are declared as basilisk fields
scalar rho[];
vector q[],u[];We declare a function that will give the pressure, the user should define this function in his code
double P0( double x);\sigma_{ij} the stress tensor given by \displaystyle \sigma_{ij} = - p_0 \delta_{ij} + \sigma_{ij}^{(v)} + \sigma_{ij}^{(K)} where \sigma_{ij}^{(v)} has the Newtonian form with vanishing compression viscosity \displaystyle \sigma_{ij}^{(v)} = \mu ( \partial_i u_j + \partial_j u_i - (2/3) (\nabla \cdot \mathbf{u})\delta_{ij} ) where \mu is the viscosity. We shall thus need
scalar mu[];and \displaystyle \sigma_{ij}^{(K)} = \lambda \left[ \left( \frac 12 \vert \nabla \rho \vert^2 + \rho \nabla^2 \rho \right) \delta_{ij} - \partial_i \rho \partial_j \rho \right] is the Korteweg stress tensor, which can be rewritten \displaystyle \sigma_{ij}^{(K)} = \lambda (\rho \nabla^2 \rho \, \delta_{ij} - T_{ij} ) (notice the sign difference with NZ95) where the tensor T is expressed in 2D as \displaystyle T = \left( \begin{array}{cc} (\partial_x \rho)^2/2 - (\partial_y \rho)^2/2 & \partial_x \rho \, \partial_y \rho \\ \partial_x \rho \, \partial_y \rho & - (\partial_x \rho)^2/2 + (\partial_y \rho)^2/2 \end{array} \right) thus we define
(const) double lambda;Mac Cormack scheme
We shall use the predictor corrector scheme in NZ95 written as \displaystyle \partial_t \mathbf{f} + \partial_x\mathbf {F}_x[\mathbf{f}] + \partial_y\mathbf {F}_y[\mathbf {f}] = 0, where \displaystyle \mathbf{f} = \left(\begin{array}{c} \rho \\ \rho u_x \\ \rho u_y \end{array}\right) and \mathbf{F}_i has functional dependence on \mathbf{f} through \displaystyle \mathbf{F}_x[\mathbf{f}] = \mathbf{F}_x(\mathbf{f}, \partial_x \mathbf{f}, \partial_y \mathbf{f}, \nabla^2 \rho) (minor typo in NZ95: \nabla \rightarrow \nabla^2) and \displaystyle \mathbf{F}_y[\mathbf{f}] = \mathbf{F}_y(\mathbf{f}, \partial_y \mathbf{f}, \partial_x \mathbf{f}, \nabla^2 \rho). One has \displaystyle \mathbf{F}_x = \left(\begin{array}{c} \rho u_x\\ \rho u_x^2 - \sigma_{xx} \\ \rho u_x u_y - \sigma_{xy} \end{array}\right) and \mathbf{F}_y is obtained analogously by exchanging the 2nd and 3rd components and exchanging x and y. (We expect basilisk to perform this exchange automatically). Explicitly, \displaystyle \mathbf{F}_x = \left(\begin{array}{c} \rho u_x\\ \rho u_x^2 + p_0(\rho) \\ \rho u_x u_y \end{array} \right) - \mu \left(\begin{array}{c} 0\\ 2 \partial_x u_x - (2/3)(\partial_x u_x + \partial_y u_y) \\ \partial_x u_y + \partial_y u_x \end{array} \right) + \lambda \left(\begin{array}{c} 0\\ T_{xx} - \rho \nabla^2 \rho \\ T_{xy} \end{array} \right). To discretize the equation we define a time step \tau, a spatial step h and forward, backwards and centered difference operators \displaystyle \partial^+_x v = \frac{v(x+h) - v(x)}h \displaystyle \partial^-_x v = \frac{v(x) - v(x-h)}h \displaystyle \partial^c_x v = \frac{v(x+h) - v(x-h)}{2h} The basilisk expressions are obvious. In addition we will need the centered Laplacian \nabla^{2,c} = \partial_x^+\partial_x^- + \partial_y^+\partial_y^- that is defined and written as in here. The predictor step is (see NZ95)
\displaystyle \mathbf{f}^* = \mathbf{f}^n - \tau \partial_x^- \mathbf{F}_x(\mathbf{f}, \partial_x^+ \mathbf{f}, \partial_y^c \mathbf{f}, \nabla^{2,c} \rho) - \tau \partial_y^- \mathbf{F}_y(\mathbf{f}, \partial_y^+ \mathbf{f}, \partial_x^c \mathbf{f}, \nabla^{2,c} \rho),
where we have written \mathbf{f} for \mathbf{f}^n for simplicity and the corrector step is
\displaystyle \mathbf{f}^{n+1} = \frac12 (\mathbf{f}^n + \mathbf{f}^*) - \frac12\tau \partial_x^+ \mathbf{F}_x(\mathbf{f^*}, \partial_x^- \mathbf{f^*}, \partial_y^c \mathbf{f^*}, \nabla^{2,c} \rho^*) - \frac12\tau \partial_y^+ \mathbf{F}_y(\mathbf{f^*}, \partial_y^- \mathbf{f^*}, \partial_x^c \mathbf{f^*}, \nabla^{2,c} \rho^*).
Initial conditions
What follows is shamelessly pumped from the centered Navier-Stokes solver
After user initialisation, we initialise the fluid properties.
We update fluid properties.
event ("properties");We set the initial timestep (this is useful only when restoring from a previous run).
dtmax = DT;
event ("stability");
}Time integration
The timestep for this iteration is controlled by the CFL condition, applied to the centered velocity field \mathbf{u}; and the timing of upcoming events.
event set_dtmax (i++,last) dtmax = DT;
event stability (i++,last) {
foreach()
foreach_dimension()
u.x[] = q.x[]/rho[];
dt = dtnext (timestep (u, dtmax));
}Mac Cormack predictor_corrector event
event predictor_corrector(i++,last)
{Laplacian
We take the Laplacian of \rho without the 1/\Delta^2 factor, the division by the cell size \Delta will be added below. We declare it as an automatic variable:
scalar laprho[],rhop[];
vector qp[];
trash({rhop});
foreach()
{
laprho[]= rho[1,0] + rho[-1,0] + rho[0,1] + rho[0,-1] - 4*rho[];Predictor step
We compute the predictor for the density: \displaystyle \rho^* = \rho - \tau [ \partial_x^- (\rho u_x) + \partial_y^- ( \rho u_y )]
rhop[] = rho[] - (dt/Delta)*(q.x[] - q.x[-1,0] + q.y[] - q.y[0,-1]);
foreach_dimension()
u.x[] = q.x[]/rho[];
}Predictor for velocity \displaystyle q_x^* = q_x - \tau \{ \partial_x^- (\rho u_x^2 + p_0 ) - (2/3) \partial_x^- [\mu( 2 \partial_x^+ u_x - \partial_y^c u_y )] - (\lambda/2) \partial_x^- [ 2\rho \nabla^{2,c} \rho - (\partial_x^+ \rho)^2 + (\partial_y^c \rho)^2 ] \displaystyle + \partial_y^- (\rho u_x u_y) - \partial_y^- [\mu( \partial_y^+ u_x - \partial_x^c u_y )] + \lambda \partial_y^- (\partial_y^+ \rho \, \partial_x^c \rho) \}
which after minor manipulation is \displaystyle u_x^* = u_x - \tau \{ \partial_x^- (\rho u_x^2 + p_0 ) + \partial_y^- (\rho u_x u_y) \displaystyle - (2/3) \partial_x^- [\mu( 2 \partial_x^+ u_x - \partial_y^c u_y )] - \partial_y^- [\mu( \partial_y^+ u_x + \partial_x^c u_y )] This is a code (not compiled) for the last line
#if naught
- (2/3)*( mu[ 0,0] * (2*(u.x[1,0]-u.x[ 0,0])/Delta - 0.5*(u.y[ 0,1]-u.y[ 0,-1]))/Delta
- mu[-1,0] * (2*(u.x[0,0]-u.x[-1,0])/Delta - 0.5*(u.y[-1,1]-u.y[-1,-1]))/Delta
)/Delta
- ( mu[ 0,0] * ( (u.x[0,1] - u.x[0, 0])/Delta + 0.5*(u.y[1, 0] - u.y[-1, 0])/Delta )
- mu[0,-1] * ( (u.x[0,0] - u.x[0,-1])/Delta + 0.5*(u.y[1,-1] - u.y[-1,-1])/Delta )
)/Delta
#endif\displaystyle + (\lambda/2) [ \partial_x^- ( (\partial_x^+ \rho)^2 - (\partial_y^c \rho)^2 - 2\rho \nabla^{2,c} \rho ) + 2 \partial_y^- (\partial_y^+ \rho \partial_x^c \rho)] \} Here is again a code (not compiled) for the above line
#if naught
+ 0.5*lambda*( sq((rho[1,0]-rho[ 0,0])/Delta) - 0.25*sq((rho[ 0,1]-rho[0 ,-1])/Delta) - 2*rho[ 0,0]*laprho[ 0,0]
-sq((rho[0,0]-rho[-1,0])/Delta) + 0.25*sq((rho[-1,1]-rho[-1,-1])/Delta) + 2*rho[-1,0]*laprho[-1,0]
+ 2*( ((rho[0,1]-rho[ 0,0])/Delta)*(rho[1, 0]-rho[-1, 0])/(2.*Delta)
-((rho[0,0]-rho[0,-1])/Delta)*(rho[1,-1]-rho[-1,-1])/(2.*Delta) )/Delta
)
#endifAfter some simplifications here is the actual code
trash({qp});
foreach()
foreach_dimension()
qp.x[] = q.x[] - (dt/Delta)*
( rho[ 0,0]*sq(u.x[ 0,0]) + P0(rho[ 0,0])
- rho[-1,0]*sq(u.x[-1,0]) - P0(rho[-1,0])
+ rho[]*u.x[]*u.y[] - rho[0,-1]*u.x[0,-1]*u.y[0,-1]
- ( (2./3.)*( mu[] * (2*(u.x[1,0] - u.x[] ) - 0.5*(u.y[0,1] + u.y[0,-1]) )
-mu[-1,0] * (2*(u.x[] - u.x[-1,0]) - 0.5*(u.y[-1,1] + u.y[-1,-1]) ) )
- mu[] * ( u.x[0,1] - u.x[] + 0.5*(u.y[1,0] - u.y[-1,0] ) )
+ mu[0,-1] * ( u.x[0,0] - u.x[0,-1] + 0.5*(u.y[1,-1] - u.y[-1,-1]) )
- 0.5*lambda*( sq(rho[1,0]-rho[] ) - 0.25*sq(rho[0,1] -rho[0,-1] ) - 2*rho[ 0,0]*laprho[]
-sq(rho[0,0]-rho[-1,0]) + 0.25*sq(rho[-1,1]-rho[-1,-1]) + 2*rho[-1,0]*laprho[-1,0]
+ (rho[0,1]-rho[0, 0]) * (rho[1, 0]-rho[-1, 0])
- (rho[0,0]-rho[0,-1]) * (rho[1,-1]-rho[-1,-1])
)
)/Delta
);Corrector step
Predicted Laplacian
foreach()
{
laprho[]= rhop[1,0] + rhop[-1,0] + rhop[0,1] + rhop[0,-1] - 4*rhop[];Corrector for density \displaystyle \rho^{n+1} = \frac 12 [ \rho^n + \rho^* ] - \frac \tau 2 [ \partial_x^+ (\rho^* u^*_x) + \partial_y^+ ( \rho^* u^*_y )]
rho[] = 0.5*(rhop[]+rho[]) - 0.5*(dt/Delta)*(qp.x[1,0] - qp.x[] + qp.y[0,1] - qp.y[]);
foreach_dimension()
u.x[] = qp.x[]/rhop[];
}Corrector for velocity \displaystyle q_x^{n+1} = \frac 12 (q_x^* + q_x^n ) - \frac \tau 2 \{ \partial_x^+ (\rho u_x^{* 2} + p_0^* ) + \partial_y^+ (\rho u_x^* u_y^*) \displaystyle - (2/3) \partial_x^+ [\mu( 2 \partial_x^- u_x^* - \partial_y^c u_y^* )] - \partial_y^+ [\mu( \partial_y^- u_x^* + \partial_x^c u_y^* )] \displaystyle + (\lambda/2) [ \partial_x^+ ( (\partial_x^- \rho^*)^2 - (\partial_y^c \rho^*)^2 - 2\rho^* \nabla^{2,c} \rho^* ) + 2 \partial_y^+ (\partial_y^- \rho \partial_x^c \rho^*)] \} We transcribe this, noticing that for any field f \displaystyle \partial_x^+ \partial_x^- f = \partial_x^- \partial_x^+ f and \displaystyle \partial_x^+ (\partial_x^- f)^2 = \partial_x^- (\partial_x^+ f)^2 .
foreach()
foreach_dimension()
q.x[] = 0.5*(qp.x[] + q.x[]) - 0.5*(dt/Delta)*
( rhop[1,0]*sq(u.x[1,0]) + P0(rhop[1,0])
- rhop[0,0]*sq(u.x[0,0]) - P0(rhop[0,0])
+ rhop[0,1]*u.x[0,1]*u.y[0,1]
- rhop[0,0]*u.x[0,0]*u.y[0,0]
- ( (2./3.)*( mu[1,0] * (2*(u.x[1,0] - u.x[ 0,0]) - 0.5*(u.y[1,1] + u.y[1,-1]) )
-mu[0,0] * (2*(u.x[0,0] - u.x[-1,0]) - 0.5*(u.y[0,1] + u.y[0,-1]) ) )
- mu[0,1] * ( u.x[0,1] - u.x[0, 0] + 0.5*(u.y[1,1] - u.y[-1,1] ) )
+ mu[0,0] * ( u.x[0,0] - u.x[0,-1] + 0.5*(u.y[1,0] - u.y[-1,0]) )
- 0.5*lambda*( sq(rhop[1,0]-rhop[] ) - 0.25*sq(rhop[0,1] -rhop[0,-1] ) - 2*rhop[1,0]*laprho[1,0]
-sq(rhop[0,0]-rhop[-1,0]) + 0.25*sq(rhop[-1,1]-rhop[-1,-1]) + 2*rhop[0,0]*laprho[0,0]
+ (rhop[0,1]-rhop[] ) * (rhop[1,0]-rhop[-1,0]) - (rhop[0,0]-rhop[0,-1])*(rhop[1,-1]-rhop[-1,-1])
)
)/Delta
);
}To do
- Work out the surface tension as a function of the equation of state.
- Create examples and tests as in NZ95.
