********************************************************************
* *
* The equations of displacements, stresses and strains in solids *
* *
********************************************************************
by: D B Spalding
Contents
--------
1. Nomenclature
2. Strain related to stress
3. Normal stress related to strain
4. Balance of forces
5. Differential equations for the displacements
6. The differential equations for the displacements; 2D
7. The differential equations for the displacements; 1D
8. Cylindrical geometry
9. The representation of the differential equations in PHOENICS
10. Concluding remarks
--------------------------------------------------------------------
1. Nomenclature
---------------
* 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
* Linear thermal expansion: et
* Normal stresses in x, y, z directions, divided by Young's modulus,
presumed uniform: nx, ny, nz
* Shear stresses in xy, xz, yz planes, divided by Young's modulus:
sxy, sxz, syz
* Displacements in x, y and z directions: u, v, w
* Poisson's ratio: R
* Various constant functions of Poisson's ratio: R1, R2, R3, R4 ...
namely:
R1 = (1 - R) / [(1 + R) (1 - 2 R)]
R2 = R / [(1 + R) (1 - 2 R)]
R3 = 1 / (1 - 2 R)
R4 = 1 / [ 2 (1 + R)]
* derived functions:
R5 = R1 - R2 R2 / R1 = 1 / (1 + R)
R6 = R2 - R2 R2 / R1 = R / [(1 + R) (1 - R)]
R7 = R3 - R3 R2 / R1 = 1 / (1 - R)
R8 = R1 - 2 R2 R2 / ( R1 + R2) ]
R9 = R3 [ 1 - 2 R2 R2 / ( R1 + R2) ]
R10 = 1 / [ 2 ( 1 + R ) ( 1 - 2 R) ]
R11 = 1 / R4 = 2 (1 + R)
R12 = 2 ( 1 + R ) / ( 1 - R )
R13 = 2 ( 1 + R ) (1 - R )
R14 = 2 ( 1 - R )
Note that, if R = 0, then R1 = R3 = R5 = R8 = R11 = 1.0,
R2 = R6 = 0.0 ,
R4 = R10 = 0.5 and
R12= 2.0
* First-order differentiation operators
with respect to x, y, z, a: ddx., ddy., etc
* Second-order differentiation operators with respect to
x & x y & y z & z x & y x & z and y & z:
d2dxx. d2dyy. d2dzz. d2dxy. d2dxz. d2dyz.
* The divergence-of-gradient operator:
div_grad.= d2dxx. + d2dyy. + d2dzz
* The multiplication operator is implied by any space in an equation
which is not adjacent to one of the operators: + - / = or a
differential operator.
* External forces per unit volume, including gravity, aceleration,
etc, in the x, y and z directions: fx, fy , fz
2. Strain related to stress
---------------------------
The following are expressions of Hooke's and Poisson's
experimentally-based laws:
2.1 Normal strains related to normal stresses
---------------------------------------------
ex = nx - R ( ny + nz ) + et (2.1-1)
ey = ny - R ( nx + nz ) + et (2.1-2)
ez = nz - R ( nx + ny ) + et (2.1-3)
[ These correspond to the top part of matrix D on page 58 of Huebner &
Thornton (H&T form now on), except that et (thermal strain) is not
mentioned by them.
In the absence of any stresses, they correspond to ex=ey=ez=et, which
seems reasonable. ]
2.2 Shear stresses related to shear strains
-------------------------------------------
sxy = R4 ( ddy.u + ddx.v ) = syx (2.2-1)
sxz = R4 ( ddz.u + ddx.w ) = sxz (2.2 2)
szy = R4 ( ddy.w + ddz.v ) = szy (2.2-3)
[ The origin of these equations is not clear. They appear to correspond
to the bottom part of the C matrix of H&T. But the coefficient s there
are (1-2R)/2, whereas R4 = 1/{2(1+R)} ]
3. Normal stress related to strain normal strain
------------------------------------------------
3.1 General 3D equations
------------------------
The following are derived from the above by algebraic
re-arrangement:
nx = R1 ex + R2 ( ey + ez ) - R3 et (3.1-1)
ny = R1 ey + R2 ( ex + ez ) - R3 et (3.1-2)
nz = R1 ez + R2 ( ex + ey ) - R3 et (3.1-3)
The details of the derivation, for nx only, are as follows:
Addition of R times equations 2.1-2 and 2.1-3 to
(1 - R) times equation 2.1-1 leads to:
(1 - R) ex + R ey + R ez
=
(1 - R) nx - (1 - R) R (ny + nz) + (1 - R) et
+ R ny - R R (nx + nz) + R et
+ R nz - R R (nx + ny) + R et (3.1-4)
in which the multipliers of both ny and nz sum to zero.
Hence,
(1 - R) ex + R ey + R ez =
(1 - R - 2 R R) nx + (1 + R) et (3.1-5)
ie
(1 - R) ex + R ey + R ez =
(1 + R) (1 - 2 R) nx + (1 + R) et (3.1-6)
ie
nx = [ (1 - R)/ {(1 + R) (1 - 2 R)} ] ex +
[ R / {(1 + R) (1 - 2 R)}] ( ey + ez) -
[ 1 / (1 - 2 R)} ] et (3.1-7)
ie
nx = R1 ex + R2 ( ey + ez ) - R3 et (3.1-1)
[ Note that the shear-stress equations (1.2-1.2.3) have not been used
here.]
3.2 2D equations
----------------
If ez (for example) is zero everywhere, which implies that the
process is two-dimensional. movement in the third dimension being
prevented, then the equations for nx and ny become:
nx = R1 ex + R2 ey - R3 et (3.2-1)
ny = R1 ey + R2 ex - R3 et (3.2-2)
nz = R2 ( ex + ey ) - R3 et (3.2-3)
while the equation for nz remains as above.
If however there is no constraint in the z direction, ez may be
non-zero, while nz obeys:
nz = 0 . (3.2-3)
Then ez may be eliminated from the equations of section 3.1, with
the result:
nx = R5 ex + R6 ey - R7 et (3.2-4)
ny = R5 ey + R6 ex - R7 et (3.2-5)
The derivation, for nx only, is as follows:
From 3.1-3, with nz = 0, leads to:
ez = - (R2/R1) ( ex + ey ) + (R3/R1) et (3.2-6)
Substitution of ez, from 3.2-6 into 3.1-1, yields:
nx = R1 ex + R2 ( ey - (R2/R1) (ex + ey))
+ (R2 R3 /R1) et - R3 et
ie
nx = (R1 - R2 R2/R1) ex + R2 (1 - R2/R1) ey
- R3 (1 - R2/R1) et (3.2-7)
ie, from the definitions of R5, R6 and R7:
nx = R5 ex + R6 ey - R7 et (3.2-4)
Expressed in terms of R, this is:
nx = [ 1 / (1 + R) ] ex +
{ R / [(1 + R) (1 - R)] } ey -
[ 1 / (1 - R) ] et (3.2-4)
3.3 1D equations
----------------
If ey and ez (for example) are zero everywhere, the equations for
nx, ny and nz become:
nx = R1 ex - R3 et (3.3-1)
ny = R2 ex - R3 et (3.3-2)
nz = R2 ex - R3 et (3.3-3)
If however there are no constraints in the y and z direction, ey and
ez may be non-zero, while ny and nz obey:
ny = 0 and
nz = 0 ;
then ey and ez may be eliminated from the equations of section 3.2,
with the result:
nx = R8 ex - R9 et (3.3-4)
The derivation is as follows. Equations 3.1-1. 3.1-2 and 3.1-3
become:
nx = R1 ex + R2 ( ey + ez ) - R3 et (3.3-5)
0 = R1 ey + R2 ( ex + ez ) - R3 et (3.3-6)
0 = R1 ez + R2 ( ex + ey ) - R3 et (3.3-7)
Since there is nothing to distinguish ey from ez, they may be taken
as equal, with the results:
nx = R1 ex + 2 R2 ey - R3 et (3.3-5)
0 = ( R1 + R2 ) ey + R2 ex - R3 et (3.3-6)
Substitution for ey from 3.3-6 in 3.3-5 yields:
nx = R1 ex + [ 2 R2 / (R1 + R2) ] [ - R2 ex + R3 et] -
R3 et (3.3-7)
ie
nx = [ R1 - 2 R2 R2 / ( R1 + R2) ] ex -
R3 [ 1 - 2 R2 R2 / ( R1 + R2) ] et (3.3-8)
In terms of R, this is:
4. Balance of forces
--------------------
The balances of the forces acting on unit volume in the three
coordinate directions lead to:
ddx. nx + ddy. sxy + ddz.sxz + fx = 0 (4.-1)
ddy. ny + ddx. sxy + ddz.syz + fy = 0 (4.-2)
ddz. nz + ddx. sxz + ddy.syz + fz = 0 (4.-3)
5. The differential equations for the displacements; 3D
-------------------------------------------------------
5.1 The x-direction displacement u
----------------------------------
Differentiation of the equation 4.-1 leads,
after substitution from the equations of section 3.1, to:
ddx. [ R1 ex + R2 ( ey + ez ) - R3 et ] +
ddy. [ R4 ( ddy.u + ddx.v) ] +
ddz. [ R4 ( ddz.u + ddx.w) ] +
fx = 0.0
Since, by definition:
ex = ddx. u
ey = ddy. v
ez = ddz. w
the strains ex, ey, ez can be eliminated from the above differential
equation, with the result:
ddx. [ R1 ddx. u + R2 ( ddy. v + ddz. w ) - R3 et ] +
ddy. [ R4 ( ddy. u + ddx. v) ] +
ddz. [ R4 ( ddz. u + ddx. w) ] +
fx = 0.0
If Poisson's ratio R is taken as constant, implying that all Rs are
also constant, this equation can be written in terms of second-order
operators, as follows:
ddx. [ R1 ddx. u + R2 ( ddy. v + ddz. w ) - R3 et ] +
R4 d2dyy. u + R4 ddx. ddy. v +
R4 d2dzz. u + R4 ddx. ddz. w +
fx = 0.0
with further rearrangement to:
ddx. [ ( R1 - R4 ) ddx. u + ( R2 + R4 ) ( ddy. v + ddz. w ) ] -
ddx. [ R3 et ] +
R4 ( d2dxx. u + d2dyy. u + d2dzz. u ) +
fx = 0.0
But, as a consequence of the above definitions,
R1 - R4 = R2 + R4 = 1 / [2 (1 + R) (1 - 2 R)] = R10
Therefore, the equation can be reduced to:
R10 ddx. dil - R3 ddx et + R4 div_grad. u + fx = 0
where
dil = ddx. u + ddy. v + ddz. w , the dilatation.
Re-arrangement so as to give prominence to the div_grad term, leads
to:
div_grad. u + ddx. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ]
+ fx / R4 = 0.0 (5.1-1)
[ This should be compared with equation C.15 of H&T. In that equation,
the multiplier of dil is 1/(1-2R).
R10/R4 = 2(1+R)/{2(1+R)(1-2R)} i.e. 1/(1-2R) . So they equations agree.
The multiplier of fx in H&T is 2(1+R)/E, which is again in agreement.
However, H&T are silent about et .]
5.2 The y- and z-direction displacements v and w
------------------------------------------------
Replacing x successively by y and z, and u by v and w, leads to the
corresponding equations:
div_grad. v + ddy. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ]
+ fy / R4 = 0.0
div_grad. w + ddz. [ ( R10 / R4 ) dil - ( R3 / R4 ) et ]
+ fz / R4 = 0.0
Note that, as a consequence of the definitions,
R10 / R4 = 1 / ( 1.0 - 2 R ) = R3
and
R3 / R4 = 2 ( 1 + R ) / ( 1 - 2 R )
The three equations can thus be written as:
--------------------------------------------------------------------
|u| |x| |x|
div_grad. |v| + R3 dd|y|. ( dil - R11 et ) + R11 f|y| = 0.0
|w| |z| |z|
(5.2-1)
--------------------------------------------------------------------
where it is to be noted that:
R3 = 1 /( 1 - 2 R)
R11 = 1 / R4 = 2 ( 1 + R )
6. The differential equations for the displacements; 2D
-------------------------------------------------------
6.1 When displacements in the third direction are zero
------------------------------------------------------
Let the zero-displacement direction be z. Then the equation of
section 5.2 still holds, for directions x and y, with the proviso
that ddz. terms are absent from both the div_grad and dil
definitions.
6.2 When stressess in the third direction are zero
--------------------------------------------------
When nz, say, is zero, the analysis must start from
nx = R5 ex + R6 ey - R7 et
rather than from
nx = R1 ex + R2 ey - R3 et .
Thereafter, apart from the substitutions of Rn+4 for Rn, and from
the absence of z-direction derivatives, the analysis proceeds as in
section 6. The resulting equations are:
--------------------------------------------------------------------
|u| |x| |x| |
div_grad. |v| + R3 dd|y|. ( dil - R11 et ) + R11 f|y| = 0.0
(6.2-1)
--------------------------------------------------------------------
whereby it should be noted that R11 remains in place, but R13 has
replaced R12.
7. The differential equations for the displacements; 1D
-------------------------------------------------------
7.1 When displacements in the second and third directions are zero
------------------------------------------------------------------
Let the zero-displacement directions be y and z. Then the equation of
section 5.2 still holds, for directions x and y, with the proviso
that ddz. terms are absent from both the div_grad and dil
definitions.
7.2 When stressess in the second and third directions are zero
--------------------------------------------------------------
When ny and nz, say, are zero, the analysis must start from:
nx = R8 ex + R9 et .
Thereafter, the analysis proceeds as in section 6. The resulting
equation is:
------------------------------------------------------------------
|
| div_grad. u + R11 ddx. ( dil + R14 et ) + fx / R4 = 0.0
| (7.2-1)
------------------------------------------------------------------
8. Cylindrical geometry
-----------------------
8.1 Nomenclature
----------------
When the grid is cylindrical polar, and x, y and z are replaced by
a (ie angle), r (ie radial distance) and z (now axial distance),
* the strains become: ea, er, ez
* the normal stresses become: na, nr, nz
* the displacements become: u, v, w
8.2 Strain related to stress
----------------------------
The following relations replace those of section 2 above:
ea = na - R ( nr + nz) + et
er = nr - R ( na + nz) + et
ez = nz - R ( na + nr) + et
sar = R4 ( ddr.u - u / r + dda.v / r) = sra
saz = R4 ( ddz.u + dda.w / r) = sza
szr = R4 ( ddr.w + ddz.v ) = srz
8.3 Normal stress related to strain, general 3D
-----------------------------------------------
Because the e-to-n relations are the same as before (ie in section
2), so are the n-to-e relations (ie as in section 3).
8.4 Balance of forces
---------------------
In the angular direction, the balance of forces becomes:
dda. na / r + ddr. sar + ddz. saz + fa +
( sra + sar ) / r = 0.0 .
In the radial direction, the balance of forces becomes:
ddr. ( nr r ) / r + dda. sar / r + ddz. szr - na / r
+ fr = 0.0 ,
which may also be expressed as:
ddr. nr + dda. sar / r + ddr. szr + ( nr - na ) / r
+ fr = 0.0 .
In the axial direction, the balance of forces becomes:
ddz. nz + ddr. ( srz r) / r + dda. saz / r
+ fz = 0.0 ,
which may also be written as:
ddz. nz + ddr. srz + dda. saz / r + srz / r
+ fz = 0.0 .
Thus, as compared with the corresponding equations appearing section
3, the following extra terms appear:
- in the equation for the angular (formerly x) direction: none;
- in the equation for the radial (formerly y) direction:
+ ( nr - na ) / r ;
- in the equation for the axial (formerly, and still, z) direction:
+ srz / r .
8.5 The differential equations for the displacements
----------------------------------------------------
Since the differences in the differential equations from their
cartesian-coordinate counterparts all derive from the above-
mentioned extra terms, it can be deduced that the differential
equations in cylindrical-polar coordinates can be expressed as:
div_grad. u + R11 dda. ( dil + R12 et ) / r + fa / R4 = 0.0
div_grad. v + R11 ddr. ( dil + R12 et )
+ [fr + ( nr - na ) / r ] / R4 = 0.0
div_grad. w + R11 ddz. ( dil + R12 et )
+ (fz + srz / r ) / R4 = 0.0
9. The representation of the differential equations in PHOENICS
--------------------------------------------------------------
9.1 Comparison of the equations for displacements and velocities
----------------------------------------------------------------
The momentum equations which are solved by PHOENICS, when the
viscosity is set to unity and the convection terms are neglected,
can be represented as:
----------------------------------------------------------
|u| |x| |x|
div_grad. |v| - dd|y|. ( p ) + mom_so|y| = 0.0
|w| |z| |z| (9.1-1)
----------------------------------------------------------
where p stands for pressure, and mom_so stands for momentum source
per unit volume.
Comparison of these equations for those for the displacements shows
that they can be made identical if the pressure is related to the
dilatation and to the thermal expansion by:
- p = R3 ( dil - R11 et )
ie - p = [1 / (1 - 2 R)] [ dil - 2 ( 1 + R ) et ]
(9.1-2)
and if:
|x| = |x|
mom_so|y| = R11 f|y|
|z| = |z|
ie:
|x| = |x|
mom_so|y| = 2 ( 1 + R ) f|y| (9.1-3)
|z| = |z|
9.2 The dilatation in PHOENICS
------------------------------
In PHOENICS, the dilatation, ie ddx. u + ddy. v + ddz. w ,
appears as a volumetric source in the mass-conservation equation.
Therefore the required relationship can be contrived by providing,
for the part of the domain which is occupied by solids, a volumetric
source term calculated as:
vol_sor = dil = - (1 / R3) p + R11 et
ie:
vol_sor = dil = - (1 - 2 R) p + 2 (1 + R ) et (9.2-1)
9.3 Normal-stress boundary conditions
-------------------------------------
The expression for nx in equation 3.1-1 can be re-written in terms
of the dilatation as follows:
nx = ( R1 - R2 ) ex + R2 ( ex + ey + ez ) - R3 et
ie
nx = ( R1 - R2 ) ex + R2 dil - R3 et (9.3-1)
Replacement of dil by p, from equation 9.1-2, then leads to:
nx = ( R1 - R2 ) ex +
R2 [- (1 - 2 R) p + 2 (1 + R ) et ] -
R3 et (9.3-2)
This can be rearranged with ex on the left-hand side and replaced by
its equivalent, ddx.u, as:
ddx. u = [ 1 / (R1 - R2) ] nx +
[ R2 (1 - 2 R) / (R1 - R2 )] p +
[ R3 - 2 R2 (1 + R) / ( R1 - R2 ) ] et (9.3-3)
which becomes, after simplification:
ddx. u = ( 1 + R) nx + R p + ( 1 + R ) et (9.3-4)
Since ddx.u ie the gradient of the displacement u, this equation
can be used as the normal-stress boundary condition of the
differential equation for u.
It can be represented by way of a source term in the computational
cell within the solid and adjacent to the fluid, as illustrated
below.
. | . | . /|
. | . | . /|
. | . | . /|
-----*------>------*------>-----*---->=====> normal stress
. | . | . /|
. | solid | solid /| fluid
. | . | . /|
Corresponding equations can be derived for the v and w boundary
conditions, namely:
ddy. v = ( 1 + R) ny + R p + ( 1 + R ) et (9.3-5)
ddz. w = ( 1 + R) nz + R p + ( 1 + R ) et (9.3-6)
A consequence is that, when a direct stress S is to be exerted at a
boundary, the RATCH and COVAL statements in the RHOENICS Q1 file
take the form:
PATCH(name,north/east/high,,,,,,)
COVAL(name,v1/u1/w1,fixflu,S*(1+R))
A further consequence is that, within EARTH, sources per unit area
are applied to cells adjacent to solid-fluid boundaries which are
equal to:
R p + ( 1 + R ) et .
10. Concluding remarks
----------------------
Further topics to be dealt with in the continuation of the current
document include:-
* allowance for non-uniform material properties;
* the relationship of the dilatation to the PHOENICS pressure
variable;
* the solution algorithm;
* a series of exact solutions of the equations, which can be used as
tests;
* a series of practically-interesting problems, which can be used so
as to interest engineers in the use of the technique.
It should also be remarked that all the above derivations require to
be checked carefully and independently.