김경종
24 KiB
arbitrary point because the solution of (7.44) and (7.45) can be determined only after one value of \phi is fixed.
The solution of (7.44) with the boundary conditions is analogous to the solution of a heat transfer problem. Once the potential function \phi has been evaluated, we can use Bernoulli's equation to calculate the pressure distribution in the fluid. Figure 7.3 illustrates a finite element analysis of the flow around an object in a channel.
text_image
A' ∂φ/n = ∂φ/∂y = 0 L1/L2 << 1 n B' vp ∂φ/n = -vp n φ100 = 1.0 ∂φ/n = 0 at boundary L1 vp ∂φ/n = vp n y x A n B ∂φ/n = -∂φ/∂y = 0 L2
Figure 7.3 Analysis of flow in a channel with an island, v_{p} is the prescribed velocity. Note that on the boundary A-A' the inflow condition requires that a negative gradient of \phi be imposed, whereas on the boundary B-B' a positive gradient is prescribed. (However, using a heat transfer program, the boundary conditions would be q^{s} = v_{p} on A-A' and q^{s} = -v_{p} on B-B' for a flow to the right because the flow is calculated as proportional to minus the gradient of the potential.)
7.3.3 Torsion
With the introduction of a stress function \phi , the elastic torsional behavior of a shaft is governed by the equation (see, for example, Y. C. Fung [A])
\frac {\partial^ {2} \phi}{\partial x ^ {2}} + \frac {\partial^ {2} \phi}{\partial y ^ {2}} + 2 G \theta = 0 \tag {7.46}
where \theta is the angle of twist per unit length and G is the shear modulus of the shaft material. The shear stress components at any point can be calculated using
\tau_ {z x} = \frac {\partial \phi}{\partial y}; \quad \tau_ {z y} = - \frac {\partial \phi}{\partial x} \tag {7.47}
and the applied torque is given by
T = 2 \int_ {A} \phi d A \tag {7.48}
where A is the cross-sectional area of the shaft. The boundary condition on \phi is that \phi is zero on the boundary of the shaft. Hence, the heat transfer equations in (7.1) and (7.2) also govern the torsional behavior of a shaft provided the appropriate field variables are used.
EXAMPLE 7.7: Evaluate the torsional rigidity of a square shaft using the two different finite element meshes shown in Fig. E7.7.
text_image
2 cm 2 cm s r 2 cm 2 cm 2 cm y 1 2 3 4 x
Mesh (a)
text_image
2 cm 2 cm 2 cm 2 cm 2 cm y 3 3 1 2 1 5 4 x
Mesh (b)
Figure E7.7 Finite element meshes used in calculation of torsional rigidity of a square shaft
Using the analogy between this torsional problem and a heat transfer problem for which the governing finite element equations have been stated in Section 7.2.3, we now want to solve
\left[ \sum_ {m} \int_ {V ^ {(m)}} \mathbf {B} ^ {(m) ^ {T}} \mathbf {k} ^ {* (m)} \mathbf {B} ^ {(m)} d V ^ {(m)} \right] \boldsymbol {\phi} = \sum_ {m} \int_ {V ^ {(m)}} \mathbf {H} ^ {(m) ^ {T}} \theta d V ^ {(m)} \tag {a}
where
\phi^ {(m)} = \mathbf {H} ^ {(m)} \phi
\boldsymbol {\phi} ^ {\prime (m)} = \mathbf {B} ^ {(m)} \boldsymbol {\phi}
\boldsymbol {\phi} ^ {T} = \left[ \begin{array}{l l l l} \phi_ {1} & \phi_ {2} & \dots & \phi_ {n} \end{array} \right]
\mathbf {k} ^ {* (m)} = \left[ \begin{array}{c c} \frac {1}{2 G} & 0 \\ 0 & \frac {1}{2 G} \end{array} \right]
\begin{array}{l} \boldsymbol {\phi} ^ {\prime (m)} = \left[ \begin{array}{c} - \tau_ {z y} \\ \tau_ {z x} \end{array} \right] \\ = \mathbf {B} ^ {(m)} \phi \\ \end{array}
T = \sum_ {m} \int_ {A ^ {(m)}} 2 \phi^ {(m)} d A ^ {(m)} \tag {b}
In each of the analysis cases we need to consider only one element because of symmetry conditions. Using mesh (a), we have for element 1,
\phi = \frac {1}{4} (1 + r) (1 + s) \phi_ {1}
and
\left[ \begin{array}{l} \frac {\partial \phi}{\partial r} \\ \frac {\partial \phi}{\partial s} \end{array} \right] = \frac {1}{4} \left[ \begin{array}{l} (1 + s) \\ (1 + r) \end{array} \right] \phi_ {1}
Hence the equations in (a) reduce to (considering a unit length of shaft)
\begin{array}{l} \left\{4 \int_ {- 1} ^ {+ 1} \int_ {- 1} ^ {+ 1} \frac {1}{4} [ (1 + s) (1 + r) ] \frac {1}{2 G} \frac {1}{4} \left[ \begin{array}{l} (1 + s) \\ (1 + r) \end{array} \right] \det \mathbf {J} d r d s \right\} \phi_ {1} \\ = 4 \theta \int_ {- 1} ^ {+ 1} \int_ {- 1} ^ {+ 1} \frac {1}{4} (1 + r) (1 + s) \det \mathbf {J} d r d s; \quad \det \mathbf {J} = 1 \\ \end{array}
or \frac{1}{3G}\phi_{1} = \theta
Hence, \phi_{1} = 3G\theta
Using (b), we thus obtain
T = 4 \int_ {- 1} ^ {+ 1} \int_ {- 1} ^ {+ 1} \frac {1}{2} (1 + r) (1 + s) (3 G \theta) d r d s = 2 4 G \theta
so that \frac{T}{\theta} = 24G
Considering next mesh (b), we recognize that the \phi values on the boundary are zero and that for element 1 we have \phi_{4} = \phi_{5} . Hence, we have only two unknowns \phi_{1} and \phi_{4} to be calculated. The interpolation functions for the eight-node element are given in Fig. 5.5. Proceeding in the same way as in the analysis with mesh (a), we obtain
\phi_ {1} = 2. 1 5 7 G \theta
\phi_ {4} = 1. 9 2 1 G \theta
T = 3 5. 2 G \theta
so that \frac{T}{\theta} = 35.2G
The exact solution of (7.46) for T/\theta is 36.1G. Hence, the analysis with mesh (a) gives an error of 33.5 percent, whereas the analysis with mesh (b) yields a result of only 2.5 percent error. ^{1}
7.3.4 Acoustic Fluid
Consider an inviscid isentropic fluid with the fluid particles undergoing only small displacements. Not including body force effects, the equations governing the response of the fluid are the momentum equations (see, for example, F. M. White [A]),
\rho \dot {\mathbf {v}} + \nabla p = \mathbf {0} \tag {7.49}
and the constitutive equation
\boldsymbol {\beta} \boldsymbol {\nabla} \cdot \mathbf {v} + \dot {p} = 0 \tag {7.50}
where v is the velocity of the fluid particles, p is the pressure, and \beta is the bulk modulus. The boundary conditions are as follows.
On the boundary S_{v} a prescribed velocity v_{n}^{S} in the direction of the unit normal vector n (pointing outward) to the fluid boundary,
\mathbf {v} \cdot \mathbf {n} \mid_ {S _ {0}} = v _ {n} ^ {s} \tag {7.51}
On the boundary S_f a prescribed pressure p^s ,
p \mid_ {s _ {f}} = p ^ {s} \tag {7.52}
To solve for the fluid motion, it is convenient to introduce the velocity potential \phi , where
\mathbf {v} = \nabla \phi ; \quad p = - \rho \dot {\phi} \tag {7.53}
With this definition (7.49) is identically satisfied [neglecting changes in density in (7.49)], and (7.50) becomes the acoustic equation
\nabla^ {2} \phi = \frac {1}{c ^ {2}} \ddot {\phi} \tag {7.54}
with the wave speed c = \sqrt{\beta/\rho} . The boundary conditions are now of course
\left. \frac {\partial \phi}{\partial n} \right| _ {s _ {n}} = v _ {n} ^ {s} \tag {7.55}
and -\rho \dot{\phi}|_{s_f} = p^s (7.56)
On comparing the governing equations (7.54) to (7.56) with the heat transfer governing equations, we recognize that a strong analogy exists. However, in the analysis of the fluid, the second time derivative of the solution variable is taken instead of the first time derivative in heat transfer analysis. Hence, for example, when using a heat transfer program to calculate the frequencies of an acoustic fluid, the frequencies sought are obtained by taking the square roots of the calculated frequencies. We demonstrate this observation in the following example.
EXAMPLE 7.8: Consider an acoustic fluid in a closed rigid cavity (see Fig. E7.8). Use a 2 \times 2 mesh of eight-node elements to model the fluid and estimate the lowest frequency of vibration of the fluid.
We use the governing fluid equation (7.54) with the boundary condition (7.55). Hence, by analogy to the development of the principle of virtual temperatures (see Example 7.1), the appropriate variational equation is
\int_ {V} \delta \phi \frac {1}{c ^ {2}} \ddot {\phi} d V + \int_ {V} (\boldsymbol {\nabla} \delta \phi) \cdot (\boldsymbol {\nabla} \phi) d V = 0 \tag {a}
Substituting the finite element interpolations into (a), we obtain
\mathbf {M} \ddot {\boldsymbol {\phi}} + \mathbf {K} \boldsymbol {\phi} = \mathbf {0} \tag {b}
text_image
12 in 20 in Fluid of bulk modulus β, density ρ Unit depth n n Acoustic fluid, β = 3.16 × 10⁵ psi, ρ = 9.35 × 10⁻⁵ lb·sec²/in⁴ State variable φ Boundary condition ∂φ/∂n = 0
Figure E7.8 Problem and finite element discretization of fluid in rigid cavity
where the Mφ term corresponds to the strain energy and the Kφ term to the kinetic energy. The corresponding eigenvalue problem is
\mathbf {K} \boldsymbol {\phi} = \omega^ {2} \mathbf {M} \boldsymbol {\phi} \tag {c}
The matrices K and M are calculated as in heat transfer analysis [the heat conduction matrix K (with unit material conductivities in all directions) corresponds to the matrix K in (b), and the heat capacity matrix C (with pc|_{heat} replaced by 1/c^{2}|_{acoustic} corresponds to the matrix M in (b)].
The problem (c) can therefore be solved directly with a heat transfer program that calculates the eigenvalues of the problem K\theta = \lambda C\theta (see Section 9.6). For a 2 \times 2 finite element discretization of eight-node elements, the results are \lambda_{1} = 0 , \lambda_{2} = (9166)^{2} , \lambda_{3} = (15,277)^{2} , and hence, the lowest nonzero calculated frequency of the problem in (c) is \omega_{2} = 9166 . We should note that the zero frequency corresponds to \phi = constant in the cavity. The analytical solution to the problem gives \omega_{2} = 9132 .
Equations (7.54) to (7.56) have been used to develop an efficient finite element formulation to model the interaction between acoustic fluids and structures (see G. C. Everstine [A], L. G. Olson and K. J. Bathe [A], and Example 7.9). Also, the analysis procedure has been generalized to consider the fluid to undergo very large particle motions (see C. Nitikitpaiboon and K. J. Bathe [B]).
EXAMPLE 7.9: Consider the problem shown in Fig. E7.9 and the model indicated. Evaluate the matrices for the fluid-structure interaction problem and calculate the lowest frequency of vibration.
The analysis of the problem requires the coupling of the fluid response with the response of the spring.
The principle of virtual work for the piston/spring gives
\overline {{{u}}} m \ddot {u} + \overline {{{u}}} k u = \overline {{{u}}} f ^ {F} + \overline {{{u}}} R (t) \tag {a}
where f^{F} corresponds to the force exerted by the fluid on the piston/spring. The “principle of
Unit cross-sectional area
Bulk modulus \beta = 2.1 \times 10^{9} Pa
Mass density \rho = 1000 \, kg/m^{3}
text_image
Rigid piston of mass m = 10³ k = 10⁷ R(t) Frictionless rollers 10 Spring element k = 10⁷ φ₁ φ₃ φ₂ 10 One 3-node one-dimensional element to represent fluid
Figure E7.9 Acoustic fluid in a cavity with piston
virtual potentials" of the fluid is
\int_ {v _ {f}} \overline {{{\phi}}} \frac {1}{c ^ {2}} \ddot {\phi} d V _ {f} + \int_ {v _ {f}} \nabla \overline {{{\phi}}} \cdot \nabla \phi d V _ {f} = \int_ {I} \overline {{{\phi}}} ^ {I} \dot {u} _ {n} d I \tag {b}
where I denotes the (wetted) fluid-structure interface and \dot{u}_{n} is the velocity of that interface (of the piston). We note that (b) is derived from (7.54) as we have shown (for the principle of virtual temperatures) in Example 7.1. Also, we have \partial\phi/\partial n|_{I} = \dot{u}_{n} , with n denoting the direction of the unit normal vector on the fluid domain (pointing outward to this domain).
We now represent the fluid domain by one three-node element. Using the developments of Chapter 5 (see specifically Section 5.3), we have corresponding to (b),
\mathbf {M} _ {F} \ddot {\boldsymbol {\phi}} + \mathbf {K} _ {F} \boldsymbol {\phi} = \mathbf {R} _ {\dot {u}}
with M_{F} , K_{F} , and R_{i} defined by
\frac {1}{c ^ {2}} \left[ \begin{array}{c c c} \frac {4}{3} & - \frac {1}{3} & \frac {2}{3} \\ & \frac {4}{3} & \frac {2}{3} \\ \text {Sym.} & & \frac {1 6}{3} \end{array} \right] \left[ \begin{array}{l} \ddot {\phi} _ {1} \\ \ddot {\phi} _ {2} \\ \ddot {\phi} _ {3} \end{array} \right] + \left[ \begin{array}{c c c} \frac {7}{3 0} & \frac {1}{3 0} & \frac {- 8}{3 0} \\ & \frac {7}{3 0} & \frac {- 8}{3 0} \\ \text {Sym.} & & \frac {1 6}{3 0} \end{array} \right] \left[ \begin{array}{l} \phi_ {1} \\ \phi_ {2} \\ \phi_ {3} \end{array} \right] = \left[ \begin{array}{l} 0 \\ \dot {u} _ {n} \\ 0 \end{array} \right] \tag {c}
Next we can couple (a) and (c) by noting that \dot{u}_{n} in (c) is equal to the time derivative of the displacement u in (a) and that f^{F} in (a) is given by the pressure in the fluid,
f ^ {F} = - \rho_ {F} \dot {\phi} | _ {I} = - \rho_ {F} \dot {\phi} _ {2}
since the cavity has unit cross-sectional area.
The coupled fluid-structure equations are
\begin{array}{l} \left[ \begin{array}{c:ccc} m & 0 & 0 & 0 \\ \hline 0 & & & \\ 0 & & - \rho_ {F} \mathbf {M} _ {F} & \\ 0 & & & \\ \end{array} \right]\left[ \begin{array}{c} \ddot {u} \\ \ddot {\phi} _ {1} \\ \ddot {\phi} _ {2} \\ \ddot {\phi} _ {3} \\ \end{array} \right] + \left[ \begin{array}{c c c c} 0 & 0 & \rho_ {F} & 0 \\ 0 & 0 & 0 & 0 \\ \rho_ {F} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array} \right]\left[ \begin{array}{c} \dot {u} \\ \dot {\phi} _ {1} \\ \dot {\phi} _ {2} \\ \dot {\phi} _ {3} \\ \end{array} \right] \\ + \left[ \begin{array}{c c c c} k & 0 & 0 & 0 \\ \hline 0 & & & \\ 0 & - \rho_ {F} \mathbf {K} _ {F} & & \\ 0 & & & \end{array} \right] \left[ \begin{array}{l} u \\ \phi_ {1} \\ \phi_ {2} \\ \phi_ {3} \end{array} \right] = \left[ \begin{array}{c} R (t) \\ 0 \\ 0 \\ 0 \end{array} \right] \tag {d} \\ \left[ \begin{array}{c:ccc} m & 0 & 0 & 0 \\ \hline 0 & & & \\ 0 & & - \rho_ {F} \mathbf {M} _ {F} & \\ 0 & & & \\ \end{array} \right]\left[ \begin{array}{c} \ddot {u} \\ \ddot {\phi} _ {1} \\ \ddot {\phi} _ {2} \\ \ddot {\phi} _ {3} \\ \end{array} \right] + \left[ \begin{array}{c c c c} 0 & 0 & \rho_ {F} & 0 \\ 0 & 0 & 0 & 0 \\ \rho_ {F} & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{array} \right]\left[ \begin{array}{c} \dot {u} \\ \dot {\phi} _ {1} \\ \dot {\phi} _ {2} \\ \dot {\phi} _ {3} \\ \end{array} \right] \\ \end{array}
TABLE E7.9 Frequencies of fluid-structure model in Fig. E7.9
| Coupled system | Piston without fluid | Fluid without piston, open cavity | Fluid with stationary piston, $k = \infty$ |
| $\omega_1 = 0$ | 0 | 0 | |
| $\omega_2 = 212$ | $\omega = \sqrt{k/m} = 100$ | 229 | 502 |
| $\omega_3 = 744$ | 822 | 1122 |
We note that to obtain symmetric coefficient matrices in (d) we have multiplied both sides of (c) by -\rho_{F} . Also, note that in (d) the coefficient matrix to the first time derivatives of the nodal variables is not a damping matrix but simply a matrix that couples the fluid and structural response. There is no physical damping in this problem.
The solution to the problem would be obtained by the time integration of the dynamic response using for example the trapezoidal rule (see Section 9.2.4). However, it is also interesting to evaluate the free-vibration frequencies of the fluid-structure system and compare them with the frequencies of the fluid and the structure when acting alone. Table E7.9 lists these frequencies. We notice that in the fluid-structure system, because of the fluid, the (lowest frequency) structural vibration occurs at a 112 percent higher frequency. This frequency can be expressed as \omega_{2} = \sqrt{(k + k')/(m + m')} , where k' and m' are an increase in stiffness and mass due to the fluid. Note also that the first frequency \omega_{1} = 0 corresponds to a “rigid body mode” with u = 0 and \phi_{1} = \phi_{2} = \phi_{3} = constant .
Further details of this formulation are given by L. G. Olson and K. J. Bathe [A] and C. Nitikitpaiboon and K. J. Bathe [B].
7.3.5 Exercises
7.9. Use a computer program to solve for the seepage flow in the problem shown. Assume a very long dam. Verify that an accurate solution of the response has been obtained.
text_image
Water level Impermeable dam 30 ft 10 ft 20 ft 60 ft Soil k = 1 ft/hr 2 ft Impermeable rock
7.10. Use a computer program to solve for the flow of a fluid around the circular object shown. Assume an inviscid fluid and planar flow conditions.
text_image
20 Circular object of radius 10 20 Uniform velocity v
7.11. Use a computer program to solve for the torsional rigidity of the shaft considered in Example 7.7.
7.12. Calculate the steady-state distribution of voltage in the specimen shown. The solution of this kind of problem is of interest in monitoring crack growth (see, for example, R. O. Ritchie and K. J. Bathe [A]).
text_image
Uniform current inflow = 1 5.2 in 45° 20 in Uniform current outflow = 1 50 in 50 in Electric conductivity k = 1.0 Line of symmetry
7.13. Use a computer program to solve for the three lowest frequencies of the fluid-spring system considered in Example 7.9. Use a coarse and a fine finite element discretization and compare your results with those given in Example 7.9.
7.14. Use a computer program to solve for the frequency of the cylinder shown oscillating in a cavity of water. (Hint: Here the \phi formulation in Section 7.3.4 is effective.)
text_image
Spring k = 1240.26 N/m Rigid cylinder Mass = 31.4 kg Diameter = 0.2 m Fluid (water) ρ = 1000 kg/m³ β = 2.1 × 10⁹ Pa Unit depth
7.4 ANALYSIS OF VISCOUS INCOMPRESSIBLE FLUID FLOWS
A frequently followed approach in the analysis of fluid mechanics problems is to develop the governing differential equations for the specific flow, geometry, and boundary conditions under consideration and then solve these equations using finite difference procedures.
However, during the last decade, much progress has been made toward the analysis of general fluid flow problems using finite element procedures, and at present very complex fluid flows are solved. Inviscid and viscous, compressible and incompressible fluid flows with or without heat transfer are analyzed and also the coupling between structural and fluid response is considered (see, for example, K. J. Bathe, H. Zhang, and M. H. Wang [A]).
The objective in this section is to discuss briefly how the finite element procedures presented earlier can also be employed in the analysis of fluid flow problems and to present some important additional techniques that are required if practical problems of fluid flows are to be solved. We consider the large application area of incompressible viscous fluid flow with or without heat transfer. The methods applicable in this field of analysis are also directly used for compressible flow solutions, although additional complex analysis phenomena (notably the calculation of shock fronts) need to be addressed when compressible fluid flow is analyzed.
In order to identify the similarities and differences in the basic finite element formulations of problems in solid (see Chapter 6) and fluid mechanics, consider the summary of the governing continuum mechanics equations given in Table 7.4. In this table and in the discussion to follow we use the indicial notation employed in Chapter 6. Considering the kinematics of a viscous fluid, the fluid particles can undergo very large motions and, for a general description, it is effective to use an Eulerian formulation. The essence of this formulation is that we focus attention on a stationary control volume and that we use this volume to measure the equilibrium and mass continuity of the fluid particles. This means that in the Eulerian formulation a separate equation is written to express the mass conservation relation—a condition that is embodied in the determinant of the deformation gradient when using a Lagrangian formulation. It further means that the inertia forces involve convective terms that, in the numerical solution, result in a nonsymmetric coefficient matrix that depends on the velocities to be calculated.
An advantage of an Eulerian formulation lies in the use of simple stress and strain rate measures, namely, measures that we use in infinitesimal displacement analysis, except that velocities must be calculated instead of displacements. However, if the domain of solution changes, such as in free surface problems, a pure Eulerian formulation would require the creation of new control volumes, and it is more effective to use an arbitrary Lagrangian-Eulerian formulation, as briefly enumerated below.
In the Lagrangian formulations (as discussed in Chapter 6), the mesh moves with (is "attached to") the material particles. Hence, the same material particle is always at the same element mesh point (given in an isoparametric element by the r , s , t coordinates). In the finite element pure Eulerian formulation the mesh points are stationary and the material particles move through the finite element mesh in whichever way the flow conditions govern such movements. In an arbitrary Lagrangian-Eulerian formulation, the mesh points move but not necessarily with the material particles. In fact, the mesh movement corresponds to the nature of the problem and is imposed by the solution algorithm. While the finite element mesh spans the complete analysis domain throughout the solution and its boundaries move with the movements of free surfaces and structural (or solid) boundaries, the fluid particles move relatively to the mesh points. This approach allows the modeling of general free surfaces and interactions between fluid flows and structures (see, for example, J. Donea, S. Giuliani, and J. P. Halleux [A], A. Huerta and W. K. Liu [A], C. Nitikitpaiboon and K. J. Bathe [B], K. J. Bathe, H. Zhang, and M. H. Wang [A] and S. Rugonyi and K. J. Bathe [A]).
TABLE 7.4 Basic continuum mechanics equations used in Lagrangian and Eulerian formulations
| Lagrangian formulation | Eulerian formulation |
| Geometric representation | Geometric representation |
| 0x3, tX3 | x3 |
| 0x1, tX1 | x1 |
| tXj = 0xj + tUi | D/Dt() = ∂/∂t() + tvj ∂/∂xj() |
| Conservation of mass | Conservation of mass |
| m = ∫0v0ρ d0V = ∫tV'ρ d'V ⇒ 0ρ/ρ = det (δX) | m = ∫0v0ρ d0V = ∫tV'ρ d'V ⇒ D/Dt(ρ) + 'ρ ∂'vi/∂xi = 0 |
| Equations of motion | Equations of motion |
| ∂/∂'xj(τij) + 'fBi= 0; 'fBi= 'fBi- 'ρ 'üi | ∂/∂xj(τij) + 'fBi= 0; 'fBi= 'fBi- 'ρ D/Dt(νi) |
| Principle of virtual displacements | Principle of virtual velocities |
| ∫τijδtεij d'V = ∫tV'fBiδuid'V + ∫tS'fSδuSd'S | ∫τij'δeij dV = ∫tV'fBiδvidV + ∫s'fSδvSdS |
| Total Lagrangian formulation: | 'τij = -'p δij + 2μ 'eij; 'eij = 1/2(νi,j + 'vj,i) |
| ∫0vδSij δ0εij d0V = 'R |