Strange numerical divergence in FSI example

[zhangmuElias]

Dear everyone,

To use water as the fluid medium, I made some changes to the benchmark case TF_fsi.
In order to keep the Reynolds number still be 200, I reduced the size of the entire computing domain by a factor of 100, and the inflow velocity $U_m$ by a factor of 10. This case can run well when the plate is rigid. However, when I added the solid plate part and activated the solid solver, some strange velocity and pressure surge emerged at the FSI surface near the free and fixed end of the plate, as you can see in the following figures. Do you know why this happened and how we can ovoid this problem?

The case and mesh file are also attached.

cylinderflap_tail_head_cy2_4D_bl.zip

Figure 1 Geometry and mesh

Figure 2 Results of rigid plate

Figure 3 Results of elastic plate

from dolfin import *
import numpy as np
from os import path
from mpi4py import MPI as pyMPI

from turtleFSI.problems import *
from turtleFSI.modules import *

parameters["ghost_mode"] = "shared_facet"
#_compiler_parameters = dict(parameters["form_compiler"])


def set_problem_parameters(default_variables, **namespace):
    # Overwrite or add new variables to 'default_variables'
    default_variables.update(dict(
        # Temporal variables
        T=2,                         # End time [s]
        dt=0.001,                      # Time step [s]
        theta=0.501,                    # Temporal scheme

        # Physical constants ('FSI 3')
        fluid_properties=[],
        solid_properties=[], 
        material_model="StVenantKirchoff",
        Um=0.2,                       # Max. velocity inlet, CDF3: 2.0 [m/s]
        rho_f=1.0e3,                  # Fluid density [kg/m3]
        mu_f=1.0e-3,                     # Fluid dynamic viscosity [Pa.s]
        rho_s=1.0e3,                  # Solid density[kg/m3]
        nu_s=0.4,                     # Solid Poisson ratio [-]
        mu_s=2.0e6,                   # Shear modulus, CSM3: 0.5E6 [Pa]
        lambda_s=4e6,                 # Solid 1st Lame Coefficient [Pa]
        fluid="fluid",                             # ["fluid", "no_fluid"] Turn off fluid and only solve the solid problem
        solid="solid",                             # ["solid", "no_solid"] Turn off solid and only solve the fluid problem


        # Problem specific
        folder="CylinderFlap-n2-tail-head-withCy-bl",      # Name of the results folder
        extrapolation="biharmonic",   # No displacement to extrapolate
        extrapolation_sub_type="constrained_disp",  # Biharmonic type
        bc_ids=[11, 12, 13, 16],          # Ids of makers for the mesh extrapolation
        #bc_ids=[2, 3, 4, 6],          # Ids of makers for the mesh extrapolation

        # Domain settings
        dx_f_id=17,       # Domain id of the fluid domain
        dx_s_id=18,       # Domain id of the solid domain

        # Solver settings
        recompute=1,                  # Compute the Jacobian matrix every iteration
        save_step=1,  # Save file frequency

        # Geometric variables
        R=0.0005,                       # Radius of the circle
        H=0.0041,                       # Total height
        L=0.25,                        # Length of domain
        f_L=0.0035,                     # Length of the flag
        f_H=0.0002,                     # Height of the flag
        c_x=0.002,                      # Center of the circle x-direction
        c_y=0.002))                     # Center of the circle y-direction

    default_variables["compiler_parameters"].update({"quadrature_degree": 5})

    return default_variables
"""
1 11 "Inlet"
1 12 "Outlet"
1 13 "Wall"
1 14 "Bar"
1 15 "Barwall"
1 16 "Circle"
2 17 "Fluid"
2 18 "Solid"
"""

def get_mesh_domain_and_boundaries(R, H, L, f_L, f_H, c_x, c_y, **namespace):
    # Read mesh
    mesh_folder = path.join(path.dirname(path.abspath(__file__)), "mesh")
    mesh_path = path.join(mesh_folder, "CylinderFlap-n2-tail-head-withCy2-4D-bl.h5")     # mesh geometry

    mesh = Mesh()

    hdf = HDF5File(mesh.mpi_comm(), mesh_path, "r")
    hdf.read(mesh, "/mesh", False)
    domains = MeshFunction("size_t", mesh, mesh.geometry().dim())
    hdf.read(domains, "/domains")
    boundaries = MeshFunction("size_t", mesh, mesh.geometry().dim() - 1)
    hdf.read(boundaries, "/boundaries")


    return mesh, domains, boundaries


def initiate(c_x, c_y, R, f_L, **namespace):
    # Coordinate for sampling statistics
    coord = [c_x + R + f_L+8*R, c_y] # point at the bar right end

    # Lists to hold results
    displacement_x_list = []
    displacement_y_list = []
    drag_list = []
    lift_list = []
    time_list = []

    return dict(displacement_x_list=displacement_x_list, displacement_y_list=displacement_y_list,
                drag_list=drag_list, lift_list=lift_list, time_list=time_list, coord=coord)


class Inlet(UserExpression):
    def __init__(self, Um, H, **kwargs):
        self.Um = 1.5 * Um
        self.H = H
        self.factor = 0
        super().__init__(**kwargs)

    def update(self, t):
        if t < 0.2:
            self.factor = 0.5 * (1 - np.cos(t * np.pi / 0.2)) * self.Um
        else:
            self.factor = self.Um

    def eval(self, value, x):
        value[0] = self.factor * x[1] * (self.H - x[1]) / (self.H / 2.0)**2
        value[1] = 0

    def value_shape(self):
        return (2,)


def create_bcs(DVP, v_deg, Um, H, boundaries, extrapolation_sub_type, **namespace):
    inlet = Inlet(Um, H, degree=v_deg)

    # Fluid velocity conditions
    u_inlet = DirichletBC(DVP.sub(1), inlet, boundaries, 11)
    u_wall = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 13)
    u_circ = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 16)
    u_barwall = DirichletBC(DVP.sub(1), ((0.0, 0.0)), boundaries, 15)

    # Pressure Conditions
    p_out = DirichletBC(DVP.sub(2), 0, boundaries, 12)

    # Assemble boundary conditions
    bcs = [u_wall, u_inlet, u_circ, u_barwall, p_out]

    # Boundary conditions on the displacement / extrapolation
    if extrapolation_sub_type != "constrained_disp_vel":
        d_wall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 13)
        d_inlet = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 11)
        d_outlet = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 12)
        d_circle = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 16)
        d_barwall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 15)
        for i in [d_wall, d_inlet, d_outlet, d_circle, d_barwall]:
            bcs.append(i)

    else:
        w_wall = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 13)
        w_inlet = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 11)
        w_outlet = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 12)
        w_circle = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 16)
        w_barwall = DirichletBC(DVP.sub(3), ((0.0, 0.0)), boundaries, 15)

        d_wall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 13)
        d_inlet = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 11)
        d_outlet = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 12)
        d_circle = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 16)
        d_barwall = DirichletBC(DVP.sub(0), ((0.0, 0.0)), boundaries, 15)

        for i in [w_wall, w_inlet, w_outlet, w_circle, w_barwall,
                  d_wall, d_inlet, d_outlet, d_circle, d_barwall]:
            bcs.append(i)

    return dict(bcs=bcs, inlet=inlet)


def pre_solve(t, inlet, **namespace):
    """Update boundary conditions"""
    inlet.update(t)


################################################################################
# the function mpi4py_comm and peval are used to overcome FEniCS limitation of
# evaluating functions at a given mesh point in parallel.
# https://fenicsproject.discourse.group/t/problem-with-evaluation-at-a-point-in
# -parallel/1188


def mpi4py_comm(comm):
    '''Get mpi4py communicator'''
    try:
        return comm.tompi4py()
    except AttributeError:
        return comm


def peval(f, x):
    '''Parallel synced eval'''
    try:
        yloc = f(x)
    except RuntimeError:
        yloc = np.inf*np.ones(f.value_shape())

    comm = mpi4py_comm(f.function_space().mesh().mpi_comm())
    yglob = np.zeros_like(yloc)
    comm.Allreduce(yloc, yglob, op=pyMPI.MIN)

    return yglob
################################################################################


def post_solve(t, dvp_, coord, displacement_x_list, displacement_y_list, drag_list, lift_list, mu_f, n,
               verbose, time_list, ds, dS, **namespace):
    d = dvp_["n"].sub(0, deepcopy=True)
    v = dvp_["n"].sub(1, deepcopy=True)
    p = dvp_["n"].sub(2, deepcopy=True)

    # Compute drag and lift
    Dr = -assemble((sigma(v, p, d, mu_f)*n)[0]*ds(16))
    Li = -assemble((sigma(v, p, d, mu_f)*n)[1]*ds(16))
    Dr += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[0]*dS(14))
    Li += -assemble((sigma(v("+"), p("+"), d("+"), mu_f)*n("+"))[1]*dS(14))

    # Append results
    drag_list.append(Dr)
    lift_list.append(Li)
    time_list.append(t)
    d_eval = peval(d, coord)
    displacement_x_list.append(d_eval[0])
    displacement_y_list.append(d_eval[1])
    # displacement_x_list.append(d(coord)[0])
    # displacement_y_list.append(d(coord)[1])

    # Print
    if MPI.rank(MPI.comm_world) == 0 and verbose:
        print("Distance x: {:e}".format(displacement_x_list[-1]))
        print("Distance y: {:e}".format(displacement_y_list[-1]))
        print("Drag: {:e}".format(drag_list[-1]))
        print("Lift: {:e}".format(lift_list[-1]))


def finished(results_folder, displacement_x_list, displacement_y_list, drag_list, lift_list, time_list, **namespace):
    if MPI.rank(MPI.comm_world) == 0:
        np.savetxt(path.join(results_folder, 'Lift.txt'), lift_list, delimiter=',')
        np.savetxt(path.join(results_folder, 'Drag.txt'), drag_list, delimiter=',')
        np.savetxt(path.join(results_folder, 'Time.txt'), time_list, delimiter=',')
        np.savetxt(path.join(results_folder, 'dis_x.txt'), displacement_x_list, delimiter=',')
        np.savetxt(path.join(results_folder, 'dis_y.txt'), displacement_y_list, delimiter=',')

[keiyamamo]

Hi @zhangmuElias
you wrote that the left end of plate is fixed, but i can't seem to find such a boundary condition. In fact, I don't think I can find any bc for the plate.

[zhangmuElias]

Dear @keiyamamo, thanks for replying. In fact I still used the name barwall for the left circle bc of the plate, which has the id of 15.
Do you think it's enough?

[zhangmuElias]

Dear all:

To update, I tested other mesh models for the case file I attached above.
It is surprising that the Laplace mesh model with volume_change, constant and small_constant successfully simulates this case with water as the fluid medium. The strange velocity surge near the head and tail of the plate disappears.
Therefore at least there are no missing settings in my attached case file above.

I guess when the fluid viscosity changes from 1 [Pa.s] to 1e-3 [Pa.s], the diffusion effect decreases which might make the solver numerically unstable. The Laplace mesh model may serve itself as an extra diffusion term to stabilise the solver.

I also tried different parameters of alpha_u in Biharmonic model from 0.01 to 0.1 to 1 to 10, but all of them failed with the same problem I reported above.
There is a possibility Biharmonic model may have some difficulties in dealing with low-viscosity problems, I need to check it from articles. I am not sure about the exact reason. Did you try using the Biharmonic model in fluids like air or water before? It might help if others can share experiences.

Figure: the result of the reduced scaled case with water as fluid (Re=200)

Best regards
ZHANG

[keiyamamo]

Hi @zhangmuElias

Good to know it worked. I don’t have much experience with biharmonic model, thanks for sharing your experience with us.

[zhangmuElias]

Dear @keiyamamo, then in case of fsi problem, which mesh model do you use? Although the Laplace model can eliminate the velocity and pressure surge problem in the low viscosity situation, I also found that whenever the displacement is a little bit too large, the Newton solver will diverge. So this modification of biharmonic to Laplace can only partly solve this problem.

[keiyamamo]

I use Laplace model and it works fine for us because the displacement is not so large. One thing that I noticed now is that your mesh size changes dramatically from boundary layer to adjacent cells. Additionally, I’m not sure if it’s a good idea to have boundary layer around the moving object since boundary layer cells have skewed shape, potentially causing “crash” of cells when deformed.

[zhangmuElias]

Dear @keiyamamo, thanks for noticing me about the mesh.
I also tested the original benchmark case TF_fsi inside the file folder problems. It also diverged at time t=3.7s, and I noticed that it was the time when the plate started to vibrate. But, as shown in the article "Slyngstad, Andreas Strøm. Verification and Validation of a Monolithic Fluid-Structure Interaction Solver in FEniCS. A comparison of mesh lifting operators. MS thesis. 2017." the laplace model can at least handle the fsi2 and fsi3 benchmark case. I am not sure why the TF_fsi using laplace diverged in my computer, have you tried that?

Figure: The result of TF_fsi case at the last time step

[keiyamamo]

Hi @zhangmuElias
Unfortunately, no. I have only tried biharmonic type for TF_fsi problem.

[zhangmuElias]

Dear @keiyamamo , ok I see, thanks!

[zhangmuElias]

Dear @keiyamamo , in these days' tests, I found the parameter of the biharmonic model $\alpha_U$ should be reduced for low-viscosity fluid. In my case, when $\alpha_U$ equals 1e-8, the problem I reported above can be solved.

[keiyamamo]

Dear @zhangmuElias , thank you for sharing it with us and good to know it worked for you! Please feel free to close the issue.

[zhangmuElias]

ok, dear @keiyamamo , thanks for the discussion with me!