See discussions, stats, and author profiles for this publication at: https://www.researchgate.

net/publication/269030853

Generalization of Michell’s solution to the plane problem theory of elasticity

in polar coordinates in the event of a radially inhomogeneous body

Conference Paper in WIT Transactions on Modelling and Simulation · September 2014

DOI: 10.2495/BE370181

CITATIONS READS

3 939

2 authors:

Vladimir I. Andreev Nikita Yu Tsybin

Moscow State University of Civil Engineering Moscow State University of Civil Engineering
81 PUBLICATIONS 404 CITATIONS 22 PUBLICATIONS 136 CITATIONS

SEE PROFILE SEE PROFILE

All content following this page was uploaded by Nikita Yu Tsybin on 28 June 2018.

The user has requested enhancement of the downloaded file.

 Boundary Elements and Other Mesh Reduction Methods XXXVII 215

Generalization of Michell’s solution to the

plane problem theory of elasticity in polar
coordinates in the event of a radially
inhomogeneous body
V. I. Andreev & N. Y. Cybin
Department of Strength of Materials,
Moscow State University of Civil Engineering, Russia

Abstract
In the plane problem of elasticity of inhomogeneous bodies problems are often
found with the radial inhomogeneity occurring in the presence of axially
symmetric physical fields (temperature, radiation, etc.). At the same time the
plane problem itself can be two-dimensional. Of homogeneous bodies Michell’s
solution for Airy’s stress function is the most well known. This solution is
presented in the form of an infinite series in the trigonometric functions with
constant coefficients. This article considers the statement of the problem in
displacements, when the main unknown chosen functions are ur ,  and vr ,  .
The solution in displacements has the advantage that if the boundary conditions
are in displacements it is not necessary to integrate Cauchy relations.
Displacements are also represented in the form of series, but unlike in Michell’s
solutions the coefficients of trigonometric functions are also functions that
depend on the radius. They are also solved in the example.
Keywords: Michell’s solution, theory of elasticity, plane problem,
inhomogeneity.

1 Introduction
This paper considers the method of solving the plane problem of the theory of
elasticity of inhomogeneous bodies in polar coordinates. For the case of radial
inhomogeneity there is a method of separation of variables which is a
generalization of Michell’s solutions for a homogeneous body [1].

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
doi:10.2495/BE370181
216 Boundary Elements and Other Mesh Reduction Methods XXXVII

2 Main equations
For the case where the mechanical characteristics of the material depend only on
the radius the equations of the plane problem in polar coordinates are
1 u u    2u     2v
2
  2   
u
  
 r
2
r r r 2  r 2 2 r r (1)
  3 v    2  u 1  
 2


r

r r r 
v 
u    3
 

r
K f  R  0;  
r
  2 v 1 v v      2u   3 u
 2     
 r r r r 2  r r r 2 
(2)
  2  2 v 1   u
  
r r  
r
v
r
 3 
 v  K f    0,  
 r 
2 2
r 
where   (r ),   (r ) ,  f denotes forced (e.g. temperature) strains, R,  are
mass forces. In the plane stress problem  and  must be replaced by
2 Е 2 К Е
*   ; *   .
  2 1   2   2 31   

3 Solution for homogeneous material

If the mechanical properties of the material are constant it is known as Michell’s
solution for the Airy stress function [1]. In the notation of the coefficients
according to [1], we can write the complete solution in the form:
a1
  a0 ln r  b0 r 2  c0 r 2 ln r  d 0 r 2   a0   r sin   (b1r 3 
2
c
 a1r 1  b1r ln r ) cos   1 r cos   (d1r 3  c1r 1 
2

 
 d1r ln r ) sin    an r n  bn r n  2  an r  n  bn r  n  2 cos n 
(3)
n2

 c r 

 n
n
 d n r n  2  cn r  n  d n r  n  2 sin n.
n2

Stresses can be determined by known formulas:

1  1  2   2 1  1  2 
r   2 ;    2 ;  r  2  .
r r r  2
r r  r r
Having determined the deformation from Hooke’s law and integrating Cauchy
relations, we obtain expressions for displacements. For the case of plane stress
state they are have the form [2]:

1 
Eu  a0  2b0 1   r  c0 21   r ln r  1   r   d 0 1   r 
r

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
 Boundary Elements and Other Mesh Reduction Methods XXXVII 217

 a 1     c 1   
 1
 2

 2b  sin    1
  2



 2d1   cos   c1 ln r  d1 1  3 r 2 

 
 c1 1   r 2  d1 1   ln r sin   a1 ln r  b1 1  3 r 2  a1 1   r 2 
 b 1   ln r cos    a n1   r

1 n
n1
 bn n  2  n  2r n1  (4)
n 2

 an n1   r n1  bn n  2  n  2r n1 cos n    – n n1   r n1 



n2
 d n n  2  n  2r n1  cn n1   r n1  d n n  2  n  2r n1 sin n;

1 
Ev  4c0 r  4d 0 r ln r  a0  c1 1    2  2d1  sin  
r

 a1 1    2  2b1  cos   а1 (2 ln r  1  ) 2  b1 5   r 2 
 a1 1   r 2  b1 1    ln r   ln r sin   с1 (2 ln r  1  ) 2 
 d1 5   r 2  c1 1   r 2  d1 1    ln r   ln r  cos    (5)
 an n1   r

 n1
 bn 4  n1   r n 1
 an n1   r  n1

n2


 bn 4  n1   r n1 sin n   сn n1   r n1 d n 4  n1   r n1 


n2
 cn n1   r n1
 d n 4  n1   r n1 cos n.
In general, the solutions of equations (1) and (2) depend on the functions
 r ,  , r ,  , K r ,   and cannot be written in closed form. In the following
for the case when the mechanical properties together with their first derivatives
are continuous functions of the radius, by the method of separation of variables
equations (1) and (2) are reduced to systems of ordinary differential equations,
which, by analogy with the Michell’s generalized solution, are called
generalized equations of the plane problem in polar coordinates.

4 Generalized equations of the plane problem for a radially

inhomogeneous body
Using the analogy with the solutions (4) and (5) for a homogeneous material, we
seek the solution of equations (1) and (2) in the form
 u   0   1     
            2    sin    3    cos  
 v    0   1   2   3 
 
(6)
      sn    
  4   sin    5   cos  
 4   5 
 
n  2   sn 
 sin n   cn  cos n,
  cn  
where  0 ,... cn are the functions depending only on the radius. It is assumed
that  f r ,  , R r ,  ,  r ,   also might also be presented by series similar
to (6).

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
218 Boundary Elements and Other Mesh Reduction Methods XXXVII

Substituting (6) into (1) and (2) by equating to zero the sum of terms that do
not contain  , as well as the amounts which are multiplied by ,  sin  , and so
on, leads to an infinite system of ordinary differential equations for the functions
i , i . Below are the equations for some functions from the expressions (6)
 0    
  2  0  0
 
2 
1 
  3
1 
 r r  r r2
(7)


   20  0  1   3 K f ,0   R0  0;
r

   3    
1  1    0  0  0  
r r2  r r2  (8)
 3K
 1  r0   0    f ,1  Т 0  0;
r r
 
 
  2  1  1  21     21   1  3 K f ,1   R1  0;
r r  r
(9)

     
 1  1  21   r1  1   Т1  0; (10)
 r r  r
  n 2
  2  sn   sn
 sn  
2 2
 sn 

n cn 
  3
n cn 
 r r  r r r2 (11)


    2sn   sn  n cn   3 K f , sn   Rsn  0;
r

   3       2 2
 ncn    sn  sn  
n cn    sn n  sn 
2
r r  r r2  r2 (12)
 3K
 n cn  r sn   sn   n f , cn  Т sn  0;
r r
    n 2
  2  сn
  cn  cn    cn 

n sn 
  3
n sn 
 r 2
r  r 2 r r2 (13)

 
    2 cn   cn  n sn   3 K f , cn   Rcn  0;
r
   3       2 2
nsn    cn  cn  
n sn    cn n  cn 
2
r r  r r2  r2 (14)
 3K
 n sn  r cn   cn   n f , sn  Т cn  0.
r r
In the above equations, the prime denotes differentiation with respect to the
radius. Obtained system of equations is partially decomposed. For example,
equations (11) and (14) form a closed system with respect to a pair of functions
 sn and  cn , and equations (12) and (13) respectively to functions  сn and

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
 Boundary Elements and Other Mesh Reduction Methods XXXVII 219

 sn . The rest of the functions in (6) can be found consistently from the
corresponding equations. For example, from (9) and (10) we can find 1 and  1
after which from (7) and (8) we can find  0 and  0 .
To equations (7)–(14) shall be added boundary conditions of any type.
Formulas for stresses can be obtained from Hooke’s law with Cauchy relations
and expressions (6).

5 Two analytical solutions

5.1 Pure torsion of thin ring
Consider a thin ring fixed on the inner contour ( r  а ), and on the outer ( r  b )
loaded with uniform tangential forces with intensity q (Fig. 1). Because the
displacements v do not dependent on  , the solution is sought in the form
v   0 . Other functions  i and all functions i are equal to zero. According to
(8)  0 satisfies the equation

 1     1  
 0   0     0     0, (15)
r   r r  
and also to boundary conditions:
  
r  a,  0  0; r  b,   q   0  0   q. (16)
 r 
C2
It is easy to show that if    0  const is a solution of (15) then  0  С1r 
r
which corresponds to the members of the formula (5), and comprising constants
с and а0 . By satisfying the conditions (16) in the case of a homogeneous
material we obtain
2
( hom) qb2  r a  ( hom) b
v     ; r  q  . (17)
2 0 a  a r  r
Proceeding to the solution for the inhomogeneous material we will consider
one of the simplest dependencies    0 r / a  . In that event equation (15) is the
Euler equation, the solution of which is
 0  C1r  C2 r   1 .
Having determined the constants from (16) we find v and  r :
2
qb 2  r a  1  b
v (inh )    ; (rinh )
 q   .
  2 0 a  a r  1   r
(18)

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
220 Boundary Elements and Other Mesh Reduction Methods XXXVII

Figure 1: Design scheme.

In expressions (17) and (18) the superscripts refer respectively to the

homogeneous and inhomogeneous material. Comparing the two solutions, we
can see that the stresses in the ring do not depend on the inhomogeneity, which
can be explained as follows. This problem is statically determinate. If you cut the
ring the inner radius of which is r0  a (Fig. 2), the stresses on the inner contour
of the ring can be calculated from the condition that the torque about the axis x is
equal to zero:
 M ( x)  q    2b  b  (r0 )    2r0  r0  0  (r0 )  qb2 / r02 .
Displacements in the two cases are different. The ratio of displacement of the
outer contour points in the inhomogeneous and homogeneous rings is

v (inh) (b) 2(k   2  1)

 ,
v ( hom) (b) (  2)(k 2  1)k 

where k  b a . Analysis shows that when   0 , which corresponds to

toughening the material with increasing radius, this ratio is always less than
unity.

Figure 2: Equilibrium of the ring: M ( x)  0 .

5.2 Ring under the action of normal and shear loads

Let us consider the problem of equilibrium in the thin ring when on the outer
surface the loads are distributed:
p  p0 (1  cos 2) 2 ; q  p0 sin 2 2 , (19)

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
 Boundary Elements and Other Mesh Reduction Methods XXXVII 221

and on the inner surface there are no loads (Fig. 3). In the case a  b it is the
Kirsch’s problem about stretching of the thin plate with a small circular hole.

Figure 3: Calculation scheme.

Let the modulus of elasticity of the ring’s material vary according to the
power law

r
E r   E0   , (20)
a
and Poisson ratio   const .
To satisfy the boundary conditions
r  a,  r   r  0;
r  b,  r  p(),  r  q() (21)

in the expressions for the stresses we can restrict the terms containing the
functions 0 , 1 ,  c 2 and  s 2 :

E    1     2 s 2  
r   0   0    c 2   c 2   cos 2 ;
(1   )(1  )  r   r  
E  0  1     2 s 2   (22)
    0    c 2  c 2   cos 2;
(1  )(1  )  r   r  

E      s 2  2 c 2  
r   1  1      s 2    sin 2 .
2(1  )  r   r  

For determining functions 0 , 1 ,  c 2 and  s 2 it is necessary to consider

equations (7), (10), (12) and (13) for n = 2, taking into account  в  R    0 .
To go from a plane strain to plane stress state we will replace  with
  E /(1   2 ) .
From the condition of uniqueness, the solution must not contain  . Thus,
from the third equality (22) it follows that
1
1   0.
r

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
222 Boundary Elements and Other Mesh Reduction Methods XXXVII

By using this equality equation (10) takes the form

1  0 .
The integral of this equation is a linear function which comes is constant at
r  a, b . Taking into account that boundary condition for  r , equation (19)
does not contain a constant so we must assume 1  0 . Then equation (7) is
simplified:

*  2 0  r0  r 20   (*)  20  (r ) 0  0 .

*
(23)

Thus, the problem can be separated into two problems: determination using
equation (23) for the axisymmetric part of the solution and the determination
from equations (12) and (13) which depends on  . Substitution of (20) into
these equations leads them to the form

0 
0  1       1 20  0; (24)
r r

 9  21    c 2  6  2  (25)

s 2  1   s 2      s2    2  c 2  0;
2
r  1    r 1   r  1    r2

    (26)
с 2  1    c 2    3  2  c 2  1    s 2  3    2  s 2  0.
2
r r r r2
Solution of equation (24) is the function [3]:
1  1  
(27)
0  С1r 2  C2r 2 ,

where   1    ,   1   2  4  1 .
Integration constants appearing in (27) can be found from the boundary
conditions for the axisymmetric component of the external load:
r  a,  r  0; r  b,  r  p0 2 . (28)
The system of two ordinary differential equations (25) and (26) can be reduced
to a single fourth-order equation:

 
r 4 cIV2  6  2  r 3с2  5   2    3 r 2 с 2 
(29)
  
  2    9  9 rс 2  3 2  3  3  9 с 2  0. 
The obtained equation can be reduced to a differential equation with constant
coefficients by introducing variable t using dependence r  exp( t ) :

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
 Boundary Elements and Other Mesh Reduction Methods XXXVII 223

d 4 c 2 d 3 c 2
  d dt
2

4
 2 3
  2      10 2
c2

dt dt

  2   2  10  ddt  3
c2 2

 3  3  9 с 2  0.

The characteristic equation corresponding to the obtained equation is

  
l 4  2l 3   2    10 l 2   2   2  10 l  
2
 3  3  3  9  0.

Using substitution   l 2  l this equation can be reduced to a quadratic

equation:

2      10   3 2  3  3  9  0.

The final solution of equation (29) is represented as
4
c 2   Di elt .
i 1

Constants Di determined from the boundary conditions for the non-

axisymmetric component:
r  a,  r   r  0;
r  b,  r  p0 cos 2 2 , , (30)
 r   p0 sin 2 2 .
Function  s 2 can be found from eqn (25). Below is an example of the
calculation performed for the following initial data:   1; b a  2;   1 3 ;
Е  2 104 MPa.
Fig. 4 shows diagrams of stresses   along three radial directions. You may
see that greatest differences of stress for inhomogeneous and homogeneous
materials do not exceed 20%. Fig. 5 shows diagrams of the displacements of the
inner contour points of the ring along the angular coordinate. It can be noted that
displacement in the inhomogeneous ring is greater than in the homogeneous ring.
This is due to the fact that at   1 the modulus decreases from the inner to the
outer contour twice, and this leads to a decrease in the total ring stiffness.

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
224 Boundary Elements and Other Mesh Reduction Methods XXXVII

Figure 4: Stresses  in the ring. 1 –   0 ; 2 –    4 ; 3 –    2 ;

——— inhomogeneous material, - - - - homogeneous material.

Figure 5: Displacements of inner contour of the ring. ——— inhomogeneous

material, - - - - homogeneous material.

6 Numerical–analytical method of solution

If the dependencies Е r  and  r  are complex enough then it is not possible to

obtain an analytical solution. In this case it is necessary to use numerical
integration methods. Reduction of partial differential equations to a system of
ordinary differential equations and their subsequent numerical solution is called
the numerical–analytical method (this method is sometimes called the
semianalytic method).

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
 Boundary Elements and Other Mesh Reduction Methods XXXVII 225

The above equations (7)–(14) and the equations for the remaining functions
i and  i included in expressions (6) can be divided into groups of two
second-order equations for a pair of functions. Two second-order equations may
be reduced to the four equations of the first order, and the resulting system of
four equations can be solved numerically using one of the computer systems
such as Matlab. In the following we give an example of an appropriate
calculation.
Here is the solution of the problem discussed in Section 5.2. We reduce
equations (12) and (13) to a system of four first-order equations. We introduce
the notation
y1  c 2 , y2  c 2 , y3   s 2 , y4  s 2 (31)
Then equations (12) and (13) can be written as

   3  y y    2 2
 ny2  ny1   y 4  4  3   n y3 
2 r r2 
r r  r2 (32)
 3K
 ny1  ry 4  y3   n f , cn  Т sn  0;
r r

y1  n 2
  2  y 2  y 2 
2

2 1
y 

ny 4 
  3
ny3 
 r r  r r r2 (33)

 
    2 y 2   y1  ny3   3 K f , cn   Rcn  0;
r
Adding to these two equations the equalities
y1  y 2 ; y3  y 4 , (34)
which follow from (31), we receive the system of four first-order equations.
Boundary conditions for non-axisymmetric component (30) can be written as
y  2 y3 y  2 y1 (35)
r  a, y2   1  0 ; y4  3 0.
a a
 E  y1  2 y3  р0
 (1  )(1  )  y2    ;
r  b,  b  2 .
 (36)
 E  y  2 y1  р
y  3  0.
 2(1  )  4 b  2

Thus equations (32)–(34) with boundary conditions (35) and (36) represent the
boundary problem for the non-axisymmetric part of the solution.
Similarly, in view of 1  0 the boundary value problem can be presented for
the axisymmetric component:
 
0  1    0    1 0  0; (37)
r r2

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
226 Boundary Elements and Other Mesh Reduction Methods XXXVII


r  a, 0   0  0 ;
а
  1  р0
(38)
E 
r  b,  0   0 
(1   )(1   )  r  2
Equation (37) with boundary conditions (38) is also easily solved using any
software package.
To determine the accuracy of the numerical–analytic methods we performed
calculations of the same sample as in Section 5.2 for the same initial data.
Interval [a, b] is divided into 100 steps. Table 1 shows the comparative values of
the stresses in an inhomogeneous ring when   45 obtained by analytical and
numerical calculations.

Table 1: Stresses in the ring.

Analytical Numerical
r a calculation calculation
r r r r
1.0 0.0 0.0 0.0 0.0
1.2 0.252 -1.266 0.263 -1.296
1.4 0.377 -1.338 0.378 -1.363
1.6 0.444 -1.108 0.444 -1.121
1.8 0.480 -0.807 0.480 -0.813
2.0 2.000 -0.500 2.000 -0.500

Comparison of the results obtained by the two methods suggests a fairly high
precision of the numerical–analytic method.

7 Conclusions
The developed method of separation of variables in the plane problem in
aggregates with a numerical method for solving a system of ordinary differential
equations allows solutions to be obtained for a wide range of problems for the
radially inhomogeneous bodies. Such problems are encountered in heat and
power generation, the construction of underground facilities and other areas of
technology in the presence of temperature [4], radiation [5] and moisture [6–8].
The solution of such problems, taking into account heterogeneity, clarifies the
stress–strain state of structures, which may lead to their optimization and to
improvements in efficiency [9–11].

Acknowledgement
This work was supported by the Ministry of Education and Science of Russia
under grant number № 7.2122.2014/К.

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)
 Boundary Elements and Other Mesh Reduction Methods XXXVII 227

References
[1] Timoshenko, S. & Goodier J.N., Theory of Elasticity, McGraw-Hill Book
Company: New York, Toronto and London, p. 116, 1951.
[2] Andreev V.I. About one method of solving of plane problem of the theory
of elasticity for radial inhomogeneous body, Applied Mechanics: (Kiev),
23(4), pp. 16–23, 1987.
[3] Kamke E. Diferentialgleichungen, Leipzig, 1959.
[4] Andreev V.I., Minayeva A.S. Modeling of equal-stressed cylinder exposed
power and thermal stresses, International Journal for Computational Civil
and Structural Engineering. Volume 7, Issue 1, pp. 71–75, 2011.
[5] Andreev V.I., Dubrovskiy I.A. Тhe irradiation effect on the stress state of
a concrete shells, International Journal for Computational Civil and
Structural Engineering Volume 6, Issue 1,2 ; рp. 39–44, 2010.
[6] Andreev V.I., Avershyev A.S. Stationary Problem of Moisture-elasticity
for Inhomogeneous thick-walled Shells, Advanced Materials Research,
Vols. 671–674, рр. 571–575, Trans Tech Publications, Switzerland, 2013.
[7] Andreev V.I., Avershyev A.S. Nonstationary problem moisture elasticity
for nonhomogeneous hollow thick-walled cylinder, Transactions of
International Conference on Fluid Structure Interaction, pp. 123–132 WIT
Press, 2013.
[8] Andreev V.I., Avershyev A.S. About Influence of Moisture on Stress State
of Soil taking into account Inhomogeneity, International Journal for
Computational Civil and Structural Engineering №9, Iss. 3–4, pp. 14–20,
2013.
[9] Andreev V.I. Optimization of thick-walled shells based on solutions of
inverse problems of the elastic theory for inhomogeneous bodies.
Computer Aided Optimum Design in Engineering XII, WIT Press,
pp. 189–201, 2012.
[10] Andreev V.I. About one way of optimization of the thick-walled shells,
Applied Mechanics and Materials, Vols. 166–169, pp. 354–358, Trans
Tech Publications, Switzerland, 2012.
[11] Andreev V.I. Inverse problems of the inhomogeneous theory of elasticity
for thick-walled shells, Computational Methods and Experimental
Measurements, Vol. 2, No 2, рp. 202–216, WIT Press, 2014.

WIT Transactions on Modelling and Simulation, Vol 57, © 2014 WIT Press
www.witpress.com, ISSN 1743-355X (on-line)

View publication stats