Tutorial 11: Fluid-Structure Interaction
In this tutorial, we will learn
- How to solve a surface-coupled multi-physics problem.
- Construct FE spaces defined in different domains.
- Define interface triangulations and integrate on them.
Problem statement
Strong form
Let $\Gamma_{\rm FS}$ be the interface between a fluid domain $\Omega_{\rm F}$ and a solid domain $\Omega_{\rm S}$. We denote by $\Gamma_{\rm F,D}$ and $\Gamma_{\rm F,N}$ the fluid boundaries with Dirichlet and Neumann conditions, respectively. The Fluid-Structure Interaction (FSI) problem reads:
find $u_{\rm F}$, $p_{\rm F}$ and $u_{\rm S}$ such that
satisfying the Dirichlet and Neumann boundary conditions
and the kinematic and dynamic conditions at the fluid-solid interface
Where $\boldsymbol{\sigma}_{\rm F}(u_{\rm F},p_{\rm F})=2\mu_{\rm F}\boldsymbol{\varepsilon}(u_{\rm F}) - p_{\rm F}\mathbf{I}$ and $\boldsymbol{\sigma}_{\rm S}(u_{\rm S})=2\mu_{\rm S}\boldsymbol{\varepsilon}(u_{\rm S}) +\lambda_{\rm S}tr(\boldsymbol{\varepsilon}(u_{\rm S}))\mathbf{I}$.
Geometry and Discrete model
In this tutorial we solve the benchmark descrived in [1], consisting on a flow over an elastic flag after a cylinder. The computational domain is defined by a channel of size $\Omega \doteq (0,4.5)\times(0,0.41)$, with an embedded cylinder of radius $R=0.05$ and center at $C=(0.2,0.2)$. The associated FE triangulation is denoted by $\mathcal{T}$, the fluid and solid domain and their associated triangulations will be denoted by $\Omega_{\rm F}$, $\Omega_{\rm S}$, $\mathcal{T}_{\rm F}$ and $\mathcal{T}_{\rm S}$, respectively.
In order to load the discrete model we first setup Gridap
using Gridap
The discrete model for the elastic flag problem is generated by loading the "../models/elasticFlag.json"
file.
model = DiscreteModelFromFile("../models/elasticFlag.json")
We can inspect the loaded geometry and associated parts by printing to a vtk
file:
writevtk(model,"model")
This will produce an output in which we can identify the different parts of the domain, with the associated labels and tags.
Part | Notation | Label | Tag |
---|---|---|---|
Solid-cylinder wall | $\Gamma_{\rm S,D_{cyl}}$ | "fixed" | 1 |
Fluid-solid interface | $\Gamma_{\rm FS}$ | "interface" | 2 |
Channel inlet | $\Gamma_{\rm F,D_{in}}$ | "inlet" | 3 |
Channel outlet | $\Gamma_{\rm F,N_{out}}$ | "outlet" | 4 |
Channel walls | $\Gamma_{\rm F,D_{wall}}$ | "noSlip" | 5 |
Fluid-cylinder wall | $\Gamma_{\rm F,D_{cyl}}$ | "cylinder" | 6 |
Fluid domain | $\Omega_{\rm F}$ | "fluid" | 7 |
Solid domain | $\Omega_{\rm S}$ | "solid" | 8 |
External conditions and properties
Boundary conditions
We apply Dirichlet boundary conditions at the channel inlet, upper and lower boundaries and on the cylinder. A parabolic profile is enforced at the channel inlet, while a no-slip condition is imposed on the remaining Dirichlet boundaries.
const U = 1.0
const H = 0.41
uf_in(x) = VectorValue( 1.5 * U * x[2] * ( H - x[2] ) / ( (H/2)^2 ), 0.0 )
uf_0(x) = VectorValue( 0.0, 0.0 )
us_0(x) = VectorValue( 0.0, 0.0 )
We consider a free tranction condition at the channel outlet
hN(x) = VectorValue( 0.0, 0.0 )
p_jump(x) = 0.0
External forces
In this test, the body forces acting on the fluid an solid are zero.
f(x) = VectorValue( 0.0, 0.0 )
s(x) = VectorValue( 0.0, 0.0 )
g(x) = 0.0
Material properties
We use a linear elastic constitutive law for the elastic flag. Given the Young's modulus $E$ and the Poisson ratio $\nu$, we can compute the Lamé constants, $\lambda$ and $\mu$, using the following function:
function lame_parameters(E,ν)
λ = (E*ν)/((1+ν)*(1-2*ν))
μ = E/(2*(1+ν))
(λ, μ)
end
Then, we get the Lamé parameters for a solid with $E=1.0$ MPa and $\nu=0.33$.
const E_s = 1.0
const ν_s = 0.33
const (λ_s,μ_s) = lame_parameters(E_s,ν_s)
The Cauchy stress tensor for the solid part is defined by $\sigma_s = 2\mu\varepsilon + \lambda tr(\varepsilon)\mathbf{I}$. Note that we use the trace operator from the LinearAlgebra
package. With the macro @law
we are able to define a function whose arguments depend on the coordinates, without the need of passing such coordinates as an argument.
using LinearAlgebra: tr
@law σ_s(ε) = λ_s*tr(ε)*one(ε) + 2*μ_s*ε
For the fluid part, we only need to define the viscosity $\mu_f$, which we set to $0.01$.
const μ_f = 0.01
The Cauchy stress tensor for the fluid part is given by $\sigma_f = \sigma^{dev}_f - p\mathbf{I}$, with $\sigma^{dev}_f=2\mu_f$ the deviatoric part of the stress. Since we use a mixed form with the pressure $p$ as an unknown, the stress law will only involve the deviatoric part.
@law σ_dev_f(ε) = 2*μ_f*ε
Numerical scheme
FE Spaces
In this tutorial we use an inf-sup stable velocity-pressure pair, $P_k/P_{k-1}$ elements, with continuous velocities and pressures. We select the same velocity interpolation space for the fluid and the solid parts, defined as follows
where $T$ denotes an arbitrary cell of the FE mesh $\mathcal{T}_{\rm X}$, and $P_k(T)$ is the usual Lagrangian FE space of order $k$ defined on a mesh of triangles or tetrahedra. Note that the sub-index $(\cdot)_{\rm X}$ stands for the fluid or solid parts, $(\cdot)_{\rm F}$ or $(\cdot)_{\rm S}$, respectively. On the other hand, the space for the pressure is only defined in the fluid domain, $\Omega_{\rm F}$, and is given by
Before creating the FE spaces, we first need to define the discrete model of each of the sub-domains, which are constructed restricting the global discrete model to the corresponding part. This is done by calling the DiscreteModel
function with the desired geometrical part label, i.e. "solid"
and "fluid"
, respectively.
model_solid = DiscreteModel(model,"solid")
model_fluid = DiscreteModel(model,"fluid")
In what follows we will assume a second-order veloticty interpolation, i.e. $k=2$
k = 2
Having set up all the ingredients, we can create the different FE spaces for the test functions. For the velocity FE spaces we call the TestFESpace
function with the corresponding discrete model, using a 2-dimensional VectorValue
type, a :Lagrangian
reference FE element of order k
and conformity :H1
. Note that we assign different Dirichlet boundary labels for the two different parts, generating the variational spaces with homogeneous Dirichlet boundary conditions, $V_{\rm F,0}$ and $V_{\rm S,0}$ .
Vf = TestFESpace(
model=model_fluid,
valuetype=VectorValue{2,Float64},
reffe=:Lagrangian,
order=k,
conformity =:H1,
dirichlet_tags=["inlet", "noSlip", "cylinder"])
Vs = TestFESpace(
model=model_solid,
valuetype=VectorValue{2,Float64},
reffe=:Lagrangian,
order=k,
conformity =:H1,
dirichlet_tags=["fixed"])
For the pressure test FE space, we use the fluid discrete model, a scalar value type, a :Lagrangian
reference FE element of order k-1
and :C0
conformity.
Qf = TestFESpace(
model=model_fluid,
valuetype=Float64,
order=k-1,
reffe=:Lagrangian,
conformity=:C0)
The trial FE spaces are generated from the test FE spaces, adding the corresponding function for the various Dirichlet boundaries, leading to $U_{\rm F,g_{\rm F}}$, $U_{\rm S,g_{\rm S}}$ and $P_{\rm F}$.
Uf = TrialFESpace(Vf,[uf_in, uf_0, uf_0])
Us = TrialFESpace(Vs,[us_0])
Pf = TrialFESpace(Qf)
Finally, we glue the test and trial FE spaces together, defining a unique test and trial space for all the fields using the MultiFieldFESpace
function. That is $Y=[V_{\rm S,0}, V_{\rm F,0}, Q_{\rm F}]^T$ and $X=[U_{\rm S,g_{\rm S}}, U_{\rm F,g_{\rm F}}, P_{\rm F}]^T$
Y = MultiFieldFESpace([Vs,Vf,Qf])
X = MultiFieldFESpace([Us,Uf,Pf])
Weak form
We now introduce the solution and test function vectors as $\mathbf{x}^h=[\mathbf{u}^h_{\rm S},\mathbf{u}^h_{\rm F}, p^h_{\rm F}]^T$ and $\mathbf{y}^h=[\mathbf{v}^h_{\rm S},\mathbf{v}^h_{\rm F}, q^h_{\rm F}]^T$. The weak form of the coupled FSI problem using the Nitche's method, see for instance [2], reads: find $\mathbf{x}^h \in X$ such that
with the following definitions:
- \[a_s(\mathbf{x}^h,\mathbf{y}^h)\]is the bilinear form associated with the solid counterpart, defined as
function a_s(x,y)
us,uf,p = x
vs,vf,q = y
ε(vs) ⊙ σ_s(ε(us))
end
- \[a_f(\mathbf{x}^h,\mathbf{y}^h)\]is the bilinear form associated with the fluid counterpart, defined as
function a_f(x,y)
us,uf,p = x
vs,vf,q = y
(ε(vf) ⊙ σ_dev_f(ε(uf))) - (∇⋅vf)*p + q*(∇⋅uf)
end
- \[a_{fs}(\mathbf{x}^h,\mathbf{y}^h)\]is the bilinear form associated with the coupling between fluid and solid counterparts. To difine this form we use the well known Nitsche's method, which enforces the continuity of fluid and solid velocities as well as the continuity of the normal stresses, see for instance [2]. The final expression for this term reads:
Where $\chi$ is a parameter that can take values $1.0$ or $-1.0$ and it is used to define the symmetric or antisymmetric version of the method, respectively.
const γ = 1.0
const χ = -1.0
function a_fs(x,y)
us_Γ, uf_Γ, p_Γ = x
vs_Γ, vf_Γ, q_Γ = y
us, uf, p = us_Γ.⁻, uf_Γ.⁺, p_Γ.⁺
vs, vf, q = vs_Γ.⁻, vf_Γ.⁺, q_Γ.⁺
α*(vf-vs)⋅(uf-us) + (vf-vs)⋅(p*n_Γfs-n_Γfs⋅σ_dev_f(ε(uf))) + (q*n_Γfs-χ*n_Γfs⋅σ_dev_f(ε(vf)))⋅(uf-us)
end
- \[l_s(\mathbf{y}^h)\]is the linear form associated with the solid counterpart, defined as
function l_s(y)
vs,vf,q = y
vs⋅s
end
- \[l_f(\mathbf{y}^h)\]is the linear form associated with the fluid counterpart, defined as
function l_f(y)
vs,vf,q = y
vf⋅f + q*g
end
- \[l_{f,\Gamma_N}(\mathbf{y}^h)\]is the linear form associated with the fluid Neumann boundary condition, defined as
function l_f_Γn(y)
vs,vf,q = y
vf⋅hN
end
Numerical integration
To define the quadrature rules used in the numerical integration of the different terms, we first need to generate the domain triangulation. Here we create the triangulation of the global domain, $\mathcal{T}$.
trian = Triangulation(model)
The solid and fluid triangulations, $\mathcal{T}_{\rm F}$ and $\mathcal{T}_{\rm S}$, are constructed from the discrete models restricted to the respective parts.
trian_solid = Triangulation(model_solid)
trian_fluid = Triangulation(model_fluid)
We also generate the triangulation and associated outer normal field for the outlet boundary, $\Gamma_{\rm F,N_{out}}$, which will be used to impose a Neumann condition.
trian_Γout = BoundaryTriangulation(model,"outlet")
n_Γout = get_normal_vector(trian_Γout)
Finally, to impose the interface condition between solid and fluid, we will need the triangulation and normal field of such interface, $\Gamma_{\rm FS}$. The interface triangulation is generated by calling the InterfaceTriangulation
function specifying the two interfacing domain models. Note that the normal field will point outwards with respect to the first entry.
trian_Γfs = InterfaceTriangulation(model_fluid,model_solid)
n_Γfs = get_normal_vector(trian_Γfs)
From the interface triangulation we can obtain the interface elements length, $h$, and the penalty parameter, $\alpha=\gamma\frac{\mu_f}{h}$, used in the Nitsche's terms.
using Gridap.Arrays
using LinearAlgebra: norm
xe = get_cell_coordinates(trian_Γfs)
he = apply(x->norm(x[2]-x[1]),xe)
α = apply(h->γ*μ_f/h, he)
Once we have all the triangulations, we can generate the quadrature rules to be applied each domain.
degree = 2*k
quad_solid = CellQuadrature(trian_solid,degree)
quad_fluid = CellQuadrature(trian_fluid,degree)
Idem for the boundary part.
bdegree = 2*k
quad_Γout = CellQuadrature(trian_Γout,bdegree)
Idem for the interface part.
idegree = 2*k
quad_Γfs = CellQuadrature(trian_Γfs,idegree)
Algebraic System of Equations
After defining the weak form of the problem and the integration quadrature rules to perform the numerical integration, we are ready to assemble the linear system of equations. In this case, the system will have the following structure:
In order to construct the system we first define the diferent discrete terms using the functions AffineFETerm
, that assemble contributions on the left-hand side and the right-hand side of the system, for the fluid and solid interior terms evaluated in the respective triangulations, $\mathcal{T}_{\rm F}$ and $\mathcal{T}_{\rm S}$. The fluid Neumann boundary condition is assembled using the function FESource
, which only affects the right-hand side of the system, evaluating the corresponding form on the triangulation of the outlet boundary, $\Gamma_{\rm F,N_{out}}$. The coupling terms are evaluated at the interface $\Gamma_{\rm FS}$ using the function LinearFETerm
, assembling only to the left-hand side of the system.
t_Ω_s= AffineFETerm(a_s,l_s,trian_solid,quad_solid)
t_Ω_f = AffineFETerm(a_f,l_f,trian_fluid,quad_fluid)
t_Γfs = LinearFETerm(a_fs,trian_Γfs,quad_Γfs)
t_Γn_f = FESource(l_f_Γn,trian_Γout,quad_Γout)
The final FE operator is constructed using the function AffineFEOperator
and passing as arguments the trial and test FE spaces, $X$ and $Y$, and all the different terms previously defined.
op = AffineFEOperator(X,Y,t_Ω_s,t_Ω_f,t_Γn_f,t_Γfs)
Finally, we call solve
to obtain the solution vector of nodal values $[\mathbf{U}^h_{\rm S},\mathbf{U}^h_{\rm F},\mathbf{P}^h_{\rm F}]^T$
uhs, uhf, ph = solve(op)
Post-processing
Visualization
The solution fields $[\mathbf{U}^h_{\rm S},\mathbf{U}^h_{\rm F},\mathbf{P}^h_{\rm F}]^T$ are defined over all the domain, extended with zeros on the inactive part. Calling the function writevtk
passing the global triangulation, we will output the global fields.
writevtk(trian,"trian", cellfields=["uhs" => uhs, "uhf" => uhf, "ph" => ph])
However, we can also restrict the fields to the active part by calling the function restrict
with the field along with the respective active triangulation.
uhs_solid = restrict(uhs, trian_solid)
uhf_fluid = restrict(uhf, trian_fluid)
ph_fluid = restrict(ph, trian_fluid)
writevtk(trian_solid,"trian_solid",cellfields=["uhs"=>uhs_solid])
writevtk(trian_fluid,"trian_fluid",cellfields=["ph"=>ph_fluid,"uhf"=>uhf_fluid])
Quantities of Interest
trian_ΓS = BoundaryTriangulation(model,["cylinder","interface"])
quad_ΓS = CellQuadrature(trian_ΓS,bdegree)
n_ΓS = get_normal_vector(trian_ΓS)
uh_ΓS = restrict(uhf_fluid,trian_ΓS)
ph_ΓS = restrict(ph_fluid,trian_ΓS)
FD, FL = sum( integrate( (n_ΓS⋅σ_dev_f(ε(uh_ΓS))) - ph_ΓS*n_ΓS, trian_ΓS, quad_ΓS ) )
println("Drag force: ", FD)
println("Lift force: ", FL)
References
[1] Turek, S., Hron, J., Madlik, M., Razzaq, M., Wobker, H., & Acker, J. F. (2011).* Numerical simulation and benchmarking of a monolithic multigrid solver for fluid-structure interaction problems with application to hemodynamics. In Fluid Structure Interaction II (pp. 193-220). Springer, Berlin, Heidelberg.
[2] Burman, E., and Fernández, M. A. Stabilized explicit coupling for fluid–structure interaction using Nitsche's method. Comptes Rendus Mathematique 345.8 (2007): 467-472.
This page was generated using Literate.jl.