...

Chapter_3.pdf

by user

on
Category: Documents
3

views

Report

Comments

Transcript

Chapter_3.pdf
Chapter 3
Linear models
The formulation of linear models is introduced on the basis of linear heat conduction, linear elastostatic and viscous flow problems. Linear models are used
very frequently in engineering practice. Nevertheless they should be viewed as
special cases of nonlinear models that account for material nonlinearities, geometric nonlinearities, mechanical contact, radiation, etc. Non-linear models are
discussed in Chapter 9.
For the purposes of the following discussion a mathematical model is understood to be a statement of a mathematical problem in the form of one or
more ordinary or partial differential equations and specification of the solution
domain and boundary conditions. This is called the strong form of the model.
Beginning with Chapter 4 generalized formulations, also known as the weak
forms, will be discussed.
3.1
Notation
The Euclidean space in n dimensions will be denoted by Rn . The Cartesian1
coordinate axes will be labeled x, y, z (in cylindrical systems r, θ, z) and a vector
u in Rn will be denoted either by ~u or u. For example, ~u ≡ u ≡ {ux uy uz }T
represents a vector in R3 .
We will employ the index notation also where the Cartesian coordinate axes
are labeled x = x1 , y = x2 , z = x3 . The index notation will be introduced
gradually in parallel with the conventional vector notation so that readers not
yet familiar with it can become familiar. The basic rules of index notation are
summarized in the following:
1. A free index is understood to range from one to three (in two spatial
dimensions from one to two).
2. The position vector of a point is xi whereas in conventional notation it is
{x y z}T . A general vector ~a ≡ {ax ay az }T is written simply as ai .
1 René
Descartes (in Latin: Renatus Cartesius) 1596-1650.
75
76
CHAPTER 3. LINEAR MODELS
3. Two free indices represent a matrix. The size of the matrix depends on
the range of indices. Thus, in three dimensions:
 


axx axy axz
a11 a12 a13
aij ≡ a21 a22 a23  ≡ ayx ayy ayz  ·
a31 a32 a33
azx azy azz
The identity matrix is represented by the Kronecker2 delta δij , defined as
follows:
(
1 if i = j
δij =
0 if i 6= j.
4. Indices following a comma represent differentiation with respect to the
variables identified by the indices. For example, if u(xi ) is a scalar function
then
∂u
∂2u
u,2 ≡
,
u,23 ≡
·
∂x2
∂x2 ∂x3
The gradient of u is simply u,i .
5. Repeated indices represent summation. For example, the scalar product
of two vectors ai and bj is ai bi ≡ a1 b1 + a2 b2 + a3 b3 . If ui = ui (xk ) is a
vector function in three dimensional space then
ui,i ≡
∂u1
∂u2
∂u3
+
+
∂x1
∂x2
∂x3
is the divergence of ui . Repeated indices are also called “dummy indices”,
since summation is performed and therefore the index designation is immaterial. For example, ai bi = ak bk . The product of two matrices aij and
bij is: cij = aik bkj .
6. The transformation rules for Cartesian vectors and tensors are presented
in Appendix C.
Example 3.1.1 The divergence theorem in index notation is:
Z
Z
ui,i dV =
ui ni dS
Ω
(3.1)
∂Ω
where ui and ui,i are continuous on the domain Ω and its boundary ∂Ω, ni is
the outward unit normal vector to the boundary, dV is the differential volume
and dS is the differential surface. We will use the divergence theorem in the
derivation of generalized formulations.
Exercise 3.1.1 Outline a derivation of the divergence theorem in two dimensions. Hint: The divergence theorem follows directly from Green’s theorem3 .
2 Leopold
3 George
Kronecker 1823-1891.
Green 1793-1841.
3.2. HEAT CONDUCTION
3.2
77
Heat conduction
The problem of heat conduction is representative of an important class of problems in engineering physics, called potential flow problems, that include, for
example, the flow of viscous fluids in porous media and electrostatic phenomena. Furthermore, the mathematical formulation of potential flow problems has
analogies with some physical phenomena that are not related to potential flow
problems, such as small deflection of membranes and torsion of elastic bars.
Mathematical models of heat conduction are based on two fundamental relationships: The conservation law and Fourier’s law. These are described in the
following.
1. The conservation law states that the quantity of heat entering any
volume element of the conducting medium equals the quantity of heat exiting the volume element plus the quantity of heat retained in the volume
element. The heat retained causes a change in temperature in the volume element which is proportional to the specific heat of the conducting
medium c (in J/(kg K) units) multiplied by the density ρ (in kg/m3 units).
The temperature will be denoted by u(x, y, z, t) where t is time.
The heat flow rate across a unit area is represented by a vector quantity
called heat flux. The heat flux is measured in W/m2 units and will be
denoted by ~q = ~q(x, y, z, t) := {qx (x, y, z, t) qy (x, y, z, t) qz (x, y, z, t)}T .
In addition to heat flux entering and leaving the volume element, heat
may be generated within the volume element, for example from chemical
reactions. The heat generated per unit volume and unit time will be
denoted by Q (in W/m3 units).
Figure 3.1: Control volume and notation for heat conduction.
Applying the conservation law to the volume element shown in Fig. 3.1,
78
CHAPTER 3. LINEAR MODELS
we have:
∆t[qx ∆y∆z − (qx + ∆qx )∆y∆z + qy ∆x∆z − (qy + ∆qy )∆x∆z+
qz ∆x∆y − (qz + ∆qz )∆x∆y + Q∆x∆y∆z] = cρ∆u∆x∆y∆z.
(3.2)
Assuming that u and ~q are continuous and differentiable and neglecting
terms that go to zero faster than ∆x, ∆y, ∆z, ∆t, we have:
∆qx =
∂qx
∆x,
∂x
∆qy =
∂qy
∆y,
∂y
∆qz =
∂qz
∆z,
∂z
∆u =
∂u
∆t.
∂t
On factoring ∆x∆y∆z∆t the conservation law is obtained:
−
∂qx
∂qy
∂qz
∂u
−
−
+ Q = cρ ·
∂x
∂y
∂z
∂t
In index notation:
−qi,i + Q = cρ
∂u
·
∂t
(3.3)
(3.4)
2. Fourier’s law states that the heat flux vector is related to the temperature gradient in the following way:
µ
¶
∂u
∂u
∂u
qx = − kxx
+ kxy
+ kxz
(3.5)
∂x
∂y
∂z
µ
¶
∂u
∂u
∂u
qy = − kyx
+ kyy
+ kyz
(3.6)
∂x
∂y
∂z
µ
¶
∂u
∂u
∂u
qz = − kzx
+ kzy
+ kzz
(3.7)
∂x
∂y
∂z
where the coefficients kxx , kxy , . . . kzz are called coefficients of thermal
conductivity (measured in W/(mK) units). It is customary to write kx :=
kxx , ky := kyy , kz := kzz . The coefficients of thermal conductivity will be
assumed to be independent of the temperature u, unless otherwise stated.
Denoting the matrix of coefficients by [K], Fourier’s law of heat conduction
can be written as:
~q = −[K]gradu.
(3.8)
The matrix of coefficients [K] is symmetric and positive-definite. The
negative sign indicates that the direction of heat flow is opposite to the
direction of the temperature gradient, i.e., the direction of heat flow is
from high to low temperature. In index notation eq. (3.8) is written as:
qi = −kij u,j .
For isotropic materials kij = kδij .
(3.9)
3.2. HEAT CONDUCTION
79
Remark 3.2.1 The values of material properties k, c and ρ vary widely. For
example for pure silver k = 406; for aluminum k = 205; for steel k = 50.2, and
for fiberglass k = 0.04 in W/(m K) units. Material properties are determined
by experiments4 . Typically only average values are reported, reliable statistical
information is rarely available.
3.2.1
The differential equation
Combining equations (3.3) through (3.7), we have:
µ
¶
µ
¶
∂
∂
∂u
∂u
∂u
∂u
∂u
∂u
kx
+ kxy
+ kxz
+
kyx
+ ky
+ kyz
+
∂x
∂x
∂y
∂z
∂y
∂x
∂y
∂z
µ
¶
∂
∂u
∂u
∂u
∂u
kzx
+ kzy
+ kz
+ Q = cρ
(3.10)
∂z
∂x
∂y
∂z
∂t
which can be written in the following compact form:
div ([K]gradu) + Q = cρ
In index notation:
(kij u,j ),i + Q = cρ
∂u
·
∂t
∂u
·
∂t
(3.11)
(3.12)
In many practical problems u is independent of time, that is, the right hand
side of eq. (3.11) and hence eq. (3.12) is zero. Such problems are called stationary
or steady state problems. The solution of a stationary problem can be viewed as
the solution of some time-dependent problem with time-independent boundary
conditions at t = ∞.
In formulating eq. (3.11) we assumed that kij are differentiable functions.
In many practical problems the solution domain is comprised of subdomains
Ωi that have different material properties. In such cases eq. (3.11) is valid
on each subdomain and on the boundaries of adjoining subdomains continuous
temperature and flux are prescribed.
For complete definition of a mathematical model initial and boundary conditions have to be specified. This is discussed in the following section.
3.2.2
Boundary and initial conditions
The solution domain will be denoted by Ω and its boundary by ∂Ω. We will
consider three kinds of boundary conditions:
1. Prescribed temperature. The temperature u = û is prescribed on
boundary region ∂Ωu .
4 See, for example, Standard Test Method for Thermal Conductivity of Solids by Means
of the Guarded-Comparative-Longitudinal Heat Flow Technique, E 1225-87, Annual Book of
ASTM Standards 14.02, ASTM 1994.
80
CHAPTER 3. LINEAR MODELS
2. Prescribed flux. The flux vector component normal to the boundary,
denoted by qn , is prescribed on the boundary region ∂Ωq . By definition;
qn := ~q · ~n ≡ −([K]gradu) · ~n ≡ −kij u,j ni
(3.13)
where ~n ≡ ni is the (outward) unit normal to the boundary. The prescribed flux on ∂Ωq will be denoted by q̂n .
3. Convective heat transfer. On boundary region ∂Ωc the flux vector
component qn is proportional to the difference between the temperature
of the boundary and the temperature of a convective medium:
qn = hc (u − uc )
(x, y, z) ∈ ∂Ωc
(3.14)
where hc is coefficient of convective heat transfer in W/(m2 K) units and
uc is the (known) temperature of the convective medium.
The sets ∂Ωu , ∂Ωq and ∂Ωc are non-overlapping and collectively cover the
entire boundary. Any of the sets may be empty.
The boundary conditions are generally time-dependent. For time-dependent
problems an initial condition has to be prescribed on Ω: u(x, y, z, 0) = U (x, y, z).
It is possible to show that eq. (3.10), subject to the enumerated boundary
conditions, has a unique solution. Stationary problems also have unique solutions with the exception that when flux is prescribed over the entire boundary
of Ω then the following condition must be satisfied:
Z
Z
Q dV =
qn dS.
(3.15)
Ω
∂Ω
This is easily seen by integrating
(kij u,j ),i + Q = 0
(3.16)
on Ω and using the divergence theorem, eq. (3.1) and the definition (3.13).
Note that if ui is a solution of eq. (3.16) then ui + C is also a solution, where
C is an arbitrary constant. Therefore the solution is unique up to an arbitrary
constant.
In addition to the three types of boundary conditions discussed in this section radiation may have to be considered. When two bodies exchange heat by
radiation then the flux is proportional to the difference of the fourth power of
their absolute temperatures, therefore radiation is a non-linear boundary condition. The boundary region subject to radiation, denoted by ∂Ωr , may overlap
∂Ωc . Radiation is discussed in Section 9.1.1.
In the following it will be assumed that the coefficients of thermal conductivity, the flux prescribed on Ωq , and hc and uc prescribed on Ωc are independent
of the temperature u. In general, this assumption is justified in a narrow range
of temperature only.
Exercise 3.2.1 Discuss the physical meaning of eq. (3.15).
3.2. HEAT CONDUCTION
81
Exercise 3.2.2 Show that in cylindrical coordinates r, θ, z the conservation law
is of the form:
1 ∂(rqr ) 1 ∂qθ
∂qz
∂u
−
−
−
+ Q = cρ ·
(3.17)
r ∂r
r ∂θ
∂z
∂t
Use two methods: (a) Apply the conservation law to a differential volume element in cylindrical coordinates and (b) transform eq. (3.3) to cylindrical coordinates.
Exercise 3.2.3 Show that there are three mutually perpendicular directions
(called principal directions) such that the heat flux is proportional to the (negative) gradient vector. Hint: Consider steady state heat conduction and let:
[K]gradu = λ gradu
then show that the principal directions are defined by the normalized characteristic vectors.
Remark 3.2.2 The result of Exercise 3.2.3 implies that the general form of
matrix [K] can be obtained by rotation from orthotropic material axes.
Exercise 3.2.4 List all of the physical assumptions incorporated into the mathematical model represented by eq. (3.10) and the boundary conditions described
in this section.
3.2.3
Symmetry, antisymmetry and periodicity
A scalar function is said to be symmetric with respect to a plane if in symmetrically located points with respect to the plane the function has equal values. On
a plane of symmetry qn = 0. A function is said to be antisymmetric with respect
to a plane if in symmetrically located points with respect to the plane the function has equal absolute values but of opposite sign. On a plane of antisymmetry
u = 0.
In many instances the domain and the coefficients have one or more planes of
symmetry and the source function and boundary conditions are either symmetric
or antsymmetric with respect to the planes of symmetry. In such cases it is often
advantageous to exploit symmetry and antisymmetry in the formulation of the
problem.
When Ω, [K], Q and the boundary conditions are periodic then a periodic
sector of Ω has at least one periodic boundary segment pair denoted by ∂Ω+
p and
+
∂Ω−
.
On
corresponding
points
of
a
periodic
boundary
segment
pair,
P
∈
∂Ω+
p
p
−
−
+
−
+
and P ∈ ∂Ωp , the boundary conditions are: u(P ) = u(P ) and qn = −qn− .
Periodic, symmetric and antisymmetric boundary conditions are illustrated in
the following example.
Example 3.2.1 Figure 3.2 represents a plate-like body of constant thickness.
It can be divided into five sectors as illustrated. Let us assume that on the
cylindrical boundary represented by the inner circle (∂Ω(i) ) constant flux is
82
CHAPTER 3. LINEAR MODELS
prescribed and on the boundary represented by the outer circle (∂Ω(o) ) constant
temperature is prescribed. The planar surfaces parallel to the x y plane are
perfectly insulated. On the surfaces represented by the elliptical boundaries let
(
qn =
C sin ϕ for −π/2 ≤ ϕ ≤ π/2
0
for π/2 < ϕ < 3π/2
where C is constant and the formula is understood to be given in the local
coordinate system of each ellipse. In this case the solution is periodic and the
solution obtained for one sector can be extended to the other sectors. Periodic
−
boundary conditions are prescribed on ∂Ω+
p and ∂Ωp .
Figure 3.2: Notation for Example 3.2.1.
If qn = C sin ϕ is prescribed on the elliptical boundaries in the local coordinate system of each ellipse then the solution is symmetric with respect to the y
axis, see Fig. 3.2(b). If qn = C cos ϕ is prescribed on the elliptical boundaries
in the local coordinate system of each ellipse then the solution is antisymmetric
with respect to the y axis.
Exercise 3.2.5 Consider the domain shown in Fig. 3.2(a). Assume that u = 0
on the boundary represented by the inner circle (∂Ω(i) ) and qn = 0 on the boundaries represented by ellipses. Prescribe sinusoidal fluxes on the outer boundary
(∂Ω(o) ) that will result in (a) periodic (b) symmetric and (c) antisymmetric
solutions.
3.2. HEAT CONDUCTION
3.2.4
83
Dimensional reduction
In many important practical applications reduction of the number of dimensions
is possible without significantly affecting the data of interest. In other words, a
mathematical model in one or two dimensions may be an acceptable replacement
for the fully three-dimensional model.
Planar problems
Consider a plate-like body shown in Fig. 3.3. The thickness tz will be assumed
constant unless otherwise stated. The mid-surface is the solution domain which
will be denoted by Ω. The solution domain lies in the xy plane. The boundary
points of Ω (shown as a dotted line) will be denoted by Γ. The unit outward
normal to the boundary is denoted by ~n.
Figure 3.3: Notation for two-dimensional models.
The heat conduction problem in two dimensions is
µ
¶
µ
¶
∂
∂u
∂u
∂
∂u
∂u
∂u
kx
+ kxy
+
kyx
+ ky
+ Q̄ = cρ
∂x
∂x
∂y
∂y
∂x
∂y
∂t
(3.18)
where the meaning of Q̄ depends on the boundary conditions prescribed on the
top and bottom surfaces (z = ±tz /2) as described in the following.
Eq. (3.18) represents one of two cases:
1. The thickness is large in comparison with the other dimensions and both
the material properties and boundary conditions are independent of z, i.e.,
u(x, y, z) = u(x, y). This is equivalent to the case of finite thickness with
the top and bottom surfaces (z = ±tz /2) perfectly insulated and Q̄ = Q.
2. The thickness is small in relation to the other dimensions. In this case
the two dimensional solution is an approximation of the three dimensional
solution that can be interpreted as the first term in the expansion of the
three dimensional solution with respect to the z coordinate. The definition
of Q̄ depends on the boundary conditions on the top and bottom surfaces
as explained below.
(a) Prescribed flux: Let us denote the heat flux prescribed on the top
(resp. bottom) surface by q̂n+ (resp. q̂n− ). Note that positive q̂n is heat
84
CHAPTER 3. LINEAR MODELS
flux exiting the body. The amount of heat exiting the body over a
small area ∆A per unit time is (q̂n+ + q̂n− )∆A. Dividing by ∆Atz heat
generated per unit volume becomes
Q̄ = Q − (q̂n+ + q̂n− )
1
·
tz
(3.19)
(b) Convective heat transfer: Let us denote the coefficient of convective
−
heat transfer at z = tz /2 (resp. z = −tz /2) by h+
c (resp. hc ) and the
corresponding temperature of the convective medium by u+
c (resp.
u−
c ) then the amount of heat exiting the body over a small area ∆A
−
−
+
per unit time is [h+
c (u − uc ) + hc (u − uc )]∆A. Therefore the heat
generated per unit volume is changed to
+
−
−
Q̄ = Q − [h+
c (u − uc ) + hc (u − uc )]
1
·
tz
(3.20)
Of course, combinations of these boundary conditions are possible. For example
flux may be prescribed on the top surface and convective boundary conditions
may be prescribed on the bottom surface.
In the two-dimensional formulation u is assumed to be constant through
the thickness. Therefore prescribing a temperature has meaning only if the
temperature is the same on the top and bottom surfaces, as well as on the side
surface in which case the solution is the prescribed temperature.
Exercise 3.2.6 Using the control volume shown in Fig. 3.4, derive eq. (3.18)
from first principles.
Figure 3.4: Control volume and notation for heat conduction in 2D.
Exercise 3.2.7 Assume that tz = tz (x, y) > 0 is a continuous and differentiable
function and the maximum value of tz is small in comparison with the other
dimensions. Assume further that convective heat transfer occurs on the surfaces
−
+
−
z = ±tz /2 with h+
c = hc = hc and uc = uc = uc . Using a control volume
similar to that shown in Fig. 3.4, but accounting for variable thickness, show
that in this case the conservation law for heat conduction in two dimensions is:
−
∂
∂u
∂
(tz qx ) −
(tz qy ) − 2hc (u − uc ) + Qtz = cρtz ·
∂x
∂y
∂t
(3.21)
3.2. HEAT CONDUCTION
85
Remark 3.2.3 In Exercise 3.2.7 the thickness tz was assumed to be continuous
and differentiable on Ω. If tz is continuous and differentiable over two or more
subdomains of Ω, but discontinuous on the boundaries of the subdomains, then
eq. (3.21) is applicable on each subdomain subject to the requirement that qn tz
is continuous on the boundaries of the subdomains.
Example 3.2.2 The plan view of a conducting medium is shown in Fig. 3.5.
The thickness tz is constant. The top and bottom surfaces (z = ±tz /2) are
perfectly insulated. On the side surfaces (x = ±b, y = ±c) the temperature is
constant (u = û0 ). Let kx = ky = k, kxy = 0 and Q = Q0 given that k and
Q0 constant. The goal is to determine the stationary temperature distribution.
The mathematical problem is to solve:
µ 2
¶
∂ u ∂2u
k
+ 2 + Q0 = 0
∂x2
∂y
on the rectangular domain shown in Fig. 3.5 with u = û0 on the boundary.
Figure 3.5: Solution domain. Example 3.2.2.
The solution of this problem can be determined by classical methods5 :
¶
n−1 µ
∞
X
Q0 b2
(−1) 2
cosh nπy/2b
u = û0 + 16
1−
cos nπx/2b.
k π 3 n=1,3,5,...
n3
cosh nπc/2b
(3.22)
This infinite series converges absolutely. It is seen that the classical solution
of this seemingly simple problem is rather complicated and in fact the exact
solution can be computed only approximately, although the truncation error
can be made arbitrarily small by computing a sufficiently large number of terms
of the infinite series.
Exercise 3.2.8 Assume that the coefficients of thermal conductivity kx , kxy ,
ky are given. Show that in a Cartesian coordinate system x0 y 0 , rotated counterclockwise by the angle α relative to the xy system, the coefficients will be:



 

cos α − sin α
kx kxy
cos α sin α
kx0 kx0 y0
·


=

0
0
0
sin α cos α
kyx ky
− sin α cos α
ky x
ky
5 See, for example, Timoshenko, S. and Goodier, J. N., Theory of Elasticity, McGraw-Hill,
New York, 2nd ed. 1951, pp. 275-276.
86
CHAPTER 3. LINEAR MODELS
Hint: A scalar {a}T [K]{a}, where {a} is an arbitrary vector, is invariant under
coordinate transformation by rotation.
Axisymmetric problems
In many important practical problems the topological description, the material
properties and the boundary conditions are axially symmetric. A simple example is a straight pipe of constant wall thickness. In such cases the problem
is usually formulated in cylindrical coordinates and, since the solution is independent of the circumferential variable, the number of dimensions is reduced to
two. In the following the z axis will be the axis of symmetry and the radial
(resp. circumferential) coordinates will be denoted by r (resp. θ). Referring to
the result of Exercise 3.2.2, and letting qθ = 0, the conservation law is
1 ∂(rqr ) ∂qz
∂u
−
+ Q = cρ ·
r ∂r
∂z
∂t
Substituting the axisymmetric form of Fourier’s law:
−
∂u
∂u
,
qz = −kz
∂r
∂z
we have the formulation of the axisymmetric heat conduction problem in cylindrical coordinates:
µ
¶
µ
¶
1 ∂
∂u
∂
∂u
∂u
rkr
+
kz
+ Q = cρ ·
(3.23)
r ∂r
∂r
∂z
∂z
∂t
qr = −kr
One or more segments of the boundary may lie on the z axis. Implied in
the formulation is that the boundary condition is the zero flux condition on
those segments. Therefore it would not be meaningful to prescribe essential
boundary conditions on those boundary segments. To show this, consider an
axisymmetric problem of heat conduction, the solution of which is independent
of z. For simplicity we assume that kr = 1. In this case the problem is essentially
one-dimensional:
µ
¶
1 d
du
r
=0
ri < r < ro .
r dr
dr
Assuming that the boundary conditions u(ri ) = ûi , u(ro ) = ûo are prescribed.
the exact solution of this problem is:
u(r) =
ûo − ûi
ûi ln ro − ûo ln ri
ln r +
·
ln ro − ln ri
ln ro − ln ri
(3.24)
Consider now the solution in an arbitrary fixed point r = % where ri < % < ro
and let ri → 0:
¶
µ
ln %
ûi ln ro / ln ri − ûo
ûo − ûi
+
= ûo .
lim u(%) = lim
ri →0
ri →0 ln ro / ln ri − 1 ln ri
ln ro / ln ri − 1
Therefore the solution is independent of ûi when ri = 0. It is left to the reader in
Exercise 3.2.9 to show that du/dr → 0 as ri → 0, hence the boundary condition
at r = 0 is the zero flux boundary condition.
3.2. HEAT CONDUCTION
87
Exercise 3.2.9 Refer to the solution given by eq. (3.24). Show that
µ
du
dr
¶
→ 0 as ri → 0
r=ri
independent of ûi and ûo .
Exercise 3.2.10 Derive eq. (3.23) by considering a control volume in cylindrical coordinates and using the assumption that the temperature is independent
of the circumferential variable.
Exercise 3.2.11 Consider water flowing in a stainless steel pipe. The temperature of the water is 80 ◦ C. The outer surface of the pipe is cooled by air flow.
The temperature of the air is 20 ◦ C. The outer diameter of the pipe is 0.20 m
and its wall thickness is 0.01 m. (a) Assuming that convective heat transfer
occurs on both the inner and outer surfaces of the pipe and u is a function of
r only, formulate the mathematical model for stationary heat transfer. (b) Assume that the coefficient of thermal conductivity of stainless steel is 20 W/mK.
(a)
(w)
Using hc = 750 W/m2 K for the water and hc = 10 W/m2 K for the air,
determine the temperature of the external surface of the pipe and the rate of
heat loss per unit length.
Heat conduction in a bar
One-dimensional models of heat conduction are discussed in this section. The
solution domain is a bar one end of which is located in x = 0, the other end
in x = `, and the cross-sectional area A = A(x) > 0 is a continuous and
differentiable function.
One-dimensional formulations are justified when the solution of the threedimensional problem represented by the one-dimensional model is a function of
x only and/or the cross-section A is small. If the solution is known to be a
function of x only then the bar is perfectly insulated along its length. If the
cross-section is small then the boundary conditions on the surface of the bar are
transferred to the differential equation.
We will assume in the following that the cross-sectional area is small and
convective heat transfer may occur along the bar, as described in Section 3.2.4.
In this case the conservation law is:
−
∂(Aq)
∂u
− cb (u − ua ) + QA = cρA
∂x
∂t
where cb = cb (x) is the coefficient of convective heat transfer of the bar (in
W/(mK) units) obtained from hc by integration:
I
cb =
hc ds
88
CHAPTER 3. LINEAR MODELS
with the contour integral taken along the perimeter of the cross-section. Therefore the differential equation of heat conduction in a bar is:
µ
¶
∂
∂u
∂u
Ak
− cb (u − ua ) + QA = cρA ·
(3.25)
∂x
∂x
∂t
One of the boundary conditions described in Section 3.2 is prescribed at x = 0
and x = `.
Example 3.2.3 Consider stationary heat flow in a partially insulated bar of
length ` and constant cross-section A. The coefficients k and cb are constant
and Q = 0. Therefore eq. (3.25) can be cast in the following form:
u00 − λ2 (u − ua ) = 0,
λ2 :=
cb
·
Ak
If the temperature ua is a linear function of x, i.e., ua (x) = a + bx and the
boundary conditions are u(0) = û0 , q(`) = q̂` then the solution of this problem
is:
u = C1 cosh λx + C2 sinh λx + a + bx
where:
C1 = û0 − a,
C2 = −
1
λ cosh λ`
µ
¶
q̂`
+ (û0 − a)λ sinh λ` + b .
k
Exercise 3.2.12 Solve the problem described in Example 3.2.3 using the following boundary conditions: q(0) = q̂0 , q(`) = h` (u(`) − Ua ) where q̂0 , h` , Ua
are given data.
Exercise 3.2.13 A perfectly insulated bar of constant cross section, length `,
thermal conductivity k, density ρ and specific heat c is subject to the initial condition u(x, 0) = U0 (constant) and the boundary conditions u(0, t) = u(`, t) = 0.
Assuming that Q = 0, show that the solution of this problem is:
µ 2 2 ¶
∞
³ x´
X
1 − cos(nπ)
n π k
u = 2U0
exp − 2
t sin nπ
·
nπ
` cρ
`
n=1
Exercise 3.2.14 Solve the problem of Exercise 3.2.13 with the initial condition
u(x, 0) = U0 (x/` − x2 /`2 ).
Remark 3.2.4 As noted in Section 2.1.4, more than one physical phenomenon
can be represented by the same differential equation. For example, flow of
viscous fluids in porous media is based on the conservation law where the flux
represents the flow rate of an incompressible viscous fluid per unit area (m/s
units) and Q represents a distributed source or sink injecting or extracting fluid
(1/s units). The potential function u represents the piezometric head and the
relationship between the piezometric gradient and the flux is formally identical
3.3. THE SCALAR ELLIPTIC BOUNDARY VALUE PROBLEM
89
to Fourier’s law given by eq. (3.8), except that it is called Darcy’s law6 , and the
elements of matrix [K] are called coefficients of permeability (m/s units).
The boundary conditions are analogous to the linear boundary conditions
described in Section 3.2; the piezometric head; the flux, or a linear combination
of the flux and piezometric head may be prescribed. In viscous flow problems
there is no physical analogy with radiation. On the other hand, one of the
boundaries may be a free surface. On a free surface the piezometric head equals
the elevation head, i.e., u equals the elevation with respect to a datum plane.
Furthermore, under stationary conditions the flux vector is tangential to the
free surface. The position of the free surface is unknown a priori and must be
determined by an iterative process.
3.3
The scalar elliptic boundary value problem
In view of Remark 3.2.4, a generic treatment of diverse physical phenomena
is possible when their common mathematical basis is exploited. We will be
concerned with the following model problem:
−div ([κ]grad u) + cu = f (x, y, z)
where

κx
[κ] := κyx
κzx
κxy
κy
κzy
(x, y, z) ∈ Ω

κxz
κyz 
κz
(3.26)
(3.27)
is a positive-definite matrix and c = c(x, y, z) ≥ 0. In index notation eq. (3.26)
reads:
−(κij u,j ),i + cu = f.
(3.28)
We will be concerned with the following boundary conditions:
1. Essential or Dirichlet boundary condition: u = û is prescribed on boundary region ∂Ωu .
When û = 0 on ∂Ωu then the Dirichlet boundary
condition is said to be homogeneous.
2. Neumann boundary condition: −[κ]grad u·~n = q̂n is prescribed on boundary region ∂Ωq . In this expression ~n is the unit outward normal to the
boundary, shown in Fig. 3.3. When q̂n = 0 on ∂Ωq then the Neumann
boundary condition is said to be homogeneous.
3. Robin boundary condition: −[κ]grad u·~n = hR (u−uR ) is given on boundary segment ∂ΩR . In this expression hR > 0 and uR are given functions.
When uR = 0 on ∂ΩR then the Robin boundary condition is said to be
homogeneous.
6 Henry
Philibert Gaspard Darcy 1803-1858.
90
CHAPTER 3. LINEAR MODELS
The boundary segments ∂Ωu , ∂Ωq , ∂ΩR and ∂Ωp are non-overlapping and collectively cover the entire boundary ∂Ω. Any of the boundary segments may be
empty.
We will consider restrictions of eq. (3.26) to one and two dimensions as well.
In one dimension we will use
µ
¶
d
du
−
κ
+ cu ≡ −(κu0 )0 + cu = f (x)
(3.29)
dx
dx
on the domain I = (0, `) with boundary conditions prescribed on the endpoints.
3.4
Linear elasticity
Mathematical models of linear elastostatic and elastodynamic problems are
based on three fundamental relationships: The strain-displacement equations,
the stress-strain relationships and the equilibrium equations. The unknowns are
the components of the displacement vector:
~u := {ux (x, y, z) uy (x, y, z) uz (x, y, z)}T ≡ ui (xj ).
1. Strain-displacement relationships. We will introduce the infinitesimal strain-displacement relationships here. A detailed derivation of these
relationships is presented in Section 9.2.1. By definition, the infinitesimal
normal strain components are:
²x ≡ ²xx :=
∂ux
∂x
²y ≡ ²yy :=
∂uy
∂y
²z ≡ ²zz :=
∂uz
∂z
(3.30)
and the shear strain components are:
µ
¶
γxy
1 ∂ux
∂uy
:=
+
2
2 ∂y
∂x
µ
¶
γyz
1 ∂uy
∂uz
≡
:=
+
2
2 ∂z
∂y
µ
¶
γzx
1 ∂uz
∂ux
≡
:=
+
2
2 ∂x
∂z
²xy = ²yx ≡
²yz = ²zy
²zx = ²xz
(3.31)
where γxy , γyz , γzx are called the engineering shear strain components. In
index notation, the state of (infinitesimal) strain at a point is characterized
by the strain tensor
1
²ij := (ui,j + uj,i ) .
(3.32)
2
2. Stress-strain relationships Mechanical stress is defined as force per
unit area (N/m2 ≡ Pa). Since one Pascal (Pa) is a very small stress, the
usual unit of mechanical stress is the Megapascal (MPa) which can be
interpreted either as 106 N/m2 or as 1 N/mm2 .
3.4. LINEAR ELASTICITY
91
Figure 3.6: Notation for stress components.
The usual engineering notation for stress components is illustrated on an
infinitesimal volume element shown in Fig. 3.6. The indexing rules are
as follows: Faces to which the positive x, y, z axes are normal are called
positive faces, the opposite faces are called negative faces. The normal
stress components are denoted by σ, the shear stresses components by τ .
The normal stress components are assigned one subscript only, since the
orientation of the face and the direction of the stress component are the
same. For example, σx is the stress component acting on the faces to which
the x-axis is normal and the stress component is acting in the positive
(resp. negative) coordinate direction on the positive (resp. negative) face.
For the shear stresses, the first index refers to the coordinate direction of
the normal to the face on which the shear stress is acting. The second index
refers to the coordinate direction in which the shear stress component is
acting.
On a positive (resp. negative) face the positive stress components are
oriented in the positive (resp. negative) coordinate directions. The reason
for this is that if we subdivide a solid body into infinitesimal hexahedral
volume elements similar to the element shown in Fig. 3.6, then a negative
face will be coincident with each positive face. By the action-reaction
principle, the forces acting on coincident faces must be equal in magnitude
and opposite in sign.
An alternative notation is σxx ≡ σx , σxy ≡ τxy , etc. which corresponds
directly to the index notation σij .
The mechanical properties of isotropic elastic materials are characterized
by the modulus of elasticity E > 0, Poisson’s ratio7 ν, and the coefficient
of thermal expansion α > 0. The stress-strain relationships, known as
7 Simeon
Denis Poisson 1781-1840.
92
CHAPTER 3. LINEAR MODELS
Hooke’s law8 , are:
1
(σx − νσy − νσz ) + αT∆
E
1
²y = (−νσx + σy − νσz ) + αT∆
E
1
²z = (−νσx − νσy + σz ) + αT∆
E
2(1 + ν)
τxy
γxy ≡ 2²xy =
E
2(1 + ν)
γyz ≡ 2²yz =
τyz
E
2(1 + ν)
γzx ≡ 2²zx =
τzx
E
²x =
(3.33)
(3.34)
(3.35)
(3.36)
(3.37)
(3.38)
where T∆ = T∆ (x, y, z) is the temperature change with respect to a reference temperature at which the strain is zero. The strain components
represent the total strain; αT∆ is the thermal strain. Mechanical strain is
defined as the total strain minus the thermal strain.
In index notation Hooke’s law can be written as
²ij =
1+ν
ν
σij − σkk δij + αT∆ δij .
E
E
(3.39)
The inverse is:
σij = λ²kk δij + 2G²ij − (3λ + 2G)αT∆ δij
(3.40)
where λ and G, called the Lamé constants9 , are defined by
λ :=
Eν
,
(1 + ν)(1 − 2ν)
G :=
E
·
2(1 + ν)
(3.41)
G is also called shear modulus and modulus of rigidity. Since λ and G
are positive, the range of Poisson’s ratio is −1 < ν < 1/2. Typically
0 ≤ ν < 1/2.
The generalized Hooke’s law states that the components of the stress tensor are linearly related to the mechanical strain tensor:
σij = Cijkl (²kl − αkl T∆ )
(3.42)
where Cijkl and αij are Cartesian tensors. By symmetry considerations
the maximum number of independent elastic constants that characterize
Cijkl is 21. The symmetric tensor αij is characterized by six independent coefficients of thermal expansion. This is the most general form of
anisotropy in linear elasticity.
8 Robert
9 Gabriel
Hooke 1635-1703.
Lamé 1795-1870.
3.4. LINEAR ELASTICITY
93
According to eq. (3.42) the elastic body is stress free when the mechanical strain tensor is zero. There may be, however, an initial stress field
(0)
σij present in the reference configuration. The initial stress field, called
residual stress, must satisfy the equilibrium equations and traction-free
boundary conditions. Residual stresses may be introduced by manufacturing processes10 , forging, machining, and various types of cold-working
operations. In composite materials there is a large difference in the coefficients of thermal expansion of the fiber and matrix. Residual stresses
are introduced when a part cools following curing operations at elevated
temperatures. When residual stresses are present then we have
(0)
σij = σij + Cijkl (²kl − αkl T∆ ).
(3.43)
Accurate determination of residual stresses is generally difficult.
3. Equilibrium Considering the dynamic equilibrium of a volume element,
similar to that shown in Fig. 3.6, except that the edges are of length
∆x, ∆y, ∆z, six equations of equilibrium are written: The resultants of
the forces and moments must vanish. Assuming that the material is not
loaded by distributed moments (body moments), consideration of moment
equilibrium leads to the conclusion that τxy = τyx , τyz = τzy , τzx = τxz ,
i.e., the stress tensor is symmetric. Assuming further that the components of the stress tensor are continuous and differentiable, application of
d’Alembert’s principle11 and consideration of force equilibrium leads to
three partial differential equations:
∂σx
∂τxy
∂τxz
∂ 2 ux
+
+
+ Fx − % 2 =0
∂x
∂y
∂z
∂t
∂τxy
∂σy
∂τyz
∂ 2 uy
+
+
+ Fy − % 2 =0
∂x
∂y
∂z
∂t
∂τxz
∂τyz
∂σz
∂ 2 uz
+
+
+ Fz − % 2 =0
∂x
∂y
∂z
∂t
(3.44)
(3.45)
(3.46)
where Fx , Fy , Fz are the components of the body force vector (in N/m3
units), % is the specific density (in kg/m3 ≡ Ns2 /m4 units). These equations are called the equations of motion. In index notation:
σij,j + Fi = %
∂ 2 ui
·
∂t2
(3.47)
For elastostatic problems the time derivative is zero and the boundary
conditions are independent of time. This yields the equations of static
equilibrium:
σij,j + Fi = 0.
(3.48)
10 For example, 7050-T7451 aluminum plates are hot-rolled, quenched, over-aged and
stretched by the imposition of 1.5 to 3.0 percent strain in the rolling direction.
11 Jean Le Rond d’Alembert 1717-1783
94
CHAPTER 3. LINEAR MODELS
3.4.1
The Navier equations
The equations of motion, called the Navier equations12 , are obtained by substituting eq. (3.40) into eq. (3.47). In elastodynamics the effects of temperature
are usually negligible, hence we will assume T∆ = 0:
Gui,jj + (λ + G)uj,ji + Fi = %
∂ 2 ui
·
∂t2
(3.49)
In elastostatic problems we have:
Gui,jj + (λ + G)uj,ji + Fi = (3λ + 2G)α(T∆ ),i .
(3.50)
Exercise 3.4.1 Derive eq. (3.50) by substituting eq. (3.40) into eq. (3.48). Indicate the rules under which the indices are changed to obtain eq. (3.50). Hint:
(²kk δij ),j = (uk,k δij ),j = uk,ki = uj,ji .
Exercise 3.4.2 Derive the equilibrium equations from first principles.
Exercise 3.4.3 In deriving eq. (3.49) and (3.50) we assumed that λ and G
are constants. Formulate the analogous equations assuming that λ and G are
smooth functions of xi ∈ Ω.
Exercise 3.4.4 Assume that Ω is the union of two or more subdomains and the
material properties are constants on each subdomain but vary from subdomain
to subdomain. Formulate the elastostatic problem for this case.
3.4.2
Boundary and initial conditions
As in the case of heat conduction, we will consider three kinds of boundary
conditions: prescribed displacements, prescribed tractions and spring boundary
conditions.
1. Prescribed displacement. One or more components of the displacement vector is prescribed on all or part of the boundary. This is called a
kinematic boundary condition.
2. Prescribed traction. One or more components of the traction vector
is prescribed on all or part of the boundary. The definition of traction
vector is given in Appendix C.1.
3. Linear spring. A linear relationship is prescribed between the traction
and displacement vector components. The general form of this relationship is:
Ti = cij (dj − uj )
(3.51)
where Ti is the traction vector, cij is a positive-definite matrix that represents the distributed spring coefficients; dj is a prescribed function that
12 Claude
Louis Marie Henri Navier 1785-1836.
3.4. LINEAR ELASTICITY
95
represents displacement imposed on the spring and uj is the (unknown)
displacement vector function on the boundary. The spring coefficients cij
(in N/m3 units) may be functions of the position xk but are independent
of the displacement ui . This is called a“Winkler spring13 ”.
A schematic representation of this boundary condition on an infinitesimal
boundary surface element is shown in Fig. 3.7 under the assumption that
cij is a diagonal matrix and therefore three spring coefficients c1 := c11 ,
c2 := c22 , c3 := c33 characterize the elastic properties of the boundary
condition.
Figure 3.7: Spring boundary condition. Schematic representation.
Fig. 3.7 should be interpreted to mean that the imposed displacement di
will cause a differential force ∆Fi to act on the centroid of the surface
element. Suspending the summation rule, the magnitude of ∆Fi is
∆Fi = ci ∆A(di − ui ),
i = 1, 2, 3
where ui is the displacement of the surface element. The corresponding
traction vector is:
∆Fi
= ci (di − ui )
∆A→0 ∆A
Ti = lim
i = 1, 2, 3.
The enumerated boundary conditions may occur in any combination. For
example, the displacement vector component u1 , the traction vector component
T2 and a linear combination of T3 and u3 may be prescribed on a boundary
segment.
In engineering practice boundary conditions are most conveniently prescribed
in the normal-tangent reference frame. The normal is uniquely defined on
13 Emil
Winkler 1835-1888.
96
CHAPTER 3. LINEAR MODELS
smooth surfaces but the tangential coordinate directions are not. It is necessary to specify the tangential coordinate directions with respect to the reference frame used. The required coordinate transformations are discussed in
Appendix C.
The boundary conditions are generally time-dependent. For time-dependent
problems the initial conditions, that is, the initial displacement and velocity
fields, denoted respectively by U (x, y, z) and V (x, y, z), have to be prescribed:
µ
~ (x, y, z) and
~u(x, y, z, 0) = U
∂~u
∂t
¶
~ (x, y, z).
=V
(x,y,z,0)
Exercise 3.4.5 Assume that the following boundary conditions are given in the
normal-tangent reference frame x0i , x01 being coincident with the normal: T10 =
c01 (d01 − u01 ); T20 = T30 = 0. Using the transformation x0i = gij xj , determine the
boundary conditions in the xi coordinate system. (See the Appendix, Section
C.3.)
3.4.3
Symmetry, antisymmetry and periodicity
Symmetry and antisymmetry of a vector function with respect to a line is illustrated in Fig. 3.8 where points A and A0 are equidistant from the y axis. In
Fig. 3.8(a) vector ~v is the mirror image of vector ~u with respect to the y axis,
the axis of symmetry. Note that vx = −ux and vy = uy . In Fig. 3.8(b) vector
~v is the antisymmetric image of vector ~u with respect to the y axis, the axis of
antisymmetry. Note that vx = ux and vy = −uy .
Figure 3.8: Symmetry and antisymmetry of a vector function in two dimensions.
In three dimensions the corresponding vector components parallel to the
plane of symmetry (resp. antisymmetry) have the same absolute value and the
same (resp. opposite) sign. The corresponding vector components normal to
the plane of symmetry (resp. antisymmetry) have the same absolute value and
opposite (resp. same) sign.
In a plane of symmetry the normal displacement and the shearing traction
components are zero. In a plane of antisymmetry the normal traction is zero
and the in-plane components of the displacement vector are zero.
3.4. LINEAR ELASTICITY
97
When the solution is periodic on Ω then a periodic sector of Ω has at least
−
one periodic boundary segment pair denoted by ∂Ω+
p and ∂Ωp . On correspond+
+
ing points of a periodic boundary segment pair, P ∈ ∂Ωp and P − ∈ ∂Ω−
p the
normal component of the displacement vector and the periodic in-plane components of the displacement vector have the same value. The normal component of
the traction vector and the periodic in-plane components of the traction vector
have the same absolute value but opposite sign.
Exercise 3.4.6 A homogeneous isotropic elastic body with Poisson’s ratio zero
occupies the domain Ω = {x, y, z | |x| < a, |y| < b, |z| < c}. Define tractions on
the boundaries of Ω such that the tractions satisfy equilibrium and the plane
z = 0 is a plane of (a) symmetry, (b) antisymmetry.
Exercise 3.4.7 Consider the domain shown in Fig. 3.2. Assume that Tn = 0
and Tt = τo where τo is constant, is prescribed on ∂Ω(o) , Tn = Tt = 0 on the
elliptical boundaries and un = ut = 0 on ∂Ω(i) . Specify periodic boundary
−
conditions on ∂Ω+
p and ∂Ωp .
3.4.4
Dimensional reduction
Owing to the complexity of three-dimensional problems in elasticity, dimensional
reduction is widely used. Various kinds of dimensional reduction are possible
in elasticity, such as planar, axisymmetric, shell, plate, beam and bar models.
Each of these model types is sufficiently important to have generated a substantial technical literature. In the following models for planar and axisymmetric
problems and axially loaded bars are discussed. Models for beams, plates and
shells will be discussed separately.
Planar elastostatic models
Let us consider a prismatic body of length `. The material points occupy the
domain Ω` , defined as follows:
Ω` = {(x, y, z) | (x, y) ∈ ω, −`/2 < z < `/2, ` > 0}
(3.52)
where ω ∈ R2 is a bounded domain The lateral boundary of the body is denoted
by
Γ` = {(x, y, z) | (x, y) ∈ ∂ω, −`/2 < z < `/2, ` > 0}
(3.53)
and the faces are denoted by
γ± = {(x, y, z) | (x, y) ∈ ω, z = ±`/2}.
(3.54)
The notation is shown in Fig. 3.9. The diameter of ω will be denoted by dω .
The material properties, the volume forces and temperature change acting on
Ω` and tractions acting on Γ` will be assumed to be independent of z. Therefore
the x y plane is a plane of symmetry. It is assumed that tractions are specified
98
CHAPTER 3. LINEAR MODELS
Figure 3.9: Notation.
on the entire boundary Γ` . This is usually called the first fundamental boundary
value problem of elasticity.
When the tractions acting on γ± are zero and `/dω << 1 then such problems are usually formulated on ω as plane stress problems. When the normal
displacements and shearing tractions acting on γ± are zero then such problems
are formulated on ω as plane strain problems, independently of the `/dω ratio.
When 1 << `/dω and the tractions acting on γ± are zero then the three dimensional thermoelasticity problem on Ω` is called the generalized plane strain
problem. In a generalized plane strain problem the stress resultants corresponding to σz must be zero:
Z
Z
Z
σz dzdy = 0
σz x dzdy = 0
σz y dzdy = 0.
∂ω
∂ω
∂ω
For details we refer to [12].
Planar elastostatic models are specializations of the Navier equations. The
mid-plane is understood to be a plane of symmetry, that is, uz (x, y, 0) = 0. The
formulation is written in unabridged notation in the following.
1. The linear strain-displacement relationships are:
∂ux
∂x
∂uy
²y :=
∂y
∂ux
∂uy
γxy :=
+
·
∂y
∂x
²x :=
(3.55)
(3.56)
(3.57)
In two dimensional problems the shearing strain components γyz and γzx
and the corresponding shear stress components τyz , τzx are zero.
3.4. LINEAR ELASTICITY
99
2. Stress-strain relationships. We will be concerned with two special
cases, the case of plane stress (σz = τyz = τzx = 0) and the case of plane
strain (²z = γyz = γzx = 0). Applying the appropriate restrictions to the
three-dimensional Hooke’s law, from eq. (3.33) and eq. (3.34), we have for
plane stress:
 

 
 
1 ν
0
 σx 
 ²x  EαT 1
E ν 1

∆
0
1
σy =
(3.58)
−

 ²
1−ν  y 
  1 − ν2
1−ν  
0
τxy
γxy
0 0
2
and, letting ²z = 0 in eq. (3.35) we have for in plane strain:
  
 
 
λ + 2G
λ
0  ²x 
 σx 
1
σy =  λ
λ + 2G 0  ²y − (3λ + 2G)αT∆ 1 · (3.59)
 
 
 
τxy
0
0
G
γxy
0
3. The static equations of equilibrium are:
∂σx
∂τxy
+
+ Fx =0
∂x
∂y
∂τyx
∂σy
+
+ Fy =0
∂x
∂y
(3.60)
(3.61)
where Fx , Fy are the components of body force vector F~ (x, y) (in N/m3
units).
The Navier equations are obtained by substituting the stress-strain and straindisplacement relationships into the equilibrium equations. For plane strain;
µ
¶
µ 2
¶
∂ ∂ux
∂uy
∂ ux
∂ 2 ux
Eα ∂T∆
(λ + G)
+
+G
+
=
− Fx (3.62)
∂x ∂x
∂y
∂x2
∂y 2
1 − 2ν ∂x
µ
¶
µ 2
¶
∂ ∂ux
∂uy
∂ uy
∂ 2 uy
Eα ∂T∆
(λ + G)
+
+G
+
=
− Fy . (3.63)
∂y ∂x
∂y
∂x2
∂y 2
1 − 2ν ∂y
The boundary conditions are most conveniently written in the normal-tangent
(nt) reference frame illustrated in Fig. 3.3. The relationship between the xy
and nt components of displacements and tractions is established by the rules
of vector transformation, described in the Appendix, Section C.3. The linear
boundary described in Section 3.4.2 are directly applicable in two dimensions.
Exercise 3.4.8 Derive eq. (3.58) and eq. (3.59) from Hooke’s law.
Exercise 3.4.9 Show that for plane stress the Navier equations are:
µ
¶
µ 2
¶
E
∂ ∂ux
∂uy
∂ ux
∂ 2 ux
Eα ∂T∆
+
+G
+
=
− Fx
2
2
2(1 − ν) ∂x ∂x
∂y
∂x
∂y
1 − ν ∂x
µ
¶
µ 2
¶
E
∂ ∂ux
∂uy
∂ uy
∂ 2 uy
Eα ∂T∆
+
+G
+
=
− Fy .
2
2
2(1 − ν) ∂y ∂x
∂y
∂x
∂y
1 − ν ∂y
(3.64)
(3.65)
100
CHAPTER 3. LINEAR MODELS
Exercise 3.4.10 Denote the components of the unit normal vector to the
boundary by nx and ny . Show that
Tx =Tn nx − Tt ny
Ty =Tn ny + Tt nx
where Tn (resp. Tt ) is the normal (resp. tangential) component of the traction
vector.
Axisymmetric elastostatic models
Axial symmetry exists when the solution domain can be generated by sweeping
a plane figure around an axis, known as the axis of symmetry, and the material
properties, loading and constraints are axially symmetric. For example, pipes,
cylindrical and spherical pressure vessels are often idealized in this way.
The radial, circumferential and axial coordinates are denoted by r, θ and z
respectively and the displacement, stress, strain and traction components are
labeled with corresponding subscripts. The problem is formulated in terms of
the displacement vector components ur (r, z) and uz (r, z).
1. The linear strain-displacement relationships are [60]:
∂ur
∂r
ur
²θ :=
r
∂uz
²z :=
∂z
µ
¶
γrz
1 ∂ur
∂uz
=
:=
+
·
2
2 ∂z
∂r
(3.66)
²r :=
²rz
2. Stress-strain relationships. For isotropic materials the
relationship is obtained from eq. (3.40):
  
 
σr 
λ + 2G
λ
λ
0 



 ²r 

  
  ²θ  EαT∆
σθ
λ
λ
+
2G
λ
0


=
−
σz   λ
λ
λ + 2G 0  
²z  1 − 2ν



 

 

0
0
0
G
τrz
γrz
(3.67)
(3.68)
(3.69)
stress-strain
 
1


 

1
(3.70)
1


 
0
3. Equilibrium. The elastostatic equations of equilibrium are [60]:
1 ∂(rσr ) ∂τrz
σθ
+
−
+ Fr = 0
r ∂r
∂z
r
1 ∂(rτrz ) ∂σz
+
+ Fz = 0.
r ∂r
∂z
(3.71)
(3.72)
Exercise 3.4.11 Write down the Navier equations for the axisymmetric model.
3.5. INCOMPRESSIBLE ELASTIC MATERIALS
3.5
101
Incompressible elastic materials
When ν → 1/2 then λ → ∞ therefore the relationship represented by eq. (3.40)
breaks down. Referring to eq. (3.39), the sum of normal strain components is
related to the sum of normal stress components by:
²kk =
1 − 2ν
σkk + 3αT∆ .
E
The sum of normal strain components is called the volumetric strain and will
be denoted by ²vol := ²kk . Defining:
σ0 :=
1
σkk
3
we have
3(1 − 2ν)
σ0 + 3αT∆
(3.73)
E
where the first term on the right is the mechanical strain, the second term is the
thermal strain. For incompressible materials, that is when ν = 1/2, ²vol = 3αT∆
is independent of σ0 . Therefore σ0 cannot be computed using Hooke’s law.
Substituting eq. (3.73) into eq. (3.40) and letting ν = 1/2 we have:
²vol ≡ ²kk =
σij = σ0 δij +
2E
(²ij − αT∆ δij ).
3
(3.74)
Substituting into eq. (3.48), and assuming that E and α are constant;
(σ0 ),i +
2E
(²ij,j − α(T∆ ),i ) + Fi = 0.
3
Writing
1
1
(ui,j + uj,i ),j = (ui,jj + uj,ij )
2
2
and interchanging the order of differentiation in the second term, we have:
²ij,j =
uj,ij = uj,ji = (uj,j ),i = (²jj ),i = 3α(T∆ ),i .
Therefore, for incompressible materials,
²ij,j =
1
3
ui,jj + α(T∆ ),i .
2
2
The problem is to determine ui such that
(σ0 ),i +
E
(ui,jj + α(T∆ ),i ) + Fi = 0
3
(3.75)
subject to the condition of incompressibility, that is, the condition that volumetric strain can be caused by temperature change but not by mechanical stress,
ui,i = 3αT∆
(3.76)
102
CHAPTER 3. LINEAR MODELS
and the appropriate boundary conditions. If displacement boundary conditions
are prescribed over the entire boundary (∂Ωu = ∂Ω) then the problem is solvable
only if the prescribed displacements are consistent with the incompressibility
condition:
Z
Z
ui ni dS = 3
αT∆ dV.
∂Ω
Ω
This follows directly from integrating eq. (3.76) and applying the divergence
theorem.
Exercise 3.5.1 An incompressible bar of constant cross section and length `
is subjected uniform temperature change (i.e., T∆ is constant). The centroidal
axis of the bar is coincident with the x1 axis. The boundary conditions are:
u1 (0) = u1 (`) = 0. The body force vector is zero. Explain how (3.74) and
(3.76) are applied in this case to find that σ11 = −EαT∆ .
Exercise 3.5.2 Write the equations of incompressible elasticity in unabridged
notation.
3.6
Stokes’ flow
The flow of viscous fluids at very low Reynolds numbers14 (Re < 1) is modeled
by the Stokes15 equations. There is a very close analogy between the Stokes
equations and the equations of incompressible elasticity discussed in Sec. 3.5.
In fluid mechanics the average normal stress is the pressure: p := −σ0 . The
vector ui represents the components of the velocity vector and the shear modulus
of the incompressible elastic solid E/3 is replaced by the coefficient of dynamic
viscosity µ (measured in Ns/m2 units):
µui,jj =p,i − Fi
ui,i =0.
(3.77)
(3.78)
Exercise 3.6.1 Write the Stokes equations in unabridged notation.
Exercise 3.6.2 Assume that velocities are prescribed over the entire boundary
for a Stokes problem (i.e., ∂Ωu = ∂Ω). What condition must be satisfied by the
prescribed velocities?
3.7
Chapter summary
The formulation of linear mathematical models for problems in heat conduction
and elasticity was described. A mathematical model was understood to be one
or more partial differential equations and the prescribed initial and boundary
conditions. The notion of mathematical model will be generalized in Chapter 4.
14 Osborne
15 George
Reynolds 1842-1912.
Gabriel Stokes 1819-1903.
3.7. CHAPTER SUMMARY
103
The model of heat conduction is based on the conservation law, a fundamental law of physics, and some empirical relationship between the derivatives of
the temperature u and the flux vector which are subject to certain limitations.
For example, the coefficients of thermal conductivity generally depend on the
value of u as well as the gradient of u. Only within a certain range of u and its
gradient can these coefficients be treated as constants.
Similarly, the model for elastic bodies is based on the conservation of momentum (in statics the equations of equilibrium) and an empirical linear relationship between the stress and strain tensors. The linear relationship holds for
small strains only. It is important to bear the limitations of a particular mathematical model in mind when considering it in the context of an engineering
decision-making process. A mathematical model must never be confused with
the physical reality that it was conceived to simulate.
Dimensional reduction should be understood as a special case of the fully
three-dimensional formulation.
It was noted that more than one physical phenomena can be represented by
a partial differential equation. This allows generic mathematical treatment of
various physical phenomena.
Fly UP