Notes on the Finite Element Method 8 The Time-dependent Vorticity Equation: Kelvin-Helmholtz Instability 10 Time Dependent Navier Stokes

Chapter 9 Time-dependent Viscous Flows

9.1 Equations and Boundary Condition.

\frac{\partial\bm{u}}{\partial t}=-\bm{\nabla}p+\nu\nabla^{2}\bm{u}.

(9.1)

where $\nu$ represents the kinematic viscosity. For incompressible flow, the mass conservation equation is:

\bm{\nabla}\cdot\bm{u}=0,

(9.2)

9.1.1 Analytic Solution

Consider the channel illustrated in figure LABEL:Fi:C6S1ShowPoseuille. Initially, for $t\leq 0$ , the velocities are everywhere zero. At $t=0$ the pressure difference across the channel is set equal to one. How does the velocity change? The solution to the steady flow problem has been described in Chapter 6, where the axial velocity $u$ is show to vary parabolically across the width of the channel (equation 6.5)

Batchelor, Chap 4, pp 193-194, shows that equation 9.1 can be made homogeneous by introducing a new dependent variable:

u^{\prime}=u+\frac{1}{2\nu}\frac{\partial p}{\partial x}(1-y^{2})

(9.3)

Look for a solution where $v$ is everywhere zero, so that $u^{\prime}$ is independent of $x$ , the equation for $u^{\prime}$ is:

\frac{\partial u^{\prime}}{\partial t}=\nu\frac{\partial^{2}u^{\prime}}{% \partial y^{2}};\quad u^{\prime}(-1,t)=u^{\prime}(1,t)=0,\quad u^{\prime}(x,y,% 0)=\frac{1}{2\nu}\frac{\partial p}{\partial x}(1-y^{2}).

(9.4)

Solutions of the form

u^{\prime}=A_{n}\cos{\frac{n\,\pi\,y}{2}}\mathrm{e}^{-\lambda_{n}^{2}\nu t},% \quad n=1,3,5\cdots

(9.5)

exist so long as:

\lambda_{n}=\pm\frac{n\,\pi}{2}

(9.6)

The constants $A_{n}$ are determined by the initial condition:

A_{n}=\frac{16}{\nu\,n^{3}\pi^{3}}\frac{\partial p}{\partial x}\sin{\frac{n\,% \pi}{2}}.

(9.7)

The analytic solution is then:

u(x,y,t)=-\frac{1}{2\nu}\frac{\partial p}{\partial x}\left[(1-y^{2})+\frac{16}% {n^{3}\pi^{3}}\sum_{n}\sin{\frac{n\,\pi}{2}}\cos{\frac{n\,\pi\,y}{2}}\mathrm{e% }^{-\lambda_{n}^{2}\nu t}\right],\quad n=1,3,5\cdots

(9.8)

This solution is illustrated in figure 9.1.

Poiseuille flow initiated by a suddenly imposed pressure
gradient, — Figure 9.1: Poiseuille flow initiated by a suddenly imposed pressure gradient, $\partial p/\partial x=-2$ , with kinematic viscosity $\nu=1$ . The analytic solution (black) is indistinguishable from the FEM solution (red). The mesh is $3\times 3$ .

9.2 Finite Element Method

The weak form of equations 9.1 through 9.2 can be written

\mathcal{M}\frac{\partial\,u}{\partial t}+\nu\mathcal{K}\,u=-\mathcal{G}^{x}\,p

(9.9)

\mathcal{M}\frac{\partial\,v}{\partial t}+\nu\mathcal{K}\,v=-\mathcal{G}^{y}\,p

(9.10)

\mathcal{D}^{x}\,u+\mathcal{D}^{y}\,v=0

(9.11)

The system of constrained global equations can be written as:

9.2.1 Open BC in terms of pressure.

Weak Form

For linear TF, the stiffness matrix is given by equation 6.16. In this case the gradient matrices are the same: $\mathcal{D}^{x}=\mathcal{G}^{x}$ and $\mathcal{D}^{x}=\mathcal{G}^{x}$ .

Assembly

Since the pressure is specified at the boundaries, there is no ambiguity regarding the mean presseure nor any need for imposing a constraint. In matrix form equations 6.10 through 6.12 are

\begin{bmatrix}\mathcal{M}&0&0\\ 0&\mathcal{M}&0\\ 0&0&0\end{bmatrix}\frac{\partial\bm{s}}{\partial t}+\begin{bmatrix}\nu\mathcal% {K}&0&\mathcal{D}^{x}\\ 0&\nu\mathcal{K}&\mathcal{D}^{y}\\ \mathcal{D}^{x}&\mathcal{D}^{y}&0\\ 0&0&1&0\end{bmatrix}\bm{s}=\mathcal{P}\frac{\partial\bm{s}}{\partial t}+% \mathcal{L}\bm{s}=0

(9.12)

where $\mathcal{M}$ is the usual mass matrix. The zero in the lower right corner of the left lost matrix remains a concern, notably when the transient problem is forced by velocity.

The trapezoidal is adopted as a scheme to step the equations forward in time.

Consider the grid illustrated in figure LABEL:Fi:C4S2simplegrid. The grid has 9 nodes; there are three variables: $u, v$ and $p$ . the total number of variabless is $27$ of which $18$ are prescribed at the boundaries (six BC on both $u$ and $v$ and six more on $p$ ) and $9$ variables are unknown. The problem for the unknown variables is then:

The solution is compared with the analytical solution on figure 9.1. The agreement is excelent.

9.2.2 Open BC in terms of velocity.

If the preceeding scheme were applied to the situation where pressure is unknown everywhere, but the velocity is given on three of the four boundaries, there would be four unknown velocities (two $u$ and two $v$ ) and nine unknown $p$ . Just as was found in chapter 6, it is possible to show that if the number of unknown pressure values exceeds the total number of unknown velocities, the matrix $\widetilde{L}$ is singular. The accepted solution is to use different order polynomials for velocity $\phi_{n}$ (e.g. quadratic polynomials) and for pressure $\phi_{n}^{p}$ (linear polynomials). In this case the matrix form is:

\begin{bmatrix}\mathcal{M}&0&0&0\\ 0&\mathcal{M}&0&0\\ 0&0&0&0\\ 0&0&0&0\\ \end{bmatrix}\frac{\partial\bm{s}}{\partial t}+\begin{bmatrix}\nu\mathcal{K}&0% &\mathcal{G}^{x}&0\\ 0&\nu\mathcal{K}&\mathcal{G}^{y}&0\\ \mathcal{D}^{x}&\mathcal{D}^{y}&0&1\\ 0&0&1&0\end{bmatrix}\bm{s}=\mathcal{P}\frac{\partial\bm{s}}{\partial t}+% \mathcal{L}\bm{s}=0

(9.13)

where now the $\mathcal{G}$ s and $\mathcal{D}s$ are no longer the same, and a constraint is included on the mean pressure.

Assembly

Consider now the reverse of the problem described in 9.2.1: Poiseuille flow is built up by increasing the velocity at the $x=-1$ boundary according to equation 9.8. The computed pressure constant shoul be constant in time.

Just refer to Chap 6 figure2.

Poiseuille in a simple grid consisting of two elements each
oulined by black lines. Black squares and numbers identify nodes
where the pressure is evaluated, Velovity is evaluated at all
notes (black and red). Velocity is prescribed at all nodes except
6 and 9. — Figure 9.2: Poiseuille in a simple grid consisting of two elements each oulined by black lines. Black squares and numbers identify nodes where the pressure is evaluated, Velovity is evaluated at all notes (black and red). Velocity is prescribed at all nodes except 6 and 9.

Any element has three vertices, where pressure is evaluared and a total of six nodes where the velocity is evaluated: the three vertices supplemented by three mid-point along the element boundaries.

Even with this scheme in place, the pressure is only determined to within a constant which would still cause the matrix to be singular. This problem has been encountered in Chapter 4. It can be dealt with by adding a Lagrange multiplier to the problem, as discussed in Chapter 4, section 4.3.1. The global problem can now be written

\begin{bmatrix}\mathcal{K}&0&\mathcal{G}^{x}&0\\ 0&\mathcal{K}&\mathcal{G}^{y}&0\\ \mathcal{D}^{x}&\mathcal{D}^{y}&0&1\\ 0&0&1&0\end{bmatrix}\begin{bmatrix}u\\ v\\ p\\ \lambda\end{bmatrix}=\mathcal{L}\bm{s}=0,\quad\mbox{where}\quad\bm{s}=\begin{% bmatrix}u\\ v\\ p\\ \lambda\end{bmatrix}

(9.14)

9.2.3 Discretization in Time.

Equation LABEL:sec:C9S2:DBC is discretized through the trapezoidal rule:

\left(\widetilde{\mathcal{P}}+\frac{\delta t}{2}\widetilde{\mathcal{L}}\right)% \widetilde{\bm{s}}^{t+1}=\left(\widetilde{\mathcal{P}}-\frac{\delta t}{2}% \widetilde{\mathcal{L}}\right)\widetilde{\bm{s}}^{t}-\delta t\left(\widehat{% \mathcal{P}}\frac{\partial\widehat{\bm{s}}}{\partial t}+\widehat{\mathcal{L}}% \,\widehat{\bm{s}}\right),

(9.15)

9.2.4 Assembly

\begin{bmatrix}\mathcal{D}^{x}&\mathcal{D}^{y}&0\\ \mathcal{M}\frac{\partial}{\partial t}+\mathcal{K}&0&\mathcal{G}^{x}\\ 0&\mathcal{M}\frac{\partial}{\partial t}+\mathcal{K}&\mathcal{G}^{y}\end{% bmatrix}\begin{bmatrix}u\\ v\\ p\end{bmatrix}=0

(9.16)

9.3 Transient Poiseuille Flow

Define response time.This response time determined by speed of sound.

9.3.1 Stiff Problems

Read: https://en.wikipedia.org/wiki/Stiff_equation

Consider a damped spring-mass system:

m\frac{d^{2}y}{dt^{2}}+r\frac{d\,y}{dt}+k\,y=0\quad y(0)=0,\quad\frac{d\,y}{dt% }(0)=0

(9.17)

where the spring constant $k=101$ and the damping coefficient $r=100$ are lage and nearly equal. Maxima gives the solution as:

\frac{100\,e^{-t}}{99}-\frac{e^{-100\,t}}{99}

(9.18)

To step equation 9.17, rewrite it as a system of two first order ODE’s:

\frac{d}{dt}\bm{s}=\begin{bmatrix}-101&-100\\ 1&0\end{bmatrix}\bm{s}=\mathcal{K}\bm{s}\quad\mbox{where}\quad\bm{s}=\begin{% bmatrix}g\\ y\end{bmatrix}

(9.19)

If the trapezoidal rule is used to discretize in time, The system of equations is

\left(\mathcal{I}-\frac{\delta t}{2}\mathcal{K}\right)\bm{s}^{t+1}=\left(% \mathcal{I}+\frac{\delta t}{2}\mathcal{K}\right)\bm{s}^{t}

(9.20)

The result is illustrated id figure xx (blue line). The time step, $\delta t=0.1$ is selcted on the basis that except for a very small time initially, the solution is expected to decay as $exp(-t)$ . The initial large oscillations of the derivative $g=dy/dt$ are sometimes referred to as ringing. This behavior is symtomatic of stiff¹¹ 1 no relation to the stiffness matris described elsewhere problems. To understand what is happening, consider the one step growth factor called $\bm{G}$ in Strang’s book. Need $\bm{G}<1$ For Trapezoidal, if $du/dt=a\,u$

\bm{G}=\frac{1+a\,\delta t/2}{1-a\,\delta t/2}

(9.21)