Update time-dependent Dirichlet BC values for FSI Problem

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.

I still haven’t figured out how to update the DirichletBC values for the moving boundary velocity. I’ve improved the MWE so it is easily readable. So far I have tthe following attempt to update the BC:

from dolfin import *
from mshr import *
import numpy as np

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)

channel = Rectangle(Point(0,0), Point(4, 2))
mesh = generate_mesh(channel, 40)

V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

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

noslip  = DirichletBC(V, Constant((0, 0)), bottom)

Right now I’m trying to update the values using a UserExpression such as the following:

class update_bc_value(UserExpression):
'''
Update the velocity value of u_moving = Function(V) for each node by simply comparing whether the coordinates are in the moving boundary.
'''
    def eval(self, value, x):
        if near(x[0], bmesh.coordinates()[i,0]) and near(x[1], bmesh.coordinates()[i,1]):
            # print("x comparison: ", (x[0], bmesh.coordinates()[i,0]))
            # print("y comparison: ", x[1], bmesh.coordinates()[i,1])
            value[0] = fsv_x
            value[1] = fsv_y

    def value_shape(self):
        return (2,)

u_moving= Function(V)
u_top.interpolate(Constant((0.0, 0.0))) # initialize with zeros
topBC = DirichletBC(V, u_moving, top)

bcu = [noslip, topBC]

I call the expression to update the values as follows:

t = 0

for n in range(10):
    for i in range(bmesh.num_vertices()):
        if (bmf[i] == 1):
            old_x = bmesh.coordinates()[i,0]
            old_y = bmesh.coordinates()[i,1]

            bmesh.coordinates()[i,1] -= 0.001*np.sin((np.pi/L)*bmesh.coordinates()[i,0])

            # compute new velocity in moving boundary 
            fsv_x = (bmesh.coordinates()[i,0] - old_x)/dt
            fsv_y = (bmesh.coordinates()[i,1] - old_y)/dt

            # update boundary values
            u_moving.interpolate(update_bc_value(element = V.ufl_element()))


   # Move mesh
    ALE.move(mesh, bmesh)
    mesh.smooth()

    # Update boundary condition
    topBC.set_value(u_moving)

From the code snippet above I’m calculating the boundary velocity of a single point that belongs to the moving boundary and then using the UserExpression previously defined to update the value of the specific corrdinate. I’m not sure if I can define a u_moving the way I did for my objective. I also didn’t figure out a proper way to set the values of topBC in a proper way.

Any thoughts are welcome.

I have solved this problem in a non-efficient procedure for my purposes (by adapting the solution provided in this link): Given the updated coordinates for each calculated point, I define the DirichletBC in a pointwise fashion (which means I cant define noslip for the top boundary line and then update it as pointwise values - at least in my understanding). So far I figured out that If I apply this approach, I must set the BC for the moving wall in a pointwise manner in every boundary configuration (whether its steady or moving). The code snipped inside the loop is as follows.

### Define boundary condition in pointwise manner
# Bottom boundary
noslipb  = DirichletBC(V, Constant((0, 0)), bottom)
bcu = [noslipb]

# Top boundary - no slip before moving
topBC = []
for i in range(bmesh.num_vertices()):
    if bmf[i]==1:
        xbc = bmesh.coordinates()[i,0]
        ybc = bmesh.coordinates()[i,1]
        valuex = 0
        valuey = 0
        topBC.append(DirichletBC(V, Constant((valuex, valuey)), "near(x[0],"+str(xbc)+") && near(x[1],"+str(ybc)+")", "pointwise"))

bcU = bcu + topBC

Still, any comments regarding this issue is welcome!