Hello everyone,
I’m working on a free surface flow problem, and this problem requires that surface stress is balanced with bulk stress, as in
In which T is the Cauchy stress tensor, K is the surface mean curvature and n is the outwards directed normal vector. I impose the stress boundary condition naturally by applying integration by parts on the stress term of the weak form. The weak form is like this:
# rate-of-deformation tensor
def DD(u):
return sym(nabla_grad(u))
# Cauchy stress tensor
def TT(u, p):
return 2*mu*DD(u) - p*I
# Interface stress tensor
def IST(sigma, kappa, n):
return sigma*kappa*n
F1 = rho*dot((u - u_n) / k, v)*dx \
+ rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(TT(U, p_n), DD(v))*dx \
- dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx
# Boundary integral term - free surface is marked with flag #1 -> ds(1)
# stress balance on interface # pressure term
b_int = inner(IST(sigma, kappa, n), v)*ds(1) - dot(p_n*n, v)*ds
F1 = F1 - b_int
a1 = lhs(F1)
L1 = rhs(F1)
I’ll be updating the surface position using the ALE method, in which the free surface is one of the boundaries of my domain (a square domain). The domain is as follows.
As you can see, the top boundary is defined as the free surface. At last, to make this application work, I need to calculate the curvature of the free surface boundary. And here is where the problem begins: I’m using the solution for this topic to calculate the curvature of all the boundaries, as I thought the approach would work for my case. But it’s not. I obtain a proper value for the curvature of a full closed circle (2pi), as used on the real example. But when I try to implement the code into my application, when I calculate the mean curvature of the surface bundary with
kappa, n = curvature(mesh, bmesh) # based on the topic mentioned
print(assemble(kappa*ds(1)))
mean curvature is zero. I don’t know if this is the right way to calculate the mean curvature of exernal boundaries like in my case. I was counting on the fact that the curvature could be obtained by
kappa = - div(n)
The code below is not a MWE, as I dont have any. Please, share some light.
from dolfin import * # FEniCs library
from mshr import *
import numpy as np
import matplotlib.pyplot as plt
from utils import curvature # third-party library based on the post described previously
def p_pointwise(FS, bflag, bmesh, bmf):
'''
FS -> VectorFunctionSpace
bflag = flag index related to the boundary
bmesh -> boundary mesh
bmf -> boundary marker mesh function
'''
topBC = []
for i in range(bmesh.num_vertices()):
if bmf[i] == bflag:
xbc = bmesh.coordinates()[i,0]
ybc = bmesh.coordinates()[i,1]
# Pointwise DirichletBC syntax
topBC.append(DirichletBC(FS, Constant(0.0),
"near(x[0],"+str(xbc)+") && near(x[1],"+str(ybc)+")", "pointwise"))
return topBC
def init_surface(mesh, bmesh, bmf, bflag):
'''
This function moves the nodes from respective boundary
'''
for i in range(bmesh.num_vertices()):
if (bmf[i] == bflag):
bmesh.coordinates()[i,1] -= 0.05*np.sin(4*(np.pi/L)*bmesh.coordinates()[i,0])
ALE.move(mesh, bmesh)
mesh.bounding_box_tree().build(mesh)
def update_surface(mesh, bmesh, bmf, bflag, u, Dt):
'''
This function updates interface position based on its
current velocity and current time step. The mesh is moved according to
'''
for i in range(bmesh.num_vertices()):
if (bmf[i] == bflag):
# get current boundary coordinates
x_ = bmesh.coordinates()[i,0]
y_ = bmesh.coordinates()[i,1]
# get displacemt from velocity at node
dx_ = u(x_,y_)[0]/Dt
dy_ = u(x_,y_)[1]/Dt
# update boundary coordinates
bmesh.coordinates()[i,0] += dx_
bmesh.coordinates()[i,1] += dy_
ALE.move(mesh, bmesh)
mesh.bounding_box_tree().build(mesh)
class Left(SubDomain):
def __init__(self):
SubDomain.__init__(self)
def inside(self, x, on_boundary):
return near(x[0], 0.0)# and on_boundary
class Right(SubDomain):
def __init__(self):
SubDomain.__init__(self)
def inside(self, x, on_boundary):
return near(x[0], L)# and on_boundary
class Bottom(SubDomain):
def __init__(self):
SubDomain.__init__(self)
def inside(self, x, on_boundary):
return near(x[1], 0.0)# and on_boundary
class Top(SubDomain): # This will later on become a free surface
def __init__(self):
SubDomain.__init__(self)
def inside(self, x, on_boundary):
return near(x[1], H)
L = 2*np.pi # Length of channel
H = 1 # Height of channel
# Create mesh
channel = Rectangle(Point(0,0), Point(L, H))
mesh = generate_mesh(channel, 50)
# Mesh-derived quantities:
d = mesh.geometry().dim()
I = Identity(d)
h = CellDiameter(mesh)
Time = 1 # Total time
num_steps = 100 # number of time steps
dt = Time / num_steps # time step size
# Bulk parameters
mu = 1 # kinematic viscosity
rho = 1000 # Density
sigma = 1 # Interfacial tension
g = 0 # Gravity
k = Constant(dt) # time step
mu = Constant(mu) # kinematic viscosity
rho = Constant(rho) # Density
sigma = Constant(sigma) # Surface tension
f = rho*Constant((0, - g)) # Body force
# Instantiate boundaries
left = Left()
right = Right()
top = Top()
bottom = Bottom()
# Extract boundary mesh
bmesh = BoundaryMesh(mesh, "exterior", True)
# Setup markers
bmf = MeshFunction("size_t", bmesh, 0)
bmf.set_all(0)
# Mark boundaries with flag
top.mark(bmf, 1)
right.mark(bmf, 2)
bottom.mark(bmf, 3)
left.mark(bmf, 4)
# initialize surface disturbance
init_surface(mesh, bmesh, bmf, bflag = 1)
# Initialize curvature
kappa, n = curvature(mesh,bmesh)
# Define function spaces (equal order interpolation):
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
# Define functions for solutions at previous and current time steps
u_n = Function(V)
u_ = Function(V)
p_n = Function(Q)
p_ = Function(Q)
# Define expressions used in variational forms
U = 0.5*(u_n + u)
# Pressure values on boundaries
P_top = p_pointwise(Q, 1, bmesh, bmf)
bcp = P_top
# Velocity boundary condition
symmetry = DirichletBC(V.sub(1), Constant(0.0), bottom) # Symmetry on bottom
bcu = [symmetry]
# Rate-of-deformation tensor for bulk
def DD(u):
return sym(nabla_grad(u))
# Cauchy stress tensor
def TT(u, p):
return 2*mu*DD(u) - p*I
# Interface stress tensor
def IST(sigma, kappa, n):
return sigma*kappa*n
def VariationalProblem(rho, mu, f, u, u_n, p_n, v, sigma, kappa, n):
# Variational problem for momentum conservation
# Step 1
F1 = rho*dot((u - u_n) / k, v)*dx \
+ rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(TT(U, p_n), DD(v))*dx \
- dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx
# Boundary integral term
# stress balance on interface # pressure term
b_int = inner(IST(sigma, kappa, n), v)*ds(1) - dot(p_n*n, v)*ds
F1 = F1 - b_int
a1 = lhs(F1)
L1 = rhs(F1)
# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/k)*div(u_)*q*dx
# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - k*dot(nabla_grad(p_ - p_n), v)*dx
# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
return A1, A2, A3, L1, L2, L3
# Use amg preconditioner if available
prec = "amg" if has_krylov_solver_preconditioner("amg") else "default"
# Use nonzero guesses - essential for CG with non-symmetric BC
parameters['krylov_solver']['nonzero_initial_guess'] = True
A1, A2, A3, L1, L2, L3 = VariationalProblem(rho, mu, f, u, u_n, p_n, v, sigma, kappa, n)
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1)#, "bicgstab", "default")
# Step 2: Pressure correction step
b2 = assemble(L2)
[bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2)#, "bicgstab", prec)
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3)#, "bicgstab", "default")
I would expect capillary effects to generate a pressure difference due to the interface curvature (with non-zero curvature and non-zero surface tension). However, I obtain a constant pressure field. Thanks in advance.