8.7. Fluid-Structure Interaction
Having deformable solids (Section 7.10) and the Navier-Stokes equations on moving domains (Section 7.5) available, obviously, both can be combined for fluid-structure interaction scenarios. We consider a 2d channel with two leaflets that will deform by the flow and thereby change the flow as well. The problem case is rather short, just combining the Navier-Stokes equations in the liquid domain and the solid equations in the deformable leaflets:
from pyoomph import *
from pyoomph.equations.navier_stokes import *
from pyoomph.equations.ALE import *
from pyoomph.equations.solid import *
class SimpleFSIProblem(Problem):
def __init__(self):
super().__init__()
self.Pinlet=15 # Inlet pressure
self.claw=GeneralizedHookeanSolidConstitutiveLaw(E=40000,nu=0.3) # Solid dynamics
self.rho, self.mu=1000 , 1 # Liquid density and dynamic viscosity
self.max_refinement_level=2 # Adaptivity level
def define_problem(self):
# Mark individual domains in the RectangularQuadMesh
def domain_name(x,y):
if (x>1.8 and x<2 and y<1.5) or (x>2.8 and x<3 and y>0.5):
return "solid"
else:
return "liquid"
self+=RectangularQuadMesh(size=[5,2],name=domain_name,N=[50,20])
# Liquid equations
leqs=MeshFileOutput()
leqs+=NavierStokesEquations(mass_density=self.rho,dynamic_viscosity=self.mu)
leqs+=PseudoElasticMesh()
leqs+=PinMeshCoordinates()@["left","right"]
leqs+=DirichletBC(mesh_y=0)@"bottom"
leqs+=DirichletBC(mesh_y=2)@"top"
leqs+=NoSlipBC()@["top","bottom"]
leqs+=(DirichletBC(velocity_y=0)+NeumannBC(velocity_x=-self.Pinlet))@"left"
leqs+=DirichletBC(velocity_y=0)@"right"
# Solid equations
seqs=MeshFileOutput()
seqs+=DeformableSolidEquations(self.claw,mass_density=2,coordinate_space="C2",scale_for_FSI=True)
seqs+=PinMeshCoordinates()@"bottom"
seqs+=PinMeshCoordinates()@"top"
# Fluid-structure interaction at the mutual interface
leqs+=FSIConnection()@"liquid_solid"
# The FSIConnection only pulls on the liquid nodes lying on the interface. Around the tips of the
# leaflets, where the imposed motion varies quickly along the interface, the first layer of liquid
# elements has to absorb the whole mismatch to the almost stationary interior and eventually inverts.
# Stiffening that layer spreads the deformation over the elements further inside.
leqs+=InterfaceMeshStiffening(5)@"liquid_solid"
# Adaptivity. The two domains are adapted individually, but pyoomph keeps both sides of the
# mutual interface at the same refinement, so the error estimator can drive it freely. This used
# to require RefineToLevel()@"liquid_solid" on BOTH sides -- forcing the interface to maximum
# refinement everywhere so that the two sides could not disagree -- which is no longer needed.
leqs+=SpatialErrorEstimator(velocity=1)
self+=leqs@"liquid"+seqs@"solid"
with SimpleFSIProblem() as problem:
problem.run(50,outstep=0.5,temporal_error=1,spatial_adapt=1)
The real main part is the FSIConnection, which must be added to the liquid side of the mutual interface. In Section 8.4, we discussed how the enforcing of continuous velocity between two liquid domains via Lagrange multipliers actually ensures the balance of tractions. The same idea is used in the pyoomph.equations.solid.FSIConnection interface. We enforce that the liquid velocity agrees with the solid velocity and thereby also ensure the balance of the tractions at the shared interface. Moreover, the fluid mesh is moved with the solid mesh. As discussed in Section 8.1, it is important to use the same scale of the test function on both sides to balance the tractions. Therefore, one has to set scale_for_FSI=True in the DeformableSolidEquations.
As opposed to the ConnectMeshAtInterface class, which moves the nodes of the meshes on both sides, the FSIConnection only moves the nodes of the liquid mesh to match those of the solid mesh. Otherwise, the particular moving mesh dynamics of the fluid domain, which does not reflect any physics, would add additional unphysical tractions to the system.
The FSIConnection constrains only those liquid nodes that lie on the interface. Around the tips of the leaflets, the imposed motion varies strongly along the interface, so the first layer of liquid elements has to absorb the entire mismatch between the strongly moving interface nodes and the almost stationary interior. Without further measures, these elements shear and squash so badly that some of them are inverted during roughly a quarter of the run, from \(t\approx13\) on. With the InterfaceMeshStiffening added on the same interface, none of them inverts. It is an interface equation, but it is tested against gradients of the bulk mesh test function, i.e. of testfunction("mesh",domain=".."): the shape function of a node that is not on the face vanishes there, but its gradient does not, so a pure surface integral does reach the interior nodes of the attached elements. It is the surface-concentrated limit of a bulk integral over a layer of the thickness of one element, which is why its argument is simply the relative extra stiffness of that layer – the deformation is thereby handed on to the elements further inside, where more of them can share it.
Note that spatial adaptivity is driven here by a SpatialErrorEstimator on the liquid velocity alone. The solid domain is given no refinement criterion of its own, yet both sides of the mutual interface must always be refined identically – an interface element on one side has to find exactly one counterpart on the other. Since the two domains are separate meshes, oomph-lib adapts them individually, so pyoomph takes care of this itself: after every adaptation, the coarser side of a connected interface is refined until the two sides carry matching facets again. You therefore do not have to constrain the interface refinement by hand, and this works under mpirun --distribute as well, where the two domains are partitioned independently of each other.
Fluid-Structure Interaction