/
index.html
63 lines (57 loc) · 14 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
<!DOCTYPE html>
<html lang="en"><head><meta charset="UTF-8"/><meta name="viewport" content="width=device-width, initial-scale=1.0"/><title>8 Incompressible Navier-Stokes · Gridap tutorials</title><link href="https://cdnjs.cloudflare.com/ajax/libs/normalize/4.2.0/normalize.min.css" rel="stylesheet" type="text/css"/><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/4.6.3/css/font-awesome.min.css" rel="stylesheet" type="text/css"/><link href="https://cdnjs.cloudflare.com/ajax/libs/highlight.js/9.12.0/styles/default.min.css" rel="stylesheet" type="text/css"/><script>documenterBaseURL="../.."</script><script src="https://cdnjs.cloudflare.com/ajax/libs/require.js/2.2.0/require.min.js" data-main="../../assets/documenter.js"></script><script src="../../siteinfo.js"></script><script src="../../../versions.js"></script><link href="../../assets/documenter.css" rel="stylesheet" type="text/css"/></head><body><nav class="toc"><h1>Gridap tutorials</h1><select id="version-selector" onChange="window.location.href=this.value" style="visibility: hidden"></select><form class="search" id="search-form" action="../../search/"><input id="search-query" name="q" type="text" placeholder="Search docs"/></form><ul><li><a class="toctext" href="../../">Introduction</a></li><li><a class="toctext" href="../t001_poisson/">1 Poisson equation</a></li><li><a class="toctext" href="../t002_validation/">2 Code validation</a></li><li><a class="toctext" href="../t003_elasticity/">3 Linear elasticity</a></li><li><a class="toctext" href="../t0041_p_laplacian/">4 p-Laplacian</a></li><li><a class="toctext" href="../t004_hyperelasticity/">5 Hyper-elasticity</a></li><li><a class="toctext" href="../t005_dg_discretization/">6 Poisson equation (with DG)</a></li><li><a class="toctext" href="../t007_darcy/">7 Darcy equation (with RT)</a></li><li class="current"><a class="toctext" href>8 Incompressible Navier-Stokes</a><ul class="internal"><li><a class="toctext" href="#Problem-statement-1">Problem statement</a></li><li><a class="toctext" href="#Numerical-Scheme-1">Numerical Scheme</a></li><li><a class="toctext" href="#Discrete-model-1">Discrete model</a></li><li><a class="toctext" href="#FE-spaces-1">FE spaces</a></li><li><a class="toctext" href="#Nonlinear-weak-form-1">Nonlinear weak form</a></li><li><a class="toctext" href="#Nonlinear-solver-phase-1">Nonlinear solver phase</a></li><li><a class="toctext" href="#References-1">References</a></li></ul></li></ul></nav><article id="docs"><header><nav><ul><li><a href>8 Incompressible Navier-Stokes</a></li></ul><a class="edit-page" href="https://github.com/gridap/Tutorials/blob/master/src/t008_inc_navier_stokes.jl"><span class="fa"></span> Edit on GitHub</a></nav><hr/><div id="topbar"><span>8 Incompressible Navier-Stokes</span><a class="fa fa-bars" href="#"></a></div></header><h1><a class="nav-anchor" id="Tutorial-8:-Incompressible-Navier-Stokes-equations-1" href="#Tutorial-8:-Incompressible-Navier-Stokes-equations-1">Tutorial 8: Incompressible Navier-Stokes equations</a></h1><p><a href="https://mybinder.org/v2/gh/gridap/Tutorials/gh-pages?filepath=v0.6.0/notebooks/t008_inc_navier_stokes.ipynb"><img src="https://mybinder.org/badge_logo.svg" alt/></a> <a href="https://nbviewer.jupyter.org/github/gridap/Tutorials/blob/gh-pages/v0.6.0/notebooks/t008_inc_navier_stokes.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 nonlinear multi-field PDEs in Gridap</li><li>How to build FE spaces whose functions have zero mean value</li></ul><h2><a class="nav-anchor" id="Problem-statement-1" href="#Problem-statement-1">Problem statement</a></h2><p>The goal of this last tutorial is to solve a nonlinear multi-field PDE. As a model problem, we consider a well known benchmark in computational fluid dynamics: the lid-driven cavity for the incompressible Navier-Stokes equations. Formally, the PDE we want to solve is: find the velocity vector <span>$u$</span> and the pressure <span>$p$</span> such that</p><div>\[\left\lbrace
\begin{aligned}
-\Delta u + \mathit{Re}\ (u\cdot \nabla)\ u + \nabla p = 0 &\text{ in }\Omega,\\
\nabla\cdot u = 0 &\text{ in } \Omega,\\
u = g &\text{ on } \partial\Omega,
\end{aligned}
\right.\]</div><p>where the computational domain is the unit square <span>$\Omega \doteq (0,1)^d$</span>, <span>$d=2$</span>, <span>$\mathit{Re}$</span> is the Reynolds number (here, we take <span>$\mathit{Re}=10$</span>), and <span>$(w \cdot \nabla)\ u = (\nabla u)^t w$</span> is the well known convection operator. In this example, the driving force is the Dirichlet boundary velocity <span>$g$</span>, which is a non-zero horizontal velocity with a value of <span>$g = (1,0)^t$</span> on the top side of the cavity, namely the boundary <span>$(0,1)\times\{1\}$</span>, and <span>$g=0$</span> elsewhere on <span>$\partial\Omega$</span>. Since we impose Dirichlet boundary conditions on the entire boundary <span>$\partial\Omega$</span>, the mean value of the pressure is constrained to zero in order have a well posed problem,</p><div>\[\int_\Omega q \ {\rm d}\Omega = 0.\]</div><h2><a class="nav-anchor" id="Numerical-Scheme-1" href="#Numerical-Scheme-1">Numerical Scheme</a></h2><p>In order to approximate this problem we chose a formulation based on inf-sub stable <span>$Q_k/P_{k-1}$</span> elements with continuous velocities and discontinuous pressures (see, e.g., [1] for specific details). The interpolation spaces are defined as follows. The velocity interpolation space is</p><div>\[V \doteq \{ v \in [C^0(\Omega)]^d:\ v|_T\in [Q_k(T)]^d \text{ for all } T\in\mathcal{T} \},\]</div><p>where <span>$T$</span> denotes an arbitrary cell of the FE mesh <span>$\mathcal{T}$</span>, and <span>$Q_k(T)$</span> is the local polynomial space in cell <span>$T$</span> defined as the multi-variate polynomials in <span>$T$</span> of order less or equal to <span>$k$</span> in each spatial coordinate. Note that, this is the usual continuous vector-valued Lagrangian FE space of order <span>$k$</span> defined on a mesh of quadrilaterals or hexahedra. On the other hand, the space for the pressure~is</p><div>\[\begin{aligned}
Q_0 &\doteq \{ q \in Q: \ \int_\Omega q \ {\rm d}\Omega = 0\}, \text{ with}\\
Q &\doteq \{ q \in L^2(\Omega):\ q|_T\in P_{k-1}(T) \text{ for all } T\in\mathcal{T}\},
\end{aligned}\]</div><p>where <span>$P_{k-1}(T)$</span> is the polynomial space of multi-variate polynomials in <span>$T$</span> of degree less or equal to <span>$k-1$</span>. Note that functions in <span>$Q_0$</span> are strongly constrained to have zero mean value. This is achieved in the code by removing one degree of freedom from the (unconstrained) interpolation space <span>$Q$</span> and adding a constant to the computed pressure so that the resulting function has zero mean value.</p><p>The weak form associated to these interpolation spaces reads: find <span>$(u,p)\in U_g \times Q_0$</span> such that <span>$[r(u,p)](v,q)=0$</span> for all <span>$(v,q)\in V_0 \times Q_0$</span> where <span>$U_g$</span> and <span>$V_0$</span> are the set of functions in <span>$V$</span> fulfilling the Dirichlet boundary condition <span>$g$</span> and <span>$0$</span> on <span>$\partial\Omega$</span> respectively. The weak residual <span>$r$</span> evaluated at a given pair <span>$(u,p)$</span> is the linear form defined as</p><div>\[[r(u,p)](v,q) \doteq a((v,q),(u,p))+ [c(u)](v),\]</div><p>with</p><div>\[\begin{aligned}
a((v,q),(u,p)) &\doteq \int_{\Omega} \nabla v \cdot \nabla u \ {\rm d}\Omega - \int_{\Omega} (\nabla\cdot v) \ p \ {\rm d}\Omega + \int_{\Omega} q \ (\nabla \cdot u) \ {\rm d}\Omega,\\
[c(u)](v) &\doteq \int_{\Omega} v \cdot \left( (u\cdot\nabla)\ u \right)\ {\rm d}\Omega.\\
\end{aligned}\]</div><p>Note that the bilinear form <span>$a$</span> is associated with the linear part of the PDE, whereas <span>$c$</span> is the contribution to the residual resulting from the convective term.</p><p>In order to solve this nonlinear weak equation with a Newton-Raphson method, one needs to compute the Jacobian associated with the residual <span>$r$</span>. In this case, the Jacobian <span>$j$</span> evaluated at a pair <span>$(u,p)$</span> is the bilinear form defined as</p><div>\[[j(u,p)]((v,q),(\delta u, \delta p)) \doteq a((v,q),(\delta u,\delta p)) + [{\rm d}c(u)](v,\delta u),\]</div><p>where <span>${\rm d}c$</span> results from the linearization of the convective term, namely</p><div>\[[{\rm d}c(u)](v,\delta u) \doteq \int_{\Omega} v \cdot \left( (u\cdot\nabla)\ \delta u \right) \ {\rm d}\Omega + \int_{\Omega} v \cdot \left( (\delta u\cdot\nabla)\ u \right) \ {\rm d}\Omega.\]</div><p>The implementation of this numerical scheme is done in Gridap by combining the concepts previously seen for single-field nonlinear PDEs and linear multi-field problems.</p><h2><a class="nav-anchor" id="Discrete-model-1" href="#Discrete-model-1">Discrete model</a></h2><p>We start with the discretization of the computational domain. We consider a <span>$100\times100$</span> Cartesian mesh of the unit square.</p><pre><code class="language-julia">using Gridap
n = 100
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(n,n))</code></pre><p>For convenience, we create two new boundary tags, namely <code>"diri1"</code> and <code>"diri0"</code>, one for the top side of the square (where the velocity is non-zero), and another for the rest of the boundary (where the velocity is zero).</p><pre><code class="language-julia">labels = FaceLabels(model)
add_tag_from_tags!(labels,"diri1",[6,])
add_tag_from_tags!(labels,"diri0",[1,2,3,4,5,7,8])</code></pre><h2><a class="nav-anchor" id="FE-spaces-1" href="#FE-spaces-1">FE spaces</a></h2><p>For the velocities, we need to create a conventional vector-valued continuous Lagrangian FE space. In this example, we select a second order interpolation.</p><pre><code class="language-julia">D = 2
order = 2
fespace1 = FESpace(
reffe=:Lagrangian, conformity=:H1, valuetype=VectorValue{D,Float64},
model=model, labels=labels, order=order, diritags=["diri0","diri1"])</code></pre><p>The interpolation space for the pressure is build as follows</p><pre><code class="language-julia">fespace2 = FESpace(
reffe=:PLagrangian, conformity=:L2, valuetype=Float64,
model=model, order=order-1, constraint=:zeromean)</code></pre><p>With the options <code>reffe=:PLagrangian</code>, <code>valuetype=Float64</code>, and <code>order=order-1</code>, we select the local polynomial space <span>$P_{k-1}(T)$</span> on the cells <span>$T\in\mathcal{T}$</span>. With the symbol <code>:PLagrangian</code> we specifically chose a local Lagrangian interpolation of type "P". Using <code>:Lagrangian</code>, would lead to a local Lagrangian of type "Q" since this is the default for quadrilateral or hexahedral elements. On the other hand, <code>constraint=:zeromean</code> leads to a FE space, whose functions are constrained to have mean value equal to zero, which is just what we need for the pressure space. With these objects, we build the test and trial multi-field FE spaces</p><pre><code class="language-julia">uD0(x) = VectorValue(0.0,0.0)
uD1(x) = VectorValue(1.0,0.0)
V = TestFESpace(fespace1)
Q = TestFESpace(fespace2)
Y = [V, Q]
U = TrialFESpace(fespace1,[uD0,uD1])
P = TrialFESpace(fespace2)
X = [U, P]</code></pre><h2><a class="nav-anchor" id="Nonlinear-weak-form-1" href="#Nonlinear-weak-form-1">Nonlinear weak form</a></h2><p>The different terms of the nonlinear weak form for this example are defined following an approach similar to the one discussed for the <span>$p$</span>-Laplacian equation, but this time using the notation for multi-field problems.</p><pre><code class="language-julia">const Re = 10.0
@law conv(x,u,∇u) = Re*(∇u')*u
@law dconv(x,du,∇du,u,∇u) = conv(x,u,∇du)+conv(x,du,∇u)
function a(y,x)
u, p = x
v, q = y
inner(∇(v),∇(u)) - inner(div(v),p) + inner(q,div(u))
end
c(v,u) = inner(v,conv(u,∇(u)))
dc(v,du,u) = inner(v,dconv(du,∇(du),u,∇(u)))
function res(x,y)
u, p = x
v, q = y
a(y,x) + c(v,u)
end
function jac(x,y,dx)
u, p = x
v, q = y
du, dp = dx
a(y,dx)+ dc(v,du,u)
end</code></pre><p>With the functions <code>res</code>, and <code>jac</code> representing the weak residual and the Jacobian, we build the nonlinear FE problem:</p><pre><code class="language-julia">trian = Triangulation(model)
quad = CellQuadrature(trian,degree=(order-1)*2)
t_Ω = NonLinearFETerm(res,jac,trian,quad)
op = NonLinearFEOperator(Y,X,t_Ω)</code></pre><h2><a class="nav-anchor" id="Nonlinear-solver-phase-1" href="#Nonlinear-solver-phase-1">Nonlinear solver phase</a></h2><p>To finally solve the problem, we consider the same nonlinear solver as previously considered for the <span>$p$</span>-Laplacian equation.</p><pre><code class="language-julia">using LineSearches: BackTracking
nls = NLSolver(
show_trace=true, method=:newton, linesearch=BackTracking())
solver = NonLinearFESolver(nls)</code></pre><p>In this example, we solve the problem without providing an initial guess (a default one equal to zero will be generated internally)</p><pre><code class="language-julia">uh, ph = solve(solver,op)</code></pre><p>Finally, we write the results for visualization (see next figure).</p><pre><code class="language-julia">writevtk(trian,"ins-results",cellfields=["uh"=>uh,"ph"=>ph])</code></pre><p><img src="../../assets/t008_inc_navier_stokes/ins_solution.png" alt/></p><h2><a class="nav-anchor" id="References-1" href="#References-1">References</a></h2><p>[1] H. C. Elman, D. J. Silvester, and A. J. Wathen. <em>Finite elements and fast iterative solvers: with applications in incompressible fluid dynamics</em>. Oxford University Press, 2005.</p><p><em>This page was generated using <a href="https://github.com/fredrikekre/Literate.jl">Literate.jl</a>.</em></p><footer><hr/><a class="previous" href="../t007_darcy/"><span class="direction">Previous</span><span class="title">7 Darcy equation (with RT)</span></a></footer></article></body></html>