An efficient algorithm for finite-difference modeling of mixed-boundary-value elastic problems

Advances in Engineering Software 37 (2006) 41–55 www.elsevier.com/locate/advengsoft

An efficient algorithm for finite-difference modeling of mixed-boundary-value elastic problems M. Zubaer Hossain, S. Reaz Ahmed*, M. Wahhaj Uddin Department of Mechanical Engineering, Bangladesh University of Engineering and Technology, Dhaka 1000, Bangladesh Received 25 April 2004; received in revised form 22 February 2005; accepted 22 February 2005 Available online 13 May 2005

Abstract The paper describes an efficient numerical scheme for the solution of displacements and stresses in mixed-boundary-value elastic problems of solid mechanics. A new variable reduction scheme is used to develop the computational model. Solution of both two- and threedimensional problems of linear elasticity is considered in the present paper. In the present approach, the problems are formulated in terms of a potential function, defined in terms of the displacement components. Compared to the conventional computational approaches, the present scheme is capable of providing numerical solution of higher accuracy with less computational effort. Application of the present finitedifference modeling scheme is demonstrated through the solutions of a number of practical stress problems of interest, and the results are compared with those obtained by the standard method of solution. The comparison of the results establishes the rationality as well as suitability of the present variable reduction scheme. q 2005 Elsevier Ltd. All rights reserved. Keywords: Variable reduction scheme; Mathematical modeling; Stress analysis; Displacement potential function; Finite-differences

1. Introduction Stress-analysis of material bodies has now become a classical subject in the field of engineering. But, somehow, the subject is still suffering from a lot of shortcomings, and thus, it is constantly coming up in the literature [1–6]. The earlier mathematical models of elasticity were very deficient in handling the practical stress problems. Elastic problems are usually formulated either in terms of deformation parameters or stress parameters. Among the existing mathematical models of two-dimensional boundary-value stress problems, the two-displacement functions approach [1] and the stress function approach [7] are noticeable. The displacement formulation involves finding two displacement functions simultaneously from two second-order elliptic partial differential equations, which is extremely difficult, and this problem becomes more serious when the boundary conditions are mixed [1]. The shortcoming of the stress * Corresponding author. Tel.: C880 2 966 5636; fax: C880 2 861 3046. E-mail address: [email protected] (S.R. Ahmed).

0965-9978/$ - see front matter q 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.advengsoft.2005.03.010

function approach is that it accepts boundary conditions only in terms of loadings. As most of the practical problems of elasticity are of mixed boundary conditions, the approach fails to provide any explicit understanding of the state of stresses at the critical regions of supports. The difficulties involved in trying to solve practical stress problems using the existing models are clearly pointed out by the authors of this paper in their previous reports [8–14] and also by Durelli [4]. For three-dimensional elastic problems, existing mathematical models either involve finding six stress components—sxx, syy, szz, sxy, syz, szx or three displacement components—ux, uy, uz, simultaneously [7]. The shortcoming of the stress formulation is that six partial differential equations are required to be solved, simultaneously for the six stress components. Also, it accepts boundary conditions only in terms of boundary loading, as the boundary restraints specified in terms of ux, uy, uz cannot be imposed satisfactorily. The displacement formulation, on the other hand, involves finding the three displacement functions simultaneously from three second-order partial differential equations of equilibrium. But solving for three functions under variously mixed conditions on the

42

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

Nomenclature x, y, z three dimensional Cartesian coordinate system E, m, G elastic modulus, Poisson’s ratio, modulus of rigidity of material ux, uy, uz displacement strain components in the x-, yand z-direction sxx, syy, szz normal stress components in the x-, y- and z-direction sxy, syz, szx shearing stress component in the xy, yz, and zx plane sr, sq radial and tangential components of stress in polar coordinates i, j, k nodal coordinate system in x-, y- and z-directions kx, ky, kz mesh lengths in the x-, y- and z-directions x kz/kx h ky/kx [A] (n!n) coefficient matrix {B} column matrix of constants {X} column matrix of variables surfaces is well-nigh impossible. As a result, no serious attempt is noticeable in the literature even to solve model problems of three-dimensional mixed-boundaryvalue elastic problems using these formulations. Accurate and reliable prediction of critical stresses is of great importance to meet the severe demand of greater reliability as well as economic challenges in design. It is noted that these critical stresses occur most frequently at the surfaces of structural components. Stress analysis of structural components is mainly handled by numerical methods. Extensive researches have been carried out in the recent years concentrating numerical solution of both twoand three-dimensional stress problems, especially with finite-element method (FEM). Of course, the adaptations of the FEM removed the major inability of managing complex boundary conditions and shapes, but the awareness of its lack of sophistication and doubtful quality of solutions, especially for the surface stresses, remained. The uncertainties associated with the finite-element prediction of stresses at the surfaces of structural components have been pointed out by several researchers [5,15,16]. On the other hand, introducing a new boundary modeling approach for finitedifference applications in solid mechanics, Dow et al. [5] verified that the accuracy of finite-difference method (FDM) in reproducing the state of stresses along the boundary was much higher than that of FEM analysis. They have, however, noted that the computational effort of the finite-difference analysis, under the new boundary modeling system, was even somewhat greater than that of FEM analysis. In a quest to develop an efficient numerical algorithm for the reliable and accurate solutions of boundary-value elastic problems, the

n j f ai

total number of unknowns displacement potential function Lame´’s strain potential function of Ref. [19] coefficients of the derivatives of j in the expressions of displacement components for three-dimensional problem bi coefficients of the derivatives of j in the expressions of displacement components for two-dimensional problem P0 intensity of radial loading for the circular disk problem (Case-II) R radius of the circular disk w uniformly distributed loading on the twodimensional beam (Case-I) s0 intensity of applied load for the three-dimensional beam (Case-III) L, D, B beam span, beam depth, beam width I moment of inertia of beam cross-section dmax maximum theoretical deflection of the beam (Case-III)

use of a new variable reduction scheme in conjunction with finite-difference method of solution have been investigated in the present research. The present paper describes the development of a new modeling approach for the finite-difference analysis of stresses and displacements in two- and three-dimensional elastic problems of solid mechanics. A new scheme of reduction of unknowns is used to develop the computational model. More specifically, the scheme reduces the problems to finding a single potential function of space variables, defined in terms of the respective displacement components of interest. Since, the present computational approach permits the reduction of variables to be evaluated at each nodal point to one, the total number of unknowns, and hence, the total number of algebraic equations to be solved is much less than that of the usual computational methods. Therefore, a tremendous saving in computational effort is achieved through the use of the present approach. Moreover, this reduction in computation will, in turn, improve the quality of solutions as well. Furthermore, one of the major factors of impediment to quality solutions of elastic problems is the treatment of the transition in boundary conditions. At the boundary, the practical problems are invariably subjected to the mixture of both known deformation and stresses. But none of the usual approaches would allow to account fully, both of these two types of boundary conditions with equal sophistication in the region of transition where boundary conditions change from one type to the other. It should be mentioned that the present modeling approach is capable of managing the mixed mode of the conditions as well as their zones of transition very efficiently. Finally, in an attempt to

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

verify the accuracy as well as reliability of the present singlefunction computational approach, a number of classical problems of solid mechanics are solved, and the corresponding finite-difference solutions are compared with those obtained by the standard method of solution as well as analytical solutions where available.

Substitution of the above expressions of ux and uy in Eqs. (1a) and (1b) gives the two equilibrium equations in terms of the function j(x,y) as: b1

2. Development of mathematical model for plane problems One of the ways to solve the two-dimensional stress problems of isotropic, elastic bodies undergoing small deformations under loading is to solve the displacement components ux and uy from the two equilibrium equations, as follows [1]

 v2 ux 1 C m v2 uy 1 K m v 2 ux C C Z0 (1a) 2 2 vxvy vx2 vy2

 v2 uy 1 C m v2 ux 1 K m v 2 uy C C Z0 2 2 vxvy vy2 vx2

  vuy E vux Cm syy Z vx 1 K m2 vy

(2b)

  E vux vuy C sxy Z 2ð1 C mÞ vy vx

(2c)

where, E is the modulus of elasticity of the material. The problem is now reduced to solving Eqs. (1a) and (1b) with specified values of the functions sxx, syy, sxy, ux, and uy on segments of boundary with all possible combination of them from segment to segment. Now, the two-dimensional problem of elasticity in terms of two unknown functions is reduced to the determination of a single function by introducing a new scheme of reduction of unknowns. In order to attain that, the displacement components are expressed as a linear combination of the second derivatives of a potential function j(x, y) as follows ux Z b 1

v2 j v2 j v2 j C b C b 2 3 vxvy vx2 vy2

(3a)

uy Z b 4

v2 j v2 j v2 j C b C b 5 6 vxvy vx2 vy2

(3b)

where the b’s are the unknown material constants for the two-dimensional problem.


 4 v4 j 1 Cm v j C b C b 2 4 4 2 vx vx3 vy 

 4 1 Km 1 Cm vj C b3 C b1 C b5 2 2 vx2 vy2 

 4 1 Km 1 Cm vj C b2 C b6 2 2 vxvy3
 1 K m v4 j C b3 Z0 2 vy4

(4a)




b4

(1b)

where, m is the Poisson’s ratio of the material. The stress components are then obtained from the following three stress–displacement relations [7]   vuy E vux Cm sxx Z (2a) vy 1 K m2 vx

43

 

 4 1 K m v4 j 1 Cm 1 Km v j C b C b 1 5 2 2 2 vx4 vx3 vy  
 
 4 1 Cm vj 1 Cm C b2 C b3 C b5 2 2 vxvy3
 4 4 1 Km vj vj C b6 4 Z 0 ð4bÞ Cb4 C b6 2 vx2 vy2 vy

The b’s are now chosen in such a way that one of the conditions of equilibrium (4a) and (4b) is automatically satisfied under all possible circumstances. This will happen only when the coefficients of all the derivatives of j present in that equation are individually zero. Let the coefficients of the derivatives of j in Eq. (4a) equate to zero, individually. Then, the resulting conditions are: b1 Z b3 Z b5 Z 0
b2 C b4


b2

1 Cm 2

(5a)

 Z0


 1 Km 1 Cm C b6 Z0 2 2

(5b)

(5c)

Thus, to be a solution of the two-dimensional stress problem j has to satisfy Eq. (4b) only. However, the values of b’s present in Eqs. (3a) and (3b) must be known in advance. It is seen that three of the six b’s in Eqs. (3a) and (3b) vanish automatically. The remaining three unknown b’s have to be known from the two Eqs. (5b) and (5c). Because there is one more unknown than the number of equations, there is no unique set of solution that satisfy the equations. However, by assuming a value for one of the b’s, the other two can be determined. Therefore, if an arbitrary value, for example, unity, is assigned to b2, the values of the unknown b’s can readily be obtained by solving Eqs. (5b)

44

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

and (5c), as follows: b2 Z 1 b4 Z K

2 1 Cm

b6 Z K

1 Km 1 Cm

(6)

(7)

It should be noted that assigning an arbitrary value to b2 is equivalent to assigning an arbitrary rigid-body initial state to the material body. Change in the value of b2 neither change the state of stress nor that of the strain in the body. Moreover, when the equilibrium Eq. (4b) is satisfied instead of Eq. (4a), i.e. if the coefficients of the individual derivatives of j in Eq. (4b) are equated to zero, the resulting set of solution for the unknown b’s are, b2 Z b4 Z b6 Z 0

 2   2  uy v uy v2 uy v2 uz 00 v ux C CG C 2 CG G Z0 vxvy vyvz vy2 vx2 vz

(9b)

 2   2  v2 uy uz v uz v2 uz 00 v ux C CG C 2 CG G Z0 vzvx vyvz vz2 vy2 vx

(9c)

0v

Substitution of the above values of b in Eq. (4b) gives the explicit expression of the governing differential equation for the solution of two-dimensional problems in terms of the potential function, j(x,y), as v4 j v4 j v4 j C2 2 2 C 4 Z 0 4 vx vx vy vy

problems of isotropic material. The three equilibrium equations in terms of the three displacement components, ux, uy and uz, are [7]  2   2  2 v ux v2 ux v2 uz 0 v ux 00 v uy C CG C 2 CG G Z0 (9a) vxvy vxvz vx2 vy2 vz

(8)

0v

2

2

where, GZ

E Eð1 KmÞ ; G0 Z and 2ð1 CmÞ ð1CmÞð1 K2mÞ

G00 Z

E : 2ð1CmÞð1 K2mÞ

The stresses at a point in a three-dimensional material body is represented by nine dependent variables, which can be obtained from the following relations [7]: sxx Z G0

vuy vux vu CM CM z vx vy vz

(10a)

syy Z M

vuy vux vu CG0 CM z vx vy vz

(10b)

szz Z M

vuy vux vu CM CG0 z vx vy vz

(10c)

b5 Z 1 1 Km b1 Z K 1 Cm b3 Z K

2 1 Cm

It can be easily verified that the substitution of the above set of b’s into the equilibrium Eq. (4a) leads to the identical governing differential equation as that of Eq. (7), thereby leaving the state of stress and strain unaffected within the material body. It is noted here that the reliability as well as suitability of the two-dimensional governing equation has been verified repeatedly through the finite-difference solutions of a number of practical problems with regular as well as irregular boundary shapes [11–14,17,18]. It is worth mentioning that the present potential function approach allows derivation of governing equation to be independent of material constants, as in the case of Airy’s stress function formulation [7]. However, the j-formulation is free from the inability of the stress function formulation in handling the mixed boundary conditions.

3. j-Formulation for three-dimensional problems In this section, the same scheme of reduction of unknowns is used to derive the governing differential equation for the solution of three-dimensional elastic



 vux vuy C Z syx vy vx

(10d)

 vuy vuz C syz ZG Zszy vz vx

(10e)

sxy Z G 

 szx ZG

 vuz vux C Zsxz vx vz

(10f)

where, MZ

mE ð1CmÞð1 K2mÞ

Likewise, the case of two-dimensional problems, the three displacement components are expressed have in terms of a potential function of space variables, j(x,y, z), as follows: ux Z a1

v2 j v2 j v2 j v2 j Ca2 2 Ca3 2 Ca4 2 vxvy vx vy vz

Ca5

v2 j v2 j Ca6 vyvz vzvx

(11a)

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

uy Z a 7

v2 j v2 j v2 j v2 j Ca Ca Ca 8 9 10 vxvy vx2 vy2 vz2

Ca11

uz Za13

v2 j v2 j Ca12 vyvz vzvx

v2 j v2 j Ca18 vyvz vzvx

(11b)

(11c)

Here, the a’s are unknown material constants for the three-dimensional problem. The corresponding three equilibrium equations in terms of the function j(x,y,z) are obtained by substituting Eqs. (11a)–(11c) into Eqs. (9a)–(9c). They are: v4 j v4 j v4 j v4 j 0 00 CGa CGa C½G a CG a  2 3 4 7 vx4 vy4 vz4 vx3 vy 4 4 4 vj v j vj C½G0 a6 CG00 a13  3 CGa5 3 CGa5 vx vz vy vz vyvz3 v4 j v4 j C½G0 a2 CGa1 CG00 a10  2 2 C½Ga4 CG00 a8  vx vy vxvy3 v4 j C½G0 a5 CG00 a12 CG00 a16  2 vx vyvz v4 j C½G0 a3 CGa1 CG00 a18  2 2 vx vz v4 j C½Ga6 CG00 a11 CG00 a14  vxvy2 vz v4 j v4 j C½Ga4 CG00 a9 CG00 a17  C½Ga6 CG00 a15  2 vxvyvz vxvz3 v4 j C½Ga3 CGa2  2 2 Z0 (12a) vy vz

G 0 a1

v4 j v4 j v4 j v4 j CG0 a8 4 CGa9 4 C½Ga8 CG00 a1  3 4 vx vy vz vx vy 4 4 vj vj CGa12 3 C½G0 a11 CG00 a14  3 vx vz vy vz v4 j v4 j C½Ga11 CG00 a15  C½G0 a7 CGa8 CG00 a4  2 2 3 vyvz vx vy 4 vj v4 j C½Ga11 CG00 a13 CG00 a6  2 C½Ga9 CGa7  2 2 vx vyvz vx vz 4 vj v4 j 0 00 00 C½G0 a8 CG00 a2  C½G a CG a CG a  12 5 16 vxvy3 vxvy2 vz 4 4 vj vj C½Ga10 CG00 a3 CG00 a18  CGa12 2 vxvyvz vxvz3 4 v j C½G0 a9 CGa8 CG00 a17  2 2 Z0 (12b) vy vz

Ga7

v4 j v4 j v4 j v4 j 0 CGa CG a CGa 14 15 16 vx4 vy4 vz4 vx3 vy 4 vj v4 j C½Ga18 CG00 a1  3 C½Ga17 CG00 a8  3 vx vz vy vz 4 v j v4 j 00 00 C½G0 a17 CG00 a9  C½Ga CG a CG a  17 4 7 vyvz3 vx2 vyvz v4 j v4 j C½Ga14 CGa13  2 2 C½G0 a13 CGa15 CG00 a6  2 2 vx vy vx vz 4 v j v4 j C½G0 a18 CG00 a3  C½G0 a16 CG00 a5 CG00 a12  3 vxvz vxvyvz2 v4 j v4 j C½G0 a14 CGa15 CG00 a11  2 2 CGa16 vy vz vxvy3 v4 j C½Ga18 CG00 a2 CG00 a10  Z0 (12c) vxvy2 vz

Ga13

v2 j v2 j v2 j v2 j Ca14 2 Ca15 2 Ca16 2 vxvy vx vy vz

Ca17

45

In this case, the a’s are chosen in such a way that two out of the above three conditions are automatically satisfied. Thus, all the coefficients of the individual derivatives of j in Eqs. (12a) and (12c) are equated to zero. Under this circumstance, combining the resulting thirty equations with each other, and utilizing the relation, GCG00 ZG0 , a total of seventeen equations are obtained for the evaluation of eighteen unknown a’s. Now, assigning an arbitrary value to one of the unknowns (here, taking a4Z1), the solution of the above set of seventeen homogenous equations yields, ai Z0

½iZ1;2;3;5;6;10;11;12;13;14;15;16;18

a4 Za17 Z1

(13)

a7 Za9 Z2ðmK1Þ a8 Z2mK1 Now, substituting the values of a’s from Eq. (13) into Eq. (12b), the explicit expression of the governing differential equation for the solution of three-dimensional problems is obtained in terms of the potential function j(x,y,z), which is given by, v4 j v4 j v4 j v4 j v4 j v4 j C 4 C 4 C2 2 2 C2 2 2 C2 2 2 Z 0 4 vx vy vz vx vy vy vz vz vx (14) The three-dimensional problem has thus been reduced to the solution of a single function j from a single differential equation of equilibrium (Eq. (14)). As mentioned in the twodimensional case, assigning an arbitrary value to a4 is equivalent to assigning a rigid-body initial state to the deformable material body. Change in the values of a4 does not change the state of stress or strain in the body. It can be mentioned here that Fung [19] presents the static equilibrium of elastic bodies in absence of body forces in terms of a different displacement potential function, f (Lame´’s strain potential), which has to satisfy the threedimensional Laplacian equation, V2fZconstant, where V2

46

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

is the Laplacian operator, V2 Z

v2 v2 v2 C 2C 2 2 vx vy vz

The relations between f, j, the stress components and the displacement components are given in Appendix A.

The physical conditions that exist on the surfaces of an elastic body are usually visualized in two different ways, namely, (a) the boundary restraints and (b) boundary loading, that is, known displacements and stresses on the boundary, respectively. Mixed-boundary-value problems are those in which the boundary conditions are specified as a mixture of boundary restraints and boundary loading, where the combination of the boundary conditions may also change from point to point on the boundary. Since, the present objective is to solve the problem in terms of the function j, all the displacement and stress components of interest are required to be expressed in terms of the function j. 4.1. Boundary conditions for two-dimensional problems The general expressions for the displacement components, ux and uy, in terms of j(x,y) have already been expressed by Eqs. (3a) and (3b) respectively. The explicit expressions of the components of displacement in terms of the function are obtained when the values of b’s from Eqs. (5a) and (6) are substituted in Eqs. (3a) and (3b), which are,

 2  1 v j v2 j uy ðx; yÞ Z K 2 2 C ð1 K mÞ 2 1 Cm vx vy

(15)

(16)

Combining Eqs. (2a)–(2c), (3a), (3b), (5a) and (6), the explicit expressions of the three stress components in terms of the potential function are obtained as follows:  E v3 j v3 j sxx ðx; yÞ Z b1 3 C ðb2 C mb4 Þ 2 2 1 Km vx vx vy  v3 j v3 j C mb C ðb3 C mb5 Þ 6 vxvy2 vy3  3  E vj v3 j Z Km 3 ð1 C mÞ2 vx2 vy vy

(17)

(19)

4.2. Boundary conditions for three-dimensional problems For the three-dimensional problems, the expressions of the three displacement components in terms of the function j(x,y,z) are obtained by combining the Eqs. (11a)–(11c) and (13), which are, ux ðx;y;zÞ Z

v2 j vxvy

(20)

  v2 j v2 j v2 j uy ðx;y;zÞ ZK 2ð1KmÞ 2 Cð1K2mÞ 2 C2ð1KmÞ 2 vx vy vz (21)

uz ðx;y;zÞ Z

2

vj vxvy

 E v3 j v3 j b4 3 C ðb1 C b5 Þ 2 2ð1 C mÞ vx vx vy  v3 j v3 j C b C ðb2 C b6 Þ 3 vxvy2 vy3   E v3 j v3 j K 3 Cm Z vx vxvy2 ð1 C mÞ2

(18)

sxy ðx; yÞ Z

4. Management of boundary conditions

ux ðx; yÞ Z

 E v3 j v3 j mb C ðb C mb Þ syy ðx; yÞ Z 1 4 2 1 K m2 vx3 vx2 vy  v3 j v3 j C b6 3 C ðb5 C mb3 Þ vxvy2 vy   E v3 j v3 j ð2 C mÞ ZK C ð1 C mÞ2 vx2 vy vy3

v2 j vyvz

(22)

Similarly, combining Eqs. (10a)–(10f), (11a)–(11c) and (13), the corresponding expressions of the six stress components in terms of the function j(x,y,z) are obtained as follows:   E v3 j v3 j v3 j ð1KmÞ 2 Km sxx ðx;y;zÞZ Km 3 (23) 1Cm vx vy vyvz2 vy  E v3 j v3 j ð2KmÞ 2 Cð2KmÞ syy ðx;y;zÞZK 1Cm vx vy vyvz2  v3 j Cð1KmÞ 3 vy

(24)

  E v3 j v3 j v3 j Km 2 Cð1KmÞ Km 3 szz ðx;y;zÞZ 1Cm vx vy vyvz2 vy (25)

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

  E v3 j v3 j v3 j Kð1KmÞ 3 Cm sxy ðx;y;zÞZ Kð1KmÞ 1Cm vx vxvy2 vxvz2 (26)   E v3 j v3 j v3 j Kð1KmÞ 2 Cm 2 Kð1KmÞ 3 syz ðx;y;zÞZ 1Cm vx vz vy vz vz (27) 

szx ðx;y;zÞZ

3

E vj 1Cm vxvyvz

normal normal normal normal

corresponding difference form when all the derivatives are replaced by their corresponding central difference expressions. The complete central difference expression of the two-dimensional fourth order bi-harmonic partial differential Eq. (7) is as follows: 2h2 ½ji;jC1 C jiC1;jK1 C ji1;jC1 C jiK1;j1  C ji;jC2 C ji;j2 K 4h2 ðh2 C 1Þ½jiC1;j C jiK1;j 

 (28)

One can thus realize that identical combination of derivatives would exist in each of the corresponding twoand three-dimensional expressions of different body parameters if the variations of the function with respect to z are neglected. The solution of the governing differential equations (Eqs. (7) and (14)) requires specification of two independent conditions at each nodal point of the bounding surfaces. The displacements and stresses at the boundary are usually specified in terms of their normal and tangential components. These normal and tangential components of displacement and stress will create six different combinations of boundary conditions when taken two of them at a time. Out of these six combinations, the possible four kinds of combination that may exist on the surfaces of practical stress problems are: (a) (b) (c) (d)

47

displacement and tangential displacement(s); displacement and tangential stress(es); stress and tangential stress(es); stress and tangential displacement(s).

In order to specify the physical conditions for a bounding surface of a three-dimensional body, basically one normal and two tangential components are required to be known at each nodal point on the surface. However, the two tangential components at a point can be combined to obtain their resultant component, so that each point on the surface will have two boundary conditions—one normal and the other tangential. Therefore, as far as numerical solution is concerned, it is evident from the expressions of the boundary conditions that all the boundary conditions of interest can readily be discretized in terms of the single function j by the method of finite-difference.

5. Finite-difference discretization of the governing equations The present computational scheme involves evaluation of a single function j at the mesh points of a uniform two- or three-dimensional rectangular mesh-network used to discretize the domain. The governing differential Eq. (7), used to evaluate the function j(x,y) only at the interior mesh points of a two-dimensional domain, can be expressed in its

K 4ðh2 C 1Þ½ji;jC1 C ji;j1  C h4 ½jiC2;j C jiK2;j  C ð6h4 C 8h2 C 6Þji;j Z 0

(29)

where, hZky/kx; kx, ky are the two mesh-lengths in x- and y-directions, respectively, as shown in Fig. 1. The corresponding grid structure of the difference Eq. (29) is shown in Fig. 1, which involves thirteen neighboring mesh points including the point of application. The finite-difference approximation of the three-dimensional governing Eq. (14) is obtained in a similar fashion when all the normal and cross derivatives are replaced by their corresponding central difference formulae. The complete central-difference approximation of the fourth order partial differential Eq. (14) for the solution of three-dimensional problems is as follows: h4 x4 ½jiC2;j;k CjiK2;j;k Cx4 ½ji;jC2;k Cji;jK2;k  Ch4 ½ji;j;kC2 Cji;j;kK2  K4ðh4 x4 Ch2 x4 Ch4 x2 Þ½jiC1;j;k CjiK1;j;k  K4ðx4 Ch2 x4 Ch2 x2 Þ½ji;jC1;k Cji;jK1;k  K4ðh4 Ch2 x2 Ch4 x2 Þ½ji;j;kC1 Cji;j;kK1  C2h2 x4 ½jiC1;jC1;k CjiK1;jK1;k CjiK1;jC1;k CjiK1;jK1;k  C2h2 x2 ½ji;jC1;kC1 Cji;jC1;kK1 Cji;jK1;kC1 Cji;jK1;kK1  C2h4 x2 ½jiC1;j;kC1 CjiK1;j;kC1 CjiC1;j;kK1 CjiK1;j;kK1  n o C 6ðh4 Ch4 x4 Cx4 ÞC8ðh2 x4 Ch2 x2 Ch4 x2 Þ ji;j;k Z 0 ð30Þ where, xZkz/kx; kz is the mesh length in z-direction (see Fig. 2). The corresponding three-dimensional grid structure of the governing Eq. (30) at an arbitrary internal mesh point (i, j, k) is presented in a simplified layer-wise representation of a threedimensional mesh-network, as shown in Fig. 2. It is, thus, clear that application of the governing differential Eq. (14) involves a total of twenty five neighboring mesh points including the point of application (i, j, k). One can easily realize the extension of the two-dimensional stencil to the case of threedimensional problems. The limitation and complexity associated with analytical solutions ultimately leads to the conclusion that a numerical modeling for the practical problems of stress analysis is the plausible approach. The finite-difference technique is used here to transfer the governing partial differential equations

48

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

Physical boundary

y j

j

i kx

i

values of the potential function, j at the mesh points of the domain concerned can be obtained from a system of linear algebraic equations resulting from the discretization of the governing equation and the associated boundary conditions. It is noteworthy that, in the present approach, each mesh point of both the two- and threedimensional computational domain has a single algebraic equation for the evaluation of the corresponding unknown. Finally, all the parameters of interest in the solution of elastic problems, namely, displacement and stress components, can readily be obtained as soon as j is known at the mesh points of the domain, as all the components are expressed as summation of different derivatives of the function j.

ky x

6. Reduction of the computational work Point of application (i, j)

Fig. 1. Application of the finite-difference stencil of the governing Eq. (7) at two typical internal nodal points of a two-dimensional body.

and also the second and third order partial differential equations associated with the boundary conditions, into their corresponding algebraic forms. In order to avoid the inclusion of points exterior to those involved by the application of the governing equilibrium equation at the interior nodal points, different finite difference forms has been used for each boundary conditions, which includes different combination of forward, backward and central-difference schemes. The computer program is then organized in such a fashion that it will automatically select the most appropriate form of the boundary conditions to be satisfied at a point on the boundary depending on the actual position of the point with respect to the global reference coordinate system. The discrete

This section explains, how the computational work is reduced due to the reduction in the number of parameters to be evaluated at the nodal points. In case of a numerical approach, the problem is reduced to the solution of a set of algebraic equations. The unknowns are the values of parameters at the nodes. The equation is ½An;n fXgn Z fBgn

(31)

To solve for n unknowns of Eq. (31), required computations are as follows: (a) Division: nC

nK1 X

ið3 C iÞ

iZ1

(b) Subtraction: nK1 X

ið3 C iÞ

iZ1

For a two-dimensional problem, the usual computational methods involve finding two parameters (like, ux, uy) at each nodal point. Here, in the present approach, only a single parameter, j has to be evaluated at each nodal point. So, the reduction of the unknowns between the usual and the present approach is from 2n to n, respectively. The equation for solving 2n unknowns, as required for the case of the usual approaches is: ½A2n;2n fXg2n Z fBg2n

(32)

The corresponding computations for finding the 2n unknowns of Eq. (32) are, (a) Division: Fig. 2. Application of the finite-difference stencil of the governing Eq. (14) at an arbitrary internal nodal point (i, j, k) of a three-dimensional body.

2n C

2X nK1 iZ1

ið3 C iÞ

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

P2nK1

(b) Subtraction: 2nK1 X

ið3 C iÞ

The increase in computation due to the additional n unknowns is )

2nK1 X

(

nK1 X 2n C ið3 C iÞ K n C ið3 C iÞ iZ1 iZ1 ( ) ( ) 2nK1 nK1 X X ið3 C iÞ K ið3 C iÞ iZ1

) ½Division ½Subtraction

iZ1

Therefore, the percentage saving in the number of division

2n C

Z

PnK1 ið3 C iÞ K iZ1 ið3 C iÞ !100 z87% P2nK1 iZ1 ið3 C iÞ

It can be shown in a similar fashion that the percentage saving in the computational work is about 96% for the case of a three-dimensional problem, where the usual computational methods use at least three parameters (for example, ux, uy, uz) at each nodal point of the computational domain. Therefore, the present j-formulation in conjunction with the finite-difference method of solution would save a significant amount of computational effort for the solution of a problem, in comparison with that required for the usual computational approaches.

iZ1

(

iZ1

Z

49

7. Examples of application of the single function computational approach

  P2nK1 PnK1 iZ1 ið3 C iÞ K n C iZ1 ið3 C iÞ

 !100 P 2n C 2nK1 iZ1 ið3 C iÞ

Using the present computational approach, two classical plane problems and a three-dimensional problem of solid mechanics have been solved numerically and the results are compared with those obtained by the standard method of solution as well as analytical solutions where available. The two plane problems (Case-I and II) considered are shown in Fig. 3(a) and (b). Case-I is an example of a deep cantilever

z87% and the percentage saving in the number of subtraction

Y (a)

(b) Uniformly distributed load, w y Po D

R O

L x

(c)

z

Fixed support

Uniformly distributed pressure, σ0

D y B L x Fig. 3. Loading and geometry of (a) two-dimensional deep cantilever beam (Case-I), (b) Solid circular disk subjected to uniform radial pressure (Case-II), (c) three-dimensional both ends fixed deep beam (Case-III).

50

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

beam of rectangular cross section with aspect ratio, L/DZ2, loaded uniformly over its whole span. Case-II is an example of a classical stress problem with curved boundary, namely, a solid circular disk, subjected to uniform radial pressure on the circular boundary. The three-dimensional problem (Case-III) is a both ends fixed deep beam (L/DZ2, L/BZ 2), loaded uniformly over the whole region of its upper surface, as shown in Fig. 3(c). The finite-difference mesh-

networks used for the solution of Case-I and II are shown in Fig. 4(a) and (b), respectively, where the exterior field nodal points are involved by the application of the governing equation at the internal nodal points in the neighborhood of the physical boundary. For the solution of the threedimensional problem, a three-dimensional rectangular mesh-network, similar to that shown in Fig. 2 has been used. In order to obtain, the numerical results of the

Physical boundary

(a)

Exterior boundary

y, j

kx

ky

x, i

Physical Boundary Exterior Nodal Point Interior Nodal Point

(b)

15.75 14.00 12.25

Y-coordinate

10.50 8.75 7.00 5.25 3.50 1.75 0.00 0

2

4

6

8

10

12

14

16

X-coordinate Fig. 4. Finite-difference discretization of (a) the cantilever beam, (b) the circular disk (field nodal points do not coincide with the boundary nodal points).

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

problems, values for Young’s modulus and Poisson’s ratio of 209 GPa and 0.3, respectively, are used. 7.1. Results of the two-dimensional cantilever beam problem (Case-I) Numerical solution for the stresses has been obtained using a relatively small number of mesh points (23!23), so that the solutions can be readily compared with the FEM solutions of identical mesh points. Here, the stresses are

(b) 12 Elementary theory FEM FDM

20 10 0 –10 –20 –30

normalized with respect to the uniform load-parameter, w. The present numerical solutions are compared with the solutions of elementary beam theory [20] and also with the finite-element solutions. The FEM solutions were obtained by using the standard finite-element code, ANSYS. Two hundred, 8-noded, plane stress, isoparametric elements (20!10) were used in obtaining the FEM solutions. Fig. 5(a) shows the bending stress distributions in the beam at the fixed end, y/LZ0, and Fig. 5(b) at a section, y/LZ0.1. The three solutions are presented in the same

0.0

0.2

0.4

0.6

0.8

Normalized bending stress, σyy / w

Normalized bending stress, σyy / w

(a) 30

4 0 –4 –8 –12

1.0

(d) 0 –2 –4 –6 Elementary theory FEM FDM

–8

0.0

0.2

0.4

0.6

0.8

2 0 –2 –4 –6

0.0

0.2

0.4

0.6

0.8

1.0

Normalized distance along depth, x / D

Normalized normal stress, σxx / w

Normalized normal stress, σxx / w

4

0.6

0.8

1.0

Elementary theory FEM FDM

0 –1 –2

0.0

0.2

0.4

0.6

0.8

1.0

Normalized distance along depth, x / D (f)

FEM FDM

0.4

1

–3

1.0

6

0.2

2

Normalized distance along depth, x / D (e)

0.0

Normalized distance along depth, x / D

Normalized shear stress, σxy / w

Normalized shear stress, σxy / w

(c)

Elementary theory FEM FDM

8

Normalized distance along depth, x / D

–10

51

2 1 0 –1 –2 –3

FEM FDM

0.0

0.2

0.4

0.6

0.8

1.0

Normalized distance along depth, x / D

Fig. 5. Comparison of different solutions at different sections of the cantilever beam (Case-I): (a) bending stresses at the fixed end, (b) bending stresses at section, y/L=0.1, (c) shearing stresses at the fixed end, (d) shearing stresses at section, y/L=0.1, (e) normal stresses at the fixed end, (f) normal stresses at section, y/L=0.1.

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

7.2. Results of the circular disk problem (Case-II) A classical axi-symmetric stressed body with curved boundary, whose exact analytical solution is known, has been chosen as the second example to demonstrate the validity of the present computational approach based on j-formulation. It also demonstrates the ability of the j-formulation to solve the problems of curved boundary by finite-difference method. A (31!35) uniform rectangular mesh-network is used to solve the disk problem. Since the physical boundary of the circular disk does not pass through the field nodal points, the principle of interpolation [17,18] is used to transfer the values of parameters specified at the boundary points, to the neighboring four field nodal points. The technique works very well for the management of irregular boundary shapes as demonstrated here through

the solution of the disk problem subjected to uniform radial compression along its circumference. The present j-solutions are compared with the available analytical solutions of the problem. The analytical solutions of the problem is symmetrical with respect to the axis of the disk and both the radial stress (sr) and the circumferential stress (sq) are equal to P0; the disk is thus in a condition of uniform compression in all the directions [7]. Also, the shear stress is zero all over the plate. The radial stress, sr along a line OY from the center O to the circumference of the disk (Fig. 3(b)) as obtained from the j-formulation, is shown and compared with the analytical results in Fig. 6(a). From the figure, it is evident that the present numerical results are in very good agreement with the exact solutions. The solution for the shearing stresses along the same radial line is compared in Fig. 6(b). The result also conforms to the analytical solutions very well. 7.3. Results of the three-dimensional beam problem (Case-III) A three-dimensional uniform rectangular mesh network (23!13!13) is used to obtain the numerical solution for the present three-dimensional both ends fixed beam (a) –1.2 Normalized radial stress, σr / P0

graph for the convenience of comparison. The comparison shows that the bending stress solution of FDM agree in general, reasonably well with the FEM solution, and also with those of elementary flexure theory. At the fixed support, the numerical solutions differ substantially, as expected, from the predictions of linear theory. The two numerical solutions, however, differ mainly at the corner region of the fixed support; the maximum bending stress predicted by the FDM approach is found to be slightly higher than that predicted by FEM. The non-linearity in the distribution as well as the discrepancy between the solutions is found to decrease with increasing distance from the fixed support, as seen in Fig. 5(b). The distribution of shearing stress at the fixed support and at section, y/LZ0.1, of the cantilever beam is shown in Fig. 5(c) and (d), respectively. The parabolic distributions of the elementary solution and the FEM solutions are included in the figures. The solutions of elementary theory differ significantly from the numerical solutions, especially at the fixed support. Moreover, the comparative analysis shows that the finite element prediction of shear stress in the neighborhood of both top and bottom surfaces are not very reliable; the discrepancy is found to remain quite significant even for sections away from the fixed end [see Fig. 5(d)]. It is noted that the present FDM method can correctly predict zero shear stress at the surfaces, including the corner points of the beam. Fig. 5(e) and (f) shows the comparison of solutions for the normal stress component sxx at sections, y/LZ0.0, and 0.1, respectively. Here, only the two numerical solutions are compared, as the elementary theory cannot predict this normal stress in the beam. It is seen from the comparison that, away from the surfaces, two solutions are identical. However, in the neighborhood of both the surfaces, finiteelement prediction of the normal stress is unreliable. The discrepancy is found to become more pronounced towards the fixed end, the most critical section in the beam. Note that the present method of solution reproduces the normal stress exactly both at the top as well as bottom surfaces of the beam.

FDM Analytical –1.1

–1.0

–0.9

–0.8 1.0

1.2

1.4

1.6

1.8

2.0

2.2

Normalized radial distance, Y / R (b) Normalized shear stress, σxy / P0

52

0.4 FDM Analytical 0.2

0.0

–0.2

–0.4 1.0

1.2

1.4

1.6

1.8

2.0

2.2

Normalized radial distance, Y / R Fig. 6. Comparison of solutions for the circular disk (Case-II): (a) radial stress, (b) shearing stress.

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

problem. Some of the results are presented here for verifying the soundness as well as reliability of the present three-dimensional j-formulation. Displacements are normalized with respect to the maximum theoretical deflection of the beam as predicted by the elementary beam theory [20], which is given by

(b)

1 z / D=1.0 z / D=0.8 z / D=0.6 z / D=0.4 z / D=0.2 z / D=0.0

0 –1 –2 –3 –4 –5 –6

0.0

0.2

0.4

0.6

0.8

–1 –2 –3 –4 –5

Normalized bending stress, σyy / σ0

Normalized bending stress, σyy / σ0

2

0

–2 z / D=0.2 z / D=0.1 z / D=0.0

–4

0.0

0.2

0.4

0.6

0.8

0.6

0.8

1.0

z / D=1.0 z / D=0.9 z / D=0.8 z / D=0.5

0

–2 z / D=0.2 z / D=0.1 z / D=0.0

0.0

0.2

0.4

0.6

0.8

1.0

Normalized distance along length, y / L (f) 0.30

x / B=0.1 x / B=0.1 x / B=0.2 x / B=0.3 x / B=0.4 x / B=0.5

0.20 0.10

Normalized shear stress, σxy / σ0

Normalized shear stress, σxy / σ0

0.4

2

Normalized distance along length, y / L

0.00 –0.10 x / B=0.6 x / B=0.7 x / B=0.8 x / B=0.9 x / B=1.0

–0.20 –0.30

0.2

4

–4

1.0

(e) 0.30

0.0

Normalized distance along length, y / L (d)

z / D=1.0 z / D=0.9 z / D=0.8 z / D=0.5

z / D=1.0 z / D=0.8 z / D=0.6 z / D=0.4 z / D=0.2 z / D=0.0

0

–6

1.0

4

(33)

1

Normalized distance along length, y / L (c)

s0 L4 384 EI

where, s0 is the uniform pressure acting on the top surface of the beam, I the moment of inertia of the beam crosssection, and L stands for the beam span. The stress

Normalized displacement, uz / δmax

Normalized displacement, uz / δmax

(a)

dmax Z

53

0.0

0.2

0.4

0.6

0.8

1.0

Normalized distance along length, y / L

x / B=0.1 x / B=0.1 x / B=0.2 x / B=0.3 x / B=0.4 x / B=0.5

0.20 0.10 0.00 –0.10

x / B=0.6 x / B=0.7 x / B=0.8 x / B=0.9 x / B=1.0

–0.20 –0.30

0.0

0.2

0.4

0.6

0.8

1.0

Normalized distance along length, y / L

Fig. 7. Solutions of the three-dimensional fixed beam (Case-III): (a) FDM solution of the displacement component, uz at different vertical sections of the surface, x/B=1, (b) FEM solution of the displacement component, uz at different vertical sections of the surface, x/B=1, (c) FDM solution of the bending stress component, syy at different vertical sections of the surface, x/B=0, (d) FEM solution of the bending stress component, syy at different vertical sections of the surface, x/B=0, (e) FDM solution of the shearing stress component, sxy at different lateral sections of the surface, z/D=1, (f) FEM solution of the shearing stress component, sxy at different lateral sections of the surface, z/D=1.

54

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

components are normalized directly with the load-parameter, s0. The present j-solutions are compared with the corresponding finite-element solutions obtained by the standard code, ANSYS. In obtaining the FEM solutions, 8-noded, three-dimensional, isoparametric solid elements (20!10!10) are used. Fig. 7(a) shows the distribution of the vertical displacement component, uz at different vertical sections of the surface, x/BZ1. As expected, the displacements at the fixed surfaces are zero and the maximum deflection occurs at the mid-section of the surface. Moreover, the deflection at the loaded side is found to be higher than that at the bottom surface, thereby conforming to the physical characteristics of the beam appropriately. Another important thing to be noted here is that the maximum vertical deflection of the beam is around 4.5 times higher than that predicted by the elementary beam theory. This may be due to the fact that the elementary beam theory is simply inadequate for the analysis of such a three-dimensional deep beam. The corresponding FEM solution is presented in the similar fashion in Fig. 7(b); one can easily compare the two solutions both qualitatively and quantitatively. It is observed that both the FEM and FDM solutions of displacements at different sections of the beam are found to be almost identical in nature and magnitude. The solutions for the bending stress component, syy at different vertical sections of the surface, x/BZ0, are presented in Fig. 7(c). The distributions clearly reflect the fact that the bending stresses at the top and bottom surfaces are the maximum, but are of opposite sense. The distributions of bending stress are observed to be almost linear except those at the fixed supports, when analyzed in the perspective of beam depth. Fig. 7(d) shows the corresponding finite element solutions; the two solutions compare well with each other, showing slight discrepancy near the two edges of the fixed support. More specifically, the two solutions are found to differ mainly around the neighbourhood of the corner zones of fixed supports, which are, in general, the zones of singularity. The obtained j-solutions of the shearing stress component, sxy for the loaded surface, z/DZ1, are presented in Fig. 7(e). The distribution of the stress component is found to be symmetrical about both the horizontal and lateral axes of the plane. Although the magnitude of the shearing stress component is not so significant compared to that of the bending stress, the solutions obtained by the two different methods are almost identical in shape and magnitude, as seen in Fig. 7(e) and (f). Slight discrepancies are encountered mainly at the edges of the restrained surfaces. Both the qualitative and quantitative results of the example problems, obtained through the present potential function approach, establish the soundness as well as appropriateness of the present variable reduction scheme. Moreover, the comparative analysis with the available

solutions verifies that the present j-formulation in conjunction with finite-difference method of solution is capable of producing accurate and reliable solutions at any critical section, either at or far from the bounding surfaces.

8. Conclusions The development of an efficient algorithm for the numerical solution of both two- and three-dimensional elastic problems has been described in the present paper. The philosophy of the present approach is, the problems are formulated in terms of a single potential function, defined in terms of the displacement components. It has been shown that the present computational approach of stress analysis would save the computational effort significantly in comparison to that taken by the usual computational methods. In addition, the quality of the solutions will be greatly improved because of the fact that the number of unknowns, and hence, the total number of algebraic equations to be solved is much less than that of existing approaches. The special feature of the present approach over the existing models is its ability in satisfying the boundary conditions exactly, whether they are specified in terms of loadings or physical restraints or any combination thereof. It should be remembered here that presently the major inability to solve engineering problems arises from the fact that round-off errors in the computational process of solving linear algebraic equations makes the results spurious when the number of unknowns in the algebraic equations becomes very large. There is no doubt that the present potential function approach would be the most promising scheme in dealing with the mixedboundary-value elastic problems of solid mechanics as well as in other engineering applications. The saving in the computational effort through the use of the potential function and the method of finite-difference can be extended to other fields of science and engineering where solutions are sought through computational methods.

Appendix A The relationships between the Lame´’s strain potential f [19], the displacement potential j of the present study, the stress and the displacement components of elastic bodies are all contained in this appendix for ready reference of interested readers. These relations are as follows:

ux Z

v2 j 1 vf Z vxvy 2G vx

(A1)

M.Z. Hossain et al. / Advances in Engineering Software 37 (2006) 41–55

  v2 j v2 j v2 j uy ZK 2ð1 K mÞ 2 C ð1 K 2mÞ 2 C 2ð1 K mÞ 2 vx vy vz 1 vf ðA2) Z 2G vy uz Z

v2 j 1 vf Z vyvz 2G vz

sxx Z

(A3)

  E v3 j v3 j v3 j v2 f ð1 K mÞ 2 K m K m Z 2 2 3 1 Cm vx vy vyvz vy vx (A4)  E v3 j v3 j ð2 K mÞ 2 C ð2 K mÞ 1 Cm vx vy vyvz2  v3 j v2 f Cð1 K mÞ 3 Z 2 vy vy

syy ZK

(A5)

  E v3 j v3 j v3 j v2 f Km 2 C ð1 K mÞ szz Z K m Z 1 Cm vx vy vyvz2 vy3 vz2 (A6)   E v3 j v3 j v3 j Kð1 K mÞ 3 C m sxy Z K ð1 K mÞ 1 Cm vx vxvy2 vxvz2 Z

v2 f vxvy

ðA7)

  E v3 j v3 j v3 j syz Z Kð1 K mÞ 2 C m 2 K ð1 K mÞ 3 1 Cm vx vz vy vz vz Z

szx Z

v2 f vyvz

 3  E vj v2 f Z 1 C m vxvyvz vzvx

ðA8)

(A9)

References [1] Uddin MW. Finite-difference solution of two-dimensional elastic problems with mixed boundary conditions. MSc Thesis, Carleton University, Canada; 1966.

55

[2] Horgan CO, Knowels JK. Recent developments concerning saint Venant’s principle. Adv Appl Mech 1983;23:179–269. [3] Durelli AJ, Ranganayakamma B. On the use of photoelasticity and some numerical methods. Photomech Speckle Metrol, SPIE 1987; 814:1–8. [4] Durelli AJ, Ranganayakamma B. Parametric solution of stresses in beams. J Eng Mech 1989;115(2):401–15. [5] Dow JO, Jones MS, Harwood SA. A new approach to boundary modeling for finite difference applications in solid mechanics. Int J Numer Methods Eng 1990;30:99–113. [6] Hardy SJ, Pipelzadeh MK. Static analysis of short beams. J Strain Anal 1991;26(1):15–29. [7] Timoshenko SP, Goodier JN. Theory of Elasticity. 3rd ed. New York: McGraw-Hill Book Company; 1979. [8] Idris ABM. A new approach to solution of mixed boundary-value elastic problems. MSc Thesis, Bangladesh University of Engineering and Technology, Bangladesh; 1993. [9] Ahmed SR. Numerical solutions of mixed boundary-value elastic problems. MSc Thesis, Bangladesh University of Engineering and Technology, Bangladesh; 1993. [10] Akanda MAS. Stress analysis of gear teeth. MSc Thesis, Bangladesh University of Engineering and Technology, Bangladesh; 1997. [11] Idris ABM, Ahmed SR, Uddin MW. Analytical solution of a 2-D elastic problem with mixed-boundary conditions. J Inst Eng (Singapore) 1996;36(6):11–17. [12] Ahmed SR, Khan MR, Islam KMS, Uddin MW. Analysis of stresses in deep beams using displacement potential function. J Inst Eng (India) 1996;77:141–7. [13] Ahmed SR, Idris ABM, Uddin MW. Numerical solution of both ends fixed deep beams. Comput Struct 1996;61(1):21–9. [14] Ahmed SR, Khan MR, Islam KMS, Uddin MW. Investigation of stresses at the fixed end of deep cantilever beams. Comput Struct 1998;69:329–38. [15] Richards TH, Daniels MJ. Enhancing finite element surface stress predictions: a semi-analytic technique for axisymmetric solids. J Strain Anal 1987;22(2):75–86. [16] Smart J. On the determination of boundary stresses in finite elements. J Strain Anal 1987;22(2):87–96. [17] Akanda MAS, Ahmed SR, Khan MR, Uddin MW. A finite-difference scheme for mixed boundary-value problems of arbitrary-shaped elastic bodies. Adv Eng Software 2000;31(3):173–84. [18] Akanda MAS, Ahmed SR, Uddin MW. Stress analysis of gear teeth using displacement potential function and finite-differences. Int J Numer Methods Eng 2001;53(7):1629–49. [19] Fung YC. Foundations of solid mechanics. Englewood Cliffs, NJ: Prentice-Hall; 1965. [20] Timoshenko SP, Maccullough GH. Elements of strengths of materials. Princeton, NJ: Van Nostrand; 1949.