## Overview Given vⁿ, pⁿ: 1. Discretize momentum → compute v* 2. Solve pressure Poisson equation: $\nabla^{2}p=\left( \frac{\rho}{\Delta t} \right) \nabla\cdot v^{\star}$ 3. Correct velocity: $v^{n+1}=v^{\star}-\left( \frac{\Delta t}{\rho} \right)\nabla p$ 4. Advance time step or iterate (steady case) ## Discretization $\frac{ \partial \rho \vec{v} }{ \partial t } + \nabla.(\rho \vec{v}\vec{v}) = -\nabla p + \mu \nabla^{2}\vec{v}+\rho \vec{g}$ $\nabla.\vec{v}=0$ ### Space Using the **finite volume method (FVM)**: - Divide the domain into control volumes. - Integrate the equations over each CV. - Apply the divergence theorem to turn volume integrals into surface integrals. Each term in the momentum equation is then discretized: - **Convective term**: needs interpolation (UDS, CDS, QUICK, etc.) - **Diffusion term**: usually central differences - **Pressure gradient**: centered or interpolated to faces - **Body forces**: added as source terms ### Time ## Pressure Correction $\frac{ \partial }{ \partial x_{i} }\left( \frac{ \partial p }{ \partial x_{i} } \right)= -\frac{ \partial }{ \partial x_{i} }\left( \frac{ \partial (\rho u_{i}u_{j}) }{ \partial x_{j} } \right)$ 1. Solve momentum for a known pressure resulting in a new velocity 2. Solve Poisson to obtain a corrected pressure 3. Correct velocity (and possibly pressure) 4. Iterate for next time step Coupled nonlinear algebraic equations: $(\rho u_{i})^{n+1}-(\rho u_{i})^{n}=\Delta t\left(- \frac{ \partial (\rho u_{i}u_{j})^{n+1} }{ \partial x_{j} }+\frac{ \partial \tau_{ij}^{n+1} }{ \partial x_{j} }-\frac{ \partial p^{n+1} }{ \partial x_{i} } \right)$ $\frac{ \partial }{ \partial x_{i} }\left( \frac{ \partial p^{n+1} }{ \partial x_{i} } \right)= -\frac{ \partial }{ \partial x_{i} }\left( -\frac{ \partial (\rho u_{i}u_{j})^{n+1} }{ \partial x_{j} } + \frac{ \partial \tau_{ij}^{n+1} }{ \partial x_{j} }\right)$ ### Poisson Equation Find pressure that enforces incompressibility using intermediate velocity correct velocity: $v^{n+1}=v^{\star}-\frac{\Delta t}{\rho}\nabla p$ plug into continuity: $\nabla \cdot \left( v^{\star}-\frac{\Delta t}{\rho}\nabla p \right)=0$ Poisson equation: $\boxed{\nabla^{2}p=\frac{\rho}{\Delta t}\nabla\cdot v^{\star}}$ ## Time Stepping ### Implicit - Use pressure from time $t^{n}$ to compute velocity at $t^{n+1}$. - Pressure correction is computed at the old time level. ### Explicit - Coupled solution: velocity and pressure are solved simultaneously. - More stable for stiff problems, especially at larger time steps. #### 2. **Implicit Time-Stepping (e.g., Backward Euler)** - Coupled solution: velocity and pressure are solved simultaneously. - More stable for stiff problems, especially at larger time steps. ## Solver types ### Nonlinear Solver - Guess $u_{i}^{n}$ and employ a nonlinear solver at each time step for $u_{i}^{n+1}$ - Fixed-point, newton-raphson, or secant ### Linearization about $t_{n}$ - nonlinear term is only quadratic, so we can neglect the 2nd order ${} \Delta u_{i}\Delta u_{j}$ term ### ADI Scheme - Split momentum equations to block tri-diagonal 1D problems ## Steady problems For steady problems, you just want to end up at the right solution, so an implicit scheme with a large timestep works well $A^{u_{i}^{n+1}}u_{i}^{n+1} = b_{u_{i}}^{n+1}-\frac{ \partial p^{n+1} }{ \partial x_{i} }$ ### Implicit Methods #### Outer Iterations Nonlinear solve iterating on $m\star$ and updates $A^{u_{i}^{m\star}}$ and $u_{i}^{m\star}$ for new velocity given the old pressure $A^{u_{i}^{m\star}}u_{i}^{m\star} = b_{u_{i}^{m\star}}^{m-1}-\frac{ \partial p^{m-1} }{ \partial x_{i} }$ enforce continuity and get a new pressure #### Inner Iterations Linear solve for final corrected $u_{i}^{n+1}$ $A^{u_{i}^{m\star}}u_{i}^{m} = b_{u_{i}^{m\star}}^{m-1}-\frac{ \partial p^{m} }{ \partial x_{i} }$ ### Projection Methods Construct a velocity field that doesn't satisfy continuity, then correct it using a pressure gradient 1. solve for intermediate velocity field $u^{\star}$ without considering pressure 2. Find $p^{n+1}$ using Poisson equation 3. Correct velocity using pressure gradient to get $u^{n+1}$ - Substitute $u_{i}^{m}=u^{m\star}_{i}+u'$ and $p^{m}=p^{m-1}+p'$ - SIMPLE: neglects contributions of $u'$ - SIMPLEC: approximates $u'$ in the pressure equation as a function of $p'$ - SIMPLER and PISO methods: iterate to obtain $u'$ ## Unsteady Problems ### Projection Methods $P u_{i}^{\star^{n+1}} = (I-\mathrm{Divergence})u_{i}^{\star^{n+1}}$ ### Helmholtz Theorem For sufficiently smooth vector fields $\vec{f}$ $\vec{f} = \vec{\nabla} \phi + \vec{\nabla} \times \vec{a}$ $\vec{\nabla}\phi$: curl free component ($\vec{\nabla}\times \vec{\nabla}\phi=0$ $\vec{\nabla}.(\vec{\nabla} \times \vec{a})=0$ ### Non-Incremental Projection Method Ignore pressure, then correct 1. Time Mom Equation $\rho u_{i}^{m+1}-\rho u_{i}^{m} = \Delta t ( A_\text{advect}^{m/m+1}+D_\text{iff}^{m/m+1}) - \Delta t \frac{ \partial P^{m+1} }{ \partial x_{i} } + BC$ 2. Prediction Equation $\rho u_{i}^{\star^{m+1}}-\rho u_{i}^{m} = \Delta t ( A_\text{advect}^{\star^{m/m+1}}+D_\text{iff}^{\star^{m/m+1}}) - 0 + BC$ Subtract 1 - 2 $\rho u_{i}^{n+1}-\rho u_{i}^{\star^{m+1}} = \Delta t(A^{m/m+1}-A^{\star^{m/m+1}} + D^{m/m+1} - D^{\star^{m/m+1}})- \Delta t \frac{ \partial p }{ \partial x_{i} }$ If A and $D$ are explicit at time $n$, $\Delta t(A-A+D-D)=0$ $\rho u_{i}^{n+1}-\rho u_{i}^{\star^{m+1}} = - \Delta t \frac{ \partial p }{ \partial x_{i} }$ eqn 3 ^ Correction equation for velocity Impose $\vec{\nabla}.\rho \vec{u}=0$ $\frac{ \partial }{ \partial x_{i} }(\rho u_{i}^{\star})^{n+1} = \Delta t \frac{ \partial }{ \partial x_{i} }\left( \frac{ \partial p^{n+1} }{ \partial x_{i} } \right)$ ^ PPE BC: use NS - at the boundary - along the normal direction Take eqn 3, $\cdot \vec{n}$ at boundary $\underbrace{ \left.(\rho \vec{u}^{n+1}.\vec{n})\right|_{\text{bnd}} - \left.(\rho \vec{u}^{\star}.\vec{n})\right|_{\text{bnd}} }_{ \text{should = 0} } = -\left.(\Delta t \vec{\nabla}P^{n+1} \cdot \vec{n})\right|_{\text{bnd}}$ $0 = -\left.(\Delta t \vec{\nabla}P^{n+1} \cdot \vec{n})\right|_{\text{bnd}}$ ### Incremental Methods Use pressure from previous time step to predict velocity, then solve for the correction pressure $\delta q$ and update pressure with $p^{n+1}=p^{n}+\delta q$ #### Rotational Incremental Desired scheme: (i) $\rho u^{n+1}_{i}=\rho u^{n}_{i}+\Delta t(A^{n}+D^{n+1})-\Delta t \frac{ \partial p^{n+1} }{ \partial x_{i} }$ use diffusion implicitly to help with stability condition Velocity predictor scheme: (ii) $\rho u^{\star^{n+1}}_{i}=\rho u^{n}_{i}+\Delta t(A^{n}+D^{\star^{n+1}})-\Delta t \frac{ \partial p^{n} }{ \partial x_{i} }$ ${} D^{\star^{n+1}}=\frac{ \partial }{ \partial x_{j} }\left( \mu \frac{ \partial u_{i}^{\star^{n+1}} }{ \partial x_{j} } \right) {}$ Velocity corrector scheme: $\rho u^{n+1}_{i}=\rho u_{i}^{\star^{n+1}}-\Delta t \frac{ \partial \delta q^{n+1} }{ \partial x_{i} }$ Pressure Poisson Equation: $\underbrace{ \frac{ \partial }{ \partial x_{i} }(\rho u^{n+1}_{i}) }_{ \text{want this to = 0} }=\underbrace{ \frac{ \partial }{ \partial x_{i} }(\rho u_{i}^{\star^{n+1}}) }_{ \neq 0 }-\Delta t\frac{ \partial }{ \partial x_{i} }\left( \frac{ \partial \delta q^{n+1} }{ \partial x_{i} } \right)$ $u_{i}^{\star^{n+1}}$ satisfies momentum equation and BC's Derive boundary condition for $\delta q^{n+1}$ Project NS momentum onto $\vec{n}$ at boundary $\left.\rho u^{n+1}\cdot \vec{n}\right|_{\delta D}=\left.\rho u^{\star^{n+1}}\cdot \vec{n}\right|_{\delta D} - \Delta t \frac{ \partial \delta q^{n+1} }{ \partial \vec{x} }\cdot \vec{n}$ $\frac{ \partial \delta q^{n+1} }{ \partial n } = 0$ Get $p^{n+1}$ by subtracting (ii) from (i) $\rho u_{i}^{n+1}-\rho u_{i}^{\star^{n+1}}= \Delta t(0 + D^{n+1}-D^{\star^{n+1}})-\Delta t \frac{ \partial }{ \partial x_{i} }(p^{n+1}-p^{n})$ $-\frac{ \partial \delta q^{n+1} }{ \partial x_{i} }= \frac{ \partial }{ \partial x_{j} } \left( \mu \frac{ \partial u_{i}^{n+1} }{ \partial x_{j} } + \mu \frac{ \partial u_{i}^{\star^{n+1}} }{ \partial x_{j} } \right) - \frac{ \partial }{ \partial x_{i} }(p^{n+1}-p^{n})$ Take momentum divergence $\frac{ \partial }{ \partial x_{i} }\left( \frac{ \partial }{ \partial x_{i} }(p^{n+1}-p^{n}-\delta q^{n+1}) \right) = \frac{ \partial }{ \partial x_{i} }\left( \frac{ \partial }{ \partial x_{j} }\left( -\mu \frac{ \partial u_{i} }{ \partial x } -\mu \frac{ \partial u_{i}^{\star^{n+1}} }{ \partial x_{j} }\right) \right)$ if $\mu=$ constant: $p^{n+1}=p^{n}+q^{n+1}-\mu \frac{ \partial u_{i}^{\star^{n+1}} }{ \partial x_{i} }$ pressure gradient at the normal boundary can change in rotational incremental | | BC normal | BC tangent | Continuity | | --------------------------- | --------- | ---------- | ---------- | | $u_{i}^{n+1}$ | yes | no | yes | | $u_{i}^{\star^{n+1}}$ | yes | yes | no | ### Fractional Step Methods Write navier stokes as $u_{i}^{n+1}=u_{i}+(C_{i}+D_{i}+P_{i})\Delta t$ solve it in three steps: $u_{i}^{\star}=u_{i}^{n}+C_{i}\Delta t$ $u_{i}^{\star\star}=u_{i}^{\star}+D_{i}\Delta t$ $u_{i}^{n+1}=u_{i}^{\star\star}+P_{i}\Delta t$ Example based on crank-nicolson: first step