Background
I would like to numerically solve the stationary Stokes equation from fluid dynamics
$\eta_i \nabla^2 \vec{u}_i - \nabla p = 0$
with the incompressibility condition
$\nabla \cdot \vec{u} = 0$
for the flow field in two immiscible liquid layers with a spatially varying surface tension $\sigma(x)$ between them. In what follows I am assuming that the surface tension has a constant gradient: $\partial_x \sigma = 1$. This may be due to a temperature or concentration gradient at the interface, which does not get modelled explicitly for the sake of simplicity.
This should lead to waves that propagate transversal to the interface:
For more details, see e.g. Fig. 9.3a in https://ocw.mit.edu/courses/mathematics/18-357-interfacial-phenomena-fall-2010/lecture-notes/MIT18_357F10_Lecture9.pdf
Problem Outline
The spatial layout of the problem is depicted as follows:
The integration domain is a two-dimensional box, which liquid 1 on top of liquid 2. Each have their respective viscosities $\eta_i$ and flow fields $\vec{u}_i$. At the top and bottom of the box are no-slip boundaries, which imply Dirichlet boundary conditions for the flow fields. At the interface the vertical components of each flow field are zero to enforce that both liquids can not mix. The flow tangential to the interface has the same magnitude due to continuity.
The surface tension at the interface adds another internal boundary condition that seems complicated at first sight:
$P_s \nabla \sigma = P_s (T_2-T_1) \vec{n}$.
Here $\vec{n}$ is the vector normal to the interface, $P_s$ projects a vector $\vec{v}$ on the interface (by removing its component normal to the interface) and the $T_i = \eta_i (\nabla\vec{u} + (\nabla\vec{u})^T)$ are strain-rate tensors. For the simple planar problem at hand, this allows us to reduce the surface tension condition to:
$\partial_x \sigma = \eta_2 \partial_y u_{x,2} - \eta_1 \partial_y u_{x,1}$
The analytical solution to this problem is:
$\vec{u}_i(x,y) = \frac{1 \pm y}{\eta_1+\eta_2}\partial_x \sigma \vec{e}_x$
Code
To verify this solution numerically, I wrote the following Mathematica code:
(** create mesh **)
boundaryMesh = ToBoundaryMesh[
"Coordinates" -> {{0, -1}, {0, 0}, {0, 1}, {2, 1}, {2, 0}, {2, -1}},
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,
5}, {5, 6}, {6, 1}, {2, 5}}]}
];
p1 = boundaryMesh["Wireframe"];
mesh = ToElementMesh[boundaryMesh, MaxCellMeasure -> 0.01];
p2 = mesh["Wireframe"];
Grid[{{p1, p2}}]
This gives me the following mesh with the required internal boundary:
The PDE and boundary conditions are defined as follows. Note that I am not defining separate flow fields for the top and bottom, but instead I tried to throw everything into one flow field with hoizontal and vertical components u, v.
\[Eta]=With[{\[Eta]1=0.5,\[Eta]2=1},If[y<0,\[Eta]1,\[Eta]2]];
(** u,v,p = horizontal speed, vertical speed, pressure **)
pde={
Inactive[Div][\[Eta] Inactive[Grad][u[x, y], {x, y}], {x, y}] + (p^(1,0))[x, y]==0,
Inactive[Div][\[Eta] Inactive[Grad][v[x, y], {x, y}], {x, y}] + (p^(0,1))[x, y]==0,
(u^(1,0))[x, y]+(v^(0,1))[x, y]==0
};
boundaryConditions = {
DirichletCondition[{u[x, y] == 0, v[x, y] == 0, p[x, y] == 0}, (y == -1 \[Or] y == 1) \[And] (0 < x < 2)],
DirichletCondition[v[x, y] == 0, y == 0]
};
(** solve flow velocities with a method that is an order higher than pressure for stable solutions **)
{xVel, yVel, pressure} = NDSolveValue[
{pde, boundaryConditions},
{u, v, p},
{x, y} \[Element] mesh,
Method -> {"FiniteElement",
"InterpolationOrder" -> {u -> 2, v -> 2, p -> 1}}
];
plot=VectorPlot[{xVel[x, y], yVel[x, y]}, {x, y} \[Element] mesh,
AspectRatio -> Automatic, StreamPoints -> 6,
StreamColorFunction -> "TemperatureMap",
StreamColorFunctionScaling -> False]
The obviously unsatisfying result is:
Questions
- The solution to my code is pretty boring (no flow), because I did not include the internal boundary condition for the surface tension:
$\partial_x \sigma = \eta_2 \partial_y u_{x,2} - \eta_1 \partial_y u_{x,1}$
with
$\partial_x \sigma = 1$.
How can I do that? It is not clear to me how to enforce an internal jump in the derivative by a prescribed amount using NeumannValue[]. Do I need to run 2 NDSolve's for each region in tandem as described here: FEM-ception: Using output of NDSolve FEM as input in another FEM problem ? And then have the solution in one domain as a boundary condition for the simulation in the other? Also note that this is more complicated than a simple jump in a spatially-dependent parameter, as solved here: Neumann boundary condition on a boundary inside the region. This problem involves a jump in a derivative by a specified amount.
- How can I set periodic boundary conditions? I tried
PeriodicBoundaryCondition[u[x, y], x == 0, TranslationTransform[{2, 0}]],PeriodicBoundaryCondition[v[x, y], x == 0, TranslationTransform[{2, 0}]]but then I got this error:NDSolveValue::femnodpbc: DirichletCondition can not be present on the target boundary of a PeriodicBoundaryConditon.
This problem is a first test for me of Mathematica's FEM capabilities, I later intend to apply the discontinuous surface tension boundary condition to a sphere geometry in 3D.
