김경종
491 lines
30 KiB
Markdown
491 lines
30 KiB
Markdown
<!-- source-page: 691 -->
|
|
|
|
To demonstrate the differences in the Lagrangian and Eulerian formulations consider the following example.
|
|
|
|
EXAMPLE 7.10: The motions of fluid particles in a duct are given by
|
|
|
|
$$
|
|
^ t x _ {1} = - 5 + \sqrt {2 5 + 1 0 ^ {0} x _ {1} + (^ {0} x _ {1}) ^ {2} + 4 t} \tag {a}
|
|
$$
|
|
|
|
(i) Calculate the velocities and accelerations of the particles. Express your results in the Lagrangian form $\dot{u}_{1}=f_{1}(^{0}x_{1},t)$ , $\ddot{u}_{1}=f_{2}(^{0}x_{1},t)$ .
|
|
(ii) Eliminate $^0 x_1$ from your expressions in (i) to obtain the spatial expressions for the velocities and accelerations.
|
|
(iii) Show that your expression for the acceleration in (ii) may also be obtained by combining the local acceleration and convective acceleration (i.e., by use of the usual Eulerian expression).
|
|
|
|
To obtain the result requested in (i) we simply differentiate the expression in (a) with respect to time and note that $d^{t}x_{1}/dt = d^{t}u_{1}/dt$ (because $d^{0}x_{1}/dt = 0$ )
|
|
|
|
$$
|
|
\frac {d ^ {\prime} u _ {1}}{d t} = \frac {2}{\left[ 2 5 + 1 0 ^ {0} x _ {1} + (^ {0} x _ {1}) ^ {2} + 4 t \right] ^ {1 / 2}} \tag {b}
|
|
$$
|
|
|
|
Similarly, $\frac{d^2 u_1}{dt^2} = \frac{-4}{[25 + 10^0 x_1 + (^0 x_1)^2 + 4t]^{3/2}}$ (c)
|
|
|
|
To express (b) and (c) in terms of $x_{1}$ we might solve from (a) for $^0 x_1$ in terms of $x_{1}$ and $t$ . However, in this simple example we note
|
|
|
|
$$
|
|
^ \prime x _ {1} + 5 = \sqrt {2 5 + 1 0 ^ {0} x _ {1} + (^ {0} x _ {1}) ^ {2} + 4 t}
|
|
$$
|
|
|
|
and hence
|
|
|
|
$$
|
|
\frac {d ^ {t} u _ {1}}{d t} = \frac {2}{^ {t} x _ {1} + 5}; \quad \frac {d ^ {2} {} ^ {t} u _ {1}}{d t ^ {2}} = \frac {- 4}{(^ {t} x _ {1} + 5) ^ {3}} \tag {d}
|
|
$$
|
|
|
|
In the Eulerian formulation, we simply write $x_{1} \equiv 'x_{1}$ , and it is implied that $x_{1}$ can be any coordinate value. Of course, the expressions in (d) are also valid for any time. We should recognize that in (d) we focus our attention on the coordinate $x_{1}$ and measure the velocity and acceleration of the particles as they pass that coordinate. In this calculation we do not use (and generally do not care to know) the initial positions of the particles.
|
|
|
|
On the other hand, if we use (b) and (c), we focus our attention on the particles given by their initial positions and we measure the velocities and accelerations of these particles at a particular given time.
|
|
|
|
We may also obtain the acceleration in (d) from the general Eulerian expression of acceleration (see Table 7.4). Here we employ (without using the superscript)
|
|
|
|
$$
|
|
\frac {D v}{D t} = \frac {\partial v}{\partial t} + \frac {\partial v}{\partial x} v; \quad v \equiv \frac {d ^ {\prime} u _ {1}}{d x _ {1}}; \quad x _ {1} \equiv x
|
|
$$
|
|
|
|
which gives
|
|
|
|
$$
|
|
\frac {D v}{D t} = \frac {\partial}{\partial t} \left(\frac {2}{x + 5}\right) + \left[ \frac {\partial}{\partial x} \left(\frac {2}{x + 5}\right) \right] \left(\frac {2}{x + 5}\right) = 0 - \frac {4}{(x + 5) ^ {3}}
|
|
$$
|
|
|
|
Hence, the local derivative of the velocity is zero and the convective part of the acceleration is $-4/(x + 5)^{3}$ .
|
|
|
|
<!-- source-page: 692 -->
|
|
|
|
# 7.4.1 Continuum Mechanics Equations
|
|
|
|
Let us summarize the continuum mechanics equations for incompressible fluid flow including heat transfer. These equations are of course developed in detail in standard textbooks on fluid mechanics (see, for example, F. M. White [A] or H. Schlichting [A]), but we summarize the equations here to state the notation and provide the basis for the derivation of the governing finite element equations. We will note that in some respects the notation below is different from the notation we used in the analysis of solid mechanics problems (see Chapter 6) because, for example, the velocity $v_{i}$ is now the basic kinematic variable to be calculated instead of the displacement solved for in the analysis of solids.
|
|
|
|
Using a stationary Cartesian reference frame $(x_{i}, i = 1, 2, 3)$ , the governing equations of incompressible fluid flow within the domain V are, using index notation, the usual summation convention, and implying that the conditions at time t are considered without use of a superscript t (which is used in Table 7.4),
|
|
|
|
momentum: $\rho \left(\frac{\partial v_i}{\partial t} + v_{i,j}v_j\right) = \tau_{ij,j} + f_i^B$ (7.57)
|
|
|
|
constitutive: $\tau_{ij} = -p\delta_{ij} + 2\mu e_{ij}$ (7.58)
|
|
|
|
continuity: $v_{i,i} = 0$ (7.59)
|
|
|
|
heat transfer: $\rho c_{p}\left[\frac{\partial\theta}{\partial t} +\theta ,_{i}v_{i}\right] = (k\theta ,_{i}),_{i} + q^{B}$ (7.60)
|
|
|
|
# Here we have
|
|
|
|
$v_{i} =$ velocity of fluid flow in direction $x_{i}$
|
|
|
|
$\rho =$ mass density
|
|
|
|
$\tau_{ij} =$ components of stress tensor
|
|
|
|
$f_{i}^{B} =$ components of body force vector
|
|
|
|
p = pressure
|
|
|
|
$\delta_{ij} =$ Kronecker delta
|
|
|
|
$\mu =$ fluid (laminar) viscosity
|
|
|
|
$e_{ij} =$ components of velocity strain tensor $= \frac{1}{2} (v_{i,j} + v_{j,i})$
|
|
|
|
$c_{p} =$ specific heat at constant pressure
|
|
|
|
$\theta =$ temperature
|
|
|
|
k = thermal conductivity
|
|
|
|
$q^{B} = \text{rate of heat generated per unit volume [this term also includes the rate of heat dissipated } (=2\mu e_{ij}e_{ij})]$
|
|
|
|
The boundary conditions corresponding to (7.57) to (7.60) are
|
|
|
|
Prescribed fluid velocities $v_{i}^{S}$ on the surface $S_{v}$
|
|
|
|
$$
|
|
v _ {i} \mid_ {S _ {v}} = v _ {i} ^ {s} \tag {7.61}
|
|
$$
|
|
|
|
<!-- source-page: 693 -->
|
|
|
|
Prescribed tractions $f_{i}^{S}$ on the surface $S_{f}$ ,
|
|
|
|
$$
|
|
\tau_ {i j} n _ {j} | _ {s _ {f}} = f _ {i} ^ {s} \tag {7.62}
|
|
$$
|
|
|
|
where the $n_{j}$ are the components of the unit normal vector n (pointing outward) to the fluid surface and the $f_{i}^{S}$ are the components of the (physical) traction vector.
|
|
|
|
Prescribed temperatures $\theta^{s}$ on the surface $S_{\theta}$ ,
|
|
|
|
$$
|
|
\theta \mid_ {S _ {\theta}} = \theta^ {S} \tag {7.63}
|
|
$$
|
|
|
|
Prescribed heat flux into the surface $S_{q}$ ,
|
|
|
|
$$
|
|
k \left. \frac {\partial \theta}{\partial n} \right| _ {s _ {q}} = q ^ {s} \tag {7.64}
|
|
$$
|
|
|
|
where $q^{s}$ is the heat flux input to the body. We note that for (7.61) to (7.64) the discussions in Sections 4.2.1 and 7.2.1 are directly applicable.
|
|
|
|
The heat flux input in $(7.64)$ comprises the effect of actually applied distributed heat flow, and the effect of convection and radiation heat transfer. These applied heat fluxes are included in the analysis as discussed in Section 7.2.
|
|
|
|
Another form of $(7.57)$ is obtained if these momentum equations are written to also include the continuity condition, see Example 7.12. This form is referred to as the conservative form (because momentum conservation is explicitly imposed) and is largely used in finite volume discretization methods, see S. V. Patankar [A].
|
|
|
|
Equations (7.57) to (7.60) are the standard Navier-Stokes equations governing the motion of a viscous, incompressible fluid in laminar flow with heat transfer. Inherent nonlinearities are due to the convective terms in (7.57) and (7.60) and the radiation boundary condition in (7.64). Additional nonlinearities arise if the viscosity coefficient depends on the temperature or on the velocity strain, if the specific heat $c_{p}$ , conductivity k, and convection and radiation coefficients depend on temperature, and of course if turbulence descriptions are included (see, for example, W. Rodi [A]).
|
|
|
|
Before proceeding to develop the finite element equations, let us assume that $\mu$ , $c_{p}$ , and k are constant and rewrite (7.57) to (7.60) into standard forms that reflect some important characteristics of fluid flows. If we substitute (7.58) and (7.59) into (7.57), we obtain
|
|
|
|
$$
|
|
\rho (\dot {v} _ {i} + v _ {i, j} v _ {j}) = - p _ {, i} + f _ {i} ^ {B} + \mu v _ {i, j j} \tag {7.65}
|
|
$$
|
|
|
|
Hence (7.65) is the momentum equation containing the incompressibility condition.
|
|
|
|
Let us now define the nondimensionalized variables
|
|
|
|
$$
|
|
x _ {i} ^ {*} = \frac {x _ {i}}{L}; \quad v _ {i} ^ {*} = \frac {v _ {i}}{v}; \quad t ^ {*} = \frac {t v}{L}
|
|
$$
|
|
|
|
$$
|
|
p ^ {*} = \frac {p}{\rho v ^ {2}}; \quad f _ {i} ^ {B *} = \frac {f _ {i} ^ {B} L}{\rho v ^ {2}} \tag {7.66}
|
|
$$
|
|
|
|
$$
|
|
\theta^ {*} = \frac {\theta - \theta_ {0}}{\Delta \theta}; \quad q ^ {B *} = \frac {q ^ {B} L}{\rho c _ {p} \Delta \theta v}
|
|
$$
|
|
|
|
where $L, v, \theta_{0}$ , and $\Delta\theta$ are chosen characteristic values for length, velocity, temperature, and temperature difference of the problem considered. Using these nondimensionalized variables, we can rewrite the momentum and energy equations in the following forms.
|
|
|
|
<!-- source-page: 694 -->
|
|
|
|
$$
|
|
\dot {v} _ {i} ^ {*} + \underbrace {v _ {i , j} ^ {*} v _ {j} ^ {*}} _ {\text { convection }} = - p _ {, i} ^ {*} + f _ {i} ^ {B ^ {*}} + \underbrace {\frac {1}{\text { Re }} v _ {i , j j} ^ {*}} _ {\text { diffusion }} \tag {7.67}
|
|
$$
|
|
|
|
and $\dot{\theta}^{*} + \underbrace{\theta_{,i}^{*}v_{i}^{*}}_{\text{convection}} = q^{B^{*}} + \underbrace{\frac{1}{\mathrm{Pe}}\theta_{,ii}^{*}}_{\text{diffusion}}$ (7.68)
|
|
|
|
where Re is the Reynolds number
|
|
|
|
$$
|
|
\mathrm{Re} = \frac {v L}{\nu}; \qquad \nu = \frac {\mu}{\rho} \tag {7.69}
|
|
$$
|
|
|
|
and Pe is the Péclet number
|
|
|
|
$$
|
|
\mathrm{Pe} = \frac {v L}{\alpha}; \quad \alpha = \frac {k}{\rho c _ {p}} \tag {7.70}
|
|
$$
|
|
|
|
in which $\nu$ and $\alpha$ are the kinematic viscosity and diffusivity of the fluid. Note that using the Prandtl number $\operatorname{Pr} = \nu / \alpha$ , we have $\operatorname{Pe} = (\operatorname{Pr})(\operatorname{Re})$ .
|
|
|
|
The relations (7.67) and (7.68) demonstrate a fundamental difficulty in the solution of fluid flow problems: as the Reynolds number increases, the fluid flow is dominated by the convection term in (7.67), and similarly, as the Péclet number increases, the heat transfer in the fluid is dominated by the convection term in (7.68). General analysis procedures must therefore be able to solve for the response governed primarily by diffusion at low Reynolds and Péclet numbers and primarily by convection at high Reynolds and Péclet numbers. We will refer to this observation again in Section 7.4.3.
|
|
|
|
# 7.4.2 Finite Element Governing Equations
|
|
|
|
The finite element solution of the continuum mechanics equations governing the fluid flow is obtained by establishing a weak form of the equations using the Galerkin procedure (see Section 3.3.4 and Examples 4.2 and 7.1). The momentum equations are weighted with the velocities, the continuity equation is weighted with pressure, and the heat transfer equation is weighted with temperature. Integrating over the domain of interest V and using the divergence theorem—to lower the order of the derivatives in the expressions and to incorporate the natural boundary conditions as forcing terms—gives the variational equations to be discretized by finite element interpolations:
|
|
|
|
momentum:
|
|
|
|
$$
|
|
\int_ {V} \overline {{v}} _ {i} \rho (\dot {v} _ {i} + v _ {i, j} v _ {j}) d V + \int_ {V} \overline {{e}} _ {i j} \tau_ {i j} d V = \int_ {V} \overline {{v}} _ {i} f _ {i} ^ {B} d V + \int_ {S _ {f}} \overline {{v}} _ {i} ^ {S} f _ {i} ^ {S} d S \tag {7.71}
|
|
$$
|
|
|
|
continuity:
|
|
|
|
$$
|
|
\int_ {V} \overline {{p}} v _ {i, i} d V = 0 \tag {7.72}
|
|
$$
|
|
|
|
<!-- source-page: 695 -->
|
|
|
|
heat transfer:
|
|
|
|
$$
|
|
\int_ {V} \overline {{\theta}} \rho c _ {p} (\dot {\theta} + \theta_ {, i} v _ {i}) d V + \int_ {V} k \overline {{\theta}} _ {, i} \theta_ {, i} d V = \int_ {V} \overline {{\theta}} q ^ {B} d V + \int_ {S _ {q}} \overline {{\theta}} ^ {S} q ^ {S} d S \tag {7.73}
|
|
$$
|
|
|
|
where the overbar denotes a virtual quantity.
|
|
|
|
We may refer to the relation in (7.71) as the principle of virtual velocities, and of course (7.73) represents the principle of virtual temperatures [see (7.7)]. The relations in (7.71) to (7.73) are very similar to the equations we used in Sections 4.2 and 7.2 for the stress and temperature analysis of solids, but the Eulerian formulations in (7.71) and (7.73) include the convective terms. Also, incompressible conditions are considered, and therefore, the finite element formulations and discussions in Sections 4.4.3 and 4.5 are now important. Indeed, we can now directly apply the mixed formulations in Section 4.4.3 but using velocity (instead of displacement) and pressure as variables.
|
|
|
|
Assume that the finite element discretization of $(7.71)$ to $(7.73)$ is performed with any one of the elements presented in Section 4.4.3, with velocity and pressure as variables and temperature as an additional variable at all velocity nodes (in practice, we would of course use the isoparametric generalizations; see Section 5.3.5). The governing matrix equations for a single element are then
|
|
|
|
$$
|
|
\mathbf {M} _ {v} \dot {\hat {\mathbf {v}}} + \left(\mathbf {K} _ {\mu v v} + \mathbf {K} _ {v v}\right) \hat {\mathbf {v}} + \mathbf {K} _ {v p} \hat {\mathbf {p}} = \mathbf {R} _ {B} + \mathbf {R} _ {S} \tag {7.74}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {v p} ^ {T} \hat {\mathbf {v}} = \mathbf {0} \tag {7.75}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {C} \dot {\boldsymbol {\theta}} + (\mathbf {K} _ {v \theta} + \mathbf {K} _ {\theta \theta}) \hat {\boldsymbol {\theta}} = \mathbf {Q} _ {B} + \mathbf {Q} _ {S} \tag {7.76}
|
|
$$
|
|
|
|
where $\hat{v}$ , $\hat{\theta}$ , and $\hat{p}$ are, respectively, the unknown nodal point velocities and temperatures, and nodal point or element internal pressure variables.
|
|
|
|
For example, in two-dimensional planar flow analysis, we have for an element in the $x_{2}, x_{3}$ plane,
|
|
|
|
$$
|
|
\begin{array}{l} \left[ \begin{array}{c c c c} \mathbf {M} _ {v _ {2}} & & & \\ & \mathbf {M} _ {v _ {3}} & & \text {Zeros} \\ & & \mathbf {0} & \\ & \text {Sym.} & & \mathbf {C} \end{array} \right] \left[ \begin{array}{c} \dot {\hat {\mathbf {v}}} _ {2} \\ \dot {\hat {\mathbf {v}}} _ {3} \\ \dot {\hat {\mathbf {p}}} \\ \dot {\hat {\boldsymbol {\theta}}} \end{array} \right] + \left[ \begin{array}{c c c c} \mathbf {K} _ {\mu v _ {2} v _ {2}} + \mathbf {K} _ {v v _ {2}} & \mathbf {K} _ {\mu v _ {2} v _ {3}} & \mathbf {K} _ {v _ {2} p} & \mathbf {0} \\ \mathbf {K} _ {\mu v _ {2} v _ {3}} ^ {T} & \mathbf {K} _ {\mu v _ {3} v _ {3}} + \mathbf {K} _ {v v _ {3}} & \mathbf {K} _ {v _ {3} p} & \mathbf {0} \\ \mathbf {K} _ {v _ {2} p} ^ {T} & \mathbf {K} _ {v _ {3} p} ^ {T} & \mathbf {0} & \mathbf {0} \\ \mathbf {0} & \mathbf {0} & \mathbf {0} & \mathbf {K} _ {v \theta} + \mathbf {K} _ {\theta \theta} \end{array} \right] \left[ \begin{array}{c} \hat {\mathbf {v}} _ {2} \\ \hat {\mathbf {v}} _ {3} \\ \hat {\mathbf {p}} \\ \hat {\boldsymbol {\theta}} \end{array} \right] \\ = \left[ \begin{array}{c} \mathbf {R} _ {B _ {2}} + \mathbf {R} _ {S _ {2}} \\ \mathbf {R} _ {B _ {3}} + \mathbf {R} _ {S _ {3}} \\ \mathbf {0} \\ \mathbf {Q} _ {B} + \mathbf {Q} _ {S} \end{array} \right] \tag {7.77} \\ \end{array}
|
|
$$
|
|
|
|
Let H and $\tilde{H}$ contain the interpolation functions corresponding to the velocities and temperature, and pressure, respectively; then
|
|
|
|
$$
|
|
\mathbf {M} _ {v _ {2}} = \mathbf {M} _ {v _ {3}} = \rho \int_ {V} \mathbf {H} ^ {T} \mathbf {H} d V \tag {7.78}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {\mu v _ {2} v _ {2}} = \int_ {V} \left(2 \mu \mathbf {H} _ {, x _ {2}} ^ {T} \mathbf {H} _ {, x _ {2}} + \mu \mathbf {H} _ {, x _ {3}} ^ {T} \mathbf {H} _ {, x _ {3}}\right) d V
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {\mu \omega_ {2} v _ {3}} = \int_ {V} \left(\mu \mathbf {H} _ {, x _ {3}} ^ {T} \mathbf {H} _ {, x _ {2}}\right) d V \tag {7.79}
|
|
$$
|
|
|
|
<!-- source-page: 696 -->
|
|
|
|
$$
|
|
\mathbf {K} _ {\mu v _ {3} v _ {3}} = \int_ {V} \left(2 \mu \mathbf {H} _ {, x _ {3}} ^ {T} \mathbf {H} _ {, x _ {3}} + \mu \mathbf {H} _ {, x _ {2}} ^ {T} \mathbf {H} _ {, x _ {2}}\right) d V
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {v v _ {2}} = \mathbf {K} _ {v v _ {3}} = \rho \int_ {V} \left(\mathbf {H} ^ {T} \mathbf {H} \hat {\mathbf {v}} _ {2} \mathbf {H} _ {, x _ {2}} + \mathbf {H} ^ {T} \mathbf {H} \hat {\mathbf {v}} _ {3} \mathbf {H} _ {, x _ {3}}\right) d V \tag {7.80}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {v _ {2} p} = - \int_ {V} \mathbf {H} _ {, x _ {2}} ^ {T} \tilde {\mathbf {H}} d V \tag {7.81}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {v _ {3} p} = - \int_ {V} \mathbf {H} _ {, x _ {3}} ^ {T} \tilde {\mathbf {H}} d V \tag {7.82}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {R} _ {B} = \int_ {V} \mathbf {H} ^ {T} \mathbf {f} ^ {B} d V \tag {7.83}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {R} _ {s} = \int_ {s _ {f}} \mathbf {H} ^ {s T} \mathbf {f} ^ {s} d S \tag {7.84}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {C} = \rho \int_ {V} c _ {p} \mathbf {H} ^ {T} \mathbf {H} d V \tag {7.85}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {\theta \theta} = \int_ {V} k \left(\mathbf {H} _ {, x _ {2}} ^ {T} \mathbf {H} _ {, x _ {2}} + \mathbf {H} _ {, x _ {3}} ^ {T} \mathbf {H} _ {, x _ {3}}\right) d V \tag {7.86}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {K} _ {v \theta} = \rho \int_ {V} c _ {p} \mathbf {H} ^ {T} \mathbf {H} \hat {\mathbf {v}} _ {2} \mathbf {H} _ {, x _ {2}} d V + \rho \int_ {V} c _ {p} \mathbf {H} ^ {T} \mathbf {H} \hat {\mathbf {v}} _ {3} \mathbf {H} _ {, x _ {3}} d V \tag {7.87}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {Q} _ {B} = \int_ {V} \mathbf {H} ^ {T} q ^ {B} d V \tag {7.88}
|
|
$$
|
|
|
|
$$
|
|
\mathbf {Q} _ {s} = \int_ {s _ {q}} \mathbf {H} ^ {s T} q ^ {s} d S \tag {7.89}
|
|
$$
|
|
|
|
Here we should note that in considering a straight boundary the components of $f^{s}$ are
|
|
|
|
$$
|
|
f _ {n} = - p + 2 \mu \frac {\partial v _ {n}}{\partial n} \tag {7.90}
|
|
$$
|
|
|
|
$$
|
|
f _ {t} = \mu \left(\frac {\partial v _ {t}}{\partial n} + \frac {\partial v _ {n}}{\partial t}\right) \tag {7.91}
|
|
$$
|
|
|
|
where n and t denote the coordinate axes in the normal and tangential directions on the boundary and $v_{n}$ and $v_{t}$ are the normal and tangential boundary velocities.
|
|
|
|
We note that since totally incompressible conditions are considered, the diagonal elements corresponding to the pressure variables are zero. Hence, even if the pressure variables stored in $\hat{p}$ correspond to element internal variables and not nodal point variables (as in the u/p formulation; see Section 4.4.3), the pressure variables cannot be statically condensed out at the element level. For such a procedure to be used, we must consider an almost incompressible condition, meaning that (7.72) would need to be replaced by (see Sections 4.4.3 and 4.5)
|
|
|
|
$$
|
|
\int_ {V} \overline {{p}} \left(\frac {\dot {p}}{\kappa} + v _ {i, i}\right) d V = 0 \tag {7.92}
|
|
$$
|
|
|
|
in which $\kappa$ (being very large) is the bulk modulus.
|
|
|
|
<!-- source-page: 697 -->
|
|
|
|
We can now refer to all the procedures discussed in Chapters 4 (notably Sections 4.4 and 4.5) and 5. They are directly applicable to obtain appropriate finite element discretizations of the governing fluid flow equations (see also the exercises at the end of this section). Of course, we note that the resulting finite element equations are in general highly nonlinear equations [see (7.77)] because of the convective terms and boundary radiation conditions and because in general the material properties are not constant (e.g., the viscosity $\mu$ depends significantly on the temperature $\theta$ ).
|
|
|
|
The finite element governing fluid flow equations—given for a two-dimensional element in (7.77)—can be written in the form
|
|
|
|
$$
|
|
\mathbf {R} - \mathbf {F} = \mathbf {0} \tag {7.93}
|
|
$$
|
|
|
|
This relation must hold for all times considered, and the solution can be obtained incrementally for times $\Delta t$ , $2\Delta t$ , $\ldots$ , as discussed in Chapters 6 and 9. In steady-state analysis, the terms $M\dot{v}$ and $C\dot{\theta}$ are neglected, and the incremental analysis can be pursued using a form of Newton-Raphson iteration for each time $t + \Delta t$ considered, time then denoting merely the load levels (see Chapter 6). In transient analysis, using implicit integration, a form of Newton-Raphson iteration should also be used for each time step. Alternatively, explicit integration can be performed on the velocity and temperature variables, while the incompressibility constraint demands implicit integration of the pressure equations.
|
|
|
|
A major difficulty in the solution of fluid flow problems is that the number of solution variables is generally very large (to obtain a realistic resolution of the fluid response, very fine finite element discretizations need to be employed) and that the coefficient matrices are nonsymmetric. Hence, explicit schemes that do not require the solution of a set of equations (see Sections 9.5.1 and 9.6.1) and iterative methods for the solution of equations are very attractive (see Section 8.3).
|
|
|
|
Finally, we may ask why the relations (7.57) and (7.60) were used to develop the finite element equations instead of (7.67) and (7.68). One reason is that (7.57) and (7.60) are directly applicable to flow with nonconstant material conditions. Another reason is that the momentum equation (7.57) when used in the Galerkin method leads to a surface force vector that contains the physical tractions given in (7.90) and (7.91), whereas if (7.67) is used, the resulting surface “force” vector contains nonphysical components (see Example 7.11). Therefore, the formulation based on (7.57) is frequently more general and natural, in particular when arbitrary Lagrangian-Eulerian formulations are developed that are used for the analysis of free surface flows and interactions between fluid flows and structures. Of course, the finite element matrix equations (7.74) to (7.76) are applicable to the solution of fluid flow and heat transfer problems in any set of consistent units, including the use of nondimensionalized variables.
|
|
|
|
In the following examples we consider two additional forms of the Navier-Stokes equations that can be used for finite element solutions.
|
|
|
|
EXAMPLE 7.11: Consider the solution of the Navier-Stokes equations [see (7.67)],
|
|
|
|
$$
|
|
v _ {i, j} v _ {j} = - p _ {, i} + \frac {1}{\mathrm{Re}} v _ {i, j j} \quad i, j = 1, 2
|
|
$$
|
|
|
|
Use the Galerkin procedure by weighting this equation with velocity and derive the naturally arising boundary terms. Consider two-dimensional analysis and compare these terms with the expressions of the physical tractions.
|
|
|
|
<!-- source-page: 698 -->
|
|
|
|
Using the Galerkin procedure, we obtain
|
|
|
|
$$
|
|
\int_ {V} \overline {{v}} _ {i} \left(v _ {i, j} v _ {j} + p _ {, i} - \frac {1}{\mathrm{Re}} v _ {i, j j}\right) d V = 0 \tag {a}
|
|
$$
|
|
|
|
The boundary terms are established as in Examples 4.2 and 7.1. Here we use the identities
|
|
|
|
$$
|
|
\overline {{{v}}} _ {i} p _ {, j} \delta_ {i j} = (\overline {{{v}}} _ {i} p \delta_ {i j}) _ {, j} - \overline {{{v}}} _ {i, j} p \delta_ {i j}
|
|
$$
|
|
|
|
and $\overline{v}_i v_{i,jj} = (\overline{v}_i v_{i,j}), j - \overline{v}_{i,j} v_{i,j}$
|
|
|
|
Hence, using the divergence theorem, we obtain from (a) the boundary term
|
|
|
|
$$
|
|
\int_ {S} \overline {{{v}}} _ {i} (- p \delta_ {i j} + \frac {1}{\mathrm{Re}} v _ {i, j}) n _ {j} d S
|
|
$$
|
|
|
|
However, this relation can be interpreted as
|
|
|
|
$$
|
|
\int_ {S} \overline {{v}} _ {i} \tilde {f} _ {i} d S
|
|
$$
|
|
|
|
where $\tilde{f}_i = (-p\delta_{ij} + \frac{1}{\mathbf{Re}} v_{i,j})n_j$
|
|
|
|
Considering a straight boundary, we have
|
|
|
|
$$
|
|
\tilde {f} _ {n} = - p + \frac {1}{\mathrm{Re}} \frac {\partial v _ {n}}{\partial n}
|
|
$$
|
|
|
|
$$
|
|
\tilde {f} _ {t} = \frac {1}{\mathrm{Re}} \frac {\partial v _ {t}}{\partial n} \tag {b}
|
|
$$
|
|
|
|
where $n$ and $t$ denote the coordinate axes in the directions normal and tangential to the surface boundary.
|
|
|
|
The relation (7.67) was obtained using nondimensionalized variables, and for given characteristic velocity and length we have $1/Re \sim \nu$ . Hence, we see that the expressions in (b) are not equal to the actual force expressions in (7.90) and (7.91).
|
|
|
|
EXAMPLE 7.12: Show that the momentum equations (7.57) can also be written as
|
|
|
|
$$
|
|
\rho \frac {\partial v _ {i}}{\partial t} + F _ {i j, j} = f _ {i} ^ {B} \tag {a}
|
|
$$
|
|
|
|
where $F_{ij} = \rho v_{j} v_{i} - \tau_{ij}$ (b)
|
|
|
|
Then identify the difference that will arise when (a) is used instead of (7.57) in the Galerkin procedure to obtain the governing finite element equations.
|
|
|
|
Using (b), we obtain
|
|
|
|
$$
|
|
F _ {i j, j} = \left(\rho v _ {j} v _ {i} - \tau_ {i j}\right) _ {, j} = \rho v _ {j, j} v _ {i} + \rho v _ {j} v _ {i, j} - \tau_ {i j, j} = \rho v _ {i, j} v _ {j} - \tau_ {i j, j} \tag {c}
|
|
$$
|
|
|
|
Hence, by inspection we see that (c) used in (a) gives (7.57).
|
|
|
|
The form (a) is referred to as the conservative form of the momentum equations because, by the divergence theorem, for any subdomain $V_{SD}$ of the fluid
|
|
|
|
$$
|
|
\int_ {V _ {S D}} F _ {i j, j} d V = \int_ {s _ {S D}} F _ {i j} n _ {j} d S
|
|
$$
|
|
|
|
where $S_{SD}$ is the surface area of the volume $V_{SD}$ and the $n_j$ are the components of the unit normal vector to $S_{SD}$ . Similarly, the energy equation can also be written in conservative form.
|
|
|
|
<!-- source-page: 699 -->
|
|
|
|
To identify the difference in the finite element discretization, we only need to compare the terms $\int_{V}\overline{v}_{i}(v_{j}v_{i})_{,j}dV$ and $\int_{V}\overline{v}_{i}(v_{i,j}v_{j})dV$ because the other terms are identical, and the difference is the term $\int_{V}\overline{v}_{i}(v_{i}v_{j,j})dV$ . The form of the momentum equations in (a) is typically used in finite volume methods (see S. V. Patankar [A]). Note, however, that if (a) is used in a finite element (Galerkin) formulation, the use of the divergence theorem gives on the surface $S_{f}$ the usual traction term (see Example 4.2) and an additional term involving the unknown velocities.
|
|
|
|
# 7.4.3 High Reynolds and High Péclet Number Flows
|
|
|
|
The finite element formulation of fluid flows presented in the previous section is a natural development when we consider the formulations presented earlier for the analysis of solids. The standard Galerkin procedure was used on the differential equations of motion and heat transfer, resulting in the “principle of virtual velocities” and the “principle of virtual temperatures.” The finite element discretization is appropriately obtained with finite elements that are stable and convergent with the incompressibility constraint. Hence, we used elements that satisfy the inf-sup condition, as discussed in Sections 4.4.3 and 4.5. With such elements, excellent results are obtained for low-Reynolds-number flows (in particular Stokes flows).
|
|
|
|
However, as we pointed out in Section 7.4.2, a major difference between the formulations for the analysis of fluid flows and of solids is the convective terms arising in the Eulerian fluid flow formulation. These convective terms give rise to the nonsymmetry in the finite element coefficient matrix, and when the convection is strong (as defined by the Reynolds and Péclet numbers; see discussion below), the system of equations is strongly nonsymmetric and then an additional numerical difficulty arises.
|
|
|
|
However, before we discuss this difficulty, let us recall that of course depending on the flow considered, as the Reynolds number increases and when a certain range is reached, the flow condition turns from laminar to turbulent. In theory, the turbulent flow could still be calculated by solving the general Navier-Stokes equations presented in the previous section, but such solution would require extremely fine discretization to model the details of turbulence. For practical flow conditions, the resulting finite element systems would be too large, by far, for present software and hardware capabilities. For this reason, it is common practice to solve the Navier-Stokes equations for the mean flow and express the turbulence effects by means of turbulent viscosity and heat conductivity coefficients and use wall functions to describe the near-wall behaviors.
|
|
|
|
The modeling of turbulence is a very large and important field (see W. Rodi [A]), and the finite element procedures we presented earlier for solution are in many regards directly applicable. However, one important factor in the finite element scheme must then also be that the Navier-Stokes equations corresponding to laminar flow at high Reynolds and Péclet numbers can be solved (with reasonably sized meshes). Namely, it is this solution that provides the basis for obtaining the solution of the turbulent flow.
|
|
|
|
Hence, in the following we briefly address the difficulty of solving high Reynolds (or Péclet) number conditions assuming laminar flow. For this purpose let us consider the simplest possible case that displays the difficulties that we encounter in general flow conditions. These difficulties arise from the magnitude of the convective terms when compared to the diffusive terms in (7.67) and (7.68). Hence, we consider a model problem of one-dimensional flow with prescribed velocity v (see Fig. 7.4). The temperature is prescribed
|
|
|
|
<!-- source-page: 700 -->
|
|
|
|
|
|
|
|
<details>
|
|
<summary>text_image</summary>
|
|
|
|
Unit cross-
|
|
sectional area
|
|
x
|
|
v
|
|
Fluid
|
|
L
|
|
(a) Problem considered
|
|
</details>
|
|
|
|
|
|
|
|
<details>
|
|
<summary>line</summary>
|
|
|
|
| x/L | Pe = 0 | Pe = 1 | Pe = 2 | Pe = 10 | Pe = 50 |
|
|
| --- | --- | --- | --- | --- | --- |
|
|
| 0.0 | 0.0 | 0.0 | 0.0 | 0.0 | 0.0 |
|
|
| 0.2 | 0.1 | 0.05 | 0.08 | 0.03 | 0.01 |
|
|
| 0.4 | 0.3 | 0.15 | 0.2 | 0.08 | 0.02 |
|
|
| 0.6 | 0.5 | 0.3 | 0.35 | 0.15 | 0.03 |
|
|
| 0.8 | 0.7 | 0.5 | 0.55 | 0.3 | 0.05 |
|
|
| 1.0 | 1.0 | 0.9 | 0.95 | 0.9 | 0.6 |
|
|
</details>
|
|
|
|
(b) Analytical solutions
|
|
|
|
|
|
|
|
<details>
|
|
<summary>text_image</summary>
|
|
|
|
-|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|
|
θi-1 θi θi+1
|
|
i-1 i i+1
|
|
</details>
|
|
|
|
(c) Finite element nodes (and finite difference stations)
|
|
Figure 7.4 Heat transfer in one-dimensional flow condition; prescribed
|
|
velocity $v; q^{B} = 0; \mathrm{Pe} = \frac{vL}{\alpha}, \alpha = k / \rho c_{p}$
|
|
|
|
at two points, which we label $x = 0$ and $x = L$ , and we want to calculate the temperature for $0 < x < L$ .
|
|
|
|
The governing differential equation obtained from (7.60) is
|
|
|
|
$$
|
|
\rho c _ {p} \frac {d \theta}{d x} v = k \frac {d ^ {2} \theta}{d x ^ {2}} \tag {7.94}
|
|
$$
|
|
|
|
with the boundary conditions
|
|
|
|
$$
|
|
\theta = \theta_ {L} \quad \text { at } x = 0 \tag {7.95}
|
|
$$
|
|
|
|
$$
|
|
\theta = \theta_ {R} \quad \text { at } x = L
|
|
$$
|
|
|
|
and the left-hand side in (7.94) represents the convective terms and the right-hand side the diffusive terms.
|
|
|
|
Of course, the convective and diffusive terms appear in a similar form in the Navier-Stokes equations [see (7.67)], which can also be considered in a one-dimensional analogously simple form. However, the resulting differential equation is nonlinear in v, whereas (7.94) is linear in $\theta$ . Since the solution of (7.94) already displays the basic solution difficulties, we prefer to consider (7.94) but recognize that the basic observations are also applicable to the solution of the Navier-Stokes equations.
|