How to deal with submeshes (fluid dans solid) in an FSI ALE context

General
1

Hello,

I am new to Fenics and I am trying to deal with FSI / ALE methodology. I would like to run a similar problem as the heart valve example described by Kristoffer Selim and Anders Logg (I found the paper here: http://www.logg.org/anders/pub/papers/SelimLogg2009a.pdf).

This is a 2D problem involving an elastic solid immersed in an incompressible viscous fluid flow.

In order to create the 2D geometry and mesh, I used FreeCAD and then gmsh. As a result, I have a .msh file containing the two mesh regions and all the boundaries. I can use dolfin-convert to get a set of xml files: .xml, facet_region.xml and physical_region.xml.

My problem is that I don’t understand how to split the mesh into mesh_F and mesh_S in order to solve the fluid problem and solid problem separately and then use ALE.move() to update the mesh positions.

Can anyone help?

Thank you in advance ^^

2

Btw, I am working on ubuntu and I used confa-forge to install Fenics 2019.1.0.

3

I do not have much experience with Fluid-Structure-Interaction in Fenics, but you can also have a look at TurtleFSI. This is is a monolithic fluid-structure interaction solver written in FEniCS.

1 Like
4

Thank you for your message Tuderic.

I saw TurtleFSI and other monolithic FSI implementations but I would like to focus on ALE.

I wrote a code which is obviously not working :expressionless: , but maybe it can be a good start.

The idea is to:
1- create the mesh, sub_meshes, facets and subdomains from .xml files
2- define the function spaces, test and trial functions, and other usefull functions
3- define the periodic boundary conditions
4- define the fluid solver thanks to Chorin’s projection method
5- define the structure solver
6- run on each timestep an iteration method until convergence (the mesh displacement is smaller than a certain value):
6.1- solve the velocity and pressure for the fluid
6.2- solve the velocity and change of velocity for the structure
6.3- move the mesh and calculate the mesh velocity

The model + mesh is:

tube_tige_model
tube_tige_model1380Ă—568 53.7 KB

22 is the inflow edge
23 is the outflow edge
24 is the tube wall for fluid mesh
25 is the fluid/structure interface
26 is the tube wall for solid mesh
27 is the fluid mesh
28 is the solid mesh

The code I am working on is the following one:

#####################################

import fenics as fe
from tqdm import tqdm

TIME_STEP_LENGTH = 0.001
N_TIME_STEPS = 20
KINEMATIC_VISCOSITY = 0.001
DT_write = 1
DT_SAVE = 10 * DT_write
d = 2
coef_poisson = 0.35
module_young = 10000
lame_mu = module_young/(2 * (1+coef_poisson))
lame_lambda = module_young * coef_poisson/((1+coef_poisson) * (1-2*coef_poisson))
F = fe.Constant((0.,0.))
smooth_param = 50

marker_SF_border = 25
marker_F = 27
marker_S = 28

maxiter = 50
TOLERANCE = 0.1

ongoing_simulation = False

def main():

if ongoing_simulation:
    xml_mesh_file = "saved_mesh.xml"
    xml_facet_file = "saved_facet_mesh.xml"
    xml_region_file = "saved_region_mesh.xml"
else:
    xml_mesh_file = "tube_tige_2D.xml"
    xml_facet_file = "tube_tige_2D_facet_region.xml"
    xml_region_file = "tube_tige_2D_physical_region.xml"
    
mesh = fe.Mesh(xml_mesh_file)

fd = fe.MeshFunction("size_t", mesh, xml_facet_file)
sd = fe.MeshFunction("size_t", mesh, xml_region_file)
    
dx_ = fe.Measure('dx', domain=mesh, subdomain_data=sd)
ds_ = fe.Measure('ds', domain=mesh, subdomain_data=fd)    

mesh_F = fe.SubMesh(mesh, sd, marker_F)
mesh_S = fe.SubMesh(mesh, sd, marker_S)

u_F_file_pvd = fe.File("../../results/tube_tige_2D/velocity/u_F_file.pvd")
p_F_file_pvd = fe.File("../../results/tube_tige_2D/pressure/p_F_file.pvd")
u_S_file_pvd = fe.File("../../results/tube_tige_2D/pressure/u_S_file.pvd")

# Define the function spaces
velocity_F_function_space = fe.VectorFunctionSpace(mesh_F, "Lagrange", 2)
pressure_F_function_space = fe.FunctionSpace(mesh_F, "Lagrange", 1)
velocity_S_function_space = fe.VectorFunctionSpace(mesh_S, "Lagrange", 2)

u_F_trial = fe.TrialFunction(velocity_F_function_space)
p_F_trial = fe.TrialFunction(pressure_F_function_space)
u_S_trial = fe.TrialFunction(velocity_S_function_space)

v_F_test = fe.TestFunction(velocity_F_function_space)
q_F_test = fe.TestFunction(pressure_F_function_space)
v_S_test = fe.TestFunction(velocity_S_function_space)

# Define u_F_mesh
    
u_F_mesh = fe.Function(velocity_F_function_space)
dms = [velocity_F_function_space.sub(i).dofmap().dofs() for i in range(2)]

# Define the Boundary Condition

no_slip_F = fe.DirichletBC(
    velocity_F_function_space,
    fe.Constant((0.0, 0.0)),
    fd,
    24
)

no_slip_S = fe.DirichletBC(
    velocity_S_function_space,
    fe.Constant((0.0, 0.0)),
    fd,
    26
)

in_velocity_F = fe.DirichletBC(
    velocity_F_function_space,
    fe.Constant((1.0, 0.0)),
    fd,
    22
) 

velocity_F_boundary_conditions = [no_slip_F, in_velocity_F]
pressure_F_boundary_conditions = []

velocity_S_boundary_conditions = [no_slip_S]
    
# Define the solution fields involved
u_F_prev = fe.Function(velocity_F_function_space)
u_F_tent = fe.Function(velocity_F_function_space)
u_F_next = fe.Function(velocity_F_function_space, name="displacement_F")
p_F_next = fe.Function(pressure_F_function_space, name="pressure_F")

u_S_prev = fe.Function(velocity_S_function_space)
u_S_next = fe.Function(velocity_S_function_space, name="displacement_S")

#initial conditions        
if ongoing_simulation:
    u_F_prev = fe.Function(velocity_F_function_space, "u_F_save.xml")
    p_F_next = fe.Function(pressure_F_function_space, "p_F_save.xml")
    u_S_prev = fe.Function(velocity_S_function_space, "u_S_save.xml")

### SOLVER - FLUID

# Weak form of the pressure poisson problem
pressure_poisson_weak_form_lhs_F = fe.inner(fe.grad(p_F_trial), fe.grad(q_F_test)) * dx_(marker_F)
pressure_poisson_weak_form_rhs_F = - 1.0 / TIME_STEP_LENGTH * fe.div(u_F_tent) * q_F_test * dx_(marker_F)

# Weak form of the velocity update equation
velocity_update_weak_form_lhs_F = fe.inner(u_F_trial, v_F_test) * dx_(marker_F)
velocity_update_weak_form_rhs_F = (
    fe.inner(u_F_tent, v_F_test) * dx_(marker_F)
    -
    TIME_STEP_LENGTH * fe.inner(fe.grad(p_F_next), v_F_test) * dx_(marker_F)
)

# Pre-Compute the system matrices
pressure_poisson_assembled_system_matrix_F = fe.assemble(pressure_poisson_weak_form_lhs_F)
velocity_update_assembled_system_matrix_F = fe.assemble(velocity_update_weak_form_lhs_F)

def solve_F(u_F_prev, p_F_next, u_mesh):

    momentum_weak_form_residuum_F = (
        1.0 / TIME_STEP_LENGTH * fe.inner(u_F_trial - u_F_prev, v_F_test) * dx_(marker_F)
        +
        fe.inner(fe.grad(u_F_prev) * (u_F_prev - u_mesh), v_F_test) * dx_(marker_F)
        +
        KINEMATIC_VISCOSITY * fe.inner(fe.grad(u_F_trial), fe.grad(v_F_test)) * dx_(marker_F)
        -
        fe.inner(F, v_F_test) * dx_(marker_F)

    )
    momentum_weak_form_lhs_F = fe.lhs(momentum_weak_form_residuum_F)
    momentum_weak_form_rhs_F = fe.rhs(momentum_weak_form_residuum_F)
        
    momentum_assembled_system_matrix_F = fe.assemble(momentum_weak_form_lhs_F)

    # Solve u_F_tent
    momentum_assembled_rhs_F = fe.assemble(momentum_weak_form_rhs_F)
    [bc.apply(momentum_assembled_system_matrix_F, momentum_assembled_rhs_F) for bc in velocity_F_boundary_conditions]
    fe.solve(
        momentum_assembled_system_matrix_F,
        u_F_tent.vector(),
        momentum_assembled_rhs_F,
        "gmres",
        "ilu",
    )

    # Solve p_F_next
    pressure_poisson_assembled_rhs_F = fe.assemble(pressure_poisson_weak_form_rhs_F)
    [bc.apply(pressure_poisson_assembled_system_matrix_F, pressure_poisson_assembled_rhs_F) for bc in pressure_F_boundary_conditions]
    fe.solve(
        pressure_poisson_assembled_system_matrix_F,
        p_F_next.vector(),
        pressure_poisson_assembled_rhs_F,
        "gmres",
        "amg",
    )

    # Solve u_F_next
    velocity_update_assembled_rhs_F = fe.assemble(velocity_update_weak_form_rhs_F)
    [bc.apply(momentum_assembled_system_matrix_F, momentum_assembled_rhs_F) for bc in velocity_F_boundary_conditions]
    fe.solve(
        velocity_update_assembled_system_matrix_F,
        u_F_next.vector(),
        velocity_update_assembled_rhs_F,
        "gmres",
        "ilu",
    )
    
    return u_F_next, p_F_next
   
    
### SOLVER STRUCTURE

def epsilon(u):
    return 0.5*(fe.grad(u) + fe.grad(u).T)

def sigma_F(u,p):
    return 0.2*KINEMATIC_VISCOSITY*epsilon(u)-fe.Identity(d)*p

def sigma_S(u):
    return 2 * lame_mu * epsilon(u) + lame_lambda * fe.div(u) * fe.Identity(d)

def solve_S(u_S_prev, u_F_next, p_F_next):

    n_F = fe.FacetNormal(mesh_F)

    structure_weak_form_lhs = fe.inner(epsilon(v_S_test), sigma_S(u_S_trial)) * dx_(marker_S)#fe.dx
    structure_weak_form_rhs = fe.inner(v_S_test, F) * dx_(marker_S) - fe.inner(v_S_test, sigma_F(u_F_next, p_F_next) * n_F) * ds_(marker_SF_border)
    
        
    structure_assembled_system_matrix = fe.assemble(structure_weak_form_lhs)
    
    # Solve u_S_next
    structure_assembled_rhs = fe.assemble(structure_weak_form_rhs)
    [bc.apply(structure_assembled_system_matrix, structure_assembled_rhs) for bc in velocity_S_boundary_conditions]
    fe.solve(
        structure_assembled_system_matrix,
        u_S_next.vector(),
        structure_assembled_rhs,
        "gmres",
        "ilu",
    )

    du_S_next = u_S_next - u_S_prev

    return u_S_next, du_S_next

### TIME SIMULATION


for t in tqdm(range(N_TIME_STEPS)):
    
    u_F_mesh.vector()[:] = 0.0
    
    for i in range(maxiter):

        # Initial mesh coordinates
        dof_coor_prev = velocity_F_function_space.dofmap().tabulate_all_coordinates(mesh_F)
        coor_x_prev = dof_coor_prev[dms[0]]
        coor_y_prev = dof_coor_prev[dms[1]]

        # Solve fluid problem
        u_F_next, p_F_next = solve_F(u_F_prev, p_F_next, u_F_mesh)
        
        # Solve structure problem
        u_S_next, du_S_next = solve_S(u_S_prev)
        
        # Move structure mesh
        mesh_S.move(du_S_next)
        
        # Move fluid mesh
        mesh_F.move(mesh_S)
        
        # Smooth mesh
        mesh_F.smooth(smooth_param)
        
        # Final mesh coordinates
        dof_coor_next = velocity_F_function_space.dofmap().tabulate_all_coordinates(mesh_F)
        coor_x_next = dof_coor_next[dms[0]]
        coor_y_next = dof_coor_next[dms[1]]
        
        # calculate u_F_mesh
        
        u_F_mesh.vector()[dms[0]] = (coor_x_next - coor_x_prev)/TIME_STEP_LENGTH
        u_F_mesh.vector()[dms[1]] = (coor_y_next - coor_y_prev)/TIME_STEP_LENGTH
        
        u_F_mesh.assign(u_F_mesh)

        # check convergence
        if fe.norm(du_S_next) < TOLERANCE:
            break
    
    if i == maxiter -1:
        raise RuntimeError("FSI iteration did no converge!")

    # Advance in time
    u_F_prev.assign(u_F_next)
    u_S_prev.assign(u_S_next)
    
    if t%DT_write==0:
        u_F_file_pvd << (u_F_next, float(t/DT_write))
        p_F_file_pvd << (p_F_next, float(t/DT_write))
        u_S_file_pvd << (u_S_next, float(t/DT_write))
        
    # Save mesh and future initial conditions
    if t%DT_SAVE==0:
        fe.File("saved_mesh.xml") << mesh
        fe.File("saved_facet_mesh.xml") << fd
        fe.File("saved_region_mesh.xml") << fd
        fe.File("u_F_save.xml") << u_F_next
        fe.File("p_F_save.xml") << p_F_next
        fe.File("u_S_save.xml") << u_S_next

if name == “main”:
main()

#####################################

So well, at the moment, I already have a problem with fe.DirichletBC whereas I used a similar code on simpler model involving fluid only.

I suppose I will also have some issues with the dx_ and ds_ (?)
I doubt the F/S interface condition coded in the structure solver will work correctly (?)
I also don’t know if the way I calculate the fluid mesh velocity is correct (saw it on another post).

If anyone can criticize my code, I’d be happy about it :smiley:

5

Hi,

I finally tried turtleFSI which is working really well. All is about monitoring different BCs, physical and model parameters and deal with the Newton solver. It’s not necessary to invent the wheel once again :thinking:

Moreover, all the code is available on github and it is possible to play with.

Thank you @tuderic ^^

1 Like
6

Hi sdc,

I want to work with 2 different meshes. And I think my problem is similar to yours. But my domains are disconnected, and I want to define pressure jump between 2 disconnected boundaries. You can see my topic in here. If you can help, that would be very helpful.

7

Hi @Ekrem_Ekici ,

I am not sure what you mean by jump conditions but maybe this can help you:

https://fenicsproject.org/pub/tutorial/sphinx1/._ftut1005.html

Cheers,

8

btw, once I have a mesh generated by FreeCAD + gmsh, I run a script I wrote to rewrite a .msh file without any redundant nodes at the boundary in order to get only one mesh with different domains (fluid and structure(s)), and then use dolfin-convert or another script to transform the .msh to a xdmf file. And so far, it seems to work well with turtleFSI.

1 Like