Numerical solution of plane problems in the theory of elasticity

NUMERICAL SOLUTION OF PLANE PROBLEMS IN THE PHEORY OF ELASTICITY* A. I.

KALANDIA

Tbilisi (Received

23

March

1966)

THIS paper gives in numerical form a method for solving the plane problems suggested by N.I. Muskhelishvili for simply connected regions. This method makes substantial use of the apparatus of conformal mapping. A discrete analogue of the Theodorsen method is used to construct the mapping function. 1. The problem of detersining the elastic equilibrium of a homogeneous and isotropic body in the case of plane strain or plane’state of stress reduces to seeking two functions up, qtCj1 of complex argument 5, holomorphic in the unit circle, in accordance with [ll

(1.1) Here f(a) is a complex function defined on the unit circle y, cr = eff and ~(5) denotes the contour values of the function realizing the conformal mapping of the given physical region onto the unit circle \
assume m

‘p(!d=

x

m

(P’(c)=

akGk,

i

010

l

Zh.

with

vychisl.

kakgk-',

$,(G)=

i

@(d

Then,

2

OD

the

= i

,Mat.

mat. Fiz.

%b’tik

for

151-c 1, (1.2)

0

bkok,

well-known

2

f(o)=

assumptions

7,

2, 269

452

;

Akok

on y.

regardl.zg

-

460,

the

1967.

(1.3)

convergence

of

series linear 63)

Kalandia

A.I.

290

(1.2) - (1.3) to y. problem (1.1) reduces to an infinite set of equations with respect to unknowns ak and ah’ (see cl1 Section ‘g am + 2 t

k;;kb,,,+k_t = A,

anrr+ jj &L,,+~-i 0

(m = 1, 2, . ..).

= A_,

(f.4)

0.5)

(m = 0, $2, . . *)*

The unknowns ak are found from (1.4), whilst the ak’ are computed successively Prom (1.5). The functions ~(5) sad q~(c> which are found, fully describe the state of stress together with the mapping function, sought. In particular, the sum of the prlnciDa1 stresses XX + YY = Re $P’(r;) /w’(C)). When the mapping function is known the plane problem is thus reduced It is proved that if the conditions of to the solution of system (1.4). statics are observed and the right-hand side is fairly smooth, system (1.4) gives a solution of boundary-value problem (1.1). 2. For definiteness we shall regard the domain S under examination in the plane of the variable z = x + iy as finite and containing within itself the origin of coordinates: the mapping function 2 =

satisfies

the normalization

@(5)

(24

conditions o(0) = 0,

w’(O) > 0.

Let the boundary of domain S, a simple closed the equation

(2.2) curve L, be given

(2.3)

P = P(Q), where p is a fairly

smooth function

Let the relationship established L and y by (2.1) be Q(0) = argo(@).

of the polar

by

angle @, 0 < @ < 2w.

between the polar

angles

of Points

From (2.2) the funcWe shall examine the function F(c) = In(w(<)/<). tion F(c) is holomorphic Par I
Plane

problems

in the

This signifies that the fanctfons harmonic on the circle. Therefore

In

theory

p[~(~)]

291

of elasticity

and Q(D) --6

are confueate

(2.5) This is the non-linear sfngular equation first indicated in [Zl for defining the function (D(6). After determining the solution of (2.5) the boundary values F(e’O) are found directly from formula (2.4). The method of successive (2.5) in practice.

approximations

Is usually

employed to solve

A theoretical study of the method is given in [31 (see also [41). It is established in [31 tbat for some proximity of domain S to the circle and smoothness of its boundary. sequence Q,(D), deterlnined from (2.5) br iteration (with initial function (Do(S) = D), converges together with the derivative to a unique solution of (2.5) and conformally maps domain S onto the circle. The result in [41 is slightly more general. In [51, with the same conditions of proximity to the circle, the convergence of the iteration method applied to a discrete form of equation (2.5) is established. This discrete analogue will be used below to construct a numerical mapping. 3. We shall denote by Y and @ the angles formed by the external normal to L at the point with polar angle 0 with the positive direction of the x axis and the radius vector respectively. We shall read the angle v from the x axis and p from the radius vector. For a known order of reading these angles we shall have

The function In o’(c) is regular the boundary values we have

in the circle-and

real

at 5 = 0. For (3 2)

In w‘(a) = In 1o’(0) 1 + i arg O’(o). But argo’(a)=

v-6=

where s is an arc of contour On the other

hand,

B+(l)--@,

ds

as m

lo’(u) I= = --, dt? m de

L read from an arbitrary

cos fi = pdQ, 1 ds, ds / d(P = p SBC0,

point

(3.3) on it.

and therefore

A.I.

292

f o’(a) 1 .= psec ~(~)#‘(~). shall obtain the equality

Kalandia

Substituting

In o’(0) = In {P =

these

expressionS

B(@)Q’(@)) + WV)

in (3.2)

we

+ BP) - O),

in which the pair of real functions on the right-hand side represents the contour Values on the circle of the conjugate harmonic functions regular in the cirtile. Therefore, taktng into account the second condition of (2.2). we obtain 211

21f 2-e

exp! -ia

lo'(a)l=

uww)l+Qw--r)coty---

Is

dr + In o’(0)

0

1.

(3.4)

This formula together with (3.3) enables the contour values a’(0) to be expressed in terms of the solution of (2.5). and certain known functions of a point of the contour L, the elements of this contour. We shall have o’(u)=

exp

I

1 2n

2n f 0

IX:,)

{B[~(r)l+m(~)-~)cot~dr+

+Ino’(o)-t~~B[~0(6)1+4)(6)--6) From (2.4)

Using this

we shall

equality

for

0[F(6)]ei@@)

the contour Equalities

efficient

--

we can use the form

= exp: 2: L

values (3.5)

at q’(a)

.

in the same way as previously,

write,

M

w(a)=

1

0 s

f@(r)-

z-O T)cot 2

dz + In o’(O) + itD(6)

1

. (3.6)

of the mapping function.

and (3.6) in (1.1)

enable

us to represent

in the form

the required

co-

Plane

4.

It

is

discrete formulae for

the

problems

in the

theory

of

elasticity

293

to use the

convenient

method of successive approximations in of the high accuracy of the simple quadrature

form because the integral

2x

& sf(t)cot- r--62

g(6) =

dT.

0

As the nodes natural)

in the interval

CO, 2~1 we shall

ik=$ Then we have [6,

take the points

(IV is

m = 0, 1, . . . , 2N - 1.

71 ZN-1

gm= 2 aA_mfA

(m = 0, 1, . . . , 2N - I),

(4.2)

A-O

where 0 ak-m

for

lk--1

for

IA-

even,

(4.3)

= 1

ffA-ftm

-cot

N

2

ml

odd.

Formula (4.2) is accurate when f(e) is a trigonometric polynomial of an order not higher than 1) - 1. In [TI the convergence of (4.2) is examined for different classes of functions f( ). According

to

(4.3)

equations

(4.2)

can be resolved

into

tao groups

N-l g2m

=

z

a2(A_m)+if2h+i~

A-O N-l k?Zm+l

=

2

aZ(A_m)_lf2Ar

m

=

0, i, . . . , N - 1,

k-0 and the application of the iteration method to the system with matrix ok_m is therefore extremely convenient. From (4.2)

the discrete

form of equation

(2.5)

is

A.I.

Kalandia

2N-I

a,,, =

2

6, -

m = 0, 1, . . . , 2N - 1,

ak_m In p [@iI,

(4.4)

k=O al

=

p[@nIl =

@(ens),

For domains close to the circle should begin from values @2=6, 5. We shall

P

P(%n)l.

in the known sense the iteration

(m = 0, 1, . . . , 2N-

1).

assume

o(a)

= FI (a) + iFz(+-+I,

ol(0)

(5.1)

where F,(6)

= A(ti)cms a(@),

Fz(6) = A(6)sin Q(6),

2s

A(6) = exp

- &

[

s 0

{vi@(x)] - o,(T)}Cot-

2-e

at

2

I

and (5.2)

Q(6) = v[@((ft) + O(6) - 6. For computing coefficients

b,

2no(a) b,=-f-S -_

exp(-

2n

in (1.3)

defined

by equalities

imft)de

(m = 0, fl,

f2,

(5.3)

, . .),

0 o’(a)

we shall make use of a rectangular formula. which gives an accurate value of the integral for trigonometric polynomials of order not higher than PA’ - 1: 1 -

St We

shall

2N.-1

2s l

F(ft)erp(-

ime)&?

= k

,c;

F(ftk)exp(-

ime,).

(5.4)

I==0

0

have 3,

=

GN

2

{F,(h)+

iFz(flk)}

exp(-

imftk)

(m =

0, fl,

(5.5)

. . .),

k-0

and from (4.2)

and (5.2)

2N-I vk-@k)+i(vm+@m-em)

I

.

(5.6)

Plane

Ultimately problem (1.1)

problems

in

the sequence Is this.

the

of elasticity

theory

of operations

for

295

solving

the boundary-value

First of all we solve the non-linear equation (4.4) and find the values of the function O,(e) at points 6,. From the given @‘mvalues for the given angle v(o) we compute the required number of coefficients 5,. After this, solving the truncated system (1.4). allowing for the requlrement that the principal moment of the external forces be zero, we find the first coefficients a,,, of the expansion of function q(c), and then from (1.5) we find the first coefficients of the function v(j). The boundary values on the network of the mapping function are computed from formula (3.6); these together with (5.6) define the discrete values in the form o’(exp(ie,)) To define be used

these

functions

inside

pm exp(F,(6,)-

0

2ni

i(Dm)

(5.7)

iF2(em)

the circle

sa--T’

1

o(5)= -

=

(a) do o’(G)=

the Cauchy formula must

1 __

2ni

Y

s

o’(a)da ~

u-t

P

followed by replacement of the integrals by finite sums by one method of approximate integration or another. On the basis of (5.4) and (5.8) we have the following approximate formulae which are suitable strictly Inside the circle:

2N-_bA exp[[i(@k+&)l -1

w(S,=&r, A=0

eX@(

ifiA)

-

6

ZN-1

w =*;

pAexp[--(@A-@A)] 2 A=0

[Fi(flA)L

iFZ(*A)l(exp(h%A)-

5)’

The functions cp, p’, v, o, o’ fully characterize the displacement field and also the distribution of the sum of the principal stresses in an elastic body. To compute all the stress components separately we must know in addition an acceptably accurate expression for the second derivat ive of’(<). 6. Enarnple. We shall examine an elastic diagonal by two concentrated forces, equal direction, applied at opposite corners. We square as two and locate its centre at the Fig. 1). We shall

denote

square stretched along the in magnitude but opposite In shall take a side of the origin of coordinates (see

by A, 6, C and D the successive

corners

of the square

A.1.

296

CorresPonding

to

the

For definiteness unit

magnitude

angles

i/m, S/QK,=/m, T/&X.

we shall

applied

Kalandia

assume that

at points

FIG. In the have

(for

oase

under

the

function

ing 0 should

Note. the

that

contour

o’(a), simple

boundary

see

for

the

discussions of

the

line

the

functions

equality

given

L be a Lyapunov

signs

this

above

from

curve.

theory

(3.4))

of

conditions

to

be valid

of

the

of

defin-

we must specify

in a specific

of

conformal to

infinity

mappings or to

selected

we shall

then

network have

sense.

for

the

function

(it

is

zero

example also

at the

points

For this reason, near these points of the or small quantities when using our algo-

rithm, and this has an adverse effect on the calculation. there are corners we can try to use the scheme directly points

are

p and v we shall

in the

When L has corners

tends

corresponding to these corners. circle we must deal with large

However,

forces AC.

1.

domain S be smooth

as we know from the to

concentrated

be omitted)

For the

that

examination v(o)

the

A and c along

are

isolated

a contour

with

from the rounded

if

But even if the nodal

‘danger’

corners.

points. Furthermore,

Plane

problems

in

the

theory

of

297

elasticity

complex potentials p and q~ will possess (not clearly defined) singularities on the circle which are characteristic of concentrated contour loads. This, of course, also gives rise to certain difficulties In obtaining an acceptable solution of the problem. To solve system (1.4) we must first of all know the Fourier coefficients of the function f = fl + if2 with positive indices. The easiest way to determine these coefficients is to expand the Cauchy type integral

in the

circle

We shall

in a power series, recall

the

formula

fi + if* = i

s

(-cl + iY,)ds,

0

where X, and Y, are components of the given vector of the stress applied to the boundary of domain s, and the integral is taken over the arc s of the boundary contour. According to the relation t = o(u), the arc s is a specific function of the circular arc y or. what amounts to the u = ei*. We can therefore regard the right-hand same thing, of point side of (6.2) as a given function of U. For our example we have ji -I- if2 = exp (in 14) on ABC, Integral

(6.1)

f1-k if2 = 0 on CDA.

now becomes (6.3)

where yu is

the

image

of ABC when mapping

t = o(u).

We have

F(5) = where a = eia on the circle

and c are the y.

Expanding this symmetry assuming

exp(in/4)

images

2n

c- 5 In , a-f

of the

function In a series c=_ a, we find

for

points

kl

A and C respectively

< 1, and in view

of the

298

A.I.

Kalandia

1.

ni g y+-f&t. [

exp(W4) F(C) = --

a

3-c

Consequently, i exp(in/4)

A

A~~-------+,

4 When a. = r/4,

for

ZA= 0,

42~t

=

-

exp

i [n/4 -(2k

- l)a] _--_ .

n(Zk--1)

the

coefficients

.4gk_1 we shall

have

j’-k

&k-i

= -

(k = 1, 2, . . .).

n (2k -.- 1)

(6.4)

If, in accordance with what was said in the note to Section 6, we now take N = 4k + 2, where k is an integer, the points of the circle & = (2k - I)*/4 (k = 1, 2, 3, 4) are in the middle between two adjacent nodes of the given network. Table 1 gives the values of the function @b(B), found by the iteration method from (4.4) accurate to 0.1 per cent for iv’ = 14. TABLE 1

(I,

I 0.0002

I 0.1753

0.3456

I 0.5941

I 0.9767

I 1.2252 I 1.3954 I 1.5710

As can be seen from Table 1 the polar angles 6 of the circle vary insignificantly when the circle is transformed. large approximation we can assume

points i.e.

of the with a

tD(#) = 9. Consequently, the assumptions concerning the correspondence boundary of the circle and its image made above when computing efficients of (6.4) hold with the desired accuracy.

of the the co-

TABLE 2

yv

I 1.9747

I

I

I

2.2030

I

3.5638

I

3.6514

I 3.9636

I

I 4.5741

I 10.9568

Plane

prob

leas

in

the

theory

of

299

elasticity

In (1.3) the coefficients b,, a = 0, f 1, . ,., f 32 were computed and a truncated system with I = 16, k = 1, 2, . . . , 16 was examined. The system obtained for a, consists of 31 linear equations with the same nunber of unknowns. Leaving out the details we give in Table 2 the discrete values of the tensile stress Yy over a line segment DA computed In this manner. It should be remembered that the series of o(G) in the circle does not converge at all at the corners, whilst the first series of (1.3) converges very slowly. Furthermore. the series of (1.2) diverge at the points of application of the concentrated loads, which in our example correspond to the corner points. For this reason the coefficients b, and a, decrease fairly slowly, which is confirmed by the computations. Nevertheless qualitatively Table 2 reproduces the contour stress well. Quantitatively the degree of approximation employed is probably not entirely satisfactory. It should be noted that this complex variable method in [El used for mapping the square.

example was previously examined where polynomial approximations

Translated

by

J.

by the were

Cornish

REFERENCES 1.

MUSKHELISHVILI, N. I. Some fundamental problems of the mathematical theory of elasticity (Nekotorye osnovnye zadachi matematlcheskii teorii uprugosti) Izd-vo Akad. Nauk SSSR, Moscow - Leningrad, 1949.

2.

THEODORSEN,T. Theory of wing section Techn. Rept. 411, 1 - 13, 1931.

3.

WARSCHAWSKI,5. On Theodorsen’ 8 method of COnfOrmal mapping of nearly circular regions. Q. Appl. Math. 3, 1, 12 - 23, 1945.

4.

OSTROWSKI,A. On the convergence of Theodorsen’s and Garrlck’s method of conforms1 mapping. Construction and Appl. of Conformal Maps. Natn. Bur. Stand. Appl. Math. Ser. 18, 149 - 163, 1952.

5.

OSTROWSKI,A. On a discontinuous Garrick’s method. Construction Bur.

Stand.

Appl.

Math.

Ser.

of arbitrary

shape.

NACA,

analogue of Theodorsen’s and and APD~. of Conformal paps. Natn. 165

-

174.

300

A.I.

Kalandia

6.

SERBIN, Ii. Numerical quadrature of some improper Integrals. 12, 2, 188 - 194, 1954. Math.

7.

KORNEICHUK, A. A. Quadrature formulae for singular integrals. In numerical methods of solving differential and integral equations, and‘quadrature formulae. Izd-vo Akad. Nauk SSS~, Moscow, 64 - ‘74, 1964.

8.

GRAY, Yech.

C.

Polynomial appl.

Math.

approximations 4,

444

- 448,

In plane elastic 1951.

problems.

0. Appl.

Q. Jl.