Hello everyone,
I’m working on a 2D free surface problem using the ALE method, and I am still writing a minimal working example for this thread. I also apologize for the big thread. Considering a sharp interface approach, I solve the 2D Navier-Stokes equations using a mixed formulation along with 2 boundary conditions for the free surface, namely a Neumann stress boundary and a kinematic condition for the velocity at the interface. My results present spurious behavior for velocity and pressure near the interface, and as I use the surface velocity to displace the interface, an erroneous displacement (interface teeth) occurs. I would like to ask for help of the community regarding the variational formulation and stabilization methods.
Any comments are extremely welcome. So far I’ve tried using smaller time-steps ( O(10^{-11}) s), more refined meshes (they are built using mshr),
The domain represents a free thin liquid sheet with free flow on left/right boundaries, symmetry on bottom surface and free surface as top boundary:
I’m using a mixed formulation with Taylor-Hood elements. For a given weight function v \in V in which V is a vector fuction space defined on the mesh, weak form of the momentum equation (after integration by parts of Cauchy stress tensor) is given as:
\int_{\Omega}{\rho(\frac{\partial u}{\partial t} + u \bullet \nabla u)v dx} = - \int_{\Omega}{T : \nabla v dx} + \int_{\Gamma}{(\hat n \bullet T)v ds}
in which the interface traction condition is as follow:
\hat n \bullet T = \nabla_s \bullet \sigma_s -(P_S +\Phi) \hat n
which is weakly imposed as
\int_{\Gamma}{(\hat n \bullet T)v} ds = \int_{\Gamma}{(\nabla_s \bullet \sigma_s )v ds} - \int_{\Gamma}{((P_S+\Phi) \hat n) v ds}
for \sigma_s being the interface stress tensor. The kinematic condition is
\hat n \bullet \vec u = \hat n \bullet \vec u_s
The passive case (after order reduction and simplification) provides:
\nabla_s \bullet \sigma_s = \sigma_{\alpha\beta} \kappa \hat n
In which \hat n is the surface outwards directed normal vector, \sigma_{\alpha \beta} is a constant interfacial tension and (P_S +\Phi)\hat n is a disjoining pressure term. \kappa is the boundary mean curvature, defined as:
\kappa = - \frac{\nabla \bullet \hat n}{|| \hat n ||}
In FEniCs:
def divK(nN):
# nN is the normal vector evaluated at the mesh dofs
return inner(Identity(2), as_tensor([[nN[0].dx(0), 0],
[0, nN[1].dx(1)]]))
def curvature(mesh, ds, marker):
# Get normal vector defined on dofs
V = VectorFunctionSpace(mesh, "CG", 2)
# Projection of the normal vector on P2 space
u = TrialFunction(V)
v = TestFunction(V)
l = TestFunction(C)
n = FacetNormal(mesh)
a = inner(u, v)*ds(marker)
L = inner(n, v)*ds(marker)
# Solve system
A = assemble(a, keep_diagonal=True)
b = assemble(L)
A.ident_zeros()
nn = Function(V)
solve(A, nn.vector(), b) # Obtain normal vector evaluated at the dofs
kappa = - divK(nn/sqrt(dot(nn,nn)))
return kappa, nn
It is also possible to define \kappa in terms of the tangent vector:
\kappa = \frac{\partial \hat t}{\partial s}
The tangent evaluated on the dofs is obtained by rotating the normal vector evaluated on the dofs:
import ufl
# nn is the normal vector evaluated on boundary dofs
t = ufl.perp(nn)
In this case the interface traction for the passive case is defined as
\nabla \bullet \sigma_s = \sigma_{\alpha \beta} \kappa \hat n = \sigma_{\alpha\beta}\frac{\partial \hat t}{\partial s}
An inspection of normal (yellow) and tangent (green) vectors near x = 0 yields
The non-zero (blue) values below the interface belong to the tangent vector. Maybe this could lead to the occurrence of spurious currents (??).
The implementation in FEniCs is as below.
def DD(u):
return sym(nabla_grad(u))
def TT(u, p, mu, I):
return 2*mu*DD(u) - p*I
def IST_curv(sigma, kappa, n):
return sigma*kappa*n
def IST_tang(sigma, kappa, n):
return sigma*dot(grad(t), t)
# Momentum conservation
a1 = rho*dot((u-u0)/k,v)*dx() + \
( rho*dot(dot(u ,grad(u) ),v) + \
inner(TT(u,p,mu, I),DD(v)) )*dx()
# Boundary integral -
b_int = inner(IST_curv(sigma, kappa, n_surf), v)*ds(surface_marker) - \
VdW(A, x)*dot(n, v)*ds(surface_marker)
# Continuity equation
a2 = (q*div(u))*dx()
L2 = 0
# Residual
F = (a1 + a2) - (b_int + L2)
Then I solve the problem with
J = derivative(F,w,dw)
problem = NonlinearVariationalProblem(F,w,bcu,J)
solver = NonlinearVariationalSolver(problemU)
# Configure the solver parameters (abs tolerance, rel tolerance, max iter...)
(no_iterations,converged) = solver.solve()
The solution is presented below is obtained using the curvature approach for the interface stress. However, the results are identical using the tangent approach, apart from the occurrence of more spurious currents near x = 0. The initial condition for velocity is zero. The disjoining pressure is inversely proportional to the cube of the interface height, so we expect a greater pressure near the minimal height of the surface that . At the first iteration, the velocity and pressure fields are as expected, but on the next steps, the occurrence of spurious currents near the interface leads the simulation to end.
Pressure at first iteration:
X component of velocity at first iteration:
Y component of velocity at first iteration:
A zoom near the cuved region around x = ± 0.00005:
The pressure around x = ± 0.00005 on the second iteration (in paraview):
Velocity around x = ± 0.00005 on the second iteration (in paraview):
Any thoughts are extremely welcome.