/
index.html
100 lines (99 loc) · 32.7 KB
/
index.html
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>11 Fluid-Structure Interaction · Gridap tutorials</title><link href="https://fonts.googleapis.com/css?family=Lato|Roboto+Mono" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.11.2/css/fontawesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.11.2/css/solid.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.11.2/css/brands.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/KaTeX/0.11.1/katex.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.3.6/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-dark.css" data-theme-name="documenter-dark"/><link class="docs-theme-link" rel="stylesheet" type="text/css" href="../../assets/themes/documenter-light.css" data-theme-name="documenter-light" data-theme-primary/><script src="../../assets/themeswap.js"></script></head><body><div id="documenter"><nav class="docs-sidebar"><a class="docs-logo" href="../../"><img src="../../assets/logo.png" alt="Gridap tutorials logo"/></a><div class="docs-package-name"><span class="docs-autofit">Gridap tutorials</span></div><form class="docs-search" action="../../search/"><input class="docs-search-query" id="documenter-search-query" name="q" type="text" placeholder="Search docs"/></form><ul class="docs-menu"><li><a class="tocitem" href="../../">Introduction</a></li><li><a class="tocitem" href="../t001_poisson/">1 Poisson equation</a></li><li><a class="tocitem" href="../t002_validation/">2 Code validation</a></li><li><a class="tocitem" href="../t003_elasticity/">3 Linear elasticity</a></li><li><a class="tocitem" href="../t004_p_laplacian/">4 p-Laplacian</a></li><li><a class="tocitem" href="../t005_hyperelasticity/">5 Hyper-elasticity</a></li><li><a class="tocitem" href="../t006_dg_discretization/">6 Poisson equation (with DG)</a></li><li><a class="tocitem" href="../t007_darcy/">7 Darcy equation (with RT)</a></li><li><a class="tocitem" href="../t008_inc_navier_stokes/">8 Incompressible Navier-Stokes</a></li><li><a class="tocitem" href="../t009_stokes/">9 Stokes equation</a></li><li><a class="tocitem" href="../t010_isotropic_damage/">10 Isotropic damage model</a></li><li class="is-active"><a class="tocitem" href>11 Fluid-Structure Interaction</a><ul class="internal"><li><a class="tocitem" href="#Problem-statement-1"><span>Problem statement</span></a></li><li><a class="tocitem" href="#Numerical-scheme-1"><span>Numerical scheme</span></a></li><li><a class="tocitem" href="#Post-processing-1"><span>Post-processing</span></a></li><li><a class="tocitem" href="#References-1"><span>References</span></a></li></ul></li></ul><div class="docs-version-selector field has-addons"><div class="control"><span class="docs-label button is-static is-size-7">Version</span></div><div class="docs-selector control is-expanded"><div class="select is-fullwidth is-size-7"><select id="documenter-version-selector"></select></div></div></div></nav><div class="docs-main"><header class="docs-navbar"><nav class="breadcrumb"><ul class="is-hidden-mobile"><li class="is-active"><a href>11 Fluid-Structure Interaction</a></li></ul><ul class="is-hidden-tablet"><li class="is-active"><a href>11 Fluid-Structure Interaction</a></li></ul></nav><div class="docs-right"><a class="docs-edit-link" href="https://github.com/gridap/Tutorials/blob/master/src/fsi_tutorial.jl" title="Edit on GitHub"><span class="docs-icon fab"></span><span class="docs-label is-hidden-touch">Edit on GitHub</span></a><a class="docs-settings-button fas fa-cog" id="documenter-settings-button" href="#" title="Settings"></a><a class="docs-sidebar-button fa fa-bars is-hidden-desktop" id="documenter-sidebar-button" href="#"></a></div></header><article class="content" id="documenter-page"><h1 id="Tutorial-11:-Fluid-Structure-Interaction-1"><a class="docs-heading-anchor" href="#Tutorial-11:-Fluid-Structure-Interaction-1">Tutorial 11: Fluid-Structure Interaction</a><a class="docs-heading-anchor-permalink" href="#Tutorial-11:-Fluid-Structure-Interaction-1" title="Permalink"></a></h1><p><a href="https://mybinder.org/v2/gh/gridap/Tutorials/gh-pages?filepath=v0.15.0/notebooks/t011_fsi_tutorial.ipynb"><img src="https://mybinder.org/badge_logo.svg" alt/></a> <a href="https://nbviewer.jupyter.org/github/gridap/Tutorials/blob/gh-pages/v0.15.0/notebooks/t011_fsi_tutorial.ipynb"><img src="https://img.shields.io/badge/show-nbviewer-579ACA.svg" alt/></a></p><p>In this tutorial, we will learn</p><ul><li>How to solve a surface-coupled multi-physics problem.</li><li>Construct FE spaces defined in different domains.</li><li>Define interface triangulations and integrate on them.</li></ul><ol><li><a href="#probStat">Problem Statement</a><ul><li><a href="#strongForm">Strong form</a></li><li><a href="#geometry">Geometry and Discrete model</a></li><li><a href="#conditions">Boundary conditions and properties</a></li></ul></li><li><a href="#numericalScheme">Numerical scheme</a><ul><li><a href="#feSpace">FE spaces</a></li><li><a href="#weakForm">Weak form</a></li><li><a href="#integration">Numerical integration</a></li><li><a href="#algebraic">Algebraic system of equations</a></li></ul></li><li><a href="#postprocess">Post-processing</a><ul><li><a href="#visualization">Visualization</a></li><li><a href="#QOIs">Quantities of interest</a></li></ul></li></ol><h2 id="Problem-statement-1"><a class="docs-heading-anchor" href="#Problem-statement-1">Problem statement</a><a class="docs-heading-anchor-permalink" href="#Problem-statement-1" title="Permalink"></a></h2><h3 id="Strong-form-1"><a class="docs-heading-anchor" href="#Strong-form-1">Strong form</a><a class="docs-heading-anchor-permalink" href="#Strong-form-1" title="Permalink"></a></h3><p>Let <span>$\Gamma_{\rm FS}$</span> be the interface between a fluid domain <span>$\Omega_{\rm F}$</span> and a solid domain <span>$\Omega_{\rm S}$</span>. We denote by <span>$\Gamma_{\rm F,D}$</span> and <span>$\Gamma_{\rm F,N}$</span> the fluid boundaries with Dirichlet and Neumann conditions, respectively. The Fluid-Structure Interaction (FSI) problem reads:</p><p>find <span>$u_{\rm F}$</span>, <span>$p_{\rm F}$</span> and <span>$u_{\rm S}$</span> such that</p><div>\[\left\lbrace
\begin{aligned}
-\nabla\cdot\boldsymbol{\sigma}_{\rm F} = f &\text{ in }\Omega_{\rm F},\\
\nabla\cdot u_{\rm F} = 0 &\text{ in } \Omega_{\rm F},\\
-\nabla\cdot\boldsymbol{\sigma}_{\rm S} = s &\text{ in }\Omega_{\rm S},\\
\end{aligned}
\right.\]</div><p>satisfying the Dirichlet and Neumann boundary conditions</p><div>\[\left\lbrace
\begin{aligned}
u_{\rm F} = g_{\rm F} &\text{ on } \Gamma_{\rm F,D},\\
n_{\rm F} \cdot \boldsymbol{\sigma}_{\rm F} = 0 &\text{ on } \Gamma_{\rm F,N},\\
u_{\rm S} = g_{\rm S} &\text{ on } \Gamma_{\rm S,D},\\
\end{aligned}
\right.\]</div><p>and the kinematic and dynamic conditions at the fluid-solid interface</p><div>\[\left\lbrace
\begin{aligned}
u_{\rm F} = u_{\rm S} &\text{ on } \Gamma_{\rm FS},\\
n_{\rm F} \cdot \boldsymbol{\sigma}_{\rm F} + n_{\rm S} \cdot \boldsymbol{\sigma}_{\rm S} = 0 &\text{ on } \Gamma_{\rm FS}.\\
\end{aligned}
\right.\]</div><p>Where <span>$\boldsymbol{\sigma}_{\rm F}(u_{\rm F},p_{\rm F})=2\mu_{\rm F}\boldsymbol{\varepsilon}(u_{\rm F}) - p_{\rm F}\mathbf{I}$</span> and <span>$\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}$</span>.</p><h3 id="Geometry-and-Discrete-model-1"><a class="docs-heading-anchor" href="#Geometry-and-Discrete-model-1">Geometry and Discrete model</a><a class="docs-heading-anchor-permalink" href="#Geometry-and-Discrete-model-1" title="Permalink"></a></h3><p>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 <span>$\Omega \doteq (0,4.5)\times(0,0.41)$</span>, with an embedded cylinder of radius <span>$R=0.05$</span> and center at <span>$C=(0.2,0.2)$</span>. The associated FE triangulation is denoted by <span>$\mathcal{T}$</span>, the fluid and solid domain and their associated triangulations will be denoted by <span>$\Omega_{\rm F}$</span>, <span>$\Omega_{\rm S}$</span>, <span>$\mathcal{T}_{\rm F}$</span> and <span>$\mathcal{T}_{\rm S}$</span>, respectively.</p><p>In order to load the discrete model we first setup Gridap</p><pre><code class="language-julia">using Gridap</code></pre><p>The discrete model for the elastic flag problem is generated by loading the <code>"../models/elasticFlag.json"</code> file.</p><pre><code class="language-julia">model = DiscreteModelFromFile("../models/elasticFlag.json")</code></pre><p>We can inspect the loaded geometry and associated parts by printing to a <code>vtk</code> file:</p><pre><code class="language-julia">writevtk(model,"model")</code></pre><p>This will produce an output in which we can identify the different parts of the domain, with the associated labels and tags.</p><table><tr><th style="text-align: center">Part</th><th style="text-align: center">Notation</th><th style="text-align: center">Label</th><th style="text-align: center">Tag</th></tr><tr><td style="text-align: center">Solid-cylinder wall</td><td style="text-align: center"><span>$\Gamma_{\rm S,D_{cyl}}$</span></td><td style="text-align: center">"fixed"</td><td style="text-align: center">1</td></tr><tr><td style="text-align: center">Fluid-solid interface</td><td style="text-align: center"><span>$\Gamma_{\rm FS}$</span></td><td style="text-align: center">"interface"</td><td style="text-align: center">2</td></tr><tr><td style="text-align: center">Channel inlet</td><td style="text-align: center"><span>$\Gamma_{\rm F,D_{in}}$</span></td><td style="text-align: center">"inlet"</td><td style="text-align: center">3</td></tr><tr><td style="text-align: center">Channel outlet</td><td style="text-align: center"><span>$\Gamma_{\rm F,N_{out}}$</span></td><td style="text-align: center">"outlet"</td><td style="text-align: center">4</td></tr><tr><td style="text-align: center">Channel walls</td><td style="text-align: center"><span>$\Gamma_{\rm F,D_{wall}}$</span></td><td style="text-align: center">"noSlip"</td><td style="text-align: center">5</td></tr><tr><td style="text-align: center">Fluid-cylinder wall</td><td style="text-align: center"><span>$\Gamma_{\rm F,D_{cyl}}$</span></td><td style="text-align: center">"cylinder"</td><td style="text-align: center">6</td></tr><tr><td style="text-align: center">Fluid domain</td><td style="text-align: center"><span>$\Omega_{\rm F}$</span></td><td style="text-align: center">"fluid"</td><td style="text-align: center">7</td></tr><tr><td style="text-align: center">Solid domain</td><td style="text-align: center"><span>$\Omega_{\rm S}$</span></td><td style="text-align: center">"solid"</td><td style="text-align: center">8</td></tr></table><p><img src="../../assets/fsi/tags.png" alt/></p><h3 id="External-conditions-and-properties-1"><a class="docs-heading-anchor" href="#External-conditions-and-properties-1">External conditions and properties</a><a class="docs-heading-anchor-permalink" href="#External-conditions-and-properties-1" title="Permalink"></a></h3><h4 id="Boundary-conditions-1"><a class="docs-heading-anchor" href="#Boundary-conditions-1">Boundary conditions</a><a class="docs-heading-anchor-permalink" href="#Boundary-conditions-1" title="Permalink"></a></h4><p>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.</p><div>\[\left\lbrace
\begin{aligned}
u_{\rm F,in}(x,y) = 1.5U\frac{y(H −y)}{\left(\frac{H}{2}\right)^2}\quad&\textrm{on }\Gamma_{\rm F,D_{in}},\\
u_{\rm F,0}(x,y) = [0, 0]\quad&\textrm{on }\Gamma_{\rm F,D_{wall}}\cup\Gamma_{\rm F,D_{cyl}},\\
u_{\rm S,0}(x,y) = [0, 0]\quad&\textrm{on }\Gamma_{\rm S,D_{cyl}}.\\
\end{aligned}
\right.\]</div><pre><code class="language-julia">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 )</code></pre><p>We consider a free tranction condition at the channel outlet</p><div>\[n_{\rm F} \cdot \boldsymbol{\sigma}_{\rm F} = \mathbf{0}\quad\textrm{on }\Gamma_{\rm F,N}\]</div><pre><code class="language-julia">hN(x) = VectorValue( 0.0, 0.0 )
p_jump(x) = 0.0</code></pre><h4 id="External-forces-1"><a class="docs-heading-anchor" href="#External-forces-1">External forces</a><a class="docs-heading-anchor-permalink" href="#External-forces-1" title="Permalink"></a></h4><p>In this test, the body forces acting on the fluid an solid are zero.</p><pre><code class="language-julia">f(x) = VectorValue( 0.0, 0.0 )
s(x) = VectorValue( 0.0, 0.0 )
g(x) = 0.0</code></pre><h4 id="Material-properties-1"><a class="docs-heading-anchor" href="#Material-properties-1">Material properties</a><a class="docs-heading-anchor-permalink" href="#Material-properties-1" title="Permalink"></a></h4><p>We use a linear elastic constitutive law for the elastic flag. Given the Young's modulus <span>$E$</span> and the Poisson ratio <span>$\nu$</span>, we can compute the Lamé constants, <span>$\lambda$</span> and <span>$\mu$</span>, using the following function:</p><pre><code class="language-julia">function lame_parameters(E,ν)
λ = (E*ν)/((1+ν)*(1-2*ν))
μ = E/(2*(1+ν))
(λ, μ)
end</code></pre><p>Then, we get the Lamé parameters for a solid with <span>$E=1.0$</span> MPa and <span>$\nu=0.33$</span>.</p><pre><code class="language-julia">const E_s = 1.0
const ν_s = 0.33
const (λ_s,μ_s) = lame_parameters(E_s,ν_s)</code></pre><p>The Cauchy stress tensor for the solid part is defined by <span>$\sigma_s = 2\mu\varepsilon + \lambda tr(\varepsilon)\mathbf{I}$</span>. Note that we use the trace operator from the <code>LinearAlgebra</code> package. Note that this function will be used as a composition (<code>∘</code>), using as argument a function whose arguments depend on the coordinates, without the need of passing such coordinates as an argument. That is <code>σ_s(ε(u)) = σ_s ∘ ε(u)</code>.</p><pre><code class="language-julia">using LinearAlgebra: tr
σₛ(ε) = λ_s*tr(ε)*one(ε) + 2*μ_s*ε</code></pre><p>For the fluid part, we only need to define the viscosity <span>$\mu_f$</span>, which we set to <span>$0.01$</span>.</p><pre><code class="language-julia">const μ_f = 0.01</code></pre><p>The Cauchy stress tensor for the fluid part is given by <span>$\sigma_f = \sigma^{dev}_f - p\mathbf{I}$</span>, with <span>$\sigma^{dev}_f=2\mu_f$</span> the deviatoric part of the stress. Since we use a mixed form with the pressure <span>$p$</span> as an unknown, the stress law will only involve the deviatoric part.</p><pre><code class="language-julia">σ_dev_f(ε) = 2*μ_f*ε</code></pre><h2 id="Numerical-scheme-1"><a class="docs-heading-anchor" href="#Numerical-scheme-1">Numerical scheme</a><a class="docs-heading-anchor-permalink" href="#Numerical-scheme-1" title="Permalink"></a></h2><h3 id="FE-Spaces-1"><a class="docs-heading-anchor" href="#FE-Spaces-1">FE Spaces</a><a class="docs-heading-anchor-permalink" href="#FE-Spaces-1" title="Permalink"></a></h3><p>In this tutorial we use an <em>inf-sup</em> stable velocity-pressure pair, <span>$P_k/P_{k-1}$</span> elements, with continuous velocities and pressures. We select the same velocity interpolation space for the fluid and the solid parts, defined as follows</p><div>\[V_{\rm X} \doteq \{ v \in [H^1(\Omega_{\rm X})]^d:\ v|_T\in [P_k(T)]^d \text{ for all } T\in\mathcal{T}_{\rm X} \},\]</div><p>where <span>$T$</span> denotes an arbitrary cell of the FE mesh <span>$\mathcal{T}_{\rm X}$</span>, and <span>$P_k(T)$</span> is the usual Lagrangian FE space of order <span>$k$</span> defined on a mesh of triangles or tetrahedra. Note that the sub-index <span>$(\cdot)_{\rm X}$</span> stands for the fluid or solid parts, <span>$(\cdot)_{\rm F}$</span> or <span>$(\cdot)_{\rm S}$</span>, respectively. On the other hand, the space for the pressure is only defined in the fluid domain, <span>$\Omega_{\rm F}$</span>, and is given by</p><div>\[Q \doteq \{ q \in C^0(\Omega):\ q|_T\in P_{k-1}(T) \text{ for all } T\in\mathcal{T}_{\rm F}\}.\]</div><p>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 <code>DiscreteModel</code> function with the desired geometrical part label, i.e. <code>"solid"</code> and <code>"fluid"</code>, respectively.</p><pre><code class="language-julia">model_solid = DiscreteModel(model,tags="solid")
model_fluid = DiscreteModel(model,tags="fluid")</code></pre><p>In what follows we will assume a second-order veloticty interpolation, i.e. <span>$k=2$</span></p><pre><code class="language-julia">k = 2</code></pre><p>Now we define the reference FE for the velocity and pressure fields. The velocity field reference FE, both for fluid and solid domains, will be defined by a 2-dimensional <code>VectorValue</code> type <code>:Lagrangian</code> reference FE element of order <code>k</code>. The reference FE for the pressure will be given by a scalar value type <code>:Lagrangian</code> reference FE element of order <code>k-1</code>.</p><pre><code class="language-julia">reffeᵤ = ReferenceFE(lagrangian,VectorValue{2,Float64},k)
reffeₚ = ReferenceFE(lagrangian,Float64,k-1)</code></pre><p>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 <code>TestFESpace</code> function with the corresponding discrete model, using the velocity reference FE <code>reffeᵤ</code> and conformity <code>:H1</code>. Note that we assign different Dirichlet boundary labels for the two different parts, generating the variational spaces with homogeneous Dirichlet boundary conditions, <span>$V_{\rm F,0}$</span> and <span>$V_{\rm S,0}$</span> .</p><pre><code class="language-julia">Vf = TestFESpace(
model_fluid,
reffeᵤ,
conformity=:H1,
dirichlet_tags=["inlet", "noSlip", "cylinder"])
Vs = TestFESpace(
model_solid,
reffeᵤ,
conformity=:H1,
dirichlet_tags=["fixed"])</code></pre><p>For the pressure test FE space, we use the fluid discrete model, the pressure reference FE <code>reffeₚ</code> and <code>:C0</code> conformity.</p><pre><code class="language-julia">Qf = TestFESpace(
model_fluid,
reffeₚ,
conformity=:C0)</code></pre><p>The trial FE spaces are generated from the test FE spaces, adding the corresponding function for the various Dirichlet boundaries, leading to <span>$U_{\rm F,g_{\rm F}}$</span>, <span>$U_{\rm S,g_{\rm S}}$</span> and <span>$P_{\rm F}$</span>.</p><pre><code class="language-julia">Uf = TrialFESpace(Vf,[uf_in, uf_0, uf_0])
Us = TrialFESpace(Vs,[us_0])
Pf = TrialFESpace(Qf)</code></pre><p>Finally, we glue the test and trial FE spaces together, defining a unique test and trial space for all the fields using the <code>MultiFieldFESpace</code> function. That is <span>$Y=[V_{\rm S,0}, V_{\rm F,0}, Q_{\rm F}]^T$</span> and <span>$X=[U_{\rm S,g_{\rm S}}, U_{\rm F,g_{\rm F}}, P_{\rm F}]^T$</span></p><pre><code class="language-julia">Y = MultiFieldFESpace([Vs,Vf,Qf])
X = MultiFieldFESpace([Us,Uf,Pf])</code></pre><h3 id="Numerical-integration-1"><a class="docs-heading-anchor" href="#Numerical-integration-1">Numerical integration</a><a class="docs-heading-anchor-permalink" href="#Numerical-integration-1" title="Permalink"></a></h3><p>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, <span>$\mathcal{T}$</span>.</p><pre><code class="language-julia">Ω = Triangulation(model)</code></pre><p>The solid and fluid triangulations, <span>$\mathcal{T}_{\rm F}$</span> and <span>$\mathcal{T}_{\rm S}$</span>, are constructed from the discrete models restricted to the respective parts.</p><pre><code class="language-julia">Ω_s = Triangulation(model_solid)
Ω_f = Triangulation(model_fluid)</code></pre><p>We also generate the triangulation and associated outer normal field for the outlet boundary, <span>$\Gamma_{\rm F,N_{out}}$</span>, which will be used to impose a Neumann condition.</p><pre><code class="language-julia">Γ_out = BoundaryTriangulation(model,tags="outlet")
n_Γout = get_normal_vector(Γ_out)</code></pre><p>Finally, to impose the interface condition between solid and fluid, we will need the triangulation and normal field of such interface, <span>$\Gamma_{\rm FS}$</span>. The interface triangulation is generated by calling the <code>InterfaceTriangulation</code> function specifying the two interfacing domain models. Note that the normal field will point outwards with respect to the first entry.</p><pre><code class="language-julia">Γ_fs = InterfaceTriangulation(model_fluid,model_solid)
n_Γfs = get_normal_vector(Γ_fs)</code></pre><p>Once we have all the triangulations, we can generate the quadrature rules to be applied each domain. This will be generated by calling the <code>Measure</code> function, that given a triangulation and an integration degree, it returns the Lebesgue integral measure <span>$d\Omega$</span>.</p><pre><code class="language-julia">degree = 2*k
dΩₛ = Measure(Ω_s,degree)
dΩ_f = Measure(Ω_f,degree)</code></pre><p>Idem for the boundary part.</p><pre><code class="language-julia">bdegree = 2*k
dΓ_out = Measure(Γ_out,bdegree)</code></pre><p>Idem for the interface part.</p><pre><code class="language-julia">idegree = 2*k
dΓ_fs = Measure(Γ_fs,idegree)</code></pre><h3 id="Weak-form-1"><a class="docs-heading-anchor" href="#Weak-form-1">Weak form</a><a class="docs-heading-anchor-permalink" href="#Weak-form-1" title="Permalink"></a></h3><p>We now introduce the solution and test function vectors as <span>$[\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm s}, p^h_{\rm f}]^T$</span> and <span>$[\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}]^T$</span>. The weak form of the coupled FSI problem using the Nitche's method, see for instance [2], reads:</p><p>find <span>$[\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm s}, p^h_{\rm f}]^T \in X$</span> such that</p><div>\[a([\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm s}, p^h_{\rm f}],[\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}])=l([\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}])\qquad\forall[\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}]^T\in Y,\]</div><p>where</p><div>\[\begin{aligned}
&a([\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm s}, p^h_{\rm f}],[\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}])\doteq a_s(\mathbf{u}^h_{\rm s},\mathbf{v}^h_{\rm s}) + a_f((\mathbf{u}^h_{\rm f},p^h_{\rm f}),(\mathbf{v}^h_{\rm f},q^h_{\rm f})) + a_{fs}((\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm f}, p^h_{\rm f}),(\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}))\\
&l([\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}])\doteq l_s(\mathbf{v}^h_{\rm s}) + l_f(\mathbf{v}^h_{\rm f}) + l_{f,\Gamma_N}(\mathbf{v}^h_{\rm f}).
\end{aligned}\]</div><pre><code class="language-julia">a((us,uf,p),(vs,vf,q)) = a_s(us,vs) + a_f((uf,p),(vf,q)) + a_fs((us,uf,p),(vs,vf,q))
l((vs,vf,q)) = l_s(vs) + l_f((vf,q)) + l_f_Γn(vf)</code></pre><p>with the following definitions:</p><ul><li>The bilinear form associated with the solid counterpart, <span>$a_s(\mathbf{u}^h_{\rm s},\mathbf{v}^h_{\rm s})$</span>, defined as</li></ul><div>\[a_s(\mathbf{u},\mathbf{v})=\int\varepsilon(\mathbf{v}):\boldsymbol{\sigma}_s(\epsilon(\mathbf{u}))\ d\Omega_{\rm S}.\]</div><pre><code class="language-julia">a_s(u,v) = ∫( ε(v) ⊙ (σₛ ∘ ε(u)) )dΩₛ</code></pre><ul><li>The bilinear form associated with the fluid counterpart, <span>$a_f((\mathbf{u}^h_{\rm f},p^h_{\rm f}),(\mathbf{v}^h_{\rm f},q^h_{\rm f}))$</span>, is defined as</li></ul><div>\[a_f((\mathbf{u},p),(\mathbf{v},q))=\int\left(\varepsilon(\mathbf{v}):\boldsymbol{\sigma}^{dev}_f(\mathbf{u}) - (\nabla\cdot\mathbf{v})\ p + q\ (\nabla\cdot\mathbf{u})\right) d\Omega_{\rm F}.\]</div><pre><code class="language-julia">a_f((u,p),(v,q)) = ∫( ε(v) ⊙ (σ_dev_f ∘ ε(u)) - (∇⋅v)*p + q*(∇⋅u) )dΩ_f</code></pre><ul><li>The bilinear form associated with the coupling between fluid and solid counterparts, <span>$a_{fs}((\mathbf{u}^h_{\rm s},\mathbf{u}^h_{\rm f}, p^h_{\rm f}),(\mathbf{v}^h_{\rm s},\mathbf{v}^h_{\rm f}, q^h_{\rm f}))$</span>, is given by:</li></ul><div>\[\begin{aligned}
a_{fs}((\mathbf{u}_{\rm s},\mathbf{u}_{\rm f}, p),((\mathbf{v}_{\rm s},\mathbf{v}_{\rm f}, q))=&\int\left(\gamma\frac{\mu_f}{h}(\mathbf{v}_{\rm f}-\mathbf{v}_{\rm s})\cdot(\mathbf{u}_{\rm f}-\mathbf{u}_{\rm s})\right.\\
&- (\mathbf{v}_{\rm f}-\mathbf{v}_{\rm s})\cdot(\mathbf{n}_{\rm fs}\cdot\boldsymbol{\sigma}^{dev}_f(\mathbf{u}_{\rm f}))
+ (\mathbf{v}_{\rm f}-\mathbf{v}_{\rm S})\cdot(p_{\rm f}\mathbf{n}_{\rm fs})\\
&- \chi(\mathbf{n}_{\rm fs}\cdot\boldsymbol{\sigma}^{dev}_f(\mathbf{v}_{\rm f})\cdot(\mathbf{u}^{\rm f}-\mathbf{u}_{\rm s})
+ \left.(q_{\rm f}\mathbf{n}_{\rm fs})\cdot(\mathbf{u}_{\rm f}-\mathbf{u}_{\rm s})\ \right) d\Gamma_{\rm FS}
\end{aligned}\]</div><pre><code class="language-julia">a_fs((us,uf,p),(vs,vf,q)) = ∫( α*(jump_Γ(vf,vs)⋅jump_Γ(uf,us)) +
jump_Γ(vf,vs)⋅t_Γfs(1,p,uf,n_Γfs) +
t_Γfs(χ,q,vf,n_Γfs)⋅jump_Γ(uf,us) )dΓ_fs</code></pre><p>Where <span>$\chi$</span> is a parameter that can take values <span>$1.0$</span> or <span>$-1.0$</span> and it is used to define the symmetric or antisymmetric version of the method, respectively. To difine this form we used 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].</p><pre><code class="language-julia">const γ = 1.0
const χ = -1.0
jump_Γ(wf,ws) = wf.⁺-ws.⁻
t_Γfs(χ,q,w,n) = q.⁺*n.⁺ - χ*n.⁺⋅(σ_dev_f∘ε(w.⁺))</code></pre><p>From the interface triangulation we can obtain the interface elements length, <span>$h$</span>, and the penalty parameter, <span>$\alpha=\gamma\frac{\mu_f}{h}$</span>, used in the Nitsche's terms.</p><pre><code class="language-julia">using Gridap.CellData
dim = num_cell_dims(Γ_fs)
h_Γfs = get_array(∫(1)dΓ_fs)
α = CellField( lazy_map(h->γ*μ_f/(h.^dim),h_Γfs), Γ_fs)</code></pre><ul><li>The linear form associated with the solid counterpart, <span>$l_s(\mathbf{v}^h_{\rm s})$</span>, is defined as</li></ul><div>\[l_s(\mathbf{v})=\int\mathbf{v}\cdot s\ d\Omega_{\rm S}.\]</div><pre><code class="language-julia">l_s(v) = ∫( v⋅s )dΩₛ</code></pre><ul><li>The linear form associated with the fluid counterpart, <span>$l_f(\mathbf{v}^h_{\rm f})$</span> is defined as</li></ul><div>\[l_f(\mathbf{v},q)=\int\mathbf{v}\cdot f + qg\ d\Omega_{\rm F}.\]</div><pre><code class="language-julia">l_f((v,q)) = ∫( v⋅f + q*g )dΩ_f</code></pre><ul><li>The linear form associated with the fluid Neumann boundary condition, <span>$l_{f,\Gamma_N}(\mathbf{v}^h_{\rm f})$</span>, is defined as</li></ul><div>\[l_{f,\Gamma_N}(\mathbf{v})=\int\mathbf{v}\cdot h\ d\Gamma_{\rm N_{out}}.\]</div><pre><code class="language-julia">l_f_Γn(v) = ∫( v⋅hN )dΓ_out</code></pre><p>The final bilinear and linear forms will be given by</p><h3 id="Algebraic-System-of-Equations-1"><a class="docs-heading-anchor" href="#Algebraic-System-of-Equations-1">Algebraic System of Equations</a><a class="docs-heading-anchor-permalink" href="#Algebraic-System-of-Equations-1" title="Permalink"></a></h3><p>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:</p><div>\[\begin{bmatrix}
\mathbf{A}_{u,\rm S}&\mathbf{A}_{u,\rm SF}&0\\
\mathbf{A}_{u,\rm FS}&\mathbf{A}_{u,\rm F}&\mathbf{A}_{up,\rm F}\\
0&\mathbf{A}_{up,\rm F}&\mathbf{A}_{p,\rm F}\\
\end{bmatrix}\begin{bmatrix}
\mathbf{U}^h_{\rm S}\\
\mathbf{U}^h_{\rm F}\\
\mathbf{P}^h_{\rm F}\\
\end{bmatrix} = \begin{bmatrix}
\mathbf{F}^h_{u,\rm S}\\
\mathbf{F}^h_{u,\rm F}\\
\mathbf{F}^h_{p,\rm F}\\
\end{bmatrix}\]</div><p>In order to construct the system we define the final FE operator, constructed using the function <code>AffineFEOperator</code> passing as arguments the bilinear and linear forms, <span>$a$</span> and <span>$l$</span>, together with the trial and test FE spaces, <span>$X$</span> and <span>$Y$</span>.</p><pre><code class="language-julia">op = AffineFEOperator(a,l,X,Y)</code></pre><p>Finally, we call <code>solve</code> to obtain the solution vector of nodal values <span>$[\mathbf{U}^h_{\rm S},\mathbf{U}^h_{\rm F},\mathbf{P}^h_{\rm F}]^T$</span></p><pre><code class="language-julia">uhs, uhf, ph = solve(op)</code></pre><h2 id="Post-processing-1"><a class="docs-heading-anchor" href="#Post-processing-1">Post-processing</a><a class="docs-heading-anchor-permalink" href="#Post-processing-1" title="Permalink"></a></h2><h3 id="Visualization-1"><a class="docs-heading-anchor" href="#Visualization-1">Visualization</a><a class="docs-heading-anchor-permalink" href="#Visualization-1" title="Permalink"></a></h3><p>The solution fields <span>$[\mathbf{U}^h_{\rm S},\mathbf{U}^h_{\rm F},\mathbf{P}^h_{\rm F}]^T$</span> are defined over all the domain, extended with zeros on the inactive part. Calling the function <code>writevtk</code> passing the global triangulation, we will output the global fields.</p><pre><code class="language-julia">writevtk(Ω,"trian", cellfields=["uhs" => uhs, "uhf" => uhf, "ph" => ph])</code></pre><p><img src="../../assets/fsi/Global_solution.png" alt/></p><p>However, we can also restrict the fields to the active part by calling the function <code>restrict</code> with the field along with the respective active triangulation.</p><pre><code class="language-julia">writevtk(Ω_s,"trian_solid",cellfields=["uhs"=>uhs])
writevtk(Ω_f,"trian_fluid",cellfields=["ph"=>ph,"uhf"=>uhf])</code></pre><p><img src="../../assets/fsi/Local_solution.png" alt/></p><h3 id="Quantities-of-Interest-1"><a class="docs-heading-anchor" href="#Quantities-of-Interest-1">Quantities of Interest</a><a class="docs-heading-anchor-permalink" href="#Quantities-of-Interest-1" title="Permalink"></a></h3><pre><code class="language-julia">Γ_S = BoundaryTriangulation(model,tags=["cylinder","interface"])
dΓₛ = Measure(Γ_S,bdegree)
n_ΓS = get_normal_vector(Γ_S)
FD, FL = sum( ∫( (n_ΓS⋅σ_dev_f(ε(uhf))) - ph*n_ΓS )*dΓₛ )
println("Drag force: ", FD)
println("Lift force: ", FL)</code></pre><h2 id="References-1"><a class="docs-heading-anchor" href="#References-1">References</a><a class="docs-heading-anchor-permalink" href="#References-1" title="Permalink"></a></h2><p>[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<em>. In Fluid Structure Interaction II (pp. 193-220). Springer, Berlin, Heidelberg.</em></p><p>[2] Burman, E., and Fernández, M. A. <em>Stabilized explicit coupling for fluid–structure interaction using Nitsche's method.</em> Comptes Rendus Mathematique 345.8 (2007): 467-472.</p><hr/><p><em>This page was generated using <a href="https://github.com/fredrikekre/Literate.jl">Literate.jl</a>.</em></p></article><nav class="docs-footer"><a class="docs-footer-prevpage" href="../t010_isotropic_damage/">« 10 Isotropic damage model</a></nav></div><div class="modal" id="documenter-settings"><div class="modal-background"></div><div class="modal-card"><header class="modal-card-head"><p class="modal-card-title">Settings</p><button class="delete"></button></header><section class="modal-card-body"><p><label class="label">Theme</label><div class="select"><select id="documenter-themepicker"><option value="documenter-light">documenter-light</option><option value="documenter-dark">documenter-dark</option></select></div></p><hr/><p>This document was generated with <a href="https://github.com/JuliaDocs/Documenter.jl">Documenter.jl</a> on <span class="colophon-date" title="Thursday 17 December 2020 09:46">Thursday 17 December 2020</span>. Using Julia version 1.5.3.</p></section><footer class="modal-card-foot"></footer></div></div></div></body></html>