Hello,
I’m learning how to update a Dirichlet BC in a time-dependent simulation. I have a moving wall displacing fluid as it moves. The full script is as follows:
from dolfin import *
from vedo.dolfin import ProgressBar, plot
from mshr import *
import numpy as np
import random
class top_wall(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], H)
class bottom_wall(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 0)
#=============== Space and time parameters
L = 2*np.pi # Length of channel
H = 1 # Height of channel
Time = 1 # Total time
num_steps = 100 # number of time steps
dt = Time / num_steps # time step size
# Fluid parameters
mu = 1 # kinematic viscosity
rho = 1 # density
#=============== Create mesh
channel = Rectangle(Point(0,0), Point(L, H))
mesh = generate_mesh(channel, 40)
top = top_wall()
bottom = bottom_wall()
bmesh = BoundaryMesh(mesh, "exterior", True)
bmf = MeshFunction("size_t", bmesh, 0)
bmf.set_all(0)
top.mark(bmf, 1) # mark moving boundary
#=============== Create variational formulation
V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)
# Apply boundary condition
noslip = DirichletBC(V, Constant((0, 0)), bottom)
moving_wall = DirichletBC(V, Constant((0, 0)), top)
bcu = [noslip, moving_wall]
# 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)
n = FacetNormal(mesh)
f = Constant((0, 0))
k = Constant(dt)
mu = Constant(mu)
rho = Constant(rho)
# Define strain-rate tensor
def epsilon(u):
return sym(nabla_grad(u))
# Define stress tensor
def sigma(u, p):
return 2*mu*epsilon(u) - p*Identity(len(u))
# Define variational problem for step 1
F1 = rho*dot((u - u_n) / k, v)*dx + \
rho*dot(dot(u_n, nabla_grad(u_n)), v)*dx \
+ inner(sigma(U, p_n), epsilon(v))*dx \
+ dot(p_n*n, v)*ds - dot(mu*nabla_grad(U)*n, v)*ds \
- dot(f, v)*dx
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)
# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
t = 0.0
for n in range(num_steps):
# Update current time
t += dt
# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1)
# Step 2: Pressure correction step
b2 = assemble(L2)
# [bc.apply(b2) for bc in bcp]
solve(A2, p_.vector(), b2)
# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3)
wall_vel = []
for i in range(bmesh.num_vertices()):
if (bmf[i] == 1): # The moving boundary is marked with #1
x_old = bmesh_coordinates()[i,0]
y_old = bmesh_coordinates()[i,1]
# Move mesh
bmesh.coordinates()[i,1] -= 0.001
x_new = bmesh_coordinates()[i,0] # it will always be the same in this MWE
y_new = bmesh_coordinates()[i,1]
# Store boundary velocity
wall_vel.append([(x_new - x_old)/dt, (y_new - y_old)/dt])
# Update the boundary velocity data -> update moving_wall with wall_vel
moving_wall.assign(np.asarray(wall_vel))
ALE.move(mesh, bmesh)
mesh.smooth()
I apologize for the not-so-small MWE, I didn’t know how to represent this issue without providing a general overview of how the boundaries are being defined and how the system is being solved. I’ve found some examples that uses Expressions to update boundary values, but I still didn’t figure out how to assign an entire array into the boundary velocity value as time goes. Also, is it good to define the boundaries as python classes for time-varying boundary conditions?
Any sugestions will be of great value.