Hello,
I’ve been looking at the dynamic mesh example in this forum post: How to achieve dynamic mesh in getfem?. I’m trying to adapt it to model a beam that’s completely surrounded by fluid, with the left side clamped.
My approach is to create two separate meshes: one for the beam and another for the fluid. The fluid mesh would have a hole where the beam is located. I think that my problem is that I have different element sizes for each mesh. While this might cause some issues, it’s necessary because I plan to implement an aerodynamic simulation later, which requires a finer mesh for the fluid than for the beam. This is my code: `# Beam Fluid-Structure interaction
############################################################################
import numpy as np
import getfem as gf
gf.util_trace_level(1)
gf.util_warning_level(1)
π = np.pi
E = 1000.
ν = 0.4
L = 10.
H_liquid = 5.
H_beam = 0.5
beam_load = 10.
μ_liquid = 0.3
v_in = 1.
p_out =0.
MESH
REGIONS NAMES:
BEAM_LEFT = 101
BEAM_RIGHT = 102
BEAM_TOP = 103
BEAM_BOTTOM = 104
BEAM_SURFACE = 201
INLET = 106
OUTLET = 107
WALLS = 108
IMPORT GEOMETRY:
m_beam = gf.Mesh(‘Import’,‘gmsh’, ‘MESH_GMSH/cantileverBeam.msh’)
m_fluid = gf.Mesh(‘Import’, ‘gmsh’,‘MESH_GMSH/FluidDomainBeam.msh’)
REGIONS
pi = np.pi
#a region is set for each of the beam’s faces
m_beam.set_region(BEAM_LEFT, m_beam.outer_faces_with_direction([-1,0],pi/6))
m_beam.set_region(BEAM_RIGHT, m_beam.outer_faces_with_direction([1,0],pi/6))
m_beam.set_region(BEAM_TOP, m_beam.outer_faces_with_direction([0,1],pi/6))
m_beam.set_region(BEAM_BOTTOM, m_beam.outer_faces_with_direction([0,-1],pi/6))
a region is set for all of the beam’s faces.
it should help with the dynamic mesh hopefully
m_beam.region_merge(BEAM_SURFACE, BEAM_LEFT)
m_beam.region_merge(BEAM_SURFACE, BEAM_RIGHT)
m_beam.region_merge(BEAM_SURFACE, BEAM_TOP)
m_beam.region_merge(BEAM_SURFACE, BEAM_BOTTOM)
wall1 = 202
wall2 = 203
m_fluid.set_region(INLET, m_fluid.outer_faces_with_direction([-1,0],pi/6))
m_fluid.set_region(OUTLET, m_fluid.outer_faces_with_direction([1,0],pi/6))
m_fluid.set_region(wall1, m_fluid.outer_faces_with_direction([0,1],pi/6))
m_fluid.set_region(wall2, m_fluid.outer_faces_with_direction([0,-1],pi/6))
m_fluid.region_merge(WALLS, wall1)
m_fluid.region_merge(WALLS, wall2)
m_fluid.set_region(BEAM_SURFACE, m_fluid.inner_faces())
INTEGRATION METHOD:
mim_fluid = gf.MeshIm(m_fluid, 5)
mim_beam = gf.MeshIm(m_beam, 5)
FEM ELEMENTS:
mfu_fluid_ = gf.MeshFem(m_fluid, 2)
mfu_fluid_.set_fem(gf.Fem(“FEM_Q2_INCOMPLETE(2)”))
mfv_fluid_ = gf.MeshFem(m_fluid, 2)
mfv_fluid_.set_classical_fem(2)
mfp_fluid = gf.MeshFem(m_fluid, 1)
mfp_fluid.set_classical_fem(1)
mfu_beam_ = gf.MeshFem(m_beam, 2)
mfu_beam_.set_fem(gf.Fem(“FEM_Q2_INCOMPLETE(2)”))
REMOVE FIXED DOFS
kept_dofs = list(
set(range(mfu_fluid_.nbdof()))
-set(mfu_fluid_.basic_dof_on_region(INLET))
-set(mfu_fluid_.basic_dof_on_region(OUTLET))
-set(mfu_fluid_.basic_dof_on_region(WALLS))
)
mfu_fluid= gf.MeshFem(‘partial’, mfu_fluid_, kept_dofs)
kept_dofs = list(
set(range(mfv_fluid_.nbdof()))
-set(mfv_fluid_.basic_dof_on_region(WALLS))
)
mfv_fluid = gf.MeshFem(‘partial’, mfv_fluid_, kept_dofs)
kept_dofs = list(
set(range(mfu_beam_.nbdof()))
-set(mfu_beam_.basic_dof_on_region(BEAM_LEFT))
)
mfu_beam = gf.MeshFem(‘partial’, mfu_beam_, kept_dofs)
MODEL
md = gf.Model(“real”)
FEM VARIABLE
md.add_fem_variable(“u_fluid”, mfu_fluid)
md.add_fem_variable(“v”, mfv_fluid)
md.add_fem_variable(“p”, mfp_fluid)
md.add_fem_variable(“u_beam”, mfu_beam)
md.add_filtered_fem_variable(“mult”, mfu_fluid, BEAM_SURFACE)
BRICKS
md.add_interpolate_transformation_from_expression(“beam”, m_fluid, m_beam, “[X(1),X(2)]”)
md.add_interpolate_transformation_from_expression(“liquid”, m_beam, m_fluid, “[X(1),X(2)]”)
md.add_initialized_data(“K”, E/(3*(1-2ν))) # bulk modulus
md.add_initialized_data(“G”, E/(2(1+ν))) #shear modulus
md.add_macro(“Dev33(A)”, “A-Id(2)/3*(Trace(A)+1)”) # Deviatoric part
md.add_macro(“F(u)”, “Id(2)+Grad(u)”) # deformation gradient
md.add_macro(“J(u)”, “Det(F(u))”) # Volume change ratio
md.add_macro(“tauD(u)”, “G*pow(J(u),-2/3)*Dev33(F(u)F(u)')") # Deviatoric stress
md.add_macro(“tauH(u)”, "Klog(J(u))”) # Hydrostatic stress
md.add_nonlinear_term(mim_fluid, “((tauH(u_fluid)*Id(2)+tauD(u_fluid))*Inv(F(u_fluid)')):Grad(Test_u_fluid)”)
md.add_nonlinear_term(mim_beam, “((tauH(u_beam)*Id(2)+tauD(u_beam))*Inv(F(u_beam)')):Grad(Test_u_beam)”)
md.add_nonlinear_term(mim_fluid, “mult.Test_u_fluid”, BEAM_SURFACE) # “load” on the top of the liquid domain
md.add_nonlinear_term(mim_beam, f"{beam_load}*Test_u_beam(2)", BEAM_TOP) # uniform load on the top of the beam
md.add_nonlinear_term(mim_fluid, “(u_fluid-Interpolate(u_beam,beam)).Test_mult”, BEAM_SURFACE) # bond liquid top surface to beam bottom surface
md.add_initialized_data(“mu”, μ_liquid)
md.add_initialized_data(“v_in”, v_in)
md.add_initialized_data(“p_out”, p_out)
md.add_nonlinear_term(mim_fluid, “mu*(Grad(v)*Inv(F(u_fluid))):(Grad(Test_v)Inv(F(u_fluid)))-pTrace(Grad(Test_v)*Inv(F(u_fluid)))” # stokes problem
“+Trace(Grad(v)Inv(F(u_fluid)))Test_p")
md.add_nonlinear_term(mim_fluid, f"100(v-[v_in,0]).Test_v", INLET)
md.add_nonlinear_term(mim_beam,"Interpolate(p,liquid)(J(u_beam)*Inv(F(u_beam))Normal).Test_u_beam", BEAM_SURFACE)
md.add_nonlinear_term(mim_fluid, "100(p-p_out)*Test_p”, OUTLET)
md.solve(“noisy”, “max_iter”, 100, “max_res”, 1e-6,
“lsearch”, “simplest”, “alpha max ratio”, 3., “alpha min”, 0.02, “alpha mult”, 0.6)
mfu_fluid.export_to_vtu(“liquid.vtu”, mfu_fluid, md.variable(“u_fluid”), “Displacements”,
mfv_fluid, md.variable(“v”), “Velocities”,
mfp_fluid, md.variable(“p”), “Pressure”)
mfu_beam.export_to_vtu(“beam.vtu”, mfu_beam, md.variable(“u_beam”), “Displacements”)
`if you need i can share with you the .msh file or the .geo file for the mesh