Continuity Equations in Rotating Parallel Plate Viscometer

Stationary Parallel Plate

Spatial-fractional derivatives for fluid flow and transport phenomena

Mohamed F. El-Amin , in Fractional-Order Modeling of Dynamic Systems with Applications in Optimization, Signal Processing and Control, 2022

3.5.1 Poiseuille flow

Poiseuille flow is fluid flow between two stationary parallel plates at the points $y=0$ and $y=h$ . (See Fig. 3.2.) Initially, the fluid is at rest, and it starts its motion suddenly under a pressure gradient. This problem can be modeled using a 1D simplified NS equation, and here three fractional models will be considered [37]:

Figure 3.2. Sketch of Poiseuille flow.

(3.27) $\frac{\partial u}{\partial t}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\nu \frac{{\partial }^{\alpha }u}{\partial y^{\alpha }}.$

Along with the initial and boundary conditions,

(3.28) $\begin{array}{c}u(y,0)=0,0<y<h,\hfill \\ u(0,t)=u(h,t)=0,t⩾0.\hfill \end{array}$

If one assumes a uniform effective pressure gradient in the x-direction, i.e., $-\partial p/\partial x=kH(t)$ , k is a constant, and $H(\cdot )$ is the Heaviside step function, the above system can be nondimensionalized to the form

(3.29) $\begin{array}{c}\frac{\partial }{\partial t}u(y,t)=1+\frac{1}{Re}\frac{{\partial }^{\alpha }}{\partial y^{\alpha }}u(y,t),\hfill \\ u(y,0)=0,0<y<1,\hfill \\ u(0,t)=u(1,t)=0,t⩾0,\hfill \end{array}$

where $Re=\sqrt{k/h\rho }\cdot h^{\alpha }/\nu$ is the generalized Reynolds number.

The fractional derivative ${\partial }^{\alpha }/\partial y^{\alpha }$ is called Riesz operator with $1<\alpha ⩽2$ , which can be defined on the finite intervals $y\in [0,h]$ , $t\in [0,T]$ as [1,38]

(3.30) $\frac{{\partial }^{\alpha }}{\partial y^{\alpha }}=-c_{\alpha }({}_{0}D_y^{\alpha }+{}_y{D}_h^{\alpha }),$

where $c_{\alpha }=\frac{1}{2\cos (\pi \alpha /2)}$ . The left-hand side and right-hand side Riemann–Liouville fractional derivatives can be, respectively, written as

(3.31) $\begin{array}{c}\hfill {}_{0}D_y^{\alpha }u(y,t)=\frac{1}{\mathrm{\Gamma }(2-\alpha )}\frac{{\partial }^2}{\partial y^2}{\int }_0^y\frac{u(\xi ,t)\mathrm{d}\xi }{{(y-\xi )}^{\alpha -1}},\hfill \\ \hfill {}_{y}D_h^{\alpha }u(y,t)=\frac{1}{\mathrm{\Gamma }(2-\alpha )}\frac{{\partial }^2}{\partial y^2}{\int }_y^h\frac{u(\xi ,t)\mathrm{d}\xi }{{(\xi -y)}^{\alpha -1}}.\hfill \end{array}$

The discrete approximation version can be written as [1,38]

(3.32) ${}_{0}D_y^{\alpha }u(y_j,t_k)\approx \frac{1}{h^{\alpha }}\sum _{l=0}^{j+1}{g}_l{u}_{j-l+1}^k\text{ and}{}_y{D}_1^{\alpha }u(y_j,t_k)\approx \frac{1}{h^{\alpha }}\sum _{l=0}^{M-j+1}{g}_l{u}_{j+l-1}^k,$

where $g_0=1$ , $g_l={(-1)}^l\alpha (\alpha -1)\cdots (\alpha -l+1)/l!$ , $l⩾1$ , with the uniform discretizations $(y_j,t_k)$ : $y_j=jh$ , $t_k=k\tau$ , $j=0,1,\dots ,M$ , $k=0,1,\dots ,N$ , such that $M,N\in {\mathbb{N}}^{+}$ , $h=1/M$ , and $\tau =T/N$ .

Therefore, the implicit difference scheme can be written as

(3.33) $\frac{u_j^k-u_j^{k-1}}{\tau }=1-\frac{1}{Re}\cdot \frac{1}{2h^{\alpha }\cos \frac{\pi \alpha }{2}}\cdot (\sum _{l=0}^{j+1}{g}_l{u}_{j-l+1}^k+\sum _{l=0}^{M-j+1}{g}_l{u}_{j+l-1}^k).$

Along with the numerical initial and boundary conditions,

(3.34) $u_j^0=0,u_0^k=u_M^k=0,j=0,1,2,\dots ,M,k=0,1,2,\dots ,N,$

which must be comparable with the exact solution of the traditional case, $\alpha =2$ ,

(3.35) $u(y,t)=\frac{y(1-y)}{2}Re-2\mathrm{Re}\sum _{n=1}^{\infty }\frac{1-{(-1)}^n}{{(n\pi )}^3}\exp (-\frac{n^2{\pi }^2}{Re}t)\sin n\pi y.$

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780323900898000088

Basic Flow Theory

W. Brian Rowe DSc, FIMechE , in Hydrostatic, Aerostatic and Hybrid Bearing Design, 2012

Pressure-Induced (Poiseuille) Flow

In Figure 2.2 , laminar flow takes place between two stationary parallel plates separated by a thin film. Flow is caused by the application of pressure at one end. The flow is resisted by shear stresses in the fluid. The velocity at the boundary is zero and increases to a maximum at the center of the thin film. The pressure across the flow is assumed to be constant—that is, ∂p/∂y  =   0. The force equilibrium on an element of width z is

Figure 2.2. Pressure-Induced (Poiseuille) Flow.

$pz\delta y-\left(p+\frac{\text{d}p}{\text{d}x}\delta x\right)z\delta y-\tau z\delta x+\left(\tau +\frac{\text{d}\tau }{\text{d}y}\delta y\right)z\delta x=0$

leading to

(2.5) $\frac{\text{d}\tau }{\text{d}y}=\frac{\text{d}p}{\text{d}x}$

Substituting τ   = η(du/dy) leads to

(2.6) $\eta \frac{{\text{d}}^{2}u}{\text{d}{y}^2}=\frac{\text{d}p}{\text{d}x}$

Integrating twice yields the velocity distribution. The boundary conditions are du/dy  =   0 and y  = h/2 and u  =   0 at y  =   0 or h:

(2.7) $u=\frac{1}{2\eta }\frac{\text{d}p}{\text{d}x}(y^2-yh)$

The integral of the velocity is the volumetric flow:

$q=\frac{z{h}^3}{12\eta }\frac{\text{d}p}{\text{d}x}\text{Volumetric}\text{flow}$

(2.8) $m=\frac{\rho z{h}^3}{12\eta }\frac{\text{d}p}{\text{d}x}\text{Mass flow}$

Equation (2.8) for mass flow applies for both liquids and gases. For gas lubrication, further integration to obtain a pressure distribution requires substitution for the variable density term whereas constant density is appropriate for incompressible liquids. For liquid flow through a thin parallel film slot of width z and length l, the pressure gradient dp/dx = (p ₁  − p ₂)/l is constant, so that

(2.9a) $q=\frac{z{h}^3}{12\eta l}(p_1-p_2)\text{Liquid flow}$

For isothermal gas flow, equation (2.4b) leads to volume flow q ₂ at absolute pressure p ₂:

(2.9b) $q_2=\frac{z{h}^3}{12\eta l}\frac{(p_1^2-p_2^2)}{2p_2}\text{Gas}\text{flow}$

In most situations, volume flow is stated in terms of free air or gas at ambient pressure. For isothermal conditions, which can be assumed in most cases, flow of free air is increased relative to flow at high pressure. The increased flow is proportional to p ₁/p ₂.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780123969941000024

Practical Guidelines for CFD Simulation and Analysis

Jiyuan Tu , ... Chaoqun Liu , in Computational Fluid Dynamics (Third Edition), 2018

7.3.2.2 Other Useful Guidelines

It is noted that for wall-attached boundary layers such as that found within a simple flow between two stationary parallel plates or the backward-facing step geometry, turbulent fluctuations are suppressed adjacent to the wall and the viscous effects become prominent in this region known as the viscous sublayer. This modified turbulent structure generally precludes the application of the two-equation models such as standard k–ɛ model, RNG k–ɛ model, and realizable k–ɛ model or even the Reynolds stress model at the near-wall region, which thereby requires special near-wall modelling procedures. Selecting an appropriate near-wall model represents another important strategy in the context of turbulence modelling. Here, the reader has to decide whether he/she adopts the so-called wall-function method, in which the near-wall region is bridged with wall functions, or a low-Reynolds-number turbulence model, in which the flow structure in the viscous sublayer is totally resolved. This decision will certainly depend on the availability of computational resources and the accuracy requirements for resolution of the boundary layer. Some useful guidance on the application of relevant low Reynolds number models that can be employed all the way through the wall is provided herein. More discussions and practical guidelines will also be given for near-wall treatments using the wall-function method in the next section below.

The k–ω model developed by Wilcox (1998), where ω is a frequency of the large eddies, has been shown to perform splendidly close to walls in boundary-layer flows. Such a model is common in the majority of commercial codes, and it works exceptionally well particularly under strong adverse pressure gradients hence its popularity in aerospace applications. Like the standard k–ɛ model, a modelled transport equation is solved for ω to determine its local distribution within the fluid flow. Nonetheless, the model is very sensitive to the free-stream value of ω, and unless great care is taken in prescribing this value, spurious results are obtained in both boundary-layer flows and free-shear flows. In general, the standard k–ɛ model is less sensitive to the free-stream values but is often inadequate under adverse pressure gradients. To overcome such problems, Menter (1994a,b) proposed to combine both the standard k–ɛ model and k–ω model, which retains the properties of k–ω close to the wall and gradually blends into the standard k–ɛ model away from the wall. This Menter's model has shown to eliminate the free-stream sensitivity problem without sacrificing the k–ω near-wall performance.

To account for strong nonequilibrium effects, the shear-stress transport (SST) variation of Menter's model (1993, 1996) leads to a significant improvement in handling nonequilibrium boundary-layer regions such as those found behind shocks and close to separation. It is therefore highly recommended for flow separation since the real flow is more likely to be much closer to separation (or more separated) than the calculations from the standard k–ɛ model suggest. Bear in mind that SST should not be viewed as a universal cure for turbulence modelling because it still inherits noticeable weaknesses. SST, for example, is less able to cope with flow recovery following flow reattachment. For this, a promising possibility is the use of a length-scale limiting device as proposed by Ince and Launder (1995). Interested readers may also wish to refer to Patel et al. (1985) for the applications of other various low Reynolds number versions of the standard k–ɛ model and the Reynolds stress model where modifications to the governing transport equations are used to deal with near-wall effects allowing these models to be deployed directly through to the wall. Alternatively, the standard k–ɛ model and the Reynolds stress model can be employed in the interior of the flow and coupled to the one-equation k-L model (Wolfshtein, 1969) that is dedicated to resolve mainly the wall region (see the review by Rodi, 1991), a so-called two-layer model. It is imperative that whatever low Reynolds number models are adopted, sufficient number of grid nodal points must be placed into a very narrow region adjacent to the wall to adequately capture the rapid variation in the flow variables.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780081011270000076

Governing Equations for CFD: Fundamentals

Jiyuan Tu , ... Chaoqun Liu , in Computational Fluid Dynamics (Third Edition), 2018

3.6 Generic Form of the Governing Equations for CFD

From the governing equations derived above for either the laminar or turbulent conditions, there are significant commonalities between these various equations. Here, we will present the three-dimensional form of the governing equations for the conservation of mass, momentum, energy, and the turbulent quantities. If we introduce a general variable ϕ and expressing all the fluid flow equations in the conservative incompressible form, the equation can be written as

(3.60) $\frac{\partial \phi }{\partial t}+\frac{\partial \left(u\phi \right)}{\partial x}+\frac{\partial \left(v\phi \right)}{\partial y}+\frac{\partial \left(w\phi \right)}{\partial z}=\frac{\partial }{\partial x}\left[\Gamma \frac{\partial \phi }{\partial x}\right]+\frac{\partial }{\partial y}\left[\Gamma \frac{\partial \phi }{\partial y}\right]+\frac{\partial }{\partial z}\left[\Gamma \frac{\partial \phi }{\partial z}\right]+S_{\phi }$

whilst in the conservative compressible form, the equation is given by

(3.61) $\frac{\partial \left(\rho \varphi \right)}{\partial t}+\frac{\partial \left(\rho u\varphi \right)}{\partial x}+\frac{\partial \left(\rho v\varphi \right)}{\partial y}+\frac{\partial \left(\rho w\varphi \right)}{\partial z}=\frac{\partial }{\partial x}\left[\Gamma \frac{\partial \varphi }{\partial x}\right]+\frac{\partial }{\partial y}\left[\Gamma \frac{\partial \varphi }{\partial y}\right]+\frac{\partial }{\partial z}\left[\Gamma \frac{\partial \varphi }{\partial z}\right]+S_{\varphi }$

Eqs (3.60) and (3.61) are the so-called the transport equations for the property ϕ. Each of them illustrates the various physical transport processes occurring in the fluid flow: the local acceleration and advection terms on the left-hand side are, respectively, equivalent to the diffusion term (Γ  =   diffusion coefficient) and the source term (S_ϕ ) on the right-hand side. Tables 3.1 and 3.2 present the governing equations for the incompressible and compressible flows in the Cartesian framework. In order to bring forth the common features, we have, of course, combined the terms that are not shared between the equations inside the source terms. It is noted that the additional source terms in the momentum equations S _u ′, S _v ′, and S _w ′ comprise of the pressure and nonpressure gradient terms and other possible sources such as gravity that influence the fluid motion, whilst the additional source term S_T in the energy equation may contain heat sources or sinks within the flow domain. It is also noted that the production P in the incompressible form of the turbulence equations can be obtained from its compressible counterpart via invoking the incompressible form of the continuity equation and division by constant density.

Table 3.1. Governing Equations for Incompressible Flow in Cartesian Coordinates

Mass conservation
(m)	$\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}=0$
Momentum equations
(Mx )	$\frac{\partial u}{\partial t}+\frac{\partial \left(uu\right)}{\partial x}+\frac{\partial \left(vu\right)}{\partial y}+\frac{\partial \left(wu\right)}{\partial z}=\frac{\partial }{\partial x}\left[\left(\nu +{\nu }_T\right)\frac{\partial u}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\nu +{\nu }_T\right)\frac{\partial u}{\partial y}\right]+\frac{\partial }{\partial z}\left[\left(\nu +{\nu }_T\right)\frac{\partial u}{\partial z}\right]+\left(S_u=-\frac{1}{\rho }\frac{\partial p}{\partial x}+S_u^{\prime }\right)$
(My )	$\frac{\partial v}{\partial t}+\frac{\partial \left(uv\right)}{\partial x}+\frac{\partial \left(vv\right)}{\partial y}+\frac{\partial \left(wv\right)}{\partial z}=\frac{\partial }{\partial x}\left[\left(\nu +{\nu }_T\right)\frac{\partial v}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\nu +{\nu }_T\right)\frac{\partial v}{\partial y}\right]+\frac{\partial }{\partial z}\left[\left(\nu +{\nu }_T\right)\frac{\partial v}{\partial z}\right]+\left(S_v=-\frac{1}{\rho }\frac{\partial p}{\partial y}+S_v^{\prime }\right)$
(Mz )	$\frac{\partial w}{\partial t}+\frac{\partial \left(uw\right)}{\partial x}+\frac{\partial \left(vw\right)}{\partial y}+\frac{\partial \left(ww\right)}{\partial z}=\frac{\partial }{\partial x}\left[\left(\nu +{\nu }_T\right)\frac{\partial w}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\nu +{\nu }_T\right)\frac{\partial w}{\partial y}\right]+\frac{\partial }{\partial z}\left[\left(\nu +{\nu }_T\right)\frac{\partial w}{\partial z}\right]+\left(S_w=-\frac{1}{\rho }\frac{\partial p}{\partial z}+S_w^{\prime }\right)$
Energy equation
(E)	$\frac{\partial T}{\partial t}+\frac{\partial \left(uT\right)}{\partial x}+\frac{\partial \left(vT\right)}{\partial y}+\frac{\partial \left(wT\right)}{\partial z}=\frac{\partial }{\partial x}\left[\left(\frac{\nu }{\mathit{Pr}}+\frac{{\nu }_T}{P{r}_T}\right)\frac{\partial T}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\frac{\nu }{\mathit{Pr}}+\frac{{\nu }_T}{P{r}_T}\right)\frac{\partial T}{\partial y}\right]+\frac{\partial }{\partial z}\left[\left(\frac{\nu }{\mathit{Pr}}+\frac{{\nu }_T}{P{r}_T}\right)\frac{\partial T}{\partial z}\right]+S_T$
Turbulence equations
(k)	$\frac{\partial k}{\partial t}+\frac{\partial \left(uk\right)}{\partial x}+\frac{\partial \left(vk\right)}{\partial y}+\frac{\partial \left(wk\right)}{\partial z}=\frac{\partial }{\partial x}\left[\frac{{\nu }_T}{{\sigma }_k}\frac{\partial k}{\partial x}\right]+\frac{\partial }{\partial y}\left[\frac{{\nu }_T}{{\sigma }_k}\frac{\partial k}{\partial y}\right]+\frac{\partial }{\partial z}\left[\frac{{\nu }_T}{{\sigma }_k}\frac{\partial k}{\partial z}\right]+\left(S_k=P-D\right)$
(ɛ)	$\frac{\partial \varepsilon }{\partial t}+\frac{\partial \left(u\varepsilon \right)}{\partial x}+\frac{\partial \left(v\varepsilon \right)}{\partial y}+\frac{\partial \left(w\varepsilon \right)}{\partial z}=\frac{\partial }{\partial x}\left[\frac{{\nu }_T}{{\sigma }_{\varepsilon }}\frac{\partial \varepsilon }{\partial x}\right]+\frac{\partial }{\partial y}\left[\frac{{\nu }_T}{{\sigma }_{\varepsilon }}\frac{\partial \varepsilon }{\partial y}\right]+\frac{\partial }{\partial z}\left[\frac{{\nu }_T}{{\sigma }_{\varepsilon }}\frac{\partial \varepsilon }{\partial z}\right]+\left(S_{\varepsilon }=\frac{\varepsilon }{k}\left(C_{\varepsilon 1}P-C_{\varepsilon 2}D\right)\right)$
where	$P=2{\nu }_T\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial z}\right)}^2\right]+{\nu }_T\left[{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)}^2\right]$ and $D=\varepsilon$

Table 3.2. Governing Equations for Compressible Flow in Cartesian Coordinates

Mass conservation
(m)	$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0$
Momentum equations
(Mx )	$\frac{\partial \left(\rho u\right)}{\partial t}+\frac{\partial \left(\rho uu\right)}{\partial x}+\frac{\partial \left(\rho vu\right)}{\partial y}+\frac{\partial \left(\rho wu\right)}{\partial z}=\frac{\partial }{\partial x}\left[\left(\mu +{\mu }_T\right)\frac{\partial u}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\mu +{\mu }_T\right)\frac{\partial u}{\partial y}\right]+\frac{\partial }{\partial z}\left[\left(\mu +{\mu }_T\right)\frac{\partial u}{\partial z}\right]+\left(S_u=-\frac{\partial p}{\partial x}+S_u^{\prime }\right)$
(My )	$\frac{\partial \left(\rho v\right)}{\partial t}+\frac{\partial \left(\rho uv\right)}{\partial x}+\frac{\partial \left(\rho vv\right)}{\partial y}+\frac{\partial \left(\rho wv\right)}{\partial z}=\frac{\partial }{\partial x}\left[\left(\mu +{\mu }_T\right)\frac{\partial v}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\mu +{\mu }_T\right)\frac{\partial v}{\partial y}\right]+\frac{\partial }{\partial z}\left[\left(\mu +{\mu }_T\right)\frac{\partial v}{\partial z}\right]+\left(S_v=-\frac{\partial p}{\partial y}+S_v^{\prime }\right)$
(Mz )	$\frac{\partial \left(\rho w\right)}{\partial t}+\frac{\partial \left(\rho uw\right)}{\partial x}+\frac{\partial \left(\rho vw\right)}{\partial y}+\frac{\partial \left(\rho ww\right)}{\partial z}=\frac{\partial }{\partial x}\left[\left(\mu +{\mu }_T\right)\frac{\partial w}{\partial x}\right]+\frac{\partial }{\partial y}\left[\left(\mu +{\mu }_T\right)\frac{\partial w}{\partial y}\right]+\frac{\partial }{\partial z}\left[\left(\mu +{\mu }_T\right)\frac{\partial w}{\partial z}\right]+\left(S_w=-\frac{\partial p}{\partial z}+S_w^{\prime }\right)$
Energy equation
(E)	$\frac{\partial \left(\rho h\right)}{\partial t}+\frac{\partial \left(\rho uh\right)}{\partial x}+\frac{\partial \left(\rho vh\right)}{\partial y}+\frac{\partial \left(\rho wh\right)}{\partial z}=\frac{\partial }{\partial x}\left[\lambda \frac{\partial T}{\partial x}\right]+\frac{\partial }{\partial y}\left[\lambda \frac{\partial T}{\partial y}\right]+\frac{\partial }{\partial z}\left[\lambda \frac{\partial T}{\partial z}\right]+\frac{\partial }{\partial x}\left[\frac{{\mu }_T}{{Pr}_T}\frac{\partial h}{\partial x}\right]+\frac{\partial }{\partial y}\left[\frac{{\mu }_T}{{Pr}_T}\frac{\partial h}{\partial y}\right]+\frac{\partial }{\partial z}\left[\frac{{\mu }_T}{{Pr}_T}\frac{\partial h}{\partial z}\right]+\frac{\partial p}{\partial t}+\Phi +S_T$
Turbulence equations
(k)	$\frac{\partial \left(\rho k\right)}{\partial t}+\frac{\partial \left(\rho uk\right)}{\partial x}+\frac{\partial \left(\rho vk\right)}{\partial y}+\frac{\partial \left(\rho wk\right)}{\partial z}=\frac{\partial }{\partial x}\left[\frac{{\mu }_T}{{\sigma }_k}\frac{\partial k}{\partial x}\right]+\frac{\partial }{\partial y}\left[\frac{{\mu }_T}{{\sigma }_k}\frac{\partial k}{\partial y}\right]+\frac{\partial }{\partial z}\left[\frac{{\mu }_T}{{\sigma }_k}\frac{\partial k}{\partial z}\right]+\left(S_k=\rho \left(P-D\right)\right)$
(ɛ)	$\frac{\partial \left(\rho \varepsilon \right)}{\partial t}+\frac{\partial \left(\rho u\varepsilon \right)}{\partial x}+\frac{\partial \left(\rho v\varepsilon \right)}{\partial y}+\frac{\partial \left(\rho w\varepsilon \right)}{\partial z}=\frac{\partial }{\partial x}\left[\frac{{\mu }_T}{{\sigma }_{\varepsilon }}\frac{\partial \varepsilon }{\partial x}\right]+\frac{\partial }{\partial y}\left[\frac{{\mu }_T}{{\sigma }_{\varepsilon }}\frac{\partial \varepsilon }{\partial y}\right]+\frac{\partial }{\partial z}\left[\frac{{\mu }_T}{{\sigma }_{\varepsilon }}\frac{\partial \varepsilon }{\partial z}\right]+\left(S_{\varepsilon }=\rho \frac{\varepsilon }{k}\left(C_{\varepsilon 1}P-C_{\varepsilon 2}D\right)\right)$
where	$P=2{\mu }_T\left[{\left(\frac{\partial u}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial z}\right)}^2\right]+{\mu }_T\left[{\left(\frac{\partial u}{\partial y}+\frac{\partial v}{\partial x}\right)}^2+{\left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)}^2+{\left(\frac{\partial w}{\partial x}+\frac{\partial u}{\partial z}\right)}^2\right]-\frac{2}{3}{\mu }_T{\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)}^2-\frac{2}{3}\rho {\mu }_{T}k\left(\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}+\frac{\partial w}{\partial z}\right)$ and $D=\varepsilon$

For compressible flow, the density and temperature can be evaluated through the equations of state, which provide the linkage between the energy equation and those of the mass and momentum equations. For a perfect gas, the following equations of state are $p=\rho RT$ where R is the gas constant and $e=C_{v}T$ where C _v is the specific heat of constant volume. For the dynamics viscosity and thermal conductivity, the variables can usually be determined via a linear or polynomial dependence on temperature.

This equation is usually used as the starting point for computational procedures in either the finite difference or finite volume methods. Algebraic expressions of this equation for the various transport properties are formulated and hereafter solved. For incompressible flow, by setting the transport property ϕ equal to 1, u, v, w, T, k, ɛ and selecting appropriate values for the diffusion coefficient Γ and source terms S_ϕ , we obtain the special forms presented in Table 3.3 for each of the partial differential equations for the conservation of mass, momentum, energy, and the turbulent quantities. In Table 3.4, by setting the transport property ϕ equal to 1, u, v, w, h, k, ɛ and selecting appropriate values for Γ and S_ϕ , we nonetheless obtain the special forms presented in Table 3.4 for the compressible form of the partial differential equations for the conservation of mass, momentum, energy, and the turbulent quantities.

Table 3.3. General Form of Governing Equations for Incompressible Flow in Cartesian Coordinates

Φ	ΓΦ	SΦ
1	0	0
u	$\nu +{\nu }_T$	$-\frac{1}{\rho }\frac{\partial p}{\partial x}+S_u^{\prime }$
v	$\nu +{\nu }_T$	$-\frac{1}{\rho }\frac{\partial p}{\partial y}+S_v^{\prime }$
w	$\nu +{\nu }_T$	$-\frac{1}{\rho }\frac{\partial p}{\partial z}+S_w^{\prime }$
T	$\frac{\nu }{\mathit{Pr}}+\frac{{\nu }_T}{P{r}_T}$	ST
k	$\frac{{\nu }_T}{{\sigma }_k}$	$P-D$
ɛ	$\frac{{\nu }_T}{{\sigma }_{\varepsilon }}$	$\frac{\varepsilon }{k}\left(C_{\varepsilon 1}P-C_{\varepsilon 2}D\right)$

Table 3.4. General Form of Governing Equations for Compressible Flow in Cartesian Coordinates

Φ	ΓΦ	SΦ
1	0	0
u	$\mu +{\mu }_T$	$-\frac{\partial p}{\partial x}+S_u^{\prime }$
v	$\mu +{\mu }_T$	$-\frac{\partial p}{\partial y}+S_v^{\prime }$
w	$\mu +{\mu }_T$	$-\frac{\partial p}{\partial z}+S_w^{\prime }$
h	$\frac{{\mu }_T}{P{r}_T}$	$\frac{\partial }{\partial x}\left[\lambda \frac{\partial T}{\partial x}\right]+\frac{\partial }{\partial y}\left[\lambda \frac{\partial T}{\partial y}\right]+\frac{\partial }{\partial z}\left[\lambda \frac{\partial T}{\partial z}\right]+\frac{\partial p}{\partial t}+\Phi +S_T$
k	$\frac{{\mu }_T}{{\sigma }_k}$	$P-D$
ɛ	$\frac{{\mu }_T}{{\sigma }_{\varepsilon }}$	$\frac{\varepsilon }{k}\left(C_{\varepsilon 1}P-C_{\varepsilon 2}D\right)$

Although we have systematically walked through the derivation of the complete set of governing equations in detail from basic conservation principles, the final general form pertaining to the fluid motion, heat transfer, etc. conforms simply to the generic form of Eq. (3.60) for incompressible flow and Eq. (3.61) for compressible flow. These equations are important generic transport equations as they can accommodate increasing complexity within the CFD model for solving more complicated problems generally found in engineering applications.

Let us focus on some typical complex engineering flow problems that are of significant interest such as multistep combustion processes of swirling turbulent reactive flows in combustors and multiphase flows involving interactions between gas bubbles and liquids in bubble columns. Solutions to these processes can easily be obtained by modelling them through additional transport equations expressed in the simple generic form of Eq. (3.60) for incompressible flow and Eq. (3.61) for compressible flow. For reactive flows, the transport of the various chemical species can be handled by the additional scalar quantities representing each of the reactive species and appropriately formulating the reaction rates in the source terms to account for the chemical reaction processes that are occurring. For bubbly flows, additional transport equations of the number density can be formulated and solved for the various gas bubble sizes that migrate alongside with the liquid in the bubble columns.

Understanding CFD is not meant to be an arduous process. On the contrary, Eqs (3.60) and (3.61), originally formulated from first principles, reinforce the adherent simplicity that is embraced for any transport property that may be required to be solved within the CFD framework.

In Sections 3.3.3 and 3.4.3, we introduced the dimensionless parameters such as the Reynolds number (Re) and Prandtl number (Pr) that may be useful to describe some similar physical phenomena of the flow and heat transfer processes. It will be demonstrated in the next example that the governing equations of mass, momentum, and energy can be nondimensionalized to reduce the number of parameters that appear in the equations. Other similar characteristics of the fluid flow are also discussed.

Example 3.9

Consider the dynamic similarity of the partial differential equations that govern a two-dimensional CFD case for a steady incompressible laminar flow between two stationary parallel plates. The channel dimensions are height H  =   0.1   m and length L  =   1   m. This is the same model geometry as previously investigated in Examples 3.7 and 3.8.

(a)

Nondimensionalize the continuity, momentum, and energy equations given by Eqs (3.12), (3.24), (3.25), and (3.35), respectively.

(b)

Using CFD, determine the flow field for both air and water with the same Reynolds number whilst adjusting the inlet velocity. Discuss the velocity profiles in the fully developed region.

(c)

Also determine the temperature field for both air and water with the same Reynolds number as in part (b). Discuss the nondimensional temperature profiles in the fully developed region with a prescribed uniform temperature of 330   K at the inlet and a wall temperature specified at 300   K.

Solution: (a) The nondimensional form of the governing equations (and boundary conditions) can be achieved by dividing all the dependent and independent flow variables by relevant and meaningful constant quantities. For lengths, the variable can be divided by a characteristic length H (which is the width of the channel), all velocities by a reference velocity u_in (which is the inlet velocity), pressure by ρu _in ² (which is twice the dynamic pressure for the channel), and temperature by a suitable temperature difference (which is T _∞  − T_s for the channel). We therefore obtain

$x⁎=\frac{x}{H},y⁎=\frac{y}{H},u⁎=\frac{u}{u_{\mathit{in}}},v⁎=\frac{v}{u_{\mathit{in}}},p⁎=\frac{p}{\rho u_{\mathit{in}}^2}\mathrm{and}T⁎=\frac{T-T_s}{T_{\infty }-T_s}$

where the asterisks denote the nondimensional variables. Introducing these variables into the governing equations of mass, momentum, and energy produces

$\frac{\partial u⁎}{\partial x⁎}+\frac{\partial v⁎}{\partial y⁎}=0-\mathrm{continuity}$

$\begin{array}{ll}\hfill & \hfill \\ u⁎\frac{\partial u⁎}{\partial x⁎}+v⁎\frac{\partial u⁎}{\partial y⁎}=\frac{1}{Re}\left(\frac{{\partial }^{2}u⁎}{\partial x{⁎}^2}+\frac{{\partial }^{2}u⁎}{\partial y{⁎}^2}\right)-\frac{\partial p⁎}{\partial x⁎}-x-\mathrm{momentum}\hfill & \hfill \\ \hfill & \hfill \\ \hfill & \hfill \end{array}$

$u⁎\frac{\partial v⁎}{\partial x⁎}+v⁎\frac{\partial v⁎}{\partial y⁎}=\frac{1}{Re}\left(\frac{{\partial }^{2}v⁎}{\partial x{⁎}^2}+\frac{{\partial }^{2}v⁎}{\partial y{⁎}^2}\right)-\frac{\partial p⁎}{\partial y⁎}-y-\mathrm{momentum}$

$u⁎\frac{\partial T⁎}{\partial x⁎}+v⁎\frac{\partial T⁎}{\partial y⁎}=\frac{1}{\mathit{RePr}}\left(\frac{{\partial }^{2}T⁎}{\partial x{⁎}^2}+\frac{{\partial }^{2}T⁎}{\partial y{⁎}^2}\right)-\mathrm{energy}$

Discussion: A major advantage of nondimensionalizing the governing equations is the significant reduction of parameters to be considered. By grouping the dimensional parameters, originally 8 (H, u_in , T _∞, T_s , k, ρ, μ, and C_p ), the nondimensionalized problem now involves only just two parameters (Re and Pr). Hence, for a given geometry, problems having the same values for the similarity parameters will have identical solutions. Another advantage of the use of similarity parameters is that results from a large number of experiments can be grouped and reported conveniently in terms of such parameters.

(b) The physical significance of the Reynolds number (Re) is investigated herein. As previously defined in Eq. (3.20), this dimensionless number requires the values of density (ρ) and dynamic viscosity (μ). For air, the values are ρ  =   1.2   kg/m³ and μ  =   2   ×   10^{−   5}  kg/m   s, whilst for water, they are ρ  =   1000   kg/m³ and μ  =   10^{−   3}  kg/m   s. If a laminar flow, for example Re  =   120, is assumed for both fluids, the inlet velocities for air and water are 0.02   m/s and 0.0012   m/s.

Based on CFD simulations, the axial dimensional and nondimensional velocity profiles for air and water in the fully developed region are given in Figs 3.9.2 and 3.9.3.

Fig. 3.9.2. Dimensional axial velocity profiles in the fully developed region.

Fig. 3.9.3. Nondimensional axial velocity profiles in the fully developed region.

Discussion: Although the two flows have different fluid properties, one being air whilst the other is water, we obtain the same flow behaviour. This is because they have the same Reynolds numbers, and the nondimensional governing equations of the x-momentum and y-momentum are identical, which leads to the same numerical results. If we consider another fluid mechanics example that is the flow of air or water over a flat plate (Fig. 3.9.4) having different lengths but with the same inlet velocities at the same Reynolds numbers, these two geometrically similar bodies have the same physical phenomena since they have the same friction coefficients (C_f ).

Fig. 3.9.4. Two geometrically similar bodies having the same friction coefficients at the same Reynolds numbers.

(c) The physical significance of the Prandtl number (Pr) is subsequently investigated herein. As previously defined in Eq. (3.35), it requires the values of the kinematic viscosity (ν) and thermal diffusivity (α). For air, the values are ν  =   1.667   ×   10^{−   5}  m²/s and α  =   1.667   ×   10^{−   5}  m²/s whilst for water ν  =   1   ×   10^{−   6}  m²/s and α  =   1.435   ×   10^{−   7}  m²/s. This yields Pr  =   1 for air and Pr  ≈   7 for water.

Based on the CFD simulations, the nondimensional temperature profiles for air and water in the fully developed region are given in Fig. 3.9.5.

Fig. 3.9.5. Nondimensional temperature profiles in the fully developed region.

Discussion: With reference to the energy diffusion term, 1/Re Pr in the nondimensional energy equation, the Reynolds number of air and water is the same; however, the Prandtl number of air is found to be much less than water leading to heat flow being diffused more in air than in water. Therefore, with less diffusion, more advection is encouraged in water than in air, which leads to a higher temperature in the fully developed region. Nevertheless, we have the same momentum diffusion for the velocity fields since the Reynolds numbers are the same for both fluids in the x-momentum and y-momentum equations; both fluids therefore have the same velocity profiles.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780081011270000039

CFD Solution Procedure: A Beginning

Jiyuan Tu , ... Chaoqun Liu , in Computational Fluid Dynamics (Third Edition), 2018

2.2.1 Creation of Geometry—Step 1

The first step in any CFD analyses is the definition and creation of geometry of the flow region, that is the computational domain for the CFD calculations. Let us consider two flow cases: a fluid flowing between two stationary parallel plates and a fluid passing through two cylinders in an open surrounding. It is important that the reader should always acknowledge the real physical flow representation of the problem that is to be solved as demonstrated by the respective physical domains in Figs 2.4 and 2.5. For the purpose of illustration, we designate the former and latter cases as Case 1 and Case 2. We shall assume that the two cylinders in Case 2 have the same length as the width W of the overall domain. We shall also assume that the width W of both flows within the three-dimensional physical domains is sufficiently large in order that the flows can be taken to be invariant along this transverse direction. Hence, Case 1 and Case 2 can simply be considered as two-dimensional computational domains for CFD calculations. These two flow cases will be repeatedly taken as illustrative examples to demonstrate the various basic steps that are involved in preprocess, solver, and postprocess stages.

Fig. 2.4. Case 1: Fluid flowing between two stationary parallel plates.

Fig. 2.5. Case 2: Fluid passing over two cylinders in an open surrounding.

There are certain distinct dissimilarities in the nature of these two flow problems. Case 1 represents an internal flow problem, whilst Case 2 is taken typically of an external flow scenario. In both cases, the fluid enters at the left boundary and exits at the right boundary of the computational domains. The main difference between these two flows is accentuated by the top and bottom boundaries, which brings about the classification of internal and external flows. In Case 1, the fluid flow is bounded within a domain of rigid walls as represented by the horizontal external walls of the two stationary parallel plates. That is not the characteristic of the fluid flow in Case 2 as the fluid can either take the inflow or outflow boundary conditions at the top and bottom boundaries.

One important aspect that the reader should always take note in the creation of the geometry for CFD calculations is to allow the flow dynamics to be sufficiently developed across the length L of these computational domains. For Case 1, we require the flow to be fully developed as it exits the domain. The physical interpretation and meaning of the concept of fully developed will be expounded in Chapter 3. For Case 2, we require to encapsulate the occurrence of complex wake-making development that persists behind the two cylinders as the flow passes over these cylinders. This phenomenon is analogous to the formation and shedding of vortices that are commonly experienced for flow past a cylinder. In this particular case, the top and bottom boundary effects may influence the flow passing over these two cylinders; the height H of the domain needs to be prescribed at a distance to sufficiently remove any of these boundary effects on the fluid flow surrounding the two cylinders but still manageable for CFD calculations. More practical guidance on this issue will be addressed in Chapter 6.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780081011270000027

Governing Equations for CFD—Fundamentals

Jiyuan Tu , ... Chaoqun Liu , in Computational Fluid Dynamics (Second Edition), 2013

3.2.1 Mass Conservation

One conservation law that is pertinent to fluid flow is matter may be neither created nor destroyed. Consider an arbitrary control volume V fixed in space and time as shown in Figure 3.1. The fluid moves through the fixed control volume, flowing across the control surface. Mass conservation requires that the rate of change of mass within the control volume is equivalent to the mass flux crossing the surface S of volume V (we define V as velocity, which is different from italic V, which represents a volume). In integral form,

Figure 3.1. Finite-control volume fixed in space.

(3.1) $\frac{d}{\mathit{dt}}\underset{V}{\int }\rho \mathit{dV}=-\underset{S}{\int }\rho V\cdot \mathbf{n}\mathit{dS}$

where n is the unit normal vector. We can apply Gauss's divergence theorem, which equates the volume integral of a divergence of a vector into an area integral over the surface that defines the volume. This is stated as

(3.2) $\underset{V}{\int }\mathit{div}\rho \text{V}\mathit{dV}=-\underset{S}{\int }\rho V\cdot \mathbf{n}\mathit{dS}$

Using the above theorem, the surface integral in Eq. (3.1) may be replaced by a volume V integral; hence the equation becomes

(3.3) $\underset{V}{\int }\left[\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \text{V}\right)\right]\mathit{dV}=0$

where ∇ ⋅ (ρ V) ≡ div ρV. Since Eq. (3.3) is valid for any size of volume V, the implication is that

(3.4) $\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho \text{V}\right)=0$

Equation (3.4) is the mass conservation. In the Cartesian coordinate system, it can be expressed as

(3.5) $\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0$

where the fluid velocity V at any point in the flow field is described by the local velocity components u, v, and w, which are, in general, functions of location (x, y, z) and time (t).

Alternatively, consider the scenario of a fluid flowing between two stationary parallel plates, as illustrated in Figure 3.2. An infinitesimally small control volume ΔxΔyΔz fixed in space (enlarged to the right of the figure) is analyzed where the mass conservation statement applies to the (u, v, w) flow field. Transport due to such motion is often referred to as advection. The conservation law requires that, for unsteady flow, the rate of increase of mass within the fluid element equals the net rate at which mass enters the control volume (inflow – outflow); in other words,

Figure 3.2. The conservation of mass in an infinitesimal control volume of a fluid flow between two stationary parallel plates.

(3.6) $\frac{\mathit{dm}}{\mathit{dt}}={\displaystyle \sum _{\mathit{in}}\dot{m}}-{\displaystyle \sum _{\mathit{out}}\dot{m}}$

The rate at which mass enters the control volume through the surface perpendicular to x may be expressed as (ρu) Δy Δz, where ρ is the local density of the fluid, and similarly through the surfaces perpendicular to y and z as (ρv) ΔxΔz and (ρw) ΔxΔy, respectively. The rate at which the mass leaves the surface at x + Δx may be expressed through Taylor expansion as

(3.7) $\left[\left(\rho u\right)+\frac{\partial \left(\rho u\right)}{\partial x}\mathrm{\Delta }x\right]\mathrm{\Delta }y\mathrm{\Delta }z+O\left(\mathrm{\Delta }x,\mathrm{\Delta }V\right)\text{where}\mathrm{\Delta }V=\mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z$

Similarly, the rate at which mass leaves the surfaces at y + Δy and z + Δz may also be expressed as

$\left[\left(\rho v\right)+\frac{\partial \left(\rho v\right)}{\partial y}\mathrm{\Delta }y\right]\mathrm{\Delta }x\mathrm{\Delta }z+O\left(\mathrm{\Delta }y,\mathrm{\Delta }V\right)$

and

(3.8) $\left[\left(\rho w\right)+\frac{\partial \left(\rho w\right)}{\partial z}\mathrm{\Delta }z\right]\mathrm{\Delta }x\mathrm{\Delta }y+O\left(\mathrm{\Delta }z,\mathrm{\Delta }V\right)$

Since the mass of the fluid element m is given by ρ ΔxΔyΔz, Eq. (3.6) becomes

(3.9) $\begin{array}{l}\frac{\partial \left(\rho \mathrm{\Delta }x\mathrm{\Delta }y\mathrm{\Delta }z\right)}{\partial t}=\left(\rho u\right)\mathrm{\Delta }y\mathrm{\Delta }z+\left(\rho v\right)\mathrm{\Delta }x\mathrm{\Delta }z+\left(\rho w\right)\mathrm{\Delta }x\mathrm{\Delta }y\\ -\left[\left(\rho u\right)+\frac{\partial \left(\rho u\right)}{\partial x}\mathrm{\Delta }x\right]\mathrm{\Delta }y\mathrm{\Delta }z-\left[\left(\rho v\right)+\frac{\partial \left(\rho v\right)}{\partial y}\mathrm{\Delta }y\right]\mathrm{\Delta }x\mathrm{\Delta }z\\ -\left[\left(\rho w\right)+\frac{\partial \left(\rho w\right)}{\partial z}\mathrm{\Delta }z\right]\mathrm{\Delta }x\mathrm{\Delta }y+\mathrm{\Delta }VO\left(\mathrm{\Delta }x,\mathrm{\Delta }y,\mathrm{\Delta }z\right)\end{array}$

In the limit, canceling terms and dividing by the constant-size ΔxΔyΔz, we obtain

(3.10) $\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho u\right)}{\partial x}+\frac{\partial \left(\rho v\right)}{\partial y}+\frac{\partial \left(\rho w\right)}{\partial z}=0$

Equation (3.10) is exactly the same form as derived in Eq. (3.5). This equation is precisely the partial differential form of the continuity equation. We have shown that the integral form in Eq. (3.1) can, after some manipulation, yield the partial differential form. This specific differential form is usually called the conservation form. Both Eq. (3.1) and Eq. (3.10) are in conservation form; the manipulation performed does not alter the situation.

In Chapter 2, a two-dimensional CFD analysis is performed for a channel flow that is described by the fluid flow between two stationary parallel plates (see Figure 3.2). This is made possible by the assumption that the dimension in the z coordinate direction is sufficiently large that the flow remains invariant along this coordinate direction. Since the fluid is taken to be incompressible, the density ρ is constant, i.e., the spatial and temporal variations in density are neglected relative to those velocity components of u, v, and w. We can obtain the continuity equation in two dimensions for an incompressible flow as

or

(3.11) $\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}=0$

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780080982434000032

Practical Guidelines for CFD Simulation and Analysis

Jiyuan Tu , ... Chaoqun Liu , in Computational Fluid Dynamics (Second Edition), 2013

6.3.4 Guidelines for Wall Boundary Conditions

Wall boundary conditions are generally employed for solid walls bounding the flow domain. Here, care should always be taken to ensure that the boundary conditions imposed on these walls are consistent with both the numerical model used and the actual physical features of the flow geometry. For fluid flowing between two stationary parallel plates ( Figure 6.17), wall boundary conditions are enforced for the channel walls. Nonetheless, these boundary conditions are also often used to bound fluid and solid regions, as was applied for the presence of the step within the flow domain in Figure 6.18 or the two cylinders inside the open surrounding flow environment in Figure 6.19.

For stationary walls, the default consideration is to assume that the no-slip condition applies, which simply means that the velocities are taken to be zero at the solid boundaries. This condition implies that the fluid flow comes to rest at the solid walls. We may also explore the possibility of modeling the boundary conditions on the solid walls as a free-slip condition, which assumes that the flow is parallel to the wall at this point. This condition corresponds to the absence of viscous effects in the continuum equations and is applied to the problem where the continuum approach breaks down as the fluid approaches the wall in viscous flow. For a general fluid-flow case, where the transport of heat takes precedence within the domain, great care must be taken to specify boundary conditions (for example, adiabatic walls or local heat fluxes) on the solid boundaries of the numerical model that properly represent the heat-transfer characteristics of the solid walls in the actual physical model.

In flow cases with moving or rotating walls, it is important that the boundary conditions that need to be specified are consistent with the motion of the solid walls. A lid-driven cavity and a rotating cylinder in a fluid environment, as shown in Figure 6.23, are some typical examples of the moving-boundary problems in CFD. These problems allow a positive or negative tangential velocity to be imposed at the top boundary of the lid-driven cavity or a clockwise or anticlockwise rotational speed to be specified on the circumferential surface of the rotating cylinder. Other, more complex CFD problems may require the use of sliding or moving meshes to better emulate the motion of a rotating impeller stirring the fluid in a tank or ocean waves hitting a ship's hull. The use of sliding or moving meshes is beyond the scope of this book. Interested readers are advised to consult the literature for more information (see for example, Ferziger and Perić, 1999).

Figure 6.23. Illustrative fluid-flow examples with moving or rotating walls.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780080982434000068

CFD Techniques—The Basics

Jiyuan Tu , ... Chaoqun Liu , in Computational Fluid Dynamics (Second Edition), 2013

4.5 Pressure–velocity coupling—"simple" scheme

The incompressible form of the conservation equations governing fluid flow are derived in Chapter 3 and are summarized in Section 3.6. Because of the incompressible assumption, the solution of the governing equations is complicated by the lack of an independent equation for pressure. In each of the momentum equations, fluid flow is driven by the contribution of the pressure gradients. With the additional equation provided by the continuity equation, this system of equations is self-contained; there are four equations for four dependents, u, v, w, and p, but no independent transport equation for pressure. The implication here is that the continuity and momentum equations are all that are required to solve for the velocity and pressure fields in an incompressible flow. For such a flow, the continuity equation is a kinematic constraint on the velocity field rather than a dynamic equation. In order to link the pressure with the velocity for an incompressible flow, one possible way is to construct the pressure field so as to guarantee conservation of the continuity equation.

In this section, we describe the basis of one of the most popular schemes of pressure–velocity coupling for an incompressible flow. It belongs to the class of iterative methods that is embodied in a scheme called SIMPLE, where the acronym stands for Semi-Implicit Method for Pressure-Linkage Equations. The SIMPLE scheme was developed for practical engineering solutions by Patankar and Spalding (1972). Ever since their pioneering work, it has found widespread application in the majority of commercial CFD codes. In the SIMPLE scheme, a guessed pressure field is used to solve the momentum equations. A pressure correction equation, deduced from the continuity equation, is then solved to obtain a pressure correction field, which in turn is used to update the velocity and pressure fields. These guessed fields are progressively improved through the iteration process until convergence is achieved for the velocity and pressure fields. The salient features of the SIMPLE scheme and the assembly of the complete iterative procedure are discussed later.

Variable arrangement on the grid: Before describing the SIMPLE scheme, the choice of arrangement on the grid requires some consideration. Among the many arrangements, the two most popular that have gained wide acceptance are the staggered and non-staggered grid arrangements.

The aim of having a staggered grid arrangement for CFD computations is to evaluate the velocity components at the control volume faces while the rest of the variables governing the flow field, such as the pressure, temperature, and turbulent quantities, are stored at the central node of the control volumes. A typical arrangement is depicted in Figure 4.10 for a structured finite-volume grid, and it can be demonstrated that the discrete values of the velocity component, u, from the x-momentum equation are evaluated and stored at the east, e, and west, w, faces of the control volume. Evaluation of the other velocity components using the y-momentum and z-momentum equations on the rest of the control volume faces allows a straightforward evaluation of the mass fluxes that are used in the pressure correction equation. This arrangement therefore provides a strong coupling between the velocities and pressure, which helps to avoid some types of convergence problems and oscillations in the pressure and velocity fields. Historically, staggered grid arrangement enjoyed its dominance within the CFD framework between the 1960s and 1980s. However, as the use of non-orthogonal grids became commonplace because of the need to handle complex geometries, alternative grid arrangements had to be explored because of some inherent difficulties in the staggered approach. In particular, if the staggered approach is used in generalized coordinates, curvature terms are required to be introduced into the equations that are usually difficult to treat numerically and may create non-conservative errors when the grid is not smooth.

Figure 4.13. Schematic representation of a multi-grid method using a V-cycle.

Figure 4.14. Schematic representation of a multi-grid method using a W-cycle.

Figure 4.10. Staggered and co-located arrangements of velocity components on a finite-volume grid (full symbols denote element vertices and open symbols at the center of the control volumes denote computational nodes for the storage of other governing variables).

Nowadays, the alternative grid arrangement that is frequently adopted in many commercial CFD codes is the non-staggered grid arrangement. Here, all the flow-field variables including the velocities are stored at the same set of nodal points. For the finite-volume grid shown in Figure 4.11, they are stored at the central node of the control volumes (open symbols). The co-located arrangement offers significant advantages in complicated domains, especially the capability of accommodating slope discontinuities or boundary conditions that may be discontinuous. Furthermore, if multi-grid methods are used, the co-located arrangement allows the ease of transfer of information between various grid levels for all the variables. This grid arrangement was out of favor for incompressible flow computation for a substantial period because of the difficulties faced in coupling the pressure with the velocity and the occurrence of oscillations in the pressure. Nevertheless, widespread use of the co-located grid arrangement became prominent once again through significant developments of the pressure–velocity coupling algorithms, such as the well-known Rhie and Chow (1983) interpolation scheme. This scheme, which has provided physically sensible solutions on structured co-located meshes, generated much interest for unstructured meshing applications. For the vast majority of general flows, this treatment ties together the pressure fields to yield smooth solutions, while only minimally affecting the mass fluxes. More details of this interpolation scheme are left to the interested reader.

Figure 4.11. Arrangement of velocity components on a control volume element of a structured grid at the central node, element faces, and element vertices.

Pressure correction equation and its solution: The SIMPLE scheme is essentially a guess-and-correct procedure for the calculation of pressure through the solution of a pressure correction equation. The method is illustrated by considering a two-dimensional steady laminar flow problem in a structured grid, as shown in Figure 4.11.

The SIMPLE scheme provides a robust method of calculating the pressure and velocities for an incompressible flow. When coupled with other governing variables, such as temperature and turbulent quantities, the calculation needs to be performed sequentially, since it is an iterative process. The sequence of operations in a typical CFD iterative process that embodies the SIMPLE scheme is given in Figure 4.12, with more details of each iterative step elaborated below.

Figure 4.12. The SIMPLE scheme.

Step 1: The iterative SIMPLE calculation process begins by guessing the pressure field, p*. During the iterative process, the discretized momentum equations are solved using the guessed pressure field. Applying the finite-volume method, the equations for the x-momentum and y-momentum that yield the velocity components, u* and v*, can be expressed in the same algebraic form as previously derived in Eq. (4.29), which can be recast into

(4.75) $a_P^u{u}_P^{*}=\sum a_{\mathit{nb}}^u{u}_{\mathit{nb}}^{*}-\frac{\partial p^{*}}{\partial x}\Delta V+b\text{'}$

(4.76) $a_p^v{v}_P^{*}=\sum a_{\mathit{nb}}^v{v}_{\mathit{nb}}^{*}-\frac{\partial p^{*}}{\partial y}\Delta V+b\text{'}$

where ΔV is the finite-control volume. Here, we simplify the above expressions by introducing a _nb to represent the presence of the neighboring coefficients and u _nb * and v _nb * to denote the neighboring nodal velocity components. The pressure gradient terms appearing in the above two equations are taken from the original source term b of the momentum equations while the other terms governing the fluid flow are left in the source term b′.

Step 2: If we define the correction p′ as the difference between the correct pressure field and the guessed pressure field p*, we get

(4.77) $p=p^{*}+p\text{'}$

Similarly, we can also define the corrections u′ and v′ to relate the correct velocities u and v to the guessed velocities u* and v*:

(4.78) $u=u^{*}+u\text{'}$

(4.79) $v=v^{*}+v\text{'}$

The algebraic form of the correct velocities u and v can also be expressed similarly to Eqs. (4.75) and (4.76), so that

(4.80) $a_P^u{u}_p=\sum a_{\mathit{nb}}^u{u}_{\mathit{nb}}-\frac{\partial p}{\partial x}\Delta V+b\text{'}$

(4.81) $a_p^v{v}_p=\sum a_{\mathit{nb}}^v{v}_{\mathit{nb}}-\frac{\partial p}{\partial y}\Delta V+b\text{'}$

Subtracting Eqs. (4.75) and (4.76) from Eqs. (4.80) and (4.81), we obtain

(4.82) $a_P^u\left(u_p-u_p^{*}\right)=\sum a_{\mathit{nb}}^u\left(u_{\mathit{nb}}-u_{\mathit{nb}}^{*}\right)-\frac{\partial \left(p-p^{*}\right)}{\partial x}\Delta V$

(4.83) $a_P^u\left(v_p-v_p^{*}\right)=\sum a_{\mathit{nb}}^u\left(v_{\mathit{nb}}-v_{\mathit{nb}}^{*}\right)-\frac{\partial \left(p-p^{*}\right)}{\partial y}\Delta V$

Using the correction formulae (4.77–4.79), the above equations may be rewritten as follows:

(4.84) $a_P^{u}u{\text{'}}_p=\sum a_{\mathit{nb}}^{u}u{\text{'}}_{\mathit{nb}}-\frac{\partial p\text{'}}{\partial x}\Delta V$

(4.85) $a_P^{v}v{\text{'}}_p=\sum a_{\mathit{nb}}^{v}v{\text{'}}_{\mathit{nb}}-\frac{\partial p\text{'}}{\partial y}\Delta V$

The SIMPLE scheme approximates Eqs. (4.84) and (4.85) by the omission of the terms ∑ a _nb ^u u′ _nb and ∑ a _nb ^v v′ _nb . The reader is reminded that this scheme is an iterative approach, and there is no reason why the formula designed to predict p′ needs to be physically correct. Hence, we are allowed to construct a formula for p′ that is simply a numerical artifice with the aim of expediting the convergence of the velocity field to a solution that satisfies the continuity equation. This is the essence of the algorithm. Once the pressure correction field is known, the correct velocities u and v can be obtained through the guessed velocities u* and v* from simplified Eqs. (4.84) and (4.85):

(4.86) $u_p=u_P^{*}-D^u\frac{\partial p\text{'}}{\partial x}$

(4.87) $v_p=v_P^{*}-D^v\frac{\partial p\text{'}}{\partial y}$

where

(4.88) $D^u=\frac{\Delta V}{a_p^u}\text{and}{D}^v=\frac{\Delta V}{a_p^v}$

Although Eqs. (4.86) through (4.88) have been developed to correct the velocities from the guessed velocities at the central node of the control volume, these correction formulae can also be generally applied to any location where the velocity components reside within the computational grid (as shown in Figure 4.9, the velocities may be located at central node P or at the control volume faces or at the vertices of the control volume). The general form of the velocity correction formulae, by removing the subscript P, can be expressed as

(4.89) $u=u_{}^{*}-D^u\frac{\partial p\text{'}}{\partial x}$

(4.90) $v=v_{}^{*}-D^v\frac{\partial p\text{'}}{\partial y}$

The derivation of a pressure correction equation utilizes the above two equations. Differentiating Eq. (4.89) by the Cartesian direction x and Eq. (4.90) by the Cartesian direction y and summing them together yields

(4.91)

By invoking the continuity equation, it is shown that the term represented by the source term of the right-hand side of Eq. (4.91) is zero, and Eq. (4.91) can be re-arranged as

(4.92) $\frac{\partial }{\partial x}\left(D^u\frac{\partial p\text{'}}{\partial x}\right)+\frac{\partial }{\partial y}\left(D^v\frac{\partial p\text{'}}{\partial y}\right)=\underset{massresidual}{\underbrace{\left(\frac{\partial u^{*}}{\partial x}+\frac{\partial v^{*}}{\partial y}\right)}}$

Interestingly, Eq. (4.92) behaves like a steady-state diffusion process in a two-dimensional domain. It is a Poisson equation—one of the well-known equations from classical physics and mathematics. The solution to this Poisson equation can be achieved through some efficient numerical solvers (conjugate-gradient and multi-grid methods), as previously discussed, to accelerate its convergence.

Step 3: Once the pressure correction p′ field is obtained, the pressure and velocity components are subsequently updated through the correction formulae of Eqs. (4.77), (4.89), and (4.90). If the solution concerns only a laminar CFD flow problem, the iteration process proceeds directly to check the convergence of the solution. If the solution is not converged, the process is repeated by returning to Step 1. The source term appearing in the pressure correction equation, Eq. (4.92), commonly known as the mass residual, is normally used in CFD computations as a criterion for terminating the iteration procedure. As the mass residual continues to diminish, the pressure correction p′ will be zero, thereby yielding a converged solution of p*   = p, u*   = u, and v*   = v.

Step 4: This step is executed if the CFD flow problem is turbulent or if it involves the transfer of heat or mass exchanges between different flow phases. Additional transport equations governing such a flow system need to be solved before convergence is checked. If the solution is not converged, the iterative process returns to Step 1 and repetitive calculations are carried out until convergence is reached.

The application of this SIMPLE scheme is best illustrated by solving the Chapter 2 CFD problem of a steady two-dimensional incompressible laminar flow in a channel, which is described in Example 4.5.

Example 4.5

Consider the case of a steady, two-dimensional, incompressible, laminar flow between two stationary parallel plates, as in Chapter 2. By obtaining the solution from a CFD code using the finite-volume method, track the progress of the intermediate values of u, v, p, p′ and the mass residual during the iterative process at a computational nodal point at the center of the channel (Figure 4.5.1).

Figure 4.5.1. Two-dimensional laminar flow in a channel with a monitoring point located at the center of the channel.

Solution

The problem is described as follows: To demonstrate the robustness of the SIMPLE scheme, the iterative process begins by employing the initial guesses: p*   =   0, u*   =   0 and v*   =   0. The discretized equations governing the momentum and pressure correction are solved using the default iterative solvers provided in the commercial CFD code. The inlet, outlet, and wall conditions remain the same as applied in Chapter 2.

Based on Eqs. (4.75), (4.76), (4.77), (4.89), (4.90), and (4.92), the calculated values of the pressure p, pressure correction p′, velocities u and v, and mass residual for the first iteration at the monitoring point are

The solution of the first iteration from above is subsequently used as intermediate values for the next iteration step; the second iteration yields

After repeated applications of the iterative process, the respective values for the pressure p, pressure correction p′, velocities u and v, and mass residual after 10 and 20 iterations are

and

From a theoretical viewpoint, the vertical velocity v is zero at the monitoring point and the iterative process confirms the trend of the prediction toward the zero value. It is seen during the iterative process that the intermediate values of this velocity are much smaller than the rest of the other governing variables; the convergence history plot of this velocity component is therefore omitted, since no quantitative comparison can be realized against the other variable convergence histories. The convergence histories for the rest of the governing variables, which include the pressure p, pressure correction p′, horizontal velocity u, and mass residual are illustrated in Figures 4.5.2 through 4.5.5.

Figure 4.5.2. Convergence history of the pressure p.

Figure 4.5.3. Convergence history of the pressure correction p′.

Figure 4.5.4. Convergence history of the horizontal velocity u.

Figure 4.5.5. Convergence history of the mass residual.

Discussion

From this worked example of a channel flow, the SIMPLE scheme provides an efficient iterative procedure for obtaining the velocity and pressure fields for an incompressible flow. The SIMPLE scheme is a robust method that produces rapid stabilization of the velocity and pressure, as seen by their respective convergence histories after 5 iterations. The mass residual that appears as a source term in the pressure correction equation, Eq. (4.80), continues to diminish during the iteration process, thus reaffirming conservation in the continuity equation. Subsequently, the pressure correction p′ is seen to be approaching zero. Hence, the corrections that are required to update the velocity field are also approaching zero. The trend of the convergence histories favors the likelihood of a converged steady-state solution.

The reader should be aware of other types of pressure–velocity coupling algorithms that employ a philosophy similar to the SIMPLE algorithm and that are employed by CFD users or are adopted in commercial CFD codes. These variant SIMPLE algorithms have been formulated with the aim of improving the robustness and convergence rate of the iterative process. We do not intend to provide the reader with all the details of the available algorithms, but we will briefly indicate and describe the modifications made to the original SIMPLE algorithm.

The SIMPLEC (SIMPLE-Consistent) algorithm by Van Doormal and Raithby (1984) follows the same iterative steps as in the SIMPLE algorithm. The main difference between SIMPLEC and SIMPLE is that the discretized momentum equations are manipulated so that the SIMPLEC velocity correction formulae omit terms that are less significant than those omitted in SIMPLE. Another pressure correction procedure that is also commonly employed is the PISO (Pressure Implicit with Splitting of Operators) algorithm proposed by Issa (1986). This pressure–velocity calculation procedure was originally developed for non-iterative computation of unsteady compressible flows. Nevertheless, it has been adapted successfully for the iterative solution of steady-state problems. PISO is simply an extension of SIMPLE with an additional corrector step that involves an additional pressure correction equation to enhance the convergence. The SIMPLER (SIMPLE-Revised) algorithm developed by Patankar (1980) also falls within the framework of two corrector steps, as in PISO. Here, a discretized equation for the pressure provides the intermediate pressure field before the discretized momentum equations are solved. A pressure correction is later solved where the velocities are corrected through the correction formulae, derived similarly to those in the SIMPLE algorithm.

There are other SIMPLE-like algorithms, such as SIMPLEST (SIMPLE-ShorTened) of Spalding (1980) or SIMPLEX of Van Doormal and Raithby (1985) or SIMPLEM (SIMPLE-Modified) of Archarya and Moukalled (1989), that share the same essence in their derivations. More details of all the above pressure–velocity coupling algorithms are left to interested readers.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780080982434000044

Viscous Fluid Flow

Nikolaos D. Katopodes , in Free-Surface Flow, 2019

Problems

5-1.

Derive the continuity equation for an incompressible fluid in a flow field that is density stratified, i.e. the density is a function of the spatial coordinates.

5-2.

Derive the continuity equation of an incompressible fluid using material coordinates as the independent variables. Hint: transform all variables to material coordinates, and integrate over time.

5-3.

Derive the continuity equation for a mixture of two substances in the presence of diffusion. Hint: consider the continuity equation for each substance, and then combine the results.

5-4.

Determine the conditions for the existence of the stream function in compressible flow.

5-5.

Consider the flow in a well-mixed estuary that can be roughly approximated by a unidirectional flow field described initially by a uniform velocity $u=1.0m/s$ . Then, due to some unusual meteorological conditions, the associated density field undergoes a temporal variation such that $\rho (t)={\rho }_0\cos \omega t$ . Compute the velocity field for $t>0$ , and specifically, determine the velocity at $t=2.0s$ and $x=1.5m$ . Use $\omega =1.0rad/s$ , and assume that $u(0,t)=1.0m/s$ for $t>0$ .

5-6.

A cylindrical container is filled with a fluid whose density is linearly stratified, i.e. it increases linearly with the depth. The container is suddenly dropped in a way that maintains its upright position, and thus the fluid is subjected to the gravitational acceleration, and acquires a vertical velocity $w=-gt$ . Determine the material derivative of the density during the fall.

5-7.

Consider a flow field described by the following expression for the stream function in polar coordinates

$\psi =U{r}_0^2(\frac{r^2{\sin }^2\theta }{r_0^2}-\cos \theta )$

where U and $r_0$ are constants. Determine if there is a stagnation point in the flow field, and if so, identify its coordinates.

5-8.

Why does a viscous fluid in a closed container eventually come to rest after being stirred momentarily?

5-9.

Find the viscous dissipation corresponding to laminar flow between two stationary parallel plates located distance B apart.

5-10.

Consider a steady, uniform horizontal flow between two stationary parallel plates located at $y=0$ and $y=B$ . The horizontal velocity components are described by $u=f(y)$ , and $w=0$ . Use the continuity equation to find an expression for $v(y)$ , i.e. the vertical component of the velocity field. Assume no-slip conditions apply on the plates.

5-11.

Consider an incompressible planar flow of a viscous fluid confined between two horizontal flat plates at a distance B apart. The upper plate moves parallel to itself at speed U and the lower plate is fixed. Let the x-axis lie on the lower plate. At steady state the velocity is given by $u=\frac{U}{B}y$ , where y is the vertical coordinate. Compute the components of the stress tensor.

5-12.

Consider an incompressible planar flow of a viscous fluid confined between two horizontal flat plates at a distance B apart. The upper plate moves parallel to itself at speed U and the lower plate is fixed. Let the x-axis lie on the lower plate. At steady state the flow is fully developed and independent of x. Show that the pressure distribution is hydrostatic in the vertical y direction.

5-13.

Consider an incompressible planar flow of a viscous fluid confined between two horizontal flat plates at a distance B apart. The upper plate moves parallel to itself at speed U and the lower plate is fixed. Assuming that the fluid is glycerine at ${20}^{\circ }$ , determine the average rate of viscous dissipation per unit mass when $B=0.01m$ and $U=0.01m/s$ .

5-14.

Consider an incompressible planar flow of a viscous fluid confined between two horizontal flat plates at a distance B apart. The upper plate moves parallel to itself at speed U and the lower plate is fixed.

a.

Compute the components of the strain rate tensor.

b.

Determine the principal directions of the strain rate tensor.

c.

Sketch the flow field and the new shape of an originally circular patch located half way between the two walls, as it is transported downstream.

5-15.

Consider the following flow field: $\mathbf{u}=(a{x}_2,0,0)$ , where a is a constant.

a.

Calculate the components of the rotation and strain rate tensors.

b.

Sketch a few streamlines and sketch what would happen to a fluid element over time.

c.

To what type of flow does this velocity field correspond?

5-16.

Consider the free-surface flow of a viscous fluid layer of depth H on top of a horizontal belt moving with a speed U. Assuming that the shear stress is zero at the free surface, determine the shear stress and velocity profiles. What would be the effect of applying a constant wind stress at the surface? What would be the effect if a constant pressure gradient is applied in the direction of flow?

5-17.

A two-dimensional stream of an incompressible fluid flows vertically down (parallel to the y-axis) and towards a horizontal solid boundary oriented along the x-axis. If the vertical velocity, v, is directly proportional to distance from the solid boundary, calculate and plot the distribution of the horizontal velocity, u, at various positions along the x-axis.

5-18.

A perforated cylindrical container of length L has a plunger attached to it, thus as the plunger is pushed, water is forced out of the cylinder. The container has a diameter D, the plunger has the same size, and no leakage is allowed through the plunger or the flat end of the container. The total area of the perforations is a fraction α of the cylindrical surface of the container. Using the macroscopic volume balance approach, derive an expression for the speed of the plunger as a function of the water velocity out of the perforations.

5-19.

A deep pond can be approximated as a cylindrical tank with cross-sectional area equal to $A=100m^2$ . The pond has an orifice in its bottom with cross-sectional area $a=10m^2$ . We need to empty the pond so we can dredge its bottom that contains contaminated sediments. The depth of water inside the pond, h, is a function of time, as the pond is being emptied. If initially the depth is equal to $h_0=10m$ , calculate the time needed for the pond to empty. Assume that the flow is irrotational.

5-20.

Consider an incompressible planar flow of a viscous fluid confined between two horizontal flat plates at a distance B apart. The upper plate moves parallel to itself at speed U and the lower plate is fixed. Let the x-axis lie on the lower plate. At steady state the flow is fully developed and independent of x. Compute the stream function.

5-21.

Consider a nearly spherical dust particle of diameter $D=50\mu m$ , falling in air whose temperature is $-{50}^{\circ }C$ , and whose pressure is $p_0=55kPa$ . The density of the particle is ${\rho }_s=1240kg/m^3$ . The density of the air is $\rho =0.86kg/m^3$ , and its viscosity is estimated to be $\mu =1.47\times {10}^{-5}kg/m\cdot s$ . Estimate the terminal velocity of the particle.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780128154892000058

CFD Mesh Generation: A Practical Guideline

Jiyuan Tu , ... Chaoqun Liu , in Computational Fluid Dynamics (Third Edition), 2018

4.4 Local Mesh Refinement

One local mesh refinement technique that is widely used in many CFD applications is the concept of a stretched grid in the near vicinity of domain walls. For a viscous flow bounded with solid boundaries, the need to cluster a large number of small cells within the physical boundary layer is more than just attempting to minimize the truncation error with the closely spaced grid points. Rather, it is a matter of utter importance that the actual flow physics is appropriately encapsulated.

Let us revisit the case for the fluid flowing between two stationary parallel plates as investigated in previous chapters. In a real physical flow, the existence of a developing boundary layer will grow in thickness as the fluid enters the left boundary and migrates downstream along the bottom wall of the domain such as illustrated in Fig. 4.13. By denoting the local thickness of the boundary layer as δ, where δ  = δ (x), it is evidently clear that the use of a coarse uniform mesh misses the physical boundary layer. The predicted viscous-like velocity profile shown at some point downstream in the fully developed region is simply due to the application of the no-slip boundary condition at the bottom wall. In contrast, the coarse stretched grid at the very least captures some of the essential features of the actual physical boundary layer. It is therefore not surprising that the accuracy of the computational solution is greatly influenced by the grid distribution inside the boundary layer region. We can further investigate this by imposing the concept of consistency to be discussed in Chapter 5 for the incompressible mass conservation in two dimensions where the truncation error yields the following:

Fig. 4.13. Two schematic illustrations demonstrating the need for local refinement in the near vicinity of the bottom wall to resolve the physical boundary layer.

(4.9) $\frac{\mathrm{\Delta }{x}^2}{6}\frac{{\partial }^{3}u}{\partial x^3}+\frac{\mathrm{\Delta }{y}^3}{6}\frac{{\partial }^{3}v}{\partial y^3}$

If the solution error is expected to follow the truncation error, it is imperative that the grid needs to be appropriately refined in order to sufficiently resolve the steep gradient of the velocity profile that exists within this region. This will help to minimize the solution error associated with the truncation error.

Local mesh refinement is also important to better resolve specific fluid dynamics problems such as upward stagnation flow and the backward-facing step geometry. The latter is one of the basic geometries used commonly in many engineering applications to better understand the phenomena of flow separation, flow reattachment, and free shear jet. The placement of stretched grids along both the horizontal and vertical directions to encapsulate the essential feature of the recirculation vortex is illustrated in Fig. 4.14.

Fig. 4.14. A schematic drawing for the backward-facing step geometry and the computational grid to capture the essential feature of the recirculation vortex.

Whilst applying the stretched grid described above, one useful consideration is to exercise care in avoiding sudden changes of the grid size away from the domain boundary. The mesh spacing should be continuous, and grid-size discontinuities should be removed as much as possible in regions of large flow changes, particularly when dealing with multiblock meshing of arbitrary mesh coupling, nonmatching cell faces, or extended changes of element types. Discontinuity in the grid size destabilizes the numerical procedure due to the accumulation of truncation errors in the critical flow regions. These errors usually contain the diffusive terms (second-order derivatives) where the discretization imposed on these derivatives requires very smooth grid changes. Making sure that the grid changes slowly and smoothly away from the domain boundary and within the domain interior will assist in overcoming divergence tendencies of the numerical calculations. It is also worthwhile noting that most inbuilt mesh generators in commercial codes and independent grid generation packages have the means of prescribing suitable mesh stretching or expansion ratios (rates of change of cell size for adjacent cells). The specification of these ratios should always be negotiated within the codes' requirements whilst generating the appropriate stretched grid.

By the mere construction of a stretched grid, the local refinement technique provides the possibilities of allocating additional grid nodal points to resolve the important fluid flow action and reducing or removing the grid nodal points from other regions where there is little or no action. Nevertheless, it should be noted that a stretched grid is an algebraically generated grid prescribed prior to the solution of the flow field being calculated. The question that needs to be carefully addressed is whether the generated stretched grid sufficiently captures the major fluid flow action or whether the real flow action is far away from the intended significant flow activity to be resolved by the generated stretched grid region that is not known a priori.

In the event where specific interest on the finer details of the local flow physics are required to be predicted, the use of embedded meshes can also be adopted. Here, such meshes do not necessarily have to be constrained to have matching cell faces. Fig. 4.15 illustrates an example of embedded Cartesian mesh without matching cell faces near the bottom boundary for the rectangular-type geometry. We can view this type of mesh as a special case of a local mesh refinement strategy. It is noted that the locally refined region in Fig. 4.15 could also have been achieved through the use of other types of elements such as triangular or polyhedral elements or a combination of different types of elements.

Fig. 4.15. An embedded Cartesian mesh without matching cell faces near the bottom boundary.

Read full chapter

URL:

https://www.sciencedirect.com/science/article/pii/B9780081011270000040

murphyuppound.blogspot.com

Source: https://www.sciencedirect.com/topics/engineering/stationary-parallel-plate

Share this post