navier-stokes.md 9.14 KB
 Praetorius, Simon committed Nov 13, 2020 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 Navier-Stokes equation ====================== We consider the incompressible Navier-Stokes equation in a domain $\Omega=(0,L_x)\times(0,L_y)$, math \varrho(\partial_t\mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u}) - \nabla\cdot\big(2\nu\mathbf{D}(\mathbf{u})\big) - \nabla p = \mathbf{f},\qquad\text{ in }\Omega \\ \nabla\cdot\mathbf{u} = 0  with velocity $\mathbf{u}$, pressure $p$, density $\varrho$, dynamic viscosity $\nu$, symmetric rate-of-strain tensor $\mathbf{D}(\mathbf{u})=\frac{1}{2}(\nabla\mathbf{u}+\nabla\mathbf{u}^\top)$, and external volume force $\mathbf{f}$, w.r.t. to boundary conditions $\mathbf{u}=\mathbf{g}$ on $\partial\Omega$. For an overview, see also [Wikipedia](https://en.wikipedia.org/wiki/Navier%E2%80%93Stokes_equations). Discretization -------------- ### Weak formulation  Praetorius, Simon committed Nov 13, 2020 18 In a weak variational formulation, we try to find $\mathbf{u}(t,\cdot)\in \mathbf{V}_\mathbf{g}:=\{\mathbf{v}\in H^1(\Omega)^d\,:\, \operatorname{tr}_{\partial\Omega}\mathbf{v} = \mathbf{g}\}$ and $p(t,\cdot)\in Q:=L_0^2(\Omega)$, such that  Praetorius, Simon committed Nov 13, 2020 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 math \int_\Omega \varrho(\partial_t\mathbf{u} + (\mathbf{u}\cdot\nabla)\mathbf{u})\cdot\mathbf{v} + \nu\nabla\mathbf{u}:\nabla\mathbf{v} + p\,\nabla\cdot\mathbf{v}\,\text{d}\mathbf{x} = \int_\Omega \mathbf{f}\cdot\mathbf{v}\,\text{d}\mathbf{x},\qquad\forall \mathbf{v}\in \mathbf{V}_0, \\ \int_\Omega q\,\nabla\cdot\mathbf{u}\,\text{d}\mathbf{x} = 0,\qquad \forall q\in Q.  for $t\in[0,T]$, given $\mathbf{u}(0,\cdot)\equiv\mathbf{u}^0$. The pair $\big(\mathbf{u}(t,\cdot), p(t,\cdot)\big)$ is thus in the product space $\mathbf{V}_\mathbf{g}\times Q$ where $\mathbf{V}_\mathbf{g}$ is also a product of $H^1$ spaces, i.e., $\mathbf{V}_\mathbf{g}=V_{g_0}\times V_{g_1}\times V_{g_2}\simeq [V_{g_i}]^d$. Later, this space will be approximated with a Taylor-Hood space of piecewise quadratic and linear functions. ### Time-discretization A linearized time-discretization with a semi-implicit backward Euler scheme then reads: For $k=0,1,2,\ldots, N$ find $\mathbf{u}^{k+1}\in\mathbf{V}_\mathbf{g},\, p\in Q$, s.t. math \int_\Omega \varrho\Big(\frac{1}{\tau}(\mathbf{u}^{k+1} - \mathbf{u}^k) + (\mathbf{u}^k\cdot\nabla)\mathbf{u}^{k+1}\Big)\cdot\mathbf{v} + \nu\nabla\mathbf{u}^{k+1}:\nabla\mathbf{v} + p\,\nabla\cdot\mathbf{v}\,\text{d}\mathbf{x} \\ = \int_\Omega \mathbf{f}(t^{k+1})\cdot\mathbf{v}\,\text{d}\mathbf{x},\qquad\forall \mathbf{v}\in \mathbf{V}_0, \\ \int_\Omega q\,\nabla\cdot\mathbf{u}^{k+1}\,\text{d}\mathbf{x} = 0,\qquad \forall q\in Q.  with $0=t^0( lagrange<2>(), flatInterleaved()), lagrange<1>(), flatLexicographic());  At this point, any Dirichlet boundary conditions are ignored in the definition of the bases. We  Praetorius, Simon committed Nov 13, 2020 53 build the Taylor-Hood basis as a composition of a product of$d\times\mathbb{V}_2$bases for  Praetorius, Simon committed Nov 13, 2020 54 55 56 57 58 59 60 the velocity components and a$\mathbb{V}_1$basis for the pressure. We have to use composite(...) to combine the different types for velocity and pressure space, while we can use power(...) for the d velocity component spaces of the same type. For efficiency, we use flat indexing strategies flatInterleaved and flatLexicographic to build a global continuous numbering of the basis functions.  Praetorius, Simon committed Nov 13, 2020 61 See the tutorial [Grids and Discrete Functions](../tutorials/grids-and-discretefunctions.md)  Praetorius, Simon committed Nov 13, 2020 62 63 64 65 66 67 68 69 70 71 72 73 for more details. Implementation -------------- We split the implementation into three parts, 1. the time derivative, 2. the stokes operator and 3. any external forces. For addressing the different components of the Taylor-Hood basis, we introduce the treepaths _v = Dune::Indices::_0 and _p = Dune::Indices::_1 for velocity and pressure, respectively. ### Problem framework We have an in-stationary problem consisting of a sequence of stationary equations  Praetorius, Simon committed Nov 13, 2020 74 75 for each timestep, thus we combine a [ProblemInstat](../reference/Problem.md#class-probleminstat) with a [ProblemStat](../reference/Problem.md#class-problemstat) .  Praetorius, Simon committed Nov 13, 2020 76 77 78  c++ using namespace AMDiS;  Praetorius, Simon committed Nov 13, 2020 79 ProblemStat prob{"stokes", grid, basis};  Praetorius, Simon committed Nov 13, 2020 80 81 prob.initialize(INIT_ALL);  Praetorius, Simon committed Nov 13, 2020 82 ProblemInstat probInstat{"stokes", prob};  Praetorius, Simon committed Nov 13, 2020 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 probInstat.initialize(INIT_UH_OLD);  For convenience, we use class template argument deduction here. ### Time derivative For simplicity of this example, we use a backward Euler time stepping and a linearization of the nonlinear advection term: c++ // define a constant fluid density double density = 1.0; Parameters::get("stokes->density", density); // a reference to 1/tau auto invTau = probInstat.invTau(); // <1/tau * u, v> auto opTime = makeOperator(tag::testvec_trialvec{}, density * invTau); prob.addMatrixOperator(opTime, _v, _v); // <1/tau * u^old, v> auto opTimeOld = makeOperator(tag::testvec{}, density * invTau * probInstat.oldSolution(_v)); prob.addVectorOperator(opTimeOld, _v); for (int i = 0; i < Grid::dimensionworld; ++i) { // <(u^old * nabla)u_i, v_i> auto opNonlin = makeOperator(tag::test_gradtrial{}, density * prob.solution(_v)); prob.addMatrixOperator(opNonlin, treepath(_v,i), treepath(_v,i)); }  ### Stokes operator The stokes part of the Navier-Stokes equation is a common pattern that emerges in several fluid equations. math \int_\Omega \nu\nabla\mathbf{u}:\nabla\mathbf{v} + p\,\nabla\cdot\mathbf{v} + q\,\nabla\cdot\mathbf{u}\,\text{d}\mathbf{x}  for$\mathbf{u},\mathbf{v}\in V$and$p,q\in Q$. c++ // define a constant fluid viscosity double viscosity = 1.0; Parameters::get("stokes->viscosity", viscosity); for (int i = 0; i < Grid::dimensionworld; ++i) { // auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity); prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i)); } // auto opP = makeOperator(tag::divtestvec_trial{}, 1.0); prob.addMatrixOperator(opP, _v, _p); // auto opDiv = makeOperator(tag::test_divtrialvec{}, 1.0); prob.addMatrixOperator(opDiv, _p, _v);  where we have used the identity math \nabla\mathbf{u}:\nabla\mathbf{v} = \sum_i \nabla u_i\cdot\nabla v_i  ### External volume force The force may be implemented as a vector-valued function of the global coordinates or any other expression that leads to a local force vector at the quadrature points, e.g. c++ using WorldVector = typename Grid::template Codim<0>::Geometry::GlobalCoordinate; auto opForce = makeOperator(tag::testvec{}, [](WorldVector const& x) { return WorldVector{1.0, 0.0}; }); prob.addVectorOperator(opForce, _v);  Numerical example ----------------- We consider a lid-driven cavity problem in a square domain$\Omega$with$L_x=L_y=1\$ no external force, and with Dirichlet boundary conditions math \mathbf{g}(\mathbf{x}) = \begin{pmatrix}0 \\ v_1 x_1(1 - x_1)(1 - x_0)\end{pmatrix}  Thus, the boundary value is zero on top, right, and bottom boundary and parabolic on the left boundary. The boundary condition is implemented, by first setting boundary ids for parts of the grid boundary and second, setting the Dirichlet value: c++ double vel = 1.0; Parameters::get("stokes->boundary velocity", v1); // define boundary values auto g = [vel](WorldVector const& x) { return WorldVector{0.0, v1*x[1]*(1.0 - x[1])*(1.0 - x[0])}; }; // set boundary conditions for velocity prob.boundaryManager()->setBoxBoundary({1,1,1,1}); prob.addDirichletBC(1, _v, _v, g);  ### Simulation The time-stepping process with a stationary problem in each iteration can be started using an AdaptInstationary manager class: c++ // set initial conditions prob.solution(_v).interpolate(g); // start simulation AdaptInfo adaptInfo("adapt"); AdaptInstationary adapt("adapt", prob, adaptInfo, probInstat, adaptInfo); adapt.adapt();  Complete Example ---------------- The full source code of the Navier-Stokes example can be found in the repository  Praetorius, Simon committed Nov 13, 2020 203 at [examples/navier_stokes.cc](https://gitlab.mn.tu-dresden.de/amdis/amdis-core/blob/master/examples/navier_stokes.cc).  Praetorius, Simon committed Nov 13, 2020 204 205 206 207 208 209 210 211 212 213 214 215 216  Compile with bash cmake --build build --target navier_stokes.2d  and run with bash ./build-cmake/examples/navier_stokes.2d examples/init/navier_stokes.dat.2d  The flow field with density=1, viscosity=1 and boundary velocity=10 will look like  Praetorius, Simon committed Nov 13, 2020 217