********************************************************************
* *
* The equations for displacements, stresses and strains in solids *
* *
********************************************************************
by: D B Spalding (latest update 04.04.06)
Contents
--------
1. Nomenclature
2. The differential equations for displacements
3. Boundary conditions for displacements
3.1 Forces and stresses at solid-fluid surfaces
3.2 Boundary conditions involving displacement gradients
3.3 Discussion
4. The rotation equations
4.1 Differential equations
4.2 Boundary conditions
5. Differential equations for displacements in terms of rotations
5.1 Mathematical manipulations
5.2 The resulting new equations
5.2.1 Equations in terms of rotation gradients
5.2.2 Expressing the rotation-gradients via displacement gradients
5.2.3 Introducing the shear stress
5.2.4 The complete differential equation
5.2.5 Discussion
5.2.6 Finite-volume expression
5.4 One- and two-dimensional situations
5.4.1 The displacement gradients at constrained boundaries
5.4.2 Gradients of ex, ey and ez at boundaries
6. Implementation in PHOENICS
6.1 Comparison of the differential equations
6.2 The material-dependent properties
6.3 The dilatation
6.4 Applied-force boundary conditions
6.5 Thermal-expansion boundary conditions
6.6 Rotation-gradient sources
7. Special cases
7.1 Uniform direct stress
7.2 Uniform thermal stress
7.3 Non-uniform direct stress
7.4 Non-uniform thermal stress
7.5 Uniform bending stress
7.6 Non-uniform bending stress
8. Solutions in terms of 'stress-functions'
8.1 Introduction
8.2 The torsion of straight bars of non-uniform cross-section
8.3 The bending of beams of non-uniform cross-section
---------------------------------------------------------------------
1. Nomenclature
---------------
* Literature reference:
SP Timoshenko and JN Goodier,
Theory of Elasticity (3rd edition)
McGraw Hill, 1970: T&G
BA Boley and JH Weiner
Theory of Thermal Stresses
John Wiley, 1960 B&W
* Cartesian coordinate directions: x, y, z
* Radius, axial distance and angle in
cylindrical-polar cordinates : r, a, z
* Strains (extensions) in x, y, and z directions: ex, ey, ez
* Shear-strain components in the three directions: exy, eyz, ezx
* Linear thermal expansion: et
* Normal stresses in x, y, z directions sx, sy, sz
* Shear stresses in xy, yz, zx planes sxy, syz, szx
* Displacements in x, y and z directions: u, v, w
* Poisson's ratio: P
* Dilatation D
* Young's modulus E
* Lame's constant lamda = E*P/[(1+P)*(1-2*P)] L
* Lame's constant mu = E/[2*(1+P)] G
These definitions imply that the expressions
L+G, 2*G and L+2*G , which appear in equations
below, have the significances:
L+G = E / [ 2 * (1+P) * (1-2P) ]
2*G = E / (1 + P )
L+2*G = E * (1-P) / [ (1+P) * (1-2P) ]
* Constant needed for expressing thermal-expansion
effects alpha * E / (1-2*P) H
where alpha is the longitudinal
thermal-expansion coefficient
* Externally-applied forces in x, y and Z dirctns X, Y, Z
* The temperature excess above the initial state T
which entails that et = alpha*T
* Differentiation with respect to x, y or z is indicated by
C:\WINDOWS\DESKTOP\ONLINE~1a next-line x, y or z
so that du/dx is written as u
x
likwise du/dy is written as u
y
likewise du/dz is written as u
z
* Thus: the dilatation, e, becomes e = u + v + w
x y z
* the gradients u , v and w , respectively,
x y z
are equal also to: ex , ey and ez by definition;
therefore the dilatation can also be written as:
e = ex + ey + ez
* the definitions of the shear-strain components (also called
'angular distortions') can likewise be written as:
exy = u + v ; eyz = v + w ; ezx = u + w
y x z y z x
* the definitions of the rotations and torques become:
x-wise rotation i = v - w
z y
y-wise rotation j = w - u
x z
z-wise rotation k = u - v
y x
********************************************************************
* *
* Important note: *
* Strictly speaking, i, j and k represent TWICE the rotations, as *
* properly derived by T&Goodier on page 233 and given *
* by them the symbols omega_x, omega_y and omega_z . *
* *
********************************************************************
x-wise torque I = Y - Z
z y
y-wise torque J = Z - X
x z
z-wise torque K = X - Y
y x
the torques being per unit volume
* div_grad is represented by a next-line ^,
so that div_grad u becomes u
^
2. The differential equations for displacements
-----------------------------------------------
From T&G, pages 241 and 457,
G*u + (L+G)*e - H*T + X = 0 [2.1]
^ x x
G*v + (L+G)*e - H*T + Y = 0 [2.2]
^ y y
G*w + (L+G)*e - H*T + Z = 0 [2.3]
^ z z
The corresponding normal stresses, from page 456 of the same
book, are:
sx = L*e + 2*G*u - H*T [2.4]
x
sy = L*e + 2*G*v - H*T [2.5]
y
sz = L*e + 2*G*w - H*T [2.6]
z
The corresponding shear stress, from pages 7 and 10 of the same
book, are:
sxy = G*exy = G*( u + v ) [2.7]
y x
syz = G*eyz = G*( v + w ) [2.8]
z y
szx = G*ezx = G*( w + u ) [2.9]
x z
with temperature gradients making no appearance because they
can affect rotations, e.g. u - v ,
y x
but can not afffect 'shear strains' (also called 'angular distortions')
such as: e.g. u + v .
y x
3. Boundary conditions for displacements
----------------------------------------
3.1 Forces and stresses at solid-fluid surfaces
-----------------------------------------------
At surfaces separating solid from fluids, the forces per unit area
in the cartesian directions x, y and z, namely X", Y" and Z", are
related to the stresses in the solid by the following equations,
given by T&G, p.236 :
X" = sx * l + sxy * m + szx * n [3.1-1]
Y" = sy * m + syz * n + sxy * l [3.1-2]
Z" = sz * n + szx * l + syz * m [3.1-3]
Here l, m and n are the direction cosines of the outward normal to
the surface, defined by:
l = x [3.1-4]
s
m = y [3.1-5]
s
n = z [3.1-6]
s
where s stands for the distance along the normal.
For surfaces aligned with the cell faces of a cartesian grid, the
values of the direction cosines are:
* for East faces: l = 1 ; m = 0 ; n = 0
* for North faces: l = 0 ; m = 1 ; n = 0
* for High faces: l = 0 ; m = 0 ; n = 1
3.2 Boundary conditions involving displacement gradients
--------------------------------------------------------
Substitution in 3.1-1, 3.1-2 and 3.1-3 of the expressions for
sx, sy, sz, sxy, syz, and szx of section 2 enables boundary
conditions for the displacements u, v and w to be written in terms
of their gradients and the external forces per unit area. This is
will be done for faces aligned with grid-cell surfaces, as follows:
* For East cell faces (l=1, m=0, n=0):
u = ( X" - L * e + H * T ) / (2 * G ) [3.2-1]
x
v = - u + Y" / G [3.2-2]
x y
w = - u + Z" / G [3.2-3]
x z
* For North cell faces (l=0, m=1, n=0):
v = ( Y" - L * e + H * T ) / (2 * G ) [3.2-4]
y
u = - v + X" / G [3.2-5]
y x
w = - v + Z" / G [3.2-6]
y z
* For High cell faces (l=0, m=0, n=1):
w = ( Z" - L * e + H * T ) / (2 * G ) [3.2-7]
z
u = - w + X" / G [3.2-8]
z x
v = - w + Y" / G [3.2-9]
z y
Each of the gradients on the left-hand side of these equations has,
when multiplied by (L+G), the significance of a source in the
displacement-balance equation.
3.3 Discussion
--------------
It is worth pointing out to fluid-dynamicists that the
just-mentioned sources are partly similar to and partly different
from those which apply to velocities at solid-fluid boundaries.
Thus, if the equations for u are considered:
(1) The X"s in each of the equations have their counterparts in the
corresponding equations for x-direction velocity; for x-directed
forces, whether applied normally or tangentially to the surface,
surely act as sources of x-directed momentum.
(2) The H * T terms do not appear in the equations for velocity;
nevertheless, it is understandable that they appear in precisely the
same form for each displacement, because thermal expansion has no
preferred direction.
(3) However, the L*e term represents a first big surprise for the
fluid-dynamicist; for the dilatation D contains all three
displacements; therefore gradients of v and w can provide sources of
u. Fluid-dynamics has not counterpart to this effect.
(4) The surprise may be somewhat allayed by the recognition that the
material property L is zero when the Poisson's Ratio is zero; so the
L*e term can be regarded as a Poisson's-Ratio effect, which
is a peculiarity of solids from which fluids are mercifully free.
(5) No such excuse alleviates the surprise occasioned by the
appearance of gradients v and y in the equations 3.2-5 and 3.2.8;
for they express an effect of the resistance to shear of solid
bodies which a fluid dynamicist might have thought as already
handled by the appearance of viscosity-like terms in the
differential.
(6) Finally, surprising or not, the unexpected terms MUST be
included if correct solutiions are to be obtained.
4. The rotation equations
-------------------------
4.1 Differential equations
--------------------------
It will be shown to be useful, for the solution of solid-stress
problems by means of PHOENICS, to make use of differential equations
for the rotations i, j and k, defined in section 1.
Differentiation of equation 2.1 with respect to y, and of equation
2.2 with respect to x, followed by subtraction, leads to:
G*k + K = 0.0 [4.1-1]
^
where k and K are respectively the rotation and the externally-
applied torque defined above.
Corresponding operations on the pairs of equations 2.2 with 2.3,
and 2.3 with 2.1, lead to two further equations, namely:
G*i + I = 0.0 [4.1-2]
^
and
G*j + J = 0.0 [4.1-3]
^
Here it is notable that the terms containing D and T make no
appearance, by reason of the identities:
D = D and T = T .
xy yx xy yx
These equations do not appear in T&G; but it is hardly possible that
they have escaped attention by specialists in the theory of
elasticity.
4.2 Boundary conditions
-----------------------
There being no guidance from T&G, the boundary conditions will be
derived from physics-based reasoning, attention being given to the
boundary conditions the k on north and east surfaces of a body of
which the cross-section in the x and y directions is indepedent of
the third dimension, z.
(1) As a first case, let the body be a beam with uniform thickness in
the y direction much smaller than its length in the x-direction; and
let the only external forces be equal and opposite couples at the
low-x and high-x ends, only the latter being free to rotate.
In these circumstances, the beam will bend in such a manner as to
have uniform curvature; therefore rotation k will be linear in x,
thus:
k = k1 * x/x1
where the length of the beam is x1 and the rotation at x=x1 equals
k1.
A linear k distribution satisfies the differential equation with
no sources at the north and south. Since there are no body forces
there either, the not unreasonable idea that "no forces means no
sources" gains some support.
It should be noted that the y-direction displacement is parabolic in
x, e.g. :
v = a*x**2
(2) Let us next consider the case in which the beam is loaded by a
tranverse y-direction force at the x=x1 end, which is known from
elementary beam theory to have as solution:
v = a*x**2 + b*x**3
This corresponds to a parabolic k distribution; and this satisfies
the differential equation for k only if there IS a finite source of
k at the north and south boundaries, on which however there are
still no external forces.
This contradicts the "no force means so source" hypothesis; and it
presents a dilemma: by calculating the moment exerted by the force
at x=x1 we can deduce what is the bending moment a the x=x
cross-section; but this implies some "action at a distance" notion
which is not in keeping with the approach via differential
equations; and it gives no guidance as to what is to be done when
the body is not beam-like, i.e. long and thin.
A way out of the dilemma may be the following:
- admit that there ARE force-like actions at the north and south
boundaries in BOTH cases;
- recognise these as the "pseudo forces" which enter the
u-displacement equation via equation 3.2-5 above;
????????? to be continued ???????
5. Differential equations for displacements in terms of rotations
-----------------------------------------------------------------
5.1 Mathematical manipulations
------------------------------
Gradients of the dilatation D appear in the equations 2.1, 2.2
and 2.3 . These can be written out in full as:
D = u + v + w [5.1-1]
x xx xy xz
D = u + v + w [5.1-2]
y xy yy yz
D = u + v + w [5.1-3]
z xz yz zz
It is now useful to differentiate the definitions of the rotation
components:
i with respect to x and y,
j with respect to y and z, and
k with respect to z and x, with the following results:
k = u - v [5.1-4]
x xy xx
k = u - v [5.1-5]
y yy xy
i = v - w [5.1-6]
z zz yz
i = v - w [5.1-7]
y yz yy
j = w - u [5.1-8]
x xx xz
j = w - u [5.1-9]
z xz zz
The following substitutions effect the desired simplifications:
(a) into 5.1 from 5.5 and 5.9, to yield:
D = u + u + u - k - j
x xx yy zz y z
i.e.
D = u + j - k [5.1-10]
x ^ z y
(b) into 5.2 from 5.4 and 5.6, to yield:
D = v + k - i [5.1-11]
y ^ x z
(c) into 5.3 from 5.7 and 5.8, to yield:
D = w + i - j [5.1-12]
z ^ y x
5.2 The resulting new differential equations
--------------------------------------------
5.2.1 Equations in terms of rotation gradients
----------------------------------------------
Substitution for dilatation gradients from 5.10, 5.11 and 5.2 into
the differential equations for the displacements yields three new
equations, namely:
-----------------------------------------------------------------
| |
| (L+2*G)*u + (L+G)*(j - k) - H*T + X = 0 [5.2-1] |
| ^ z y x |
| |
| (L+2*G)*v + (L+G)*(k - i) - H*T + Y = 0 [5.2-2] |
| ^ x z y |
| |
| (L+2*G)*w + (L+G)*(i - j) - H*T + Z = 0 [5.2-3] |
| ^ y x z |
| |
-----------------------------------------------------------------
These equations are easier to solve than these of section 2, from
which they were derived, because:
* the terms involving the rotation gradients can be evaluated
separately; and indeed they can often be represented by
analytical expressions in terms of x, y and z; and
* the linkages between u, v and w, formerly expressed by the
dilatation gradients, have disappeared.
It is true that some linkages remain because D appears in the
boundary-conditions listed in section 3; but these are appreciably
weaker.
5.2.2 Expressing the rotation-gradients via displacement gradients
------------------------------------------------------------------
The definitions of j and k imply that:
(j - k) = w - u - u + v
z y zx zz yy yx
which can be written as:
(j - k) = w - u - u + v
z y xz zz yy xy
i.e. as:
(j - k) = ez - u - u + ey
z y x zz yy x
i.e. as:
(j - k) = ez + ey - u - u [5.2-4]
z y x x zz yy
Similar substitutions lead to:
(k - i) = ex + ez - v - v [5.2-5]
x z y y xx zz
and:
(i - j) = ex + ey - w - w [5.2-6]
y x z z xx yy
5.2.3 Introducing the shear stress
----------------------------------
The shear stresses sxy, syz and sxz are related to terms in equations
2.2-4, 2.2-5 and 2.2-6 by the following relations:
sxy/G = u + v [2.7]
y x
syz/G = v + w [2.8]
z y
sxz/G = w + u [2.9]
x z
With the aid of these equations the terms u, u, v, v, w, w
yy zz xx zz yy xx
can be eliminated, with the results:
(j - k) = 2*(ez + ey) - sxy/G - sxz/G [5.2-7]
z y x x y z
(k - i) = 2*(ex + ez) - sxy/G - syz/G [5.2-8]
x z y y x z
(i - j) = 2*(ex + ey) - sxz/g - syz/G [5.2-9]
y x z z x y
These equations are useful only where the shear stresses are known,
for example at solid-fluid boundaries.
5.2.4 The complete differential equation
----------------------------------------
It will be seen that equation 5.2-4 contains the terms: - u - u .
yy xx
Since those two terms appear with positive signs in the u term of
^
the differential equation 5.2-1, it apppears appropriate it cancel
them out. The complete differential equation can now be written as:
(L+2*G)*u + G*( u + u) + (L+G)*(ey + ez) - H*T + X = 0 [5.2-10]
xx yy zz x x x
with corresponding equations for v and w, namely:
(L+2*G)*v + G*( v + v) + (L+G)*(ex + ez) - H*T + Y = 0 [5.2-11]
yy xx zz y y y
(L+2*G)*w + G*( w + w) + (L+G)*(ey + ez) - H*T + Z = 0 [5.2-12]
zz xx zz z z z
Let it now be observed that, the definitions of the relevant
quantities imply that H, defined above as: alpha * E / (1-2*P) ,
can also be expressed as: alpha * ( 3*L + 2*G).
With alpha abbreviated to A, the roles of L and G in the equations
can be better clarified by re-writing them (omitting the *, for
brevity, except between A and T) as:
(L+2G)u + G( u + u) + (L+G)(ey + ez) - (3L+2G)A*T + X = 0 [5.2-13]
xx yy zz x x x
(L+2G)v + G( v + v) + (L+G)(ex + ez) - (3L+2G)A*T + Y = 0 [5.2-14]
yy xx zz y y y
(L+2G)w + G( w + w) + (L+G)(ey + ez) - (3L+2G)A*T + Z = 0 [5.2-15]
zz xx zz z z z
5.2.5 Discussion
----------------
It must now be admitted that these equations could have been derived
directly from those section 2, and that the introduction of rotation
has proved to be a needless digression.
Thus, the first term of 5.2-13 can be split into two, thus:
(L+2G)u = (L+G)u + Gu
xx xx xx
Therefore the first three terms can be written as:
(L+G)u + Gu + G( u + u) + (L+G)(ey + ez)
xx xx yy zz x x
i.e. as:
G( u + u + u) + (L+G)( u + ey + ez)
xx yy zz xx x x
i.e. as:
G u + (L+G)e
^ x
These are precisely the first two terms of equation 2.1 .
Which form of equation is likely to prove most convenient for
numerical solution? Probably that of section 5.2.4 because, in a
segregated-variables solution procedure, it is possible to express
all the terms containing dependent variable (u, v or w) by way of
the finite-volume 'molecule', so that only those containing the
strain gradients need to be inserted explicitly as source terms.
Let it however be remembered that the introduction of the strain
gradients into the above expresions is not an obligatory move; for
ex can also be expressed as u .
y xy
In the following section it will appear that the second of these
has some advantages.
5.2.6 Finite-volume expression
------------------------------
Whatever the form of the differential equation and boundary
conditions, it is how they are expressed in finite-voume (i.e. FV)
form which is decisive for the satisfactory working of the computer
program. Attention will now be given to the FV equation for the
displacement u, the grid being cartesian.
The following diagram represents the grid in two dimensions; but the
existence of the z (High-Low) direction must also be borne in mind.
! ! ! ! ! !
! ! ! N ! ! !
! ! ! --> ! ! !
! ^ ! !/////////! ! !
! ! ! !/////////! ! !
!--dyn----!----k----^----kn---^----k----!---^-----!------------------
! ! ! NW/////////NE ! ! !
! ! ! W !/// P ///! E ! ! !
! v ! --> ! --> ! --> ! dyp !
! ^ ! !\\\\\\\\\! ! ! !
! ! ! !\\\\\\\\\! ! ! !
!--dys----!----k----^----ks---^----k----!---v-----!------------------
! ! ! SW\\\\\\\\\SE ! !
! ! ! !\\\ S \\\! ! !
! v ! ! --> ! ! !
! ! !<--dxp-->! ! !
! ! <--dxw--><-- dxe--> ! !
Consider the terms of equation 5.2-1 which concern displacements,
namely:
(L+2*G)*u + (L+G)*(j - k)
^ z y
Integration of the first term over the control volume surrounding P
leads ( with obvious notation) to:
(L+2G)*{ dyp*dzp*[( uE -uP) / dxe - ( uP - uW) / dxw]
+ dxp*dzp*[( uN -uP) / dyn - ( uP - uS) / dys]
+ dxp*dyp*[( uH -uP) / dyh - ( uP - uL) / dzl] }
This is a typical diffusion-type expression, requiring no comment.
Consider now the second part of the second term, involving the
gradient of the z-direction rotation component, k, the integration
of which over the control volume surrounding P yields:
- (L+G)*dxp*dzp*(kn - ks)
where kn and ks are the values of k in the control volumes shaded
////// and \\\\\\ respectively.
These quantities can be represented in discretized form as:
kn = ( uN - uP ) / dyn - ( vNE - vNW) / dxp , and
ks = ( uP - uS ) / dys - ( vSE - vSW) / dxp
Corresponding expressions can be derived for the j-gradient term,
involving uH, uL, wHE, wHW, wLE, wLW, dzh and dzl. However the
necessary conclusions can be drawn by confining attention to the
xy-plane.
All the above terms involve displacement differences. Their combined
influence can therefore be represented as:
( uE - uP ) * [ (L+2G)* dyp*dzp / dxe ]
+ ( uW - uP ) * [ (L+2G)* dyp*dzp / dxw ]
+ ( uN - uP ) * [ G* dxp*dzp / dyn ]
+ ( uS - uP ) * [ G dyp*dzp / dys ]
+ ( vNE - vNW) * [ (L+G) * dzp ]
- ( vSE - vSW) * [ (L+G) * dzp ]
with corresponding z-direction terms.
Consider now the situation in which the point P is located at the
solid-fluid interface, as shown below:
! !
! ! N
! ! -->
^ ! !/////.
! ! !/////.
dyn----!----k----^----kn ^
! NNW NW/////. !
! ! W !/// P. !
v ! --> ! --> dyp
^ ! !\\\\\. !
! ! !\\\\\. !
dys----!----k----^----ks- v-
! SSW SW\\\\\.
! ! !\\\ S.
v ! ! -->
! !<--dxp-->!
The rotations kn and ks cannot be computed as before, because
there is no direct way of calculating dv/dx. This quantity could
be computed indirectly from the shear-stress boundary condition 2.7
thus:
sxy = G*exy = G*( u + v ) [2.7]
y x
kn = 2*dy/dy - sxy
whence:
kn = 2*(uN - uP) / dyn - sxyn/G
ks = 2*(uP - uS) / dys - sxys/G
However, there does not appear to be any advantage in so doing.
Instead, it appears preferable to presume that dv/dx varies little
with x, so that its value can be extrapolated from the next cells to
the west. The variation of dv/dx with y is what is important.
5.3 Cylindrical polar coordinates
---------------------------------
5.3.1 Equations for displacements in terms of polar co-ordinates
----------------------------------------------------------------
According to B&W, p. 257, the equations for the displacements in
terms of gradients of rotation can be written in polar
co-ordinates as follows:
(L+2G)D - (L+G)(j - k) - (3L+2G)A*T + X = 0 [5.3-1]
rx z y rx
(L+2G)D - (L+G)( k - i) - (3L+2G)A*T + Y = 0 [5.2-2]
y rx z y
(L+2G)D - (L+G)(ri - j) - (3L+2G)A*T + Z = 0 [5.2-3]
z ry rx z
wherein:
* the dilatation, D = rv + u + w [5.2-4]
ry rx z
* the rotation about the z-axis, k = ru - v [5.2-5]
ry rx
* the rotation about the x-axis, i = v - w [5.2-6]
z y
* the rotation about the y-axis, j = w - y [5.2-7]
rx y
and, x is now measured in radians;
r stands for radial distance, which is identical to y,
so that dr=dy;
where r appears it has the significance of a multiplier
whereas the appearance of y on the second line signifies
differentiation, so that:
rv stands for [d(rv)/dy]/r
ry
u stands for [d(u)/dx]/r
rx
ru stands for [d(ru)/dy]/r
ry
v stands for [d(v)/dx]/r
rx
and so on.
5.4 One- and two-dimensional situations
---------------------------------------
5.4.0 Definitions and preliminary considerations
------------------------------------------------
(a) Definitions and examples
----------------------------
It frequently arises that the shape of the solid in question, and
the directions and distributions of displacement-producing
influences, whether mechanical or thermal, result in distributions
of displacement and temperature which vary in one direction only.
[ For brevity, the words 'width', 'thickness' and length will be
associated with the cartesian coordinates x, y and z respectively.
In axi-symmetrical situations, for which polar coordinates may be
used, the radial direction will always be y; and x will represent
angle around the z-axis. ]
Examples of one-dimensional distributions are:
A. A long wire of uniform diameter is pulled by equal and opposite
forces, one at each end. Its longitudinal-direction strain and
stress are uniform along its length; but the longitudinal
displacement varies linearly with longitudinal position (z).
B. A uniform-sectioned block is compressed by equal and opposite
forces applied uniformly to high-x and low-x faces. All
displacements, stresses and strains depend on the x coordinate alone;
and this is true whether its thickness and/or length are large or
small compared with its width.
C. A similar block is placed suddenly into contact at its high-x and
low-x faces with objects of uniform but unequal temperature. The
temperature of the material within the block changes with time, as a
consequence of heat conduction; but the temperature at any particular
time depends only on x.
Examples of two-dimensional distributions are:
d. The forces on the high-x and low-x faces of a block such as that
in example B. vary with its position in the width dimension, y. Then
the stresses and displacements vary with both x and y, but not with
the third (length) coordinate, z.
e. In examples B. and C., the width varies with x, but not the
length. Then the stresses, displacements and temperatures vary
with both x and y, but not with z.
f. The block of example A. is clamped firmly at the low-x face and
subjected at the other to a force which is uniform over that face
but which is directed in the y direction. The block therefore acts
as a cantilever beam, subject to bending stresses.
(b) A difference between the solid-stress and fluid-flow equations
------------------------------------------------------------------
Fluid-flow specialists are familiar with the fact that, although all
real flow phenomena are three-dimensional, much can be gained from
the study of one- and two-dimensional idealisations; however they
enjoy an advantage which the solid-stress specialist does not,
namely:
* the dimensions which are NOT the 'one', or one of the 'two', need
NOT BE CONSIDERED AT ALL.
Thus, if a fluid flows in such a manner that its properties vary
with x and y but not with z, the terms which are connected with z,
for example: u or u , are all zero.
z zz
Where displacements, strains and stresses in solids are concerned,
however, this is not so. Thus, as will be explained below (see
section 5.4. ), u may be zero without u being zero also.
v vv
Why is this? In part it is due to the phenomenon expressed by the
non-zero value of Poisson's ratio, namely: stresses in the
x-direction normally cause either strains or stresses in the y-
and z-directions as well as strains in the y-direction.
However, this is not the whole of the matter; for inter-connexions
between the three directions are implied by the fact that gradients
of the dilatation D, which involve gradients of u, v and w, appear
in the differential equations for equations each one of these
displacements.
It seems therefore wisest to warn the CFD specialist that, although
the equations governing displacements are superficially similar to
those governing velocities, there exist subtle differences which
those seeking to solve them must respect.
(c) Consequences for the equation-solving procedures
----------------------------------------------------
(i) The solid-stress-analysis literature pays little attention to
one-dimensional situations; and, when it does so, it is usually
postulated that the material is unconstrained in the other two
directions, as in the wire under tension of case A. above.
(ii) By contrast, special nomenclature is employed for two-dimensional
situations, attention being focussed on two extremes, namely, if the
x-y plane is that in which variations are to be determined:
* In extreme (1), the material is free to expand or contract in the
z-direction. This is called the 'PLANE-STRESS' situation, because
direct stresses in the z-direction are negligible.
It arises, typically, when the z-direction dimensions of the solid
in question are small compared with those in the other two
directions. However, this is far from being sufficient; for it is
possible that such a body, however thin, could be confined between
two totally inflexible planes, which prevented their z-direction
movement.
For the 'plane-stress' condition to prevail, such confining surfaces
must NOT exist.
* In extreme (2), the material is NOT free to expand or contract in
the z-direction, being constrained, it is supposed by rigid surfaces
to which it must adhere, while still being free to slide past it in
the x or y directions.
This is called the 'PLANE-STRAIN' situation, because strains in the
z-direction are negligible.
This situation arises when inflexible confining surfaces of the
just-discussed kind DO exist; but it arises also when the
z-direction dimensions of the solid are large compared with those in
the x- and y-directions; for then, if significant variations with z
did exist, they would be opposed and nullified with the
consequential shear stresses on planes of constant x and y.
(iii) For one-dimensional situations, the fact that the solid-stress
literature appears to have neglected them does not mean that they
SHOULD be ignored. It appears therefore appropriate to define two
additional situations, namely:
(3) LINEAR STRESS, in which (if z is the longitudinal direction)
stresses in the x- and y-directions are neglible. The 'wire' of
example A. illustrates this.
(4) LINEAR STRAIN, in which the solid is enclosed in a totally rigid
casing which prevents its boundaries from moving in (say) the y-
and z-directions.
(d) Particular consequences for the numerical analyst
-----------------------------------------------------
The fluid-dynamicist, when dealing with a one-dimensional situation,
will discretize space in the appropriate (say, x) direction; and the
only solved-for velocity will then be u. Velocities in the other
directions, namely v and w, will make no appearance.
Similarly, if the problem is two-dimensional (with variations in x
and y, say), the u and v velocities will be solved for; but not the
z-direction velocity, w.
If a computer code designed for computing velocities is to be employed
(as PHOENICS is) for computing displacements, the differences beteen
the equation systems, mentioned in (b) above, must be addressed; and
actions may be needed in three areas, namely:
(1) the calculation of the dilatation;
(2) the calculation of the forces at boundaries; and
(3) (when thermal effects are significant) by provision of sources
of displacement in the solved-for directions.
5.4.1 The displacement gradients at constrained boundaries
----------------------------------------------------------
Let it be supposed that a solid rectangular block is subjected to
uniform pressure on its low-x and high-x faces. The problem might be
regarded as being one-dimensional, with only equation 5.13 requiring
to be solved.
This is true; but the y and z directions cannot be completely
disregarded for the reason that the x-direction behaviour depends on
whether the low-y and high-y, and/or the low-z and high-z, faces are
held fixed, allowed to move freely, or partially constrained.
The dependence of the x-direction displacements on the y- and
z-direction constraints is expressed by way of the dilatation, i.e.
the sum of the strains ex, ey and ez; therefore ey and ez must be
given values.
Two extreme cases are common, namely:
(a) the strains are zero because the surfaces are held fixed; and
(b) the stresses are zero, because the surfaces are unconstrained.
However, although neglected by most textbooks, a third situation may
arise, namely:
(c) that of linear constraint, whereby the displacement may be
proportional to the force exerted at the boundary.
Since case(c) contains the other two as special cases, with
proportionality constant zero for case (a) and infinite for case(b),
it is this which will be considered.
Let the forces exerted on the high-x, high-y and high-z boundaries
be respectively: - bx*ex , - by*ey and - bz*ez , where ex, ey
and ez are the displacements at the constrained boundaries. Then it
follows from equations 3.1, 3.2 and 3.3, that:
ex = ( -bx*ex - L*e + H*T ) / (2*G)
i.e.
ex = ( - L*e + H*T) / ( 2*G + bx ) [5.16]
ey = ( -by*ey - L*e + H*T ) / (2*G)
i.e.
ez = ( - L*e + H*T) / ( 2*G + by ) [5.17]
ez = ( -bz*ez - L*e + H*T ) / (2*G)
i.e.
ez = ( - L*e + H*T) / ( 2*G + bz ) [5.18]
These are the equations to be used when (and ONLY when) the strain
in question has NOT been deduced from the differential equation for
the corresponding displacement.
5.4.2 Gradients of ex, ey and ez at such boundaries
---------------------------------------------------
(a) The problem
---------------
Equation 5.13 contains the term u , i.e. u + u + u ;
^ xx yy zz
and, in a one-dimensional situation, with x as the direction in
which u is allowed to vary, it might be supposed that the
u and u
yy zz
terms might be neglected. This is not the case however, at
least when thermal expansion is present.
Consider, for example, the case in which the temperature T varies
with x, and the proportional linear thermal expansion, et;
the proportionality coefficient is the material constant, alpha.
Let it also be supposed that the y- and z-direction surfaces are not
constrained in any way, so that the three strains at any point obey
the equality relation:
ex = ey = ez = et [5.19]
Then initially-plane (constant-x) surfaces will become spheres of
radius equal to: alpha * T .
x
Of course, if movement in x is allowed, but not in z, the
initially-plane surfaces will become cylinders of the same radius
about an axis normal to the xy plane; and, if movement in z is
allowed, but not in x, the initially-plane surfaces will become
cylinders of the same radius about an axis normal to the xz plane.
Finally, if movement is allowed in neither the y nor z directions,
plane surfaces will remain plane.
(b) The distribution of u
-------------------------
Consider the simple case in which et is linear in x, for example:
et = a * ( - 1 + 2*x/l) [5.20]
which corresponds to the case of a block of thickness l, which is
cooled to a lower-than-average temperature where x=0 and heated to a
correspondingly higher-than-average temperature at the opposite
face, the linear temperature being brought abour, in the steady
state, by heat conduction.
In this case, alpha * T = 2*a/l [5.21]
x
The uncontrained nature of the process entails that the linear
extensions in each direction, ex, ey and ez are all exactly equal
to et, from which the following expressions concerning first and
second gradients of displacement and rotation can be derived:-
(i) u : Differentiation of u = ex (by definition) = et
xx x
leads, via 5.20, to: u = 2*a/l [5.22]
xx
(ii) u : Integration of v = ey (by definition) = et
yy y
leads, via 5.20, to: v = a*(-1+2*x/l)*y + const [5.23]
whereafter differentiation with respect to x leads to:
v = 2*a*y/l [5.24]
x
This relation is useful because, since there can be no shear
stresses, equation 2.7 leads to:
u + v = 0,
y x
from which it follows, via differentiation of 5.24, that:
u = - v = - 2*a/l [5.25]
yy xy
(iii) u : A similar manipulation of the corresponding terms
zz involving w and z lead to:
u = - 2*a/l [5.26]
zz
(c) Consequences for one- and two-dimensional analyses
------------------------------------------------------
It follows that the supposition mooted at the start of section
5.4.2, namely that the u and u tems might be neglected is
yy zz
far from the truth. If they were neglected, u would equal 2*a/l;
^
but, when they are taken into account, u equals - 2*a/l !
^
(d) The rotation-gradient terms
-------------------------------
Equation 5.13 also contains the term: (L+G)*(j - k) .
z y
In one- or two-dimensional analyses one might be tempted to set
z- and y-direction gradients to zero here also; but that we be
equally wrong, as will now be explained.
The term j , by reason of the definition of j, can be expressed as:
z
j = w - u
z xz zz
Since w = a*(-1 + 2*x/l)z + constant, derived in a simlar manner to
5.23, differentiation leads to:
w = 2*a/l ; and so, by reason of 5.26, we deduce that
xz
j = 4*a/l [5.27]
z
A similar argument shows that k has the same value; and the upshot
y
is that, to account for the effect of the unconsidered y and z
directions, it is necessary to add to the terms in the equation:
-(L+2*G)*(4*a/l) to account for the second differential
coefficients in the first term, and
(L+G)*(8*a/l) to account for the rotation gradients, wherein
it will be recalled that 2*a/l = alpha * T/x .
6. Implementation in PHOENICS
-----------------------------
6.1 Comparison of the differential equations
--------------------------------------------
PHOENICS treats the displacements as though they were velocities;
therefore what is to be done can be deduced by comparing equation
5.13, for example, with the corresponding equation for the
x-direction velocity.
This is, for the steady state, with uniform effective viscosity emu,
and with the convection terms omitted:
emu*u + p + source_terms + X = 0 [6.1]
^ x
Comparison of this equation with 5.13 dictates that:
(a) the emu store should contain L+2*G ;
(b) the pressure-gradient term should be omitted;
(c) (L+G)*(k - i ) should appear as a source term;
z y
(d) so also should - H*T .
x
Similar conclusions can be drawn about the equations for the y and x
directions.
The equations are of course solved in finite-volume form, as far as
possible in the same way for both displacements and velocities. How
this is done will now be discussed.
6.2 The material-dependent properties
-------------------------------------
(a) The PHOENICS PROPS file
---------------------------
PHOENICS provides information about several properties of fluid and
solid materials by way of the so-called PROPS file, which is arranged
in columns, thus:
IMAT density viscosity specific thermal thermal compress
i.e. or Poisson heat conduct expnsn. -ibility or
PRPS 's ratio capacity -ivity coefff. 1/Young's Modulus
The 'or' in columns 3 and 6 refers to whether the material is fluid
or solid. Thus, in the case of solid aluminium, where the entry is:
<100> ALUMINIUM at 27 deg c
100 2700.0 0.35 896.0 204.0 2.35e-5 1.47e-11
the 0.35 is the Poisson's ratio and the 1.47e-11 is the reciprocal
of the Young's Modulus.
(b) The three-dimensionally-stored quantities
---------------------------------------------
When so instructed, PHOENICS creates three-dimensional storage for
each of the above six properties; therefore, unless other
arrangements are made, it will place the viscosity or the Poisson's
ratio in the corresponding storage location according to whether
However, within the solid it is COMPOSITE properties which are
needed, as expressed by the terms L+G, L+2G, H, etc which appear in
the differential equations and boundary conditions. Therefore, in
order to simplify the coding when PHOENICS reads the above-described
PROPS ile, it places, for solid-containing cells:
L+2*G in the viscosity store,
L+G in the compressibility store, and
H in the thermal-expansion coefficient store,
these being computed, in accordance with definitions, from:
L+2*G = E * (1-P) / [ (1+P) * ( 1-2*P)] [6.2]
L+G = E / [ 2 * ( 1+P ) * (1-2*P) ] [6.3]
H = alpha * E / (1-2*P) [6.4]
where alpha is the thermal-expansion coefficient.
These three composite properties are those appearing in the
differential equation of section 5.2. In the boundary conditions of
section 3, the L and G which are needed on their own are easily
deduced from:
G = (L+2*G) - (L+G) [6.5]
L = 2*(L+G) - (L+s*G) [6.6]
When desired, the Poisson's ratio, P, and the thermal-expansion
coefficient, alpha, can be recovered from the stored quantities via:
P = 1.0 - 0.5 * (L+2G) / (L+G) [6.7]
alpha = H / [4.0*(L+G) - (L+2G) ] [6.8]
(c) A device to prevent round-off error.
---------------------------------------
As is exemplified by the above data for aluminium, the reciprocal of
the Young's Modulus is usually of the order of 1.e-11. For this
reason, PHOENICS internally multiplies it by 1.e11, in order that
excessively small coefficients shall not be generated , so inducing
round-off error. Applied forces are correspondingly divided by 1.e11
so as to ensure that displacements, stresses and strains are
correctly reported.
6.3 The dilatation
------------------
Although the dilatation, D, does not appear in the differential
equations which are solved, it does so in the boundary conditions.
Therefore D is calculated within the solver; and its values are
stored in the locations which, for fluid regions, are filled by
values of the pressure.
For reduced-dimensionality situations, of the kind discussed
in section 5.2 above, the following practice is adopted:
* If NX = 1, so that no equation for u is solved, the strain ex is
computed from equation 5.16, the constant bx being supplied by way
of a spedat statement in the q1 file thus:
SPEDAT(BOUNDARY,XCONST,R,value of constant bx)
* If NY = 1, so that no equation for u is solved, the strain ex is
computed from equation 5.17, the constant bx being supplied by way
of a spedat statement in the q1 file thus:
SPEDAT(BOUNDARY,YCONST,R,value of constant by)
* If NX = 1, so that no equation for u is solved, the strain ex is
computed from equation 5.18, the constant bx being supplied by way
of a spedat statement in the q1 file thus:
SPEDAT(BOUNDARY,ZCONST,R,value of constant bz)
Often the values of the constants are either:
* 1.e20, i.e. a large number allowing no displacement,
or
* 0.0, implying no constraint whatever.
6.4 Applied-force boundary conditions
-------------------------------------
As shown in section 3, displacement-gradient boundary conditions
depend on D, T and G, as well as on the applied force. However, if
a force is defined by way of a patch of FIXFLU type, having FORC as
the first four characters of its name, the fourth argument of the
corresponding COVAL can be the force itself; for then the D, T and
G influences will be introduced automatically.
Forces exerted by the fluid on the solid are represented
internally, i.e. without user intervention.
6.6 Rotation-gradient sources
-----------------------------
PHOENICS provides, when necessary, cell-by-cell storage space for
three rotation-gradient sources, namely:
RGRX = (k - i ) [6.9]
z y
RGRY = (i - j ) [6.10]
x z
RGRZ = (j - k ) [6.11]
y x
Their cell-by-cell distribution may often be calculated, once for
all, at the start of the flow simulation. Alternately they may be
computed by solving the differential equations for i, j and k.
7. Special cases
----------------
7.1 Uniform direct stress
-------------------------
a. The problem
--------------
A rectangular parallelepiped, of dimensions lx, ly and lz, is held
at its low-x, low-y and low-z faces and subjected to an applied
tensile force of X" per unit area at its high-x face.
Three cases are considered:
(i) both high-y and high-z faces are free to move
(ii) both high-y and high-z faces are prevented from moving
(iii) the high-y face is free to move but the high-z face is not.
Case(i): motion allowed in y and z directions.
* Application of Hooke's law, and the definitions of Young's modulus
E and the Poisson's ratio P lead directly to:
ex = X"/E [7.1]
ey = P*X"/E [7.2]
ez = P*X"/E [7.3]
from which the dilatation may be deduced as
D = (X"/E) * ( 1 - 2P ) [7.4]
* As a consequence of 7.1, the maximum displacement will equal
ex * lx, i.e. (X"/E)*lx. [7.1a]
* It is a useful check on their derivations to ensure that the above
differential and boundary-condition equations are satisfied by
equations 7.1, 2, 3 and 4.
The differential equations 5.1, 2 and 3 reduce to
u = 0 , v = 0 and w = 0
^ ^ ^
because there are no sources from rotation gradients, temperature
gradients or applied forces within the body.
These equations are all satisfied by uniformity of ex, ey and ez,
which are the gradients of u, v and w, expressed by 7.1, 2 and 3.
The boundary conditions reduce, in the absence of temperature
effects and of Y" and Z", to
ex = ( X" - L*e ) / (2*G) [3.1a]
ey = ( - L*e ) / (2*G) [3.2a]
ez = ( - L*e ) / (2*G) [3.3a]
Summation leads to:
D = ( X" - 3*L*e ) / (2*G)
and so to
D * ( 2*G + 3*L) = X"
The definitions of G and L imply that:
2*G + 3*L = E / ( 1 - 2P )
wherefore
D * E / ( 1-2P) = X"
which is consistent with 7.4 .
Therefore both the differential and boundary-condition equations
have passed their tests.
Case (ii): motion disallowed in y and z directions.
* The solution cannot for this case be obtained by direct
application of Hooke's law, and the definitions of Young's
modulus and the Poisson's ratio; but simultaneous solution of
the general equation leads to the following results:
** The differential equations 5.1, 2 and 3 dictate constancy of
u, v and w gradients because of the absence of source terms;
therefore ex. ey and ez are all constants.
** ey and ez are both zero, by reason of the prevention of movement
in the y and z directions.
** D and ex are therefore identical, because D = ex + ey + ez by
definition.
** Boundary condition 3.1 thus reduces to:
ex = ( X" - L*ex ) / (2*G) [3.1b]
with the consequences:
ex = X" / ( L + 2*G )
and so, from the definitions of L and G,
ex = ( X"/E )*(1+P)*(1-2P)/(1-P) [7.5]
** Boundary condition 3.2 enables the y-direction stress to be
deduced, thus:
ey = 0 = ( Y" - L*e ) / (2*G) [3.2b]
whence
Y" = L * D
and so, since D equals ex, which is given by equation 7.5 ,
Y" = X" * P / (1-P) [7.6]
** Similarly, boundary condition 3.3 leads to:
Z" = X" * P / (1-P) [7.7]
Case (iii): motion disallowed in y but not in z directions.
** this case is characterised by ey=0 and sz=0. Therefore 3.3a and
3.2b are valid.
** Then the dilatation is given by:
D = ex + ey + ez
= ex - L*e/2*G
so that
D = ex * 2G / (L + 2*G)
** Substitution in 3.1 now leads to:
ex = [ X" - L * ex * 2*G / ( L + 2*G)] / 2*G
and so to
ex = (X"/E) * (1+P) * (1-P) [7.8]
** Finally,
sy = L * D
= L * (X"/E) * (1+P) * (1-P)
= (X"/E) * P / (1-2*P) [7.9]
7.2 Uniform temperature
-----------------------
a. The problem
--------------
A rectangular parallelepiped is held at its low-x, low-y and low-z
faces and is heated to a temperature T above its reference value.
Four cases are considered:
(i) the high-x, high-y and high-z faces are free to move;
(ii) the high-x face is fixed but the high-y and high-z faces
remain free;
(iii) the high-x and high-y faces are fixed but the high-z
remains free;
(ix) all faces are fixed;
Case(i): motion allowed in all three directions
Since there are no stresses in any direction, addition of equations
3.1, 3.2 and 3.3, together with the substitution of D for ex+ey+ez,
leads to:
( 2*G + 3*L )*e = H*T
which, when G, L and H are expressed via their definitions in terms
of E, P and alpha in the nomenclature, reduces to:
D = 3 * alpha*T [7.10]
and so
ex = ey = ez = et = alpha*t [7.11]
Neither poisson's ratio P, nor Young's modulus E, plays any part in
these formulae, which correspond to stress-free isotropic thermal
expansion.
Case(ii): motion prevented only in the x-direction
In this case, ex must be zero; therefore the relationship between the
dilatation D and the thermal strain et can be derived by adding
equations 3.1 and 3.2, and setting ey+ez equal to D.
The result, after substitution from the definitional relationships,
is:
D = alpha*T * 2 * (1 + P) [7.12]
ey = ez = alpha*T * (1 + P) [7.13]
Understandably, when P is zero, the dilatation is alpha*T times 2,
rather than times 3 as in case (i), because the material can expand
in only two directions. With finite P however, the x-direction
compressive stress causes additional strains in the y- and
z-directions; so the dilatation becomes larger.
The magnitude of the x-direction stress is obtained by substitution
from 7.12, and the definitional relationships, into 3.1. The result
is:
sx = - alpha*E*T [7.14]
Case(iii): motion prevented in both the x- and y-directions
When ez is the only non-zero strain, the dilatation is deduced by
setting D=ez in equation 3.3, with the result, after the usual
substitutions:
D = ez = alpha*T * ( 1+P ) / ( 1-P ) [7.15]
Once again, the dilatation increases with increase in Poisson's ratio,
and indeed twice as rapidly as before, because there are two
compressive stresses (i.e. those in the x- and y-directions)
contributing to z-direction strain.
The corresponding x- and y-direction stresses are:
sx = - (alpha*E)*T / ( 1-P ) [7.16]
Case(iv): Motion prevented in all three directions.
In this case, ex, ey, ez and therefore D also, are all zero; and the
stresses in the three directions are directly deducible from 3.1, 3.2
or 3.3 as:
sx = sy = sz = - H*T = - (alpha*E) * T / ( 1-2*P) [7.17]
This shows that the stress increases with Poisson's ratio,
corresponding to the fact that compression in one direction tends to
cause expansion in a direction at right angles which, if it is to be
prevented, requires a further compressive stress.
7.3 Non-uniform direct stress
-----------------------------
Let the force applied to the high-x face of the block in section 7.1
be removed; but let the same total force be applied uniformly through
the body, as though by a gravitational body force. We may then
expect the displacement at the end to be exactly one half of what
was before.
Let the value of this force, per unit volume, be X, which will be
equal to the previous X" divided by lx.
Once again, three cases will be considered, namely:
(i) strain is permitted in the y- and z-directions;
(ii) strain is not permitted is not permitted in either y- or
z-directions;
(iii) strain is permitted in the y- but not the z-direction .
In all cases, the direct stress per unit area which is needed to
balance the force is: X*(lx - x), having its maximum value at the
fixeb end where x=0 and its minimum at the free end where x=lx .
Case (i) Y- and z-direction movement permitted.
----------------------------------------------
The distribution of x-direction displacement will be deducible by
integrating the differential equation:
ex = u = (X/E)*( lx - x )
x
with the result:
u = (X/E)*( lx*x - 0.5*x**2 ) [7.18]
As expected, the maximum displacement, namely 0.5*(X/E)*lx**2 where
x-lx, is exactly half that given by equation 7.1a (with X"= X*lx).
The y- and z-direction strains will be - P*ex, so that the
dilatation will be:
D = ( 1.0 - P ) * (X/E)*( lx - x ) [7.19]
which has the same value at x=0 as in section 7.1, but which
diminishes to zero at the free end (x=lx))
The stresses in the y- and z-directions will of course be zero.
Case(ii) No movement permitted in y- or z-directions
----------------------------------------------------
The analysis of this case follows the same pattern was established
for case (ii) of section 7.1; indeed all the equations there
continue to be valid if X" is replaced by X*(lx - x).
The stresses, strains and dilatation of course now vary linearly
with x.
Case(iii) No movement permitted in y- or z-directions
----------------------------------------------------
What has just been written about case (ii) applies to case(iii) also.
7.4 Non-uniform temperature
------------------------------
(a) Linear temperature profile in an unconstrained block
--------------------------------------------------------
(1) The problem
---------------
* Let it now be supposed that the low-x face of the rectangular block
is held at the temperature -T0 and the high-x face at temperature
+T0. If the thermal conductivity is uniform and there are no other
heat sources, the steady-state temperature distribution will obey:
T = T0 * ( -1 + 2x/l ) , [7.20]
where l is the x-direction dimension of the block
* Let the block be free to expand in all three directions, so that
the corresponding stresses are all zero.
It is easy to see that the local strains are represented by:
ex = ey = ez = et = e/3 = alpha*T; [7.21]
so they are also linear in x; and the x-direction displacements u
can be derived by integration of the definitional relation:
u = alpha*T0*( -1 + 2x/l) [7.22]
x
with the result:
u = alpha*T0 * ( -x + x**2/l ) [7.23]
and the consequence that u has a minimum = (1/4)*alpha*T0*l at
x=l/2 and the value zero at the high-x face.
* The block does not change in overall length because the shrinkage
near its colder surface is exactly compensated by its expansion near
the hotter surface.
* The only material property which influences the solution is alpha;
for, there being no stress, neither Young's modulus nor Poisson's
ratio have any role tp play.
(2) Comparison with the equations of sections 3 and 5
-----------------------------------------------------
* The above solution has been arrived at by direct physical
arguments; but it must also be consistent with, and derivable from,
the above-derived general equations. This will now be checked.
* Equation 3.1 enables the stress in the x-direction, X", to be
deduced as:
X" = ex*2*G + L*e - H*T [3.1c]
Substitution from the solutions for ex, D and T leads to:
X" = 2*G*alpha*T + L*3*alpha*T - H*T
whereafter introduction of the definitions of G, L and H leads to:
X" = alpha*T*E * ( 1/(1+P) + 3*P/((1+P)*(1-2P)) - 1/(1-2P) )
= alpha*T*E * ( 1-2P + 3P - (1+P) ) / ((1+P)*(1-2P)) )
= alpha*T*E * ( 0 ) / ((1+P)*(1-2P)) )
= 0
as expected.
* Equation 5.13 can be checked as follows:
(L+2*G)*u + (L+G)*(j - k ) - H*T + X = 0 [5.13a]
^ z y x
As has already been explained in section 5.4.2,
(L+2*G)*u = - (L+2*G)*T*alpha
^ x
and
(L+G)*(j - k ) = 4*(L+G)*T*alpha
z y x
Since X. the externally-applied force equals zero, the equation
reduces to:
[ - (L+2*G) + 4*(L+G) - H ] * T = 0
x
Omission of the factor T , and introducing the definitions of L, G
x
and H, enable theis equation to be written as:
3*L + 2*G - H = 0
i.e.
3*P/[(1+P)*(1-2P)] + 1/(1+P) - 1/(1-2P) = 0
i.e.
3*P + 1-2*P - 1 - P = 0
i.e.
(3 - 2 - 1)*P + 1 - 1 = 0
which is evidently true.
(b) Linear temperature profile in a laterally-constrained block
---------------------------------------------------------------
To be supplied.
7.5 Non-uniform bending stress
-------------------------------------
...... to be continued ................
8. Solutions in terms of 'stress-functions'
-------------------------------------------
8.1 Introduction
----------------
Although the differential equations derived above, with the dependent
variables u, v, w, i, j, k and T, provide a complete description of
the general three-dimensional stress-strain-distribution problem, it
is sometimes easier to solve other differential equations which can
be derived from them.
Such equations are those in which the same laws are re-formulated in
terms of derived equations of which the dependent variables are
scalars, known as 'stress-functions.
In general, three such functions are required; and their
differential equations are of fourth order. However there are two
practically important cases in which a single second-order
differential equation suffices.
These cases involve the two-dimensional stress-strain distributions
caused by:
(a) the torsion of straight bars, and
(b) the bending of beams.
Both are of practical importance in engineering; and their study
also throws valuable light on how boundary conditions are to be
formulated in general three-dimensional problems when the rotation
equations have to be solved.
They are discussed in sections 8.1 and 8.2 respectively.
In both cases the treatment follows closely that which is to be
found in T&G, chapters 10 and 11, which itself is modelled on the
1855 treatise of St Venant.
St Venant employed the so-called 'semi-inverse method', which
involved assuming analytical expressions for the strain components
and proving that these conformed with the fundamental differential
equations if the forces were appropriately applied.
8.2 The torsion of straight bars of non-uniform cross-section
-------------------------------------------------------------
a. Expressions for the strain components
----------------------------------------
St Venant's assumptions were, when z measures the distance along the
length of the bar and x and y are the distances measured within its
cross-section:
(1) The u and v displacements within the cross-section are given by:
u = - theta * z * y (8.2 - 1),
v = theta * z * x (8.2 - 2),
and
(2) the w displacement is given by:
w = theta * psi(x,y) (8.2 - 3),
wherein theta is the rotation angle per unit length, and psi is a
function to be determined, which describes the 'warping of
cross-sectional surfaces.
These assumptions imply, by reason of the definitions of section 1
above:
ex = ey = ez = exy = 0 (8.2 - 4)
and
ezx = theta * ( psi - y ) (8.2 - 5)
x
eyz = theta * ( psi + x ) (8.2 - 6)
y
Expressed in terms of the components of rotation, these assumptions
imply, again by reason of the definitions of section 1:
i = theta * ( x - psi ) (8.2 - 7),
y
j = theta * ( psi + y ) (8.2 - 8),
x
k = - theta * 2 * z (8.2 - 9).
b. Corresponding expressions for the stress components
------------------------------------------------------
T&G show further that the components of stress can be
correspondingly expressed as:
sx = sy = sz = sxy = 0 (8.2 - 10)
and
szx = G * theta * ( psi - y ) (8.2 - 11)
x
syz = G * theta * ( psi + x ) (8.2 - 12)
y
At the lateral boundary of the bar, because there are no external
forces, the szx and syz components of the shear stress are linked by
the relation:
szx * l + syz * m = 0 (8.2 - 13)
where i and m are the direction cosines of the outward-directed
normal to the boundary surface.
c. (Alternative) equations to be solved
---------------------------------------
T&G show that the equations of section 2 imply that psi obeys the
differential equation:
psi + psi = 0 (8.2 - 14)
xx yy
However, they also show that the boundary conditions are simpler for
a similar equation, derived from (8.2-14), with the dependent variable
phi. This is:
phi + phi = F (8.2 - 15)
xx yy
wherein the 'stress-function' phi is related to the shear stresses
by:
szx = phi (8.2 - 16)
y
syz = - phi (8.2 - 17)
x
and
F = - 2 * G * theta (8.2 - 18)
The boundary condition for phi at the bar boundary then reduces
simply to:
phi = 0 (8.2 - 19) .
Finally it is shown that the twisting moments at the two ends of the
bar are given by:
moment = 2 * integral phi * d area (8.2 - 20)
It will be found useful, below, to express the rotation components
in terms of the phi. This can be done as follows:
* addition of 8.2-6 to 8.2-7 leads to:
i = 2 * theta * x - eyz
= 2 * theta * x - syz / G
* subtraction of 8.2-5 from 8.2-8 leads to:
j = 2 * theta * y + ezx
= 2 * theta * y + szx / G
* substitution for syz and szx from 8.2-16 and 8.2-17 then leads to:
i = 2 * theta * x - phi / G (8.2 - 21)
x
j = 2 * theta * y + phi / G (8.2 - 22)
y
d. Solution by way of PHOENICS
------------------------------
Equation (8.2-15), with boundary condition (8.2-19), is of course
easily solved by means of PHOENICS, no matter what shape the solid
posesses; and indeed the more general equation in which the material
properties can vary with x and y, namely:
( G phi ) + ( G phi ) = F (8.2 - 23)
x x y y
is equally easy.
[ Note that such a case appears on p.160 of O.Zienkiewicz's 1967
book entitled The Finite-Element Method; it will be instructive to
compare the solution given there with that computed by way of
PHOENICS ].
All that is necessary is to:
* declare a scalar variable called phi;
* de-activate its convection and time-dependent terms by way of the
TERMS command;
* impose the phi=0 condition on its boundary;
* supply a uniform source per unit volume to cells within is
boundary; and
* activate the linear-equation solver.
e. The bar of elliptic cross-section
------------------------------------
Analytical solutions are provided for several bar cross-sections
by T&G; these can also serve as tests of the accuracy of solutions
computed by PHOENICS.
That for the bar of elliptical cross-section, of which the semi-axis
lengths are a and b, is expressible as:
a**2 * b**2 * F ( x**2 y**2 )
phi = --------------- * ( ---- + ---- - 1 ) (8.2-24)
2*(a**2 + b**2) ( a**2 b**2 )
F = - 2 * moment * (a**2 + b**2) / ( pi * a**3 * b**3)
(8.2-25)
szx = - 2 * moment * y / ( pi * a * b**3 ) (8.2-26)
syz = 2 * moment * x / ( pi * a**3 * b ) (8.2-27)
Differentiation of 8.2-24 with respect first to x and then to y
leads to:
phi = b**2 * F * x / (a**2 + b**2)
x
and
phi = a**2 * F * y / (a**2 + b**2)
y
Substitution of these expressions into 8.2-21 and 8.2-22, with F
substituted from 8.2-18, leads to:
i = 2 * theta * x * a**2 / ( a**2 + b**2 ) (8.2-28)
and
j = 2 * theta * y * b**2 / ( a**2 + b**2 )
(8.2-29)
The third rotation component k, meanwhile, continues to obey the
already-presented equation, valid for all cross-sections:
k = - 2 * theta * z (8.2 - 9).
It thus appears that the distributions of i, j and k are all
linear and have magnitudes which are proportional to the distance
from the origin of coordinates.
e. Discussion
-------------
It is interesting to consider how the stresses and strains could
have been calculated WITHOUT recourse to the stress-function
technique, i.e. directly from the equations for the rotation
components and the displacements embodied in PHOENICS.
(1) Differential equation 4.1 for k is easily solved; for the source
K is zero within the bar; the boundary conditions are k=0 at z=0
and k=-2*theta*z_end at the end of the bar; and there are no forces
at the lateral boundary which can act as sources of k.
[ Here it must be remembered the k is defined as TWICE the true
rotation, and that theta represents the TWIST PER UNIT LENGTH. ]
Therefore numerical solution of the k equation will lead to the
linear distribution of k expressed by 8.2-9.
It is noteworthy that the gradient of k is finite only in the z
direction, and that this gradient does not appear as a source in any
of the differential equations for the displacements 5.13, 5.14 or
5.15.
(2) Let us assume (and check later) that gradients of i and j are
finite only in the x and y directions respectively and so also make
no appearance in the equations for the displacements.
(3) Then the equations for u and v, namely 5.13 and 5.14, are also
easily solved numerically; for their values a z=0 and z=z_end are
known, from the imposed twist; and their are no forces in their
respective directions imposed anywhere else.
Their numerical solutions will undoubtedly agree with equations
8.2-1 and 8.2-2 above.
(4) There also appear to be no forces in the z-direction which could
affect the w equation. However, the boundary shear-stress conditions
represented by 8.2-11, 8.2.12 and 8.2-13 are relevant.
These may be usefully re-written in terms of ezx, eyz and gradients
of w as:
ezx = w - theta*y ; eyz= w - theta*x ; ezx*l + eyl*m =0
x y
It thus follows that there exist easy-to-implement gradient boundary
conditions for the w equation which allows numerical solution of the
equation to proceed.
(5) With u, v and w now computed, the rotations i and j can now be
evaluated; and it will undoubtedly be found that these quantities
have gradients in the x- and y-directions respectively, so
justifying the assumption (2).
(6) It may be concluded that solution of the PHOENICS solid-stress
equation system can indeed be conducted straightforwardly for
problems in volving bars of arbitrary cross-section.
8.3 The bending of beams of non-uniform cross-section
-----------------------------------------------------