Notes on the Finite Element Method 5 Two Dimensional Waves on a Rotating Sheet of Water 7 Non-linear Advection

Chapter 6 Steady Viscous flow

6.1 Equations and Boundary Conditions

Viscous flow is governed by conservation of mass and momentum. For steady flow, the momentum equation states that shear and pressure forces are in balance.

\bm{\nabla}p=\nu\nabla^{2}\bm{u}.

(6.1)

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

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

(6.2)

Boundary conditions specify the two components $u, v$ of the velocity vector on solid boundaries. For the remaining BC there are two options: either the pressure difference alon the channel can be prescribed, or the velocities at the entrance ( $x=-1$ ) is given. These are Dirichlet conditions in either case, however the assemble method is very different dependig on which conditios are fixed along the open voundaries.

The streamfunction and vorticity are related variables. The vorticity is the curl of the velocity field. In two dimensions

\omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}

(6.3)

The stream function $\psi$ satisfies:

\nabla^{2}\psi=-\omega

(6.4)

The two dimensional Poiseuille solution describes steady viscous flow in a channel with parallel wall, when the entrance velocity varies parabolically with maximum value equal to one. Taket he wall separation to be $H^{*}$ and the channel length to be $L^{*}$ , the solution is:

u(y)=-\frac{1}{2\nu}\frac{\partial p}{\partial x}(1-y^{2}),\quad v=0,\quad p=-% 2\nu\,x+p_{0}.

(6.5)

where $x=2x^{*}/L^{*}$ and $y=2y^{*}/H^{*}$ , and $p_{0}$ is an arbitrary constant. In this solution the velocity is independent of $x$ , the abcissa parallel to the walls and the pressure in independent of the ordinate, $y$ . The vorticity and streamfunction ares also a functions of $y$ alone:

\omega=2\,y,\quad\psi=y-\frac{y^{3}}{3}+\psi_{0}.

(6.6)

The vorticity gradient is constant across the channel and the exists a vorticity gardient at each wall. $\psi_{0}$ is another arbitrary constant. The solution is illustrated in figure

Poiseuille flow in a constant width channel. The velocity
profile is parabolic, the vorticity profile is linear and the
streamfunction varies as a cubic polynomial. The pressure is
independent of

Figure 6.1: Poiseuille flow in a constant width channel. The velocity profile is parabolic, the vorticity profile is linear and the streamfunction varies as a cubic polynomial. The pressure is independent of

y

and decreases linearly from the left opening.

This analytic solution is illustrated in figure 6.1. Profile of the velocity $u$ . the vorticity $\omega$ and the streamfunction $\psi$ , variables that are independent of $x$ are illustrated in the top frame. The pressure (independent of $y$ ) is illustrated beneath.

6.2 Finite Element Method

The weak form of equations 6.2 and 7.1 is obtained in the usual way:

\begin{split}\displaystyle\nu\int\!\phi_{m}\left(\frac{\partial^{2}u}{\partial x% ^{2}}+\frac{\partial^{2}u}{\partial y^{2}}\right)\,da+\int\!\phi_{m}\frac{% \partial p}{\partial x}\,da&\displaystyle=0\\ \displaystyle\nu\int\!\phi_{m}\left(\frac{\partial^{2}v}{\partial x^{2}}+\frac% {\partial^{2}v}{\partial x^{2}}\right)\,da+\int\!\phi_{m}\frac{\partial p}{% \partial y}\,da&\displaystyle=0\\ \displaystyle\int\!\phi_{m}\left(\frac{\partial u}{\partial x}+\frac{\partial v% }{\partial y}\right)\,da&\displaystyle=0\end{split}

(6.7)

The pressure and the velocities are expandent in termd of trial functions. The following forms are chosen to allow the possibilty of different trial functions for velocity and pressure:

(u,v)(x,y)=\sum_{n=1}^{nu}\phi_{n}(x,y)(u_{n},v_{n})\quad\mbox{where}\quad\phi% _{n}(x,y)=\sum_{k=1}\!C_{k,n}x^{p(k)}y^{q(k)}

(6.8)

p=\sum_{n=1}^{np}\phi^{p}_{n}(x,y)p_{n}\quad\mbox{where}\quad\phi^{p}_{n}(x,y)% =\sum_{k=1}\!C_{k,n}^{p}x^{p(k)}y^{q(k)}

(6.9)

Introducing these expansions into equation 6.7 leads to the usual system of equations for the three dependent variables:

\int\!\phi_{m}\frac{\partial\phi_{n}^{p}}{\partial x}\,da\,p_{n}+\int\!\left(% \frac{\partial\phi_{m}}{\partial x}\frac{\partial\phi_{n}}{\partial x}+\frac{% \partial\phi_{m}}{\partial y}\frac{\partial\phi_{n}}{\partial y}\right)\,da\,u% _{n}=\int_{s}\;\phi_{m}\nabla\phi_{n}\cdot\bm{n}\,u_{n}\,ds

(6.10)

\int\!\phi_{m}\frac{\partial\phi_{n}^{p}}{\partial y}\,da\,p_{n}+\int\!\left(% \frac{\partial\phi_{m}}{\partial x}\frac{\partial\phi_{n}}{\partial x}+\frac{% \partial\phi_{m}}{\partial y}\frac{\partial\phi_{n}}{\partial y}\right)\,da\,v% _{n}=\int_{s}\;\phi_{m}\nabla\phi_{n}\cdot\bm{n}\,v_{n}\,ds

(6.11)

\int\!\phi^{p}_{m}\frac{\partial\phi_{n}}{\partial x}\,da\,u_{n}+\int\!\phi^{p% }_{m}\frac{\partial\phi_{n}}{\partial y}\,da\,v_{n}

(6.12)

Define global matrices $\mathcal{D}^{x},\mathcal{D}^{y}$ , $\mathcal{G}^{x},\mathcal{G}^{y}$ and $\mathcal{K}$ with elements:

D_{m,n}^{x}=\int\!\phi^{p}_{m}\frac{\partial\phi_{n}}{\partial x}\,da\quad D_{% m,n}^{y}=\int\!\phi^{p}_{m}\frac{\partial\phi_{n}}{\partial y}\,da

(6.13)

G_{m,n}^{x}=\int\!\phi_{m}\frac{\partial\phi_{n}^{p}}{\partial x}\,da\quad G_{% m,n}^{y}=\int\!\phi_{m}\frac{\partial\phi_{n}^{p}}{\partial y}\,da

(6.14)

\mathcal{K}_{m,n}=\int\!\left(\frac{\partial\phi_{m}}{\partial x}\frac{% \partial\phi_{n}}{\partial x}+\frac{\partial\phi_{m}}{\partial y}\frac{% \partial\phi_{n}}{\partial y}\right)\,da

(6.15)

6.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}$ ; For the canonical triangle:

\mathcal{D}e^{x}=\frac{A}{6}\begin{bmatrix}-1&1&0\\ -1&1&0\\ -1&1&0\end{bmatrix}\quad\mathcal{D}e^{y}=\frac{A}{6}\begin{bmatrix}-1&0&1\\ -1&0&1\\ -1&0&1\end{bmatrix}

(6.16)

Assembly

In matrix form equations 6.10 through 6.12 are

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

(6.17)

The zero in the lower right corner of the left lost matrix is a concern, as we will see it can cause the matrix to be singular. Fortunately this is not the case when pressure boundary conditions are imposed.

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:

\widetilde{\mathcal{L}}\widetilde{\bm{s}}=-\widehat{\mathcal{L}}\widehat{\bm{s}}

(6.18)

where $\widetilde{\mathcal{L}}$ is a $9\times 9$ matrix and $\widehat{\mathcal{L}}$ is $9\times 18$ . The solution vector has been partitioned into $\widehat{\bm{s}}$ , the known values at Dirichlet boundary points, and $\widetilde{\bm{s}}$ , the unknown values.

\widetilde{\mathcal{L}}=\frac{1}{3}\begin{bmatrix}12&-6&&&&&1&2&\\ -6&24&-6&&&&-1&0&1\\ &-6&12&&&&&-2&-1\\ &&&12&-6&&-1&1&\\ &&&-6&24&-6&-2&0&2\\ &&&&-6&12&&-1&1\\ -1&1&&1&2&&&&\\ -2&0&2&-1&0&1&&&\\ &-1&1&&-2&-1&&&\end{bmatrix},\;\mbox{and}\;\widehat{\mathcal{L}}\widehat{\bm{s% }}=\frac{1}{2}\begin{bmatrix}1\\ 2\\ 1\\ 0\\ 0\\ 0\\ 0\\ 0\\ 0\end{bmatrix}\;\mbox{with the result:}\;\begin{bmatrix}u_{6}&=&0.5\\ u_{7}&=&0.5\\ u_{9}&=&0.5\\ v_{6}&=&0\\ v_{7}&=&0\\ v_{9}&=&0\\ p_{5}&=&0\\ p_{6}&=&-0\\ p_{8}&=&0\end{bmatrix}.

(6.19)

6.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$ . It is possible to show that if the number of unknow 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). This combination of trial function is called mixed TF.

Weak Form

To achhieve this the matrix elements $C_{k,n}$ in equation 6.8is $6\times 6$ , whereas matrix elements $C_{k,n}^{p}$ in equation LABEL:eq:C6S2:exppressi remains $3\times 3$ . Now matrices $\mathcal{D}^{x}$ and $\mathcal{G}^{x}$ are no longer equal. For the canonical triangle:

\mathcal{D}e^{x}=\frac{A}{6}\begin{bmatrix}-1&1&0&1&0&-1\\ 0&-1&1&1&0&-1\\ 0&0&0&2&0&-1\end{bmatrix};\quad\mathcal{G}e^{x}=\frac{A}{6}\begin{bmatrix}0&0&% 0\\ -1&1&0\\ 0&0&0\\ -1&1&0\\ 0&0&0\\ -1&1&0\end{bmatrix}.

(6.20)

Assembly

Consider for example Poiseuille flow in a two dimensional channel, defined by four nodes (rather than nine in the example above). Velocity trial functions are taken to be quadratic polynomial and pressure trial functions are take to be linear polynomials. The pressure is determined at the four corners of the grid, whereas velocity is determined on each of the nine nodes drawn in 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. — Figure 6.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}

(6.21)

6.2.3 Dirichlet Boundary Conditions

As noted in section 1.3.3, the Dirichlet boundary points are dealt with by rewriting equation 9.14

\widetilde{\mathcal{L}}=\frac{1}{3}\begin{bmatrix}16&-4&&&-1&1&-1&1&0\\ -4&8&&&&&-1&1&0\\ &&16&-4&-1&-1&1&1&0\\ &&-4&8&&-1&&1&0\\ 1&&1&&&&&&3\\ -1&2&1&1&&&&&3\\ 1&1&-1&&&&&&3\\ -1&1&-1&-1&&&&&3\\ &&&&3&3&3&3&0\end{bmatrix},\;\mbox{and}\;\widehat{\mathcal{L}}\widehat{\bm{s}}% =\frac{1}{3}\begin{bmatrix}4\\ 0\\ 0\\ 0\\ 1\\ 1\\ 2\\ 0\\ 0\end{bmatrix}\;\mbox{with the result:}\;\begin{bmatrix}u_{6}&=&1\\ u_{9}&=&1\\ v_{6}&=&0\\ v_{9}&=&0\\ p_{1}&=&2\\ p_{2}&=&-2\\ p_{3}&=&2\\ p_{4}&=&-2\\ \lambda&=&0\end{bmatrix}.

(6.22)

in exact agreement with equation 6.5. Note in passing thatthe determinant of a modified $\widetilde{\mathcal{L}}$ modified to remove the Lagrangian contraint is zero.

6.2.4 Evaluating the solution at some arbitrary point

In many cases it may be required to evalute the solution at some location other than a grid node.

6.3 Results

6.3.1 Convergent Channel.

Figure 6.3: Viscous flow in a converging channel

6.3.2 Cavity Flow

Viscous flow in a Cavity. Along the top the velocity is
forced to unity. Left: streamfunction; the red areas in the lower
corners are regions of separated flow. Middle: vorticity, the
black lines represent the — Figure 6.4: Viscous flow in a Cavity. Along the top the velocity is forced to unity. Left: streamfunction; the red areas in the lower corners are regions of separated flow. Middle: vorticity, the black lines represent the $\omega=0$ contour; near the surface the vorticity in negative, and changes sign with depth. The vorticity is negative in the lower corners. Right: pressure, The max value is $p=591$ in the upper left corner.

6.3.3 Flow Behind a Step.

Viscous flow around a step. The velocity profile at the
letf (W) boundary is parabolic. Note the small recirculation are
near the corner at

Figure 6.5: Viscous flow around a step. The velocity profile at the letf (W) boundary is parabolic. Note the small recirculation are near the corner at

x=-1

y=-1

6.3.4 Flow Around a Cylinder

The cylinder is embedded in a flow that extends 6 units in the $x$ -direction and four units in the $y$ -direction. The cylinder is centered on $0,0$ with a radius $0.2$ . The boundary condition set the $x$ -velocity to one along the N, W, and S boundaries, with lateral velocity set to zero. The FEM solution is illustrated in figure 7.1, where only the central portion of the computed flowfield is shown. The streamfunction shows the flow streaming by the cylinder in a pattern that is summetric about the $y$ axis. The vorticity exhibits a maximum at the bottomost part of the cylinder, where the pressure gradient driven vorticity sourge islargets. A minimum in vorticity occurs at the diametrically opposite point. The pressure is maximum at the upstream stagnation point and minimum on the dowstrem stagnation point. There is a lrge pressure gradient around the body.

Viscous flow around a cylinder. The u-velocity is set to
one along the N, W and S boundaries. The solution is found in a
domain that extends from — Figure 6.6: Viscous flow around a cylinder. The u-velocity is set to one along the N, W and S boundaries. The solution is found in a domain that extends from $-3\leq x\leq 3$ and $-2\leq y\leq 2$ . Only the central portion is illustrated.