- Email: [email protected]
A general mathematical formulation for finite-difference solution of mixed-boundary-value problems of anisotropic materials
Computers and Structures 83 (2005) 35–51 www.elsevier.com/locate/compstruc
A general mathematical formulation for finite-difference solution of mixed-boundary-value problems of anisotropic materials S. Reaz Ahmed *, M. Zubaer Hossain, M. Wahhaj Uddin Department of Mechanical Engineering, Bangladesh University of Engineering and Technology, Dhaka 1000, Bangladesh Received 2 July 2003; accepted 9 August 2004
Abstract This paper presents a general mathematical model, especially suitable for finite-difference analysis of stresses and displacements of the plane elastic problems of solid mechanics. The present formulation covers the problems of anisotropic, orthotropic and isotropic materials, in which the problem is formulated in terms of a single potential function, defined in terms of the displacement components. In addition, the formulation contains a new scheme of reduction of unknowns to be solved for a particular problem. Compared to the conventional computational approaches, the present scheme gets solution of higher accuracy and in extremely shorter time. The application of the present scheme is demonstrated here through a classical problem of solid mechanics, and the results are compared with the available solutions in the literature. 2004 Elsevier Ltd. All rights reserved. Keywords: Stress analysis; Displacement potential function; Anisotropic materials; Finite-difference method
1. Introduction Analysis of stress in a material body is usually a three-dimensional problem. However, in most practical cases, the problem can reasonably be approximated to either a two-dimensional or a one-dimensional one. One-dimensional problems of stress analysis are very straightforward and included in undergraduate curricula. But the problems of two-dimensional stress analysis have a different story. In the case of very simple structures like beam, column, circular disc, regular shaped plate and thin shells, having simple geometry and boundary conditions, analytical methods are used in *
Corresponding author. E-mail address: [email protected] (S.R. Ahmed).
their stress analysis. However, the analytical approach is limited, because general closed form solution can be obtained only for very ideal cases. Moreover, the analytical results obtained for the practical problems of elasticity, usually with mixed boundary conditions, are invariably approximate as they include various idealizations. Among the existing mathematical models for twodimensional elastic problems, the stress function approach [1] and the displacement functions approach (ux, uy-formulation, [2]) are noticeable. Successful application of the stress function formulation in conjunction with the finite-difference technique has been reported for the solution of plane elastic problems where the conditions on the boundary are prescribed in terms of stresses only [3–5]. Further, Conway and Ithaca [6] extended the stress function formulation in the form of Fourier
0045-7949/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2004.08.007
36
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
Nomenclature x, y, z 1, 2, 3 E11 E22 l12 l21 G12 E, l, G
orthogonal loading axes principle material axes elastic modulus of the material in 1-direction elastic modulus of the material in 2-direction major PoissonÕs ratio minor PoissonÕs ratio shear modulus in the 1–2 plane elastic modulus, PoissonÕs ratio, modulus of rigidity of isotropic material h fiber orientation ½Q transformed material stiffness matrix Qij coefficients of the transformed stiffness matrix Qij coefficients of the stiffness matrix ux, uy, uz displacement components in the x-, y- and z-direction un, ut normal and tangential component of displacement on the boundary rxx, ryy normal stress components in the x- and ydirection rxy shearing stress component in the xy plane exx, eyy normal strain components in the x- and ydirection cxy shearing strain component in the xy plane rn, rt normal and tangential components of stress on the boundary
integrals to the case where the material is orthotropic, and obtained analytical solutions for a number of problems. Mixed-boundary-value problems are those in which the boundary conditions are specified as a mixture of boundary restraints and boundary loadings. The shortcoming of the stress function formulation is that it accepts boundary conditions in terms of loadings only. Boundary restraints specified in terms of ux and uy cannot be satisfactorily imposed on the stress function. As most of the practical problems of elasticity are of mixed boundary conditions, the approach fails to provide any explicit understanding of the stress distribution in the regions of restrained boundaries, which are, in general, the most critical zones in terms of stresses. The difficulties involved in trying to solve practical stress problems using the stress function approach are pointed out in our previous researches [7–10], and also by Durelli [11]. Two-dimensional stress analysis in practice is mainly handled by numerical methods. The major numerical methods in use are (a) the method of finite-difference and (b) the method of finite-element. Finite-difference method is an ideal numerical approach for solving par-
n
normal component of the force on the boundary segment t tangential component of the force on the boundary segment / angle between the unit positive normal to the boundary and positive x-axis X, Y components of stress on the boundary, along the direction of coordinate axes a length of boundary segment, AB i, j nodal coordinate system in x- and ydirections h, k mesh lengths in the x- and y-directions n k/h [A] coefficient matrix {B} constant matrix m order of the coefficient matrix P number of parameters at each node N total no of nodes w displacement potential function ai coefficients of the derivatives of w in the expressions of displacement components bi coefficients of the general governing differential equation of w Dn coefficient composed of elastic constants of orthotropic material w uniformly distributed loading on the beam L, D, B beam span, beam depth, beam width
tial differential equations. The difference equations that model the governing equations are simple to code and the global coefficient matrix possesses a banded structure that is amenable to efficient solution. Despite these ideal characteristics, finite-difference method has been supplanted in most solid mechanics applications in engineering by the finite-element method. This is mainly because of the limited flexibility of the method in managing the boundary conditions and shapes. However, the supremacy is brought back to the finite-difference technique here from the finite-element method through a new formulation of stress analysis of solid structures. As the new formulation establishes a priority of finite-difference technique over finite-element method, the philosophy and approach of the two methods are recapitulated here in brief.
1.1. Finite-element method Disadvantages: (a) The body is not in one piece but an assemblage of elements connected only at the nodes.
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
(b) Variations of parameters over individual elements are assumed to be very simple as in polynomials of limited order. (c) Values of parameters across element boundary are usually discontinuous. (d) Number of parameters (like ux, uy, uz) to be evaluated at each nodal point is usually more than that in finite-difference. (e) The total number of unknowns to be solved is much more than that in finite-difference method and hence the solution is likely to be crude and the computational time is more. (f) It cannot handle singularity. Advantages: (a) It works when other methods fail. (b) It is very good in managing boundaries.
complicated
1.2. Finite-difference method Disadvantages: (a) Derivatives in the governing equations are replaced by their finite-differences, resulting in some approximation. (b) Complicated boundary shapes and conditions are very difficult, almost impossible to handle. (c) It cannot handle singularity very well. Advantages: (a) The body is in one piece and thus there is no approximation on this account. (b) It permits reduction of parameters to be evaluated at the mesh points to one. (c) Total number of unknowns is less than that in finiteelement and hence, solution is more accurate and the computational time is less.
37
was much higher than that of FEM analysis. However, they noted that the computational effort of the finite-difference analysis, under the new boundary modeling system, is even somewhat greater than that of FEM analysis. The present paper describes the development of a general formulation for the finite-difference analysis of stresses and displacements in two-dimensional elastic solids. Basically, the present research is an attempt to extend the capability of our displacement potential formulation to the problems of composite materials. The formulation covers all the cases of general anisotropic, orthotropic and isotropic materials. Specifically, the formulation has been developed in such a fashion that the governing equations for the case of orthotropic and isotropic materials can be obtained as a special case of the general formulation for anisotropic materials. More importantly, the present scheme reduces the problem to finding a single parameter at the nodal points whereas the existing schemes find more than one at each node, and hence a tremendous saving in computational work is achieved. In addition, the present formulation provides the way to handle the mixed mode of boundary conditions very efficiently. Finally, the application of the present scheme is demonstrated for a classical problem of solid mechanics, and the solutions are compared with those of the existing methods.
2. Development of the general governing equations One of the ways to solve the stress problems of twodimensional elastic bodies of HookÕs material is to solve the displacement components ux and uy from the equilibrium equations. With reference to a rectangular coordinate system, in absence of body forces, the two differential equations of equilibrium for general anisotropic materials, in terms of the displacement components, are as follows [1,15]: o2 ux o2 ux o2 ux o2 uy þ Q66 2 þ Q16 2 þ 2Q16 2 ox oxoy oy ox o2 uy o2 uy þ Q12 þ Q66 þ Q26 2 ¼ 0 oxoy oy
Q11 It is noted that reliable and accurate prediction of critical stresses, mostly at the surfaces of engineering components, is of great importance to meet the severe demand for greater reliability as well as economic challenges in design. However, the uncertainties associated with the finite-element prediction of surface stresses have been clearly pointed out by several researchers [12–14]. Dow et al. [14] introduced a new boundary modeling approach for the finite-difference applications of displacement formulation of solid mechanics. In this connection, they reported that the accuracy of the finite-difference method in reproducing the state of stresses along the bounding surfaces
Q16
ð1aÞ
o2 ux o2 ux o2 ux o2 uy þ Q12 þ Q66 þ Q26 2 þ Q66 2 2 ox oxoy oy ox
þ 2Q26
o2 uy o2 uy þ Q22 2 ¼ 0 oxoy oy
ð1bÞ
where, the coefficients Qij are the elements of the transformed material stiffness matrix, ½Q, which are given by [15],
38
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
Q11 ¼ Q11 cos4 h þ 2ðQ12 þ 2Q66 Þsin2 hcos2 h þ Q22 sin4 h
gle function by introducing a scheme of reduction of unknowns. In order to attain that, the displacement components are expressed in terms of a potential function w of space variables as follows:
Q12 ¼ Q12 sin4 h þ cos4 h þ ðQ11 þ Q22 4Q66 Þsin2 hcos2 h 4
2
2
4
Q22 ¼ Q11 sin h þ 2ðQ12 þ 2Q66 Þsin hcos h þ Q22 cos h Q16 ¼ ðQ11 Q12 2Q66 Þ sin hcos3 h þ ðQ12 Q22 þ 2Q66 Þsin3 h cos h
ux ¼ a1
o2 w o2 w o2 w þ a2 þ a3 2 2 ox oxoy oy
ð4aÞ
uy ¼ a4
o2 w o2 w o2 w þ a þ a 5 6 ox2 oxoy oy 2
ð4bÞ
where, aÕs are unknown material constants. Substituting the above expressions of ux and uy in Eq. (1), we obtain the equilibrium equations in terms of the function w, which are
Q26 ¼ ðQ11 Q12 2Q66 Þsin3 h cos h þ ðQ12 Q22 þ 2Q66 Þ sin hcos3 h Q66 ¼ ðQ11 þ Q22 2Q12 2Q66 Þsin2 hcos2 h þ Q66 sin4 h þ cos4 h
E11 E211 ¼ ; 1 l12 l21 E11 l212 E22 E22 E11 E22 ¼ Q22 ¼ 1 l12 l21 E11 l212 E22 Q11 ¼
Q66 ¼ G12 ;
o4 w o4 w þ 2a3 Q16 þ a2 Q66 þ a6 Q12 þ a6 Q66 þ a5 Q26 2 2 ox oy oxoy 3 4 o w þ a3 Q66 þ a6 Q26 ¼0 ð5aÞ oy 4
l12 E22 l21 E11 l E11 E22 ¼ ¼ 12 1 l12 l21 1 l12 l21 E11 l212 E22
Q16 ¼ Q26 ¼ 0
The stress–strain relations for the general anisotropic materials are expressed through the transformed material stiffness matrix as follows [15]: 8 9 2 > Q11 < rxx > = 6 ryy ¼ 4 Q12 > > : ; rxy Q16
Q12 Q22 Q26
38 9 Q16 > < exx > = 7 Q26 5 eyy > : > ; cxy Q66
ð2Þ
In the above displacement formulation, the stresses are expressed in terms of the displacement components for their determination. The three stress–displacement relations for the plane problems of anisotropic materials, obtained from Eq. (2), are as follows: oux ouy oux ouy rxx ¼ Q11 þ Q12 þ Q16 þ ð3aÞ ox oy oy ox oux ouy oux ouy þ Q22 þ Q26 þ ryy ¼ Q12 ox oy oy ox
ð3bÞ
oux ouy oux ouy þ Q26 þ Q66 þ ox oy oy ox
ð3cÞ
rxy ¼ Q16
o4 w ox4
o4 w þ a2 Q11 þ 2a1 Q16 þ a5 Q16 þ a4 Q12 þ a4 Q66 ox3 oy þ a3 Q11 þ 2a2 Q16 þ a1 Q66 þ a6 Q16 þ a5 Q12 þ a5 Q66 þ a4 Q26
and
Q12 ¼ Q21 ¼
a1 Q11 þ a4 Q16
We now have to solve Eq. (1) with specified values of the functions rxx, ryy, rxy, ux and uy on segments of boundary with all possible combination of them. In the present approach, the two-dimensional problem of elasticity is reduced to the determination of a sin-
a1 Q16 þ a4 Q66
o4 w ox4
o4 w þ a2 Q16 þ a1 Q12 þ a1 Q66 þ a5 Q66 þ 2a4 Q26 ox3 oy þ a3 Q16 þ a2 Q12 þ a2 Q66 þ a1 Q26 þ a6 Q66 þ 2a5 Q26 þ a4 Q22 o4 w o4 w þ a3 Q12 þ a3 Q66 þ a2 Q26 þ 2a6 Q26 þ a5 Q22 2 2 ox oy oxoy 3 o4 w þ a3 Q26 þ a6 Q22 ¼0 ð5bÞ oy 4
Let us now choose aÕs in such a way that Eq. (5a) is automatically satisfied under all the circumstances. This will happen when coefficients of all the derivatives present in Eq. (5a) are individually zero. That is, when a1 Q11 þ a4 Q16 ¼ 0
ð6aÞ
a2 Q11 þ 2a1 Q16 þ a5 Q16 þ a4 Q12 þ a4 Q66 ¼ 0
ð6bÞ
a3 Q11 þ 2a2 Q16 þ a1 Q66 þ a6 Q16 þ a5 Q12 þ a5 Q66 þ a4 Q26 ¼ 0
ð6cÞ
2a3 Q16 þ a2 Q66 þ a6 Q12 þ a6 Q66 þ a5 Q26 ¼ 0
ð6dÞ
a3 Q66 þ a6 Q26 ¼ 0
ð6eÞ
Thus, for w to be a solution of the stress problem, it has to satisfy Eq. (5b) only. However, the aÕs of Eq. (4) must be known in advance. Here, we have five equations (Eq. (6)) for obtaining six unknown aÕs. We can thus assign an arbitrary value to one of these six unknowns and
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
solve for the remaining unknowns from Eq. (6). Let us assume a2 = 1 and, writing Eq. (6) in matrix form, we get 2
Q11
6 Q 6 66 6 6 2Q16 6 6 4 0
0
Q16
Q11
Q26
0
Q12 þ Q66
2Q16
0
Q66 9 8 0 > > > > > > > > > > > = < 2Q16 > ¼ Q11 > > > > > > > Q66 > > > > > ; : 0
0 Q12 þ Q66
Q16
0
Q26
0
0
y 2
38 9 a1 > > > > > > > 7> > Q16 a > 7> 3> 7< = 7 a4 0 7 > > > 7> a5 > > > Q12 þ Q66 5> > > > ; : > a 6 Q26
1
0
θ
x
ð7Þ
It is noted that, any other values of a2 could have been chosen leaving aside the homogeneous solutions. Different values of a2 will make the coefficients of the unknowns w different. The values of a2 which will lead to the diagonal dominance of the algebraic equations of unknown w at the nodal points are desirable in terms of increasing accuracy and shorter computational time. a2 = 1 leads to most simplified differential equations and hence more well-conditioned algebraic equations. The problem has now been reduced to the solution of a single function w from a single equation (Eq. (5b)). Let the equation be written as b1
39
o4 w o4 w o4 w o4 w o4 w þ b2 3 þ b3 2 2 þ b4 þ b5 4 ¼ 0 4 3 ox ox oy ox oy oxoy oy ð8Þ
then the coefficients bÕs are given by b1 ¼ a1 Q16 þ a4 Q66 b2 ¼ a1 Q12 þ Q66 þ a2 Q16 þ 2a4 Q26 þ a5 Q66
x,y,z - orthogonal loading axes
Z, 3
Fig. 1. Coordinate system for stress analysis of two-dimensional composite structure.
relations are thus obtained here as a special case of the general anisotropic material by substituting the value of h = 0. The stress–strain relations for orthotropic materials are expressed through the stiffness matrix, as follows [15]: 8 9 2 38 9 0 > Q11 Q12 > < rxx > < exx > = = 6 7 ryy ¼ 4 Q12 Q22 0 5 eyy ð9Þ > > > : : > ; ; cxy rxy 0 0 Q66 The two equilibrium equations, for the problems of orthotropic materials in terms of the displacement components, ux and uy, obtained from Eq. (1) are, Q11
o2 ux o2 uy o2 ux þ ðQ12 þ Q66 Þ þ Q66 2 ¼ 0 2 ox oxoy oy
ð10aÞ
Q66
o2 uy o2 ux o2 uy þ Q22 2 ¼ 0 þ ðQ12 þ Q66 Þ 2 ox oxoy oy
ð10bÞ
b3 ¼ a1 Q26 þ a2 Q12 þ Q66 þ a3 Q16 þ a4 Q22 þ 2a5 Q26 þ a6 Q66
b4 ¼ a2 Q26 þ a3 Q12 þ Q66 þ a5 Q22 þ 2a6 Q26 b5 ¼ a3 Q26 þ a6 Q22 Therefore, the governing differential equation (8) can be obtained explicitly once the values of aÕs of Eq. (7) are known. 3. Special cases of the general governing equation 3.1. Case-1: Orthotropic materials In this section, the governing equation of the w-formulation is deduced for the orthotropic materials where the principal material axes are aligned with the natural body axes, that is, the value of h is 0, in Fig. 1. All the
1,2,3 - principal material axes
The corresponding two equilibrium equations in terms of the function w, as obtained directly from Eq. (5) (or combining Eqs. (4) and (10)), are as follows: o4 w o4 w þ ða2 Q11 þ a4 Q12 þ a4 Q66 Þ 3 4 ox ox oy o4 w þ ða1 Q66 þ a3 Q11 þ a5 Q12 þ a5 Q66 Þ 2 2 ox oy o4 w o4 w þ ða2 Q66 þ a6 Q12 þ a6 Q66 Þ þ a3 Q66 4 ¼ 0 3 oxoy oy ð11aÞ
a1 Q11
40
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
o4 w o4 w þ ða1 Q12 þ a1 Q66 þ a5 Q66 Þ 3 4 ox ox oy o4 w þ ða2 Q12 þ a2 Q66 þ a4 Q22 þ a6 Q66 Þ 2 2 ox oy o4 w o4 w þ ða3 Q12 þ a3 Q66 þ a5 Q22 Þ þ a Q ¼0 6 22 oxoy 3 oy 4 ð11bÞ
a4 Q66
Likewise the case of general anisotropic material, assuming a2 = 1, the solutions for the remaining aÕs are obtained from the conditions required for the Eq. (11a) to be automatically satisfied. The values of aÕs thus obtained, are as follows: a1 ¼ a3 ¼ a5 ¼ 0 a2 ¼ 1 a4 ¼
a6 ¼
E211 Dn
ð12Þ
G12 ðE11 Dn
3.2. Case-2: Isotropic materials In addition to metallic materials, many of the particulate composites and also the randomly oriented discontinuous fiber composites can reasonably be treated as the isotropic materials. The corresponding equations for the problems of isotropic materials can readily be obtained from the general formulation of anisotropic materials by substituting the following conditions: 8 h ¼ 0o > > > E 11 ¼ E22 ¼ E > > : G12 ¼ G ¼ E=f2ð1 þ lÞg When the above conditions are substituted in Eq. (2), the corresponding stress–strain relations for isotropic materials can be expressed through the stiffness matrix, as follows: 8 9 2 E=1 l2 > < rxx > = 6 ryy ¼ 4 lE=1 l2 > > : ; rxy 0
l212 E22 Þ
where, Dn ¼ l12 E11 E22 þ G12 ðE11 l212 E22 Þ. Note that the solutions for aÕs can also be obtained from the matrix Eq. (7) when the condition of orthotropic material is inserted. Knowing the values of aÕs, the coefficients of different derivatives of Eq. (11b), that is, bÕs, are obtained as follows: E2 G12 b1 ¼ a4 Q66 ¼ 11 Dn
lE=1 l2 E=1 l2 0
38 9 0 > < exx > = 7 e 0 5 yy > : > ; cxy E=2ð1 þ lÞ
ð15Þ
Similarly, the two equilibrium equations in terms of the displacement components, obtained from Eq. (1) are o2 ux 1 þ l o2 uy 1 l o2 ux þ ¼0 þ ox2 oxoy oy 2 2 2
ð16aÞ
o2 uy 1 þ l o2 ux 1 l o2 uy þ þ ¼0 2 2 oy 2 oxoy ox2
ð16bÞ
b2 ¼ b4 ¼ 0
The corresponding equilibrium equations in terms of the function w, obtained from Eqs. (5) and (14), are given by
b3 ¼ a2 ðQ12 þ Q66 Þ þ a4 Q22 þ a6 Q66
a1
¼
E11 E22 ð2l12 G12 E11 Þ Dn
b5 ¼ a6 Q22 ¼
E11 E22 G12 Dn
Substituting the values of bÕs in Eq. (8) the explicit expression of the governing differential equation for the case of orthotropic materials becomes
E11 G12
o4 w o4 w þ E22 ðE11 2l12 G12 Þ 2 2 4 ox ox oy
þ E22 G12
o4 w ¼0 oy 4
ð13Þ
E o4 w 1 l2 ox4 4 E lE E ow þ a2 þ a4 þ 2 2 1l 1l 2ð1 þ lÞ ox3 oy 4 E E lE E ow þ a5 þ a3 þ a1 þ 2 2 1l 2ð1 þ lÞ 1l 2ð1 þ lÞ ox2 oy 2 4 E lE E ow þ a6 þ a2 þ 2ð1 þ lÞ 1 l2 2ð1 þ lÞ oxoy 3 þ a3
E o4 w ¼0 2ð1 þ lÞ oy 4
ð17aÞ
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51 a4
E o4 w 2ð1 þ lÞ ox4 4 lE E E ow þ a5 þ a1 þ 2 1l 2ð1 þ lÞ 2ð1 þ lÞ ox3 oy 4 lE E E E ow þ a4 þ a2 þ þ a6 2 2 1l 2ð1 þ lÞ 1l 2ð1 þ lÞ ox2 oy 2 4 lE E E ow þ a þ a3 þ 5 1 l2 2ð1 þ lÞ 1 l2 oxoy 3 E o4 w þ a6 ¼0 1 l2 oy 4
n y φ
A
(90-φ)
ut
un O C
uy ux t
Boundary segment
ð17bÞ
As before, assuming a2 = 1, the corresponding values of the aÕs for the case of isotropic materials can be obtained from Eq. (12) or directly from Eq. (17a), and are as follows:
41
x Fig. 2. Components of displacement on the boundary segment, AOC.
a1 ¼ a3 ¼ a5 ¼ 0 a2 ¼ 1
a6 ¼
2 1þl
–Y ð18Þ
y
1l 1þl
φ
σxx
When the above values of aÕs and the conditions of Eq. (14) are substituted in the general expressions of bÕs in Eq. (8), we get b1 ¼
b5 ¼
σxy
–X a
φ C
E
a sinφ
t B
σxy σyy
ð1 þ lÞ2
(90-φ)
x
b2 ¼ b4 ¼ 0 b3 ¼
n
A
a cosφ
a4 ¼
Fig. 3. Components of stress on the boundary segment, AB.
2E ð1 þ lÞ2
4. Boundary conditions
E ð1 þ lÞ2
Substitution of the above values of bÕs in Eq. (8) gives the governing differential equation for the case of isotropic material as o4 w o4 w o4 w þ2 2 2þ 4 ¼0 4 ox ox oy oy
ð19Þ
Note that the above governing bi-harmonic equation for isotropic materials, obtained as a special case of the general formulation for anisotropic materials, is found to be the identical to that used in our displacement potential function formulation for mixed boundary-value problems of elasticity. Reliability as well as suitability of that formulation has already been verified in our previous researches [7–10,16,17].
The physical conditions specified along the edge or boundary of a body are either known restraints or forces or some kind of combination of the two. Both the restraints, that is, known displacements and the forces are defined by their respective components. They are (Figs. 2 and 3), (a) (b) (c) (d)
Normal displacement (un) Tangential displacement (ut) Normal stress (rn) Tangential stress (rt)
The solution governing equation requires two known boundary conditions at any point on the boundary. Therefore, the possible boundary conditions created out of the four components are (1) (rn, rt), (2) (rn, un), (3) (rn, ut), (4) (rt, un), (5) (rt, ut), and (6) (un, ut).
42
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
However, the combinations (rn, un) and (rt, ut) do not occur in practice. The remaining four possible boundary conditions are therefore (1) (rn, rt), (2) (rn, ut), (3) (rt, un), and (4) (un, ut).
Although the known conditions are usually specified in terms of the normal and tangential components of stress or displacement of a point on the boundary, the parameters within the body are known as x- and y-components of stress and displacement. Since our objective is to solve the problem in terms of a single function w, all the boundary conditions of interest are required to be expressed in terms of the function. The components of displacement in terms of w are ux ðx; y Þ ¼ a1
o2 w o2 w o2 w þ a þ a 2 3 ox2 oxoy oy 2
ð20aÞ
uy ðx; y Þ ¼ a4
o2 w o2 w o2 w þ a5 þ a6 2 2 ox oxoy oy
ð20bÞ
Expressions for the components of stress for the interior points of an anisotropic body are obtained when the displacement components in the stress–displacement relations (Eq. (3)) are replaced by Eq. (4). Thus, we get 3
o w rxx ðx; yÞ ¼ a1 Q11 þ a4 Q16 ox3 o3 w þ a1 Q16 þ a2 Q11 þ a4 Q12 þ a5 Q16 ox2 oy o3 w þ a2 Q16 þ a3 Q11 þ a5 Q11 þ a6 Q16 oxoy 2 o3 w þ a3 Q16 þ a6 Q12 oy 3 ð21aÞ 3
o w ryy ðx; yÞ ¼ a1 Q12 þ a4 Q26 ox3 o3 w þ a1 Q26 þ a2 Q12 þ a4 Q22 þ a5 Q26 ox2 oy o3 w þ a2 Q26 þ a3 Q12 þ a5 Q22 þ a6 Q26 oxoy 2 3 o w þ a3 Q26 þ a6 Q22 oy 3 ð21bÞ
o3 w rxy ðx; yÞ ¼ a1 Q16 þ a4 Q66 ox3 o3 w þ a1 Q66 þ a2 Q16 þ a4 Q26 þ a5 Q66 ox2 oy o3 w þ a2 Q66 þ a3 Q16 þ a5 Q26 þ a6 Q66 oxoy 2 3 o w þ a3 Q66 þ a6 Q26 oy 3 ð21cÞ 4.1. Normal and tangential components of displacement In Fig. 2, the positive normal direction on the boundary is outward, positive tangential direction is clockwise and the positive /, the angle drawn from positive normal to positive x-axis, is clockwise. With these conventions, the relations between interior and boundary components of displacement can be written as ux ¼ un cos / þ ut sin / uy ¼ un sin / ut cos /
ð22Þ
Solving Eqs. (22) for un and ut, we have, un ðx; yÞ ¼ ux cos / þ uy sin / ut ðx; yÞ ¼ ux sin / uy cos /
ð23Þ
Combining Eqs. (20) and (23), the expressions for the normal and tangential components of displacement in terms of w, on the boundary of an anisotropic material, are given by un ðx; yÞ ¼ R1
o2 w o2 w o2 w þ R3 2 þ R2 2 ox oxoy oy
ð24aÞ
ut ðx; yÞ ¼ R4
o2 w o2 w o2 w þ R þ R 5 6 ox2 oxoy oy 2
ð24bÞ
where R1 ¼ a1 cos / þ a4 sin / R2 ¼ a2 cos / þ a5 sin / R3 ¼ a3 cos / þ a6 sin / R4 ¼ a1 sin / a4 cos / R5 ¼ a2 sin / a5 cos / R6 ¼ a3 sin / a6 cos /
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
4.2. Normal and tangential components of stress
ð25Þ
Y ¼ ryy sin / þ rxy cos /
Again, the normal and tangential components of stress in terms of the components, X and Y , can be written as rn ¼ X cos / þ Y sin /
Q11 Q12 sin 2/ Q16 cos 2/ 2 Q Q22 T 5 ¼ 12 sin 2/ Q26 cos 2/ 2 Q Q26 T 6 ¼ 16 sin 2/ Q66 cos 2/ 2 T4 ¼
With reference to Fig. 3, if the condition of equilibrium is applied for the forces on a segment of boundary AB, the components of stress along the direction of coordinate axes are X ¼ rxx cos / þ rxy sin /
43
As far as numerical computation is concerned, it is thus evident from the above expressions that all the boundary conditions of interest can easily be discretized in terms of the single function w in the method of finite-difference.
ð26Þ
rt ¼ X sin / Y cos /
5. Special cases of the boundary expressions
or 5.1. Case-1: Orthotropic materials
rn ¼ rxx cos2 / þ ryy sin2 / þ rxy sin 2/ r r xx yy sin 2/ rxy cos 2/ rt ¼ 2
ð27Þ
Now, if the expressions of body parameters, namely rxx, ryy, and rxy of Eq. (21), are substituted in Eq. (27), the normal and tangential components of stress in terms of w are 3
3
3
ux ðx; yÞ ¼
3
ow ow ow ow rn ðx; yÞ ¼ P 1 3 þ P 2 2 þ P 3 þ P4 3 ox ox oy oxoy 2 oy
Explicit expressions for the body parameters in terms of the function w are obtained from the general equations (20) and (21) when the respective conditions of orthotropic materials are substituted in them. The conditions are,
ð28aÞ
o2 w oxoy
uy ðx; yÞ ¼
where,
o2 w 1 o2 w E211 2 þ G12 E11 l212 E22 Dn ox oy 2
P 1 ¼ a1 T 1 þ a4 T 3
ð29bÞ E11 G12 o3 w o3 w rxx ðx; yÞ ¼ E11 2 l12 E22 3 ox oy oy Dn
P 2 ¼ a1 T 3 þ a2 T 1 þ a4 T 2 þ a5 T 3 P 3 ¼ a2 T 3 þ a3 T 1 þ a5 T 2 þ a6 T 3 P 4 ¼ a3 T 3 þ a6 T 2
ryy ðx; yÞ ¼
ð29dÞ
T 2 ¼ Q12 cos2 / þ Q22 sin2 / þ Q26 sin 2/
rxy ðx; yÞ ¼
2
T 3 ¼ Q16 cos / þ Q26 sin / þ Q66 sin 2/ and rt ðx; yÞ ¼ P 5
o3 w o3 w o3 w o3 w þ P þ P þ P 6 7 8 ox3 ox2 oy oxoy 2 oy 3
where P 5 ¼ a1 T 4 þ a4 T 6 P 6 ¼ a1 T 6 þ a2 T 4 þ a4 T 5 þ a5 T 6 P 7 ¼ a2 T 6 þ a3 T 4 þ a5 T 5 þ a6 T 6 P 8 ¼ a3 T 6 þ a6 T 5
ð29cÞ
E11 E22 o3 w o3 w ðl12 G12 E11 Þ 2 G12 3 Dn ox oy oy
T 1 ¼ Q11 cos2 / þ Q12 sin2 / þ Q16 sin 2/
2
ð29aÞ
ð28bÞ
E11 G12 o3 w o3 w E11 3 l12 E22 ox oxoy 2 Dn
ð29eÞ
Similarly, evaluating the coefficients RÕs for the orthotropic materials and substituting in Eq. (24), the expressions for the normal and tangential components of displacement are E211 sin / o2 w o2 w þ cos / 2 Dn ox oxoy 2 G12 E11 l12 E22 sin / o2 w Dn oy 2
un ðx; yÞ ¼
ð30aÞ
44
ut ðx; yÞ ¼
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
E211 cos / o2 w o2 w þ sin / 2 Dn ox oxoy G12 E11 l212 E22 cos / o2 w þ Dn oy 2
rxy ðx; yÞ ¼ ð30bÞ
When the respective values of the coefficients PÕs and TÕs are substituted in the general equation (28), the expressions for the normal and tangential components of stress for points on the boundary of an orthotropic body are,
ð32eÞ
Evaluating the corresponding values of the coefficients RÕs for the case of isotropic material and substituting them in Eq. (24), the expressions for the normal and tangential components of displacement are 2 o2 w o2 w sin / 2 þ cos / 1þl ox oxoy
un ðx; yÞ ¼
E211 G12 sin 2/ o3 w Dn ox3 E11 E11 G12 cos2 / E11 E22 sin2 / þ l12 E22 G12 sin2 / þ Dn
rn ðx; yÞ ¼
1l o2 w sin / 2 1þl oy
ut ðx; yÞ ¼
o3 w l E11 E22 G12 sin 2/ o3 w 2 þ 12 ox oy Dn oxoy 2 E11 E22 G12 l12 cos2 / þ sin2 / o3 w oy 3 Dn
3 ow o3 w þ l ox3 oxoy 2 ð1 þ l Þ2 E
ð33aÞ
2 o2 w o2 w cos / 2 þ sin / 1þl ox oxoy þ
1l o2 w cos / 2 1þl oy
ð33bÞ
ð31aÞ E2 G12 cos 2/ o3 w rt ðx; yÞ ¼ 11 Dn ox3 E11 ðE11 G12 þ E11 E22 l12 E22 G12 Þ sin 2/ þ 2Dn
o3 w l E11 E22 G12 cos 2/ o3 w 12 2 ox oy Dn oxoy 2
Similarly, evaluating the coefficients PÕs and TÕs for the case of isotropic material, the expressions for the normal and tangential components of stress for points on the boundary of an isotropic body, are obtained from Eq. (28) as follows: rn ðx; yÞ ¼
E11 E22 G12 ð1 l12 Þ sin 2/ o3 w þ 2Dn oy 3
þ ð31bÞ þ
5.2. Case-2: Isotropic materials Here, the boundary conditions for the case of isotropic materials are derived as a special case of those for the general anisotropic material. Substituting the conditions of Eq. (14) along with the respective aÕs of Eq. (18), expressions for the body parameters in terms of the function w can be obtained from the general Eqs. (20) and (21) as follows: ux ðx; yÞ ¼
o2 w oxoy
uy ðx; yÞ ¼
rxx ðx; yÞ ¼
ð32aÞ
2 1 ow o2 w 2 2 þ ð1 lÞ 2 1 þ l ox oy
E ð1 þ lÞ2
ryy ðx; yÞ ¼
o3 w o3 w l 3 2 ox oy oy
o3 w o3 w þ ð 2 þ l Þ ox2 oy oy 3 ð1 þ lÞ2 E
ð32bÞ
ð32cÞ
ð32dÞ
E ð1 þ lÞ
sin 2/
2
E ð1 þ lÞ
2
El
E
cos2 / 2sin2 / lsin2 /
sin 2/
ð1 þ lÞ2 ð1 þ lÞ2
o3 w ox3 o3 w ox2 oy
o3 w oxoy 2
lcos2 / þ sin2 /
o3 w oy 3 ð34aÞ
rt ðx; yÞ ¼
E 2
cos 2/
o3 w Eð3 þ lÞ þ sin 2/ ox3 2ð1 þ lÞ2
ð1 þ l Þ o3 w El o3 w cos 2/ 2 ox oy ð1 þ lÞ2 oxoy 2 þ
Eð1 lÞ 2ð1 þ lÞ
2
sin 2/
o3 w oy 3
ð34bÞ
One can thus easily realize that all the boundary conditions for the problems of orthotropic as well as isotropic materials can be deduced from those of the general expressions of anisotropic material. It is noted that the differential equations of the boundary conditions for the case of isotropic material are found to be identical with those used in our w-formulation [10,16,17].
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51 Stencil for both orthotropic and isotropic material
(a)
45
(b)
i-forward & j-forward j
k h
k i j j
y
i
y i
h
x
Stencil for anisotropic material
x i-forward & j-forward
(c)
h k j
y i
x Fig. 4. (a) Finite-difference grid structure for the governing equilibrium equations; (b) Finite-difference stencil for the normal and tangential components of displacements on the boundary (Eq. (37)); (c) Finite-difference stencil for the normal and tangential components of stress on the boundary (Eq. (38)).
6. Finite-difference discretization of the governing equations The present computational scheme involves evaluation of a single function w at the mesh points of a uniform rectangular mesh network over a geometrically regular or irregular region. The governing differential equation (8), which is used to evaluate the function w at the internal mesh points of an anisotropic body, can be expressed in its corresponding difference form using central difference operators. The complete finite-difference expression of the governing equation is given by,
46
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
24b1 n4 þ 16b3 n2 þ 24b5 wði; jÞ 2
2
8n ð2n b1 þ b3 Þfwði þ 1; jÞ þ wði 1; jÞg þ b2 n3 fwði þ 2; j þ 1Þ wði þ 2; j 1Þ þwði 2; j 1Þ wði 2; j þ 1Þg þ b4 nfwði þ 1; j þ 2Þ wði 1; j þ 2Þþwði 1; j 2Þ wði þ 1; j 2Þg þ 4b1 n4 fwði þ 2; jÞ þ wði 2; jÞg þ 4b5 fwði; j þ 2Þ þwði; j 2Þg þ 2b2 n3 þ 4b3 n2 2b4 n fwði þ 1; j þ 1Þ wði 1; j 1Þg þ 2b2 n3 þ 4b3 n2 þ 2b4 n
expressions (for example, i- & j-forward, i- & j-backward, i-forward & j-backward, and i-backward & j-forward), for each of the boundary conditions, are developed for different regions of boundary. For example, the i- and j-forward difference formulae for the normal and tangential components of displacement on the boundary of an anisotropic body, corresponding to Eq. (24), are written as follows: R2 fwði þ 2; j þ 2Þ 4wði þ 2; j þ 1Þ 4hk þ 16wði þ 1; j þ 1Þ 4wði þ 1; j þ 2Þ R1 þ3wði; j þ 2Þ þ 3wði þ 2; jÞg þ 2 wði 1; jÞ h R1 3R2 R3 wði þ 1; jÞ þ 2 wði; j 1Þ þ 2 hk k h R3 3R2 wði; j þ 1Þ þ 2 hk k 9R2 2R1 2R3 2 2 wði; jÞ ð37aÞ þ 4hk h k
un ðx; yÞ ¼
fwði þ 1; j 1Þ þ wði 1; j þ 1Þg 8ðn2 b3 þ 2b5 Þ fwði; j þ 1Þ þ wði; j 1Þg ¼ 0 ð35Þ
where, n = k/h (Fig. 4). The corresponding grid structure of the governing equation (35) for an arbitrary internal mesh point (i, j) is shown in Fig. 4(a). It is thus clear that, application of the governing equation involves a total of 21 neighboring mesh points, and when this point of application (i, j) becomes an immediate neighbor to the physical boundary, the equation will involve mesh points not only interior but also exterior to the physical boundary. After substitution of the respective values of bÕs for the orthotropic or isotropic materials into Eq. (35), the corresponding central difference equations can easily be obtained. The finite-difference approximation of the governing equation thus obtained for isotropic materials is, 2n2 fwði þ 1; j þ 1Þ þ wði þ 1; j 1Þ þ wði 1; j þ 1Þ þwði 1; j 1Þg þ wði; j þ 2Þ 4n2 n2 þ 1 fwði þ 1; jÞ þ wði 1; jÞg 4 n2 þ 1 fwði; j þ 1Þ þ wði; j 1Þg þ n4 fwði þ 2; jÞ þwði 2; jÞg þ wði; j 2Þ þ 6n4 þ 8n2 þ 6 wði; jÞ ¼ 0 ð36Þ
The corresponding grid structure of Eq. (36), as shown in Fig. 4(a), involves a total of 13 neighboring mesh points. From the figures, it is clear that the use of the above central difference stencils of the governing equations, especially for the points in the immediate neighborhood of the physical boundary, demands the inclusion of an imaginary (false) boundary, exterior to the physical boundary of the domain concerned. As the differential equations associated with the boundary conditions contain second- and third-order partial derivatives of the function w, the use of centraldifference expressions eventually leads to the inclusion of points exterior to the false boundary. In order to avoid the occurrence, different versions of the finite-difference expressions like forward, backward and central have been adopted in a combined form to generate the difference equations for the boundary conditions. It is noted here that the order of local truncation error has been kept the same for all the expressions developed, that is, O(h2). Four different sets of finite-difference
R5 fwði þ 2; j þ 2Þ þ 16wði þ 1; j þ 1Þ 4hk 4wði þ 2; j þ 1Þ 4wði þ 1; j þ 2Þ R4 þ3wði; j þ 2Þ þ 3wði þ 2; jÞg þ 2 wði 1; jÞ h R4 3R5 R6 wði þ 1; jÞ þ 2 wði; j 1Þ þ 2 hk k h R6 3R5 wði; j þ 1Þ þ 2 hk k 9R5 2R4 2R6 2 2 wði; jÞ ð37bÞ þ 4hk h k
ut ðx; yÞ ¼
Similarly, from Eq. (28), the corresponding i- and j-forward difference formulae for the normal and tangential components of stress on the boundary are as follows: 5P 1 3P 2 3P 3 5P 4 rn ðx; yÞ ¼ wði; jÞ þ þ þ h3 h2 k hk 2 k3 2P 2 2P 3 wði þ 1; j þ 1Þ þ þ h2 k hk 2 6P 1 3P 2 4P 3 þ 2 þ 2 wði þ 1; jÞ h3 2h k hk 4P 2 3P 3 6P 4 wði; j þ 1Þ þ þ h2 k 2hk 2 k3 3P 1 3P 2 þ wði 1; jÞ 2h3 2h2 k 3P 3 3P 4 wði; j 1Þ þ 2hk 2 2k 3 3P 1 P 3 þ 2 wði þ 2; jÞ þ h3 hk
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
P 2 3P 4 P1 wði; j þ 2Þ 3 wði þ 3; jÞ þ 2 3 hk k 2h P2 þ 2 f4wði 1; j þ 1Þ wði 1; j þ 2Þ 2h k P4 wði þ 1; j þ 2Þg 3 wði; j þ 3Þ 2k P3 f4wði þ 1; j 1Þ wði þ 2; j þ 1Þ þ 2hk 2 wði þ 2; j 1Þg
½Am;m fX gm ¼ fBgm
þ
47
ð39Þ
for m unknowns; and
½Amþ1;mþ1 fX gmþ1 ¼ fBgmþ1
for m þ 1 unknowns
ð40Þ
To reduce Eq. (40) to (39), we have to perform (i) (m + 1)(m + 1) = (m + 1)2 divisions and (ii) (m + 1)m subtractions ð38aÞ
5P 5 3P 6 3P 7 5P 8 þ 2 þ 2 þ 3 wði; jÞ h3 h k hk k 2P 6 2P 7 þ wði þ 1; j þ 1Þ þ h2 k hk 2 6P 5 3P 6 4P 7 þ 2 þ 2 wði þ 1; jÞ h3 2h k hk 4P 6 3P 7 6P 8 wði; j þ 1Þ þ þ h2 k 2hk 2 k3 3P 5 3P 6 þ wði 1; jÞ 2h3 2h2 k 3P 7 3P 8 wði; j 1Þ þ 2hk 2 2k 3 3P 5 P 7 þ þ 2 wði þ 2; jÞ h3 hk P 6 3P 8 P5 þ 2 þ 3 wði; j þ 2Þ 3 wði þ 3; jÞ hk k 2h P6 þ 2 f4wði 1; j þ 1Þ wði 1; j þ 2Þ 2h k P8 wði þ 1; j þ 2Þg 3 wði; j þ 3Þ 2k P7 f4wði þ 1; j 1Þ wði þ 2; j þ 1Þ þ 2hk 2 wði þ 2; j 1Þg ð38bÞ
rt ðx; yÞ ¼
The corresponding finite-difference stencils for the normal and tangential components of displacement and stress for points on the boundary are illustrated in Fig. 4(b) and (c). It is interesting to note that the boundary stencils for the displacement and stress components take identical form for the general anisotropic, orthotropic and isotropic materials.
When m = 0, the increase in computation due to addition of one unknown is one division. That is, when we have only one unknown, computation involved is one division. When m = 1, the increase is 4 divisions and 2 subtractions. Let us consider the total number of nodes = N; Number of parameters at each node (like ux, uy) = P. Therefore, total number of unknowns, m, in Eq. (39) is given by m = NP. For a two-dimensional problem, in both existing finite-difference and finite-element techniques, P = 2. Here, in the present finite-difference scheme, we have to find only one parameter, w, i.e., P = 1. So the reduction of the unknowns between the two techniques are from 2N of finite-element to N of the present finite-difference scheme. The increase in computation due to the additional N unknowns is
(a)
Uniformly distributed load, w
D
x
0.1L
0.8L
0.1L
y
(b)
Physical boundary
Imaginary boundary
7. Scheme for reduction of computational work x, i
This section explains how the computational work is reduced due to reduction in the number of parameters to be evaluated at the nodal points. Both in finite-element and in finite-difference methods, 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
k h
y, j
Fig. 5. (a) Loading and geometry of a simply supported deep beam; (b) Finite-difference discretization of the beam domain.
48
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
(i) (N + 1)2 + (N + 2)2 + (N + 3)2 + + (2N)2 [division], and (ii) N(N + 1) + (N + 1)(N + 2) + (N + 2)(N + 3) + + (2N1)(2N) [subtraction] Let us now evaluate the percentage reduction in the computational efforts for the individual cases of division and subtraction. (i) Division:
ffi 87:5% when N is large
(ii) Subtraction: Saving in the number of subtraction 2N ð2N 1Þð4N 1Þ 2Nð2N 1Þ þ ¼ 6 2 N ðN 1Þð2N 1Þ N ðN 1Þ þ 6 2
Saving in the number of division ¼
Percentage saving in division alone ð2N Þð2N þ 1Þð4N þ 1Þ N ðN þ 1Þð2N þ 1Þ 100 ¼ 2N ð2N þ 1Þð4N þ 1Þ N þ1 100% ¼ 1 2ð4N þ 1Þ
ð2NÞð2N þ 1Þð4N þ 1Þ N ðN þ 1Þð2N þ 1Þ 6 6
1.5
1.0 σxx / (w/B)
σxx / (w/B)
1.0
1.5 FDM FEM Simple theory
0.5 0.0 -0.5 -1.0 -0.6
-0.4
-0.2
0.0 y/D
0.2
0.4
0.0
-1.0 -0.6
0.6
-0.4
-0.2
(b)
0.0 y/D
0.2
0.4
0.6
0.2
0.4
0.6
3
2 FDM FEM Simple theory
2
σxx / (w/B)
σxx / (w/B)
0.5
-0.5
(a)
1
FDM FEM Simple theory
0
FDM FEM Simple theory
1 0 -1
-1
-2 -2 -0.6
(c)
-0.4
-0.2
0.0 y/D
0.2
0.4
-3 -0.6
0.6
(d)
-0.4
-0.2
0.0 y/D
Fig. 6. Normalized bending stress distributions at the midspan section for the beam: (a) L/D = 0.5; (b) L/D = 1.0; (c) L/D = 1.5; (d) L/D = 2.0.
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
Percentage saving in subtraction alone " # 2 2N ð4N 2 1Þ N ðN3 1Þ 3 100 ¼ 2N ð4N 2 1Þ 3
ðN 2 1Þ 100% ¼ 1 2ð4N 2 1Þ ffi 87:5% when N is large It is thus verified that the present finite-difference method of solution can save approximately 87% of computational work for the solution of a two-dimensional problem, in comparison with that required for usual finite-element solution. 8. Example of application of the w-formulation Using the present w-formulation a classical problem of solid mechanics has been solved numerically and the results are compared with the available solutions in the literature. The present isotropic problem has been
49
solved here as a special case of the general anisotropic material. The problem is a deep beam of rectangular cross section, simply supported at the ends and carries a uniformly distributed load over its whole span. The geometry and loading are illustrated in Fig. 5(a). It is noted that the analysis of the present problem is being constantly looked into by several researchers in the field of solid mechanics [3,4,11,18,19]. However, each solution possesses certain limitations, and eventually none of the solutions are found to conform to all the physical characteristics of the problem appropriately. A relatively small number of mesh points (23 · 23) are considered here for the finite-difference solution of the problem, so that the solutions can be readily compared with the FEM results at identical sections of the beams. The finite-difference mesh-network used for the present problem is shown in Fig. 5(b), where an imaginary (false) boundary, exterior to the physical boundary, is included. It is noted that the error of the present computational approach is of the order of h2. It has been
0.6
0.5 FDM FEM Simple theory
0.4
0.4 σxy / (w/B)
σxy / (w/B)
0.3 0.2 0.1
0.2
0.0
0.0
FDM FEM Simple theory
-0.1 -0.6
-0.4
-0.2
(a)
0.0 y/D
0.2
0.4
-0.2 -0.6
0.6
-0.4
-0.2
(b)
0.0 y/D
0.2
0.4
0.6
0.4
0.6
1.0
0.6
0.4
σxy / (w/B)
σxy / (w/B)
0.8
0.2 FDM FEM Simple theory
0.0
0.6 0.4 0.2
FDM FEM Simple theory
0.0
-0.2 -0.6
(c)
-0.4
-0.2
0.0 y/D
0.2
0.4
0.6
-0.6
(d)
-0.4
-0.2
0.0 y/D
0.2
Fig. 7. Normalized shearing stress distributions at section: (a) x/L = 0.25 for the beam, L/D = 0.5; (b) x/L = 0.25 for the beam, L/D = 1.0; (c) x/L = 0.25 for the beam, L/D = 1.5; (d) x/L = 0.25 for the beam, L/D = 2.0.
50
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
Table 1 Comparison of maximum bending stress at section, x/L = 0.5 Beam size
Normalized maximum bending stress Simple theory
Tension side, at y/D = 0.5 L/D = 0.5 0.15 L/D = 1.0 0.60 L/D = 1.5 1.35 L/D = 2.0 2.40 Compression side, at y/D = 0.5 L/D = 0.5 0.15 L/D = 1.0 0.60 L/D = 1.5 1.35 L/D = 2.0 2.40
FEM
FDM
1.27 1.30 1.69 2.44
1.33 1.35 1.71 2.46
0.32 0.50 1.22 2.25
0.27 0.47 1.24 2.24
Table 2 (a) Comparison of maximum shearing stress at section, x/ L = 0.25; (b) Comparison of shearing stress at the section, y/ D = 0.5, x/L = 0.25 Beam size
(a) L/D = 0.5 L/D = 1.0 L/D = 1.5 L/D = 2.0
Normalized maximum shearing stress Simple theory
FEM
FDM
0.19 0.38 0.56 0.75
0.41 0.43 0.49 0.78
0.40 0.43 0.57 0.76
Normalized shearing stress (b) L/D = 0.5 L/D = 1.0 L/D = 1.5 L/D = 2.0
0.00 0.00 0.00 0.00
0.06 0.18 0.27 0.34
0.00 0.00 0.00 0.00
0.2 0.0
σyy / (w/B)
-0.2
L / D = 0.5 L / D = 1.0 L / D = 1.5 L / D = 2.0
-0.4 -0.6 -0.8 -1.0 -1.2 -0.6
-0.4
-0.2
0.0 y/D
0.2
0.4
0.6
Fig. 8. Finite-difference solutions for the normal stress component, ryy at section, x/L = 0.25, of the beams.
however confirmed by our FDM results that no significant change in accuracy is observed with decreasing h, when it is sufficiently small, as used in the present solution. The governing differential equation is used to evaluate the function w only at the internal mesh points. The discrete values of the displacement potential function w(x, y), at the mesh points of the domain concerned are obtained from a system of linear algebraic equations, resulting from the discretization of the governing equation and the associated boundary conditions, by the use of direct method of solution (CholeskeyÕs triangular-decomposition method). Finally, all the parameters of interest in the solution, namely, displacement and stress components, are obtained as soon as w is known at the mesh points of the domain, as all the components are expressed as summation of different derivatives of the function w. In the case of present problem, stresses are normalized with respect to the uniform loading w, and beam width B, consistent with that of Ref. [19]. In order to obtain the numerical results, values for YoungÕs modulus and PoissonÕs ratio of 209 GPa and 0.3, respectively, are used. The present solutions are compared with the solutions of simple beam theory and also with the finite-element solutions [19]. The finite-element solutions were obtained by using a regular mesh of 200, 8-noded, plane stress, isoparametric elements (20 · 10). Fig. 6 presents the bending stress distributions at the mid-sections of the beams of aspect ratio, L/D = 0.5, 1.0, 1.5, and 2.0. Along with the present finite-difference solutions, the FEM results [19] and also the results of simple bending theory of beam are included for the convenience of comparison. The mid-section has been chosen here for consideration since it is the section of maximum bending moment and also it is almost free from the local influence of the loading. It is seen that, for all the beams, the bending stress curves agree reasonably well with the finite-element solutions but differ significantly, as expected, from those of the simple flexure theory. For higher L/D ratios (L/D > 1.0) the finite-element predictions are found to be identical with those of our FDM results. As expected, the nonlinearity in the distributions as well as the deviations from the simple theory is found to decrease as the L/D ratio is increased. The distribution of shearing stress at section x/ L = 0.25 is shown in Fig. 7 for different beams. The parabolic distributions based on simple theory and the corresponding FEM solutions are included in the figures. For small L/D ratios, the shear stress distributions are significantly different from that obtained from the simple theory. The maximum shear stress predicted by our FDM approach is higher for all the beams inspected than those predicted by the simple theory. However, the deviation is found to decrease as L/D ratio increases. The FEM solutions, in general, agree well with the pre-
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
sent FDM results. However, it is interesting to note that the finite-element prediction of maximum shear stress for the beam of L/D = 1.5 is found to take a lower value even than that predicted by simple theory, which is in contrast to the solutions obtained for beam of other L/ D ratios. More importantly, the FEM predictions of shear stress in the neighborhood of both top and bottom surfaces are highly unreliable, even though the section concerned is free from singularity. This discrepancy is found to remain quite significant even in the case of the beam of L/D = 2.0. Therefore, the present investigation verifies the fact that the present finite-difference method of solution is free from the inability of FEM in predicting the actual state of stresses at the surfaces. Tables 1 and 2 list some of the numerical data for the bending and shear stresses at different sections of the beams; one can easily compare the accuracy of the present solutions with those of FEM and simple beam theory. Finally, Fig. 8 presents the distributions of normal stress component ryy at section x/L = 0.25. It should be noted here that the finite-difference analysis reproduces the normal stress exactly both at the top as well as at bottom boundaries of all the beams considered. 9. Conclusions The present paper demonstrates the extension of the displacement potential formulation of isotropic materials to the case of general anisotropic materials. It has been shown that the governing equations for the isotropic as well as orthotropic materials can readily be obtained as a special case of the proposed general formulation. The present approach of stress analysis for the plane problems of solid mechanics by finite-difference technique saves us approximately 87% of computational effort in comparison to that taken by exiting FE and FD methods. In addition, the quality of solutions will be improved significantly because of the fact that the solution of any algebraic system deteriorates with increasing number of unknowns and eventually fails when the number of unknowns become very large. The present approach not only increases our horizon of managing bigger problems in stress analysis but also in other fields of science and engineering. Summarizing, we would like to point out that the use of the present approach is expected to be the most promising one to save computational effort as well as to increase the horizon of capabilities in terms of solving bigger problems with higher accuracy.
51
References [1] Timoshenko D, Goodier JN. Theory of elasticity. 3rd ed. New York, NY: McGraw-Hill Book Company; 1979. [2] Uddin MW. Finite difference solution of two-dimensional elastic problems with mixed boundary conditions. M.Sc. thesis, Carleton University, Canada, 1966. [3] Conway HD, Chow L, Morgan GW. Analysis of deep beams. J Appl Mech Trans ASME 1951;18(2):163–72. [4] Chow L, Conway HD, Winter G. Stresses in deep beams. Trans ASCE, 1952, Paper No. 2557. [5] Chapel RE, Smith HW. Finite-difference solutions for plane stresses. AIAA J 1968;6(6):1156–7. [6] Conway HD, Ithaca NY. Some problems of orthotropic plane stress. J Appl Mech Trans ASME 1953:72–6. Paper No. 52-A-4. [7] Ahmed SR, Khan MR, Islam KMS, Uddin MW. Analysis of stresses in deep beams using displacement potential function. J Inst Engrs (India) 1996;77:141–7. [8] Idris ABM, Ahmed SR, Uddin MW. Analytical solution of a 2-D elastic problem with mixed boundary conditions. J Inst Engrs (Singapore) 1996;36(6):11–7. [9] Ahmed SR, Idris ABM, Uddin MW. Numerical solution of both ends fixed deep beams. Comp Struct 1996;61(1): 21–9. [10] Ahmed SR, Khan MR, Islam KMS, Uddin MW. Investigation of stresses at the fixed end of deep cantilever beams. Comp Struct 1998;69:329–38. [11] Durelli AJ, Ranganayakamma B. Parametric solution of stresses in beams. J Eng Mech 1989;115(2):401–15. [12] 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. [13] Smart J. On the determination of boundary stresses in finite elements. J Strain Anal 1987;22(2):87–96. [14] Dow JO, Jones MS, Harwood SA. A new approach to boundary modeling for finite difference applications in solid mechanics. Int J Numer Meth Eng 1990;30:99–113. [15] Jones RM. Mechanics of composite materials. McGrawHill Book Company; 1975. [16] Akanda MAS, Ahmed SR, Khan MR, Uddin MW. A finite-difference scheme for mixed boundary-value problems of arbitrary-shaped elastic bodies. Adv Eng Softw 2000;31(3):173–84. [17] Akanda MAS, Ahmed SR, Uddin MW. Stress analysis of gear teeth using displacement potential function and finitedifferences. Int J Numer Meth Eng 2001;53(7):1629–49. [18] Suzuki S. Stress analysis of short beams. AIAA J 1986;24:1396–8. [19] Hardy SJ, Pipelzadeh MK. Static analysis of short beams. J Strain Anal 1991;26(1):15–29.
A general mathematical formulation for finite-difference solution of mixed-boundary-value problems of anisotropic materials S. Reaz Ahmed *, M. Zubaer Hossain, M. Wahhaj Uddin Department of Mechanical Engineering, Bangladesh University of Engineering and Technology, Dhaka 1000, Bangladesh Received 2 July 2003; accepted 9 August 2004
Abstract This paper presents a general mathematical model, especially suitable for finite-difference analysis of stresses and displacements of the plane elastic problems of solid mechanics. The present formulation covers the problems of anisotropic, orthotropic and isotropic materials, in which the problem is formulated in terms of a single potential function, defined in terms of the displacement components. In addition, the formulation contains a new scheme of reduction of unknowns to be solved for a particular problem. Compared to the conventional computational approaches, the present scheme gets solution of higher accuracy and in extremely shorter time. The application of the present scheme is demonstrated here through a classical problem of solid mechanics, and the results are compared with the available solutions in the literature. 2004 Elsevier Ltd. All rights reserved. Keywords: Stress analysis; Displacement potential function; Anisotropic materials; Finite-difference method
1. Introduction Analysis of stress in a material body is usually a three-dimensional problem. However, in most practical cases, the problem can reasonably be approximated to either a two-dimensional or a one-dimensional one. One-dimensional problems of stress analysis are very straightforward and included in undergraduate curricula. But the problems of two-dimensional stress analysis have a different story. In the case of very simple structures like beam, column, circular disc, regular shaped plate and thin shells, having simple geometry and boundary conditions, analytical methods are used in *
Corresponding author. E-mail address: [email protected] (S.R. Ahmed).
their stress analysis. However, the analytical approach is limited, because general closed form solution can be obtained only for very ideal cases. Moreover, the analytical results obtained for the practical problems of elasticity, usually with mixed boundary conditions, are invariably approximate as they include various idealizations. Among the existing mathematical models for twodimensional elastic problems, the stress function approach [1] and the displacement functions approach (ux, uy-formulation, [2]) are noticeable. Successful application of the stress function formulation in conjunction with the finite-difference technique has been reported for the solution of plane elastic problems where the conditions on the boundary are prescribed in terms of stresses only [3–5]. Further, Conway and Ithaca [6] extended the stress function formulation in the form of Fourier
0045-7949/$ - see front matter 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2004.08.007
36
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
Nomenclature x, y, z 1, 2, 3 E11 E22 l12 l21 G12 E, l, G
orthogonal loading axes principle material axes elastic modulus of the material in 1-direction elastic modulus of the material in 2-direction major PoissonÕs ratio minor PoissonÕs ratio shear modulus in the 1–2 plane elastic modulus, PoissonÕs ratio, modulus of rigidity of isotropic material h fiber orientation ½Q transformed material stiffness matrix Qij coefficients of the transformed stiffness matrix Qij coefficients of the stiffness matrix ux, uy, uz displacement components in the x-, y- and z-direction un, ut normal and tangential component of displacement on the boundary rxx, ryy normal stress components in the x- and ydirection rxy shearing stress component in the xy plane exx, eyy normal strain components in the x- and ydirection cxy shearing strain component in the xy plane rn, rt normal and tangential components of stress on the boundary
integrals to the case where the material is orthotropic, and obtained analytical solutions for a number of problems. Mixed-boundary-value problems are those in which the boundary conditions are specified as a mixture of boundary restraints and boundary loadings. The shortcoming of the stress function formulation is that it accepts boundary conditions in terms of loadings only. Boundary restraints specified in terms of ux and uy cannot be satisfactorily imposed on the stress function. As most of the practical problems of elasticity are of mixed boundary conditions, the approach fails to provide any explicit understanding of the stress distribution in the regions of restrained boundaries, which are, in general, the most critical zones in terms of stresses. The difficulties involved in trying to solve practical stress problems using the stress function approach are pointed out in our previous researches [7–10], and also by Durelli [11]. Two-dimensional stress analysis in practice is mainly handled by numerical methods. The major numerical methods in use are (a) the method of finite-difference and (b) the method of finite-element. Finite-difference method is an ideal numerical approach for solving par-
n
normal component of the force on the boundary segment t tangential component of the force on the boundary segment / angle between the unit positive normal to the boundary and positive x-axis X, Y components of stress on the boundary, along the direction of coordinate axes a length of boundary segment, AB i, j nodal coordinate system in x- and ydirections h, k mesh lengths in the x- and y-directions n k/h [A] coefficient matrix {B} constant matrix m order of the coefficient matrix P number of parameters at each node N total no of nodes w displacement potential function ai coefficients of the derivatives of w in the expressions of displacement components bi coefficients of the general governing differential equation of w Dn coefficient composed of elastic constants of orthotropic material w uniformly distributed loading on the beam L, D, B beam span, beam depth, beam width
tial differential equations. The difference equations that model the governing equations are simple to code and the global coefficient matrix possesses a banded structure that is amenable to efficient solution. Despite these ideal characteristics, finite-difference method has been supplanted in most solid mechanics applications in engineering by the finite-element method. This is mainly because of the limited flexibility of the method in managing the boundary conditions and shapes. However, the supremacy is brought back to the finite-difference technique here from the finite-element method through a new formulation of stress analysis of solid structures. As the new formulation establishes a priority of finite-difference technique over finite-element method, the philosophy and approach of the two methods are recapitulated here in brief.
1.1. Finite-element method Disadvantages: (a) The body is not in one piece but an assemblage of elements connected only at the nodes.
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
(b) Variations of parameters over individual elements are assumed to be very simple as in polynomials of limited order. (c) Values of parameters across element boundary are usually discontinuous. (d) Number of parameters (like ux, uy, uz) to be evaluated at each nodal point is usually more than that in finite-difference. (e) The total number of unknowns to be solved is much more than that in finite-difference method and hence the solution is likely to be crude and the computational time is more. (f) It cannot handle singularity. Advantages: (a) It works when other methods fail. (b) It is very good in managing boundaries.
complicated
1.2. Finite-difference method Disadvantages: (a) Derivatives in the governing equations are replaced by their finite-differences, resulting in some approximation. (b) Complicated boundary shapes and conditions are very difficult, almost impossible to handle. (c) It cannot handle singularity very well. Advantages: (a) The body is in one piece and thus there is no approximation on this account. (b) It permits reduction of parameters to be evaluated at the mesh points to one. (c) Total number of unknowns is less than that in finiteelement and hence, solution is more accurate and the computational time is less.
37
was much higher than that of FEM analysis. However, they noted that the computational effort of the finite-difference analysis, under the new boundary modeling system, is even somewhat greater than that of FEM analysis. The present paper describes the development of a general formulation for the finite-difference analysis of stresses and displacements in two-dimensional elastic solids. Basically, the present research is an attempt to extend the capability of our displacement potential formulation to the problems of composite materials. The formulation covers all the cases of general anisotropic, orthotropic and isotropic materials. Specifically, the formulation has been developed in such a fashion that the governing equations for the case of orthotropic and isotropic materials can be obtained as a special case of the general formulation for anisotropic materials. More importantly, the present scheme reduces the problem to finding a single parameter at the nodal points whereas the existing schemes find more than one at each node, and hence a tremendous saving in computational work is achieved. In addition, the present formulation provides the way to handle the mixed mode of boundary conditions very efficiently. Finally, the application of the present scheme is demonstrated for a classical problem of solid mechanics, and the solutions are compared with those of the existing methods.
2. Development of the general governing equations One of the ways to solve the stress problems of twodimensional elastic bodies of HookÕs material is to solve the displacement components ux and uy from the equilibrium equations. With reference to a rectangular coordinate system, in absence of body forces, the two differential equations of equilibrium for general anisotropic materials, in terms of the displacement components, are as follows [1,15]: o2 ux o2 ux o2 ux o2 uy þ Q66 2 þ Q16 2 þ 2Q16 2 ox oxoy oy ox o2 uy o2 uy þ Q12 þ Q66 þ Q26 2 ¼ 0 oxoy oy
Q11 It is noted that reliable and accurate prediction of critical stresses, mostly at the surfaces of engineering components, is of great importance to meet the severe demand for greater reliability as well as economic challenges in design. However, the uncertainties associated with the finite-element prediction of surface stresses have been clearly pointed out by several researchers [12–14]. Dow et al. [14] introduced a new boundary modeling approach for the finite-difference applications of displacement formulation of solid mechanics. In this connection, they reported that the accuracy of the finite-difference method in reproducing the state of stresses along the bounding surfaces
Q16
ð1aÞ
o2 ux o2 ux o2 ux o2 uy þ Q12 þ Q66 þ Q26 2 þ Q66 2 2 ox oxoy oy ox
þ 2Q26
o2 uy o2 uy þ Q22 2 ¼ 0 oxoy oy
ð1bÞ
where, the coefficients Qij are the elements of the transformed material stiffness matrix, ½Q, which are given by [15],
38
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
Q11 ¼ Q11 cos4 h þ 2ðQ12 þ 2Q66 Þsin2 hcos2 h þ Q22 sin4 h
gle function by introducing a scheme of reduction of unknowns. In order to attain that, the displacement components are expressed in terms of a potential function w of space variables as follows:
Q12 ¼ Q12 sin4 h þ cos4 h þ ðQ11 þ Q22 4Q66 Þsin2 hcos2 h 4
2
2
4
Q22 ¼ Q11 sin h þ 2ðQ12 þ 2Q66 Þsin hcos h þ Q22 cos h Q16 ¼ ðQ11 Q12 2Q66 Þ sin hcos3 h þ ðQ12 Q22 þ 2Q66 Þsin3 h cos h
ux ¼ a1
o2 w o2 w o2 w þ a2 þ a3 2 2 ox oxoy oy
ð4aÞ
uy ¼ a4
o2 w o2 w o2 w þ a þ a 5 6 ox2 oxoy oy 2
ð4bÞ
where, aÕs are unknown material constants. Substituting the above expressions of ux and uy in Eq. (1), we obtain the equilibrium equations in terms of the function w, which are
Q26 ¼ ðQ11 Q12 2Q66 Þsin3 h cos h þ ðQ12 Q22 þ 2Q66 Þ sin hcos3 h Q66 ¼ ðQ11 þ Q22 2Q12 2Q66 Þsin2 hcos2 h þ Q66 sin4 h þ cos4 h
E11 E211 ¼ ; 1 l12 l21 E11 l212 E22 E22 E11 E22 ¼ Q22 ¼ 1 l12 l21 E11 l212 E22 Q11 ¼
Q66 ¼ G12 ;
o4 w o4 w þ 2a3 Q16 þ a2 Q66 þ a6 Q12 þ a6 Q66 þ a5 Q26 2 2 ox oy oxoy 3 4 o w þ a3 Q66 þ a6 Q26 ¼0 ð5aÞ oy 4
l12 E22 l21 E11 l E11 E22 ¼ ¼ 12 1 l12 l21 1 l12 l21 E11 l212 E22
Q16 ¼ Q26 ¼ 0
The stress–strain relations for the general anisotropic materials are expressed through the transformed material stiffness matrix as follows [15]: 8 9 2 > Q11 < rxx > = 6 ryy ¼ 4 Q12 > > : ; rxy Q16
Q12 Q22 Q26
38 9 Q16 > < exx > = 7 Q26 5 eyy > : > ; cxy Q66
ð2Þ
In the above displacement formulation, the stresses are expressed in terms of the displacement components for their determination. The three stress–displacement relations for the plane problems of anisotropic materials, obtained from Eq. (2), are as follows: oux ouy oux ouy rxx ¼ Q11 þ Q12 þ Q16 þ ð3aÞ ox oy oy ox oux ouy oux ouy þ Q22 þ Q26 þ ryy ¼ Q12 ox oy oy ox
ð3bÞ
oux ouy oux ouy þ Q26 þ Q66 þ ox oy oy ox
ð3cÞ
rxy ¼ Q16
o4 w ox4
o4 w þ a2 Q11 þ 2a1 Q16 þ a5 Q16 þ a4 Q12 þ a4 Q66 ox3 oy þ a3 Q11 þ 2a2 Q16 þ a1 Q66 þ a6 Q16 þ a5 Q12 þ a5 Q66 þ a4 Q26
and
Q12 ¼ Q21 ¼
a1 Q11 þ a4 Q16
We now have to solve Eq. (1) with specified values of the functions rxx, ryy, rxy, ux and uy on segments of boundary with all possible combination of them. In the present approach, the two-dimensional problem of elasticity is reduced to the determination of a sin-
a1 Q16 þ a4 Q66
o4 w ox4
o4 w þ a2 Q16 þ a1 Q12 þ a1 Q66 þ a5 Q66 þ 2a4 Q26 ox3 oy þ a3 Q16 þ a2 Q12 þ a2 Q66 þ a1 Q26 þ a6 Q66 þ 2a5 Q26 þ a4 Q22 o4 w o4 w þ a3 Q12 þ a3 Q66 þ a2 Q26 þ 2a6 Q26 þ a5 Q22 2 2 ox oy oxoy 3 o4 w þ a3 Q26 þ a6 Q22 ¼0 ð5bÞ oy 4
Let us now choose aÕs in such a way that Eq. (5a) is automatically satisfied under all the circumstances. This will happen when coefficients of all the derivatives present in Eq. (5a) are individually zero. That is, when a1 Q11 þ a4 Q16 ¼ 0
ð6aÞ
a2 Q11 þ 2a1 Q16 þ a5 Q16 þ a4 Q12 þ a4 Q66 ¼ 0
ð6bÞ
a3 Q11 þ 2a2 Q16 þ a1 Q66 þ a6 Q16 þ a5 Q12 þ a5 Q66 þ a4 Q26 ¼ 0
ð6cÞ
2a3 Q16 þ a2 Q66 þ a6 Q12 þ a6 Q66 þ a5 Q26 ¼ 0
ð6dÞ
a3 Q66 þ a6 Q26 ¼ 0
ð6eÞ
Thus, for w to be a solution of the stress problem, it has to satisfy Eq. (5b) only. However, the aÕs of Eq. (4) must be known in advance. Here, we have five equations (Eq. (6)) for obtaining six unknown aÕs. We can thus assign an arbitrary value to one of these six unknowns and
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
solve for the remaining unknowns from Eq. (6). Let us assume a2 = 1 and, writing Eq. (6) in matrix form, we get 2
Q11
6 Q 6 66 6 6 2Q16 6 6 4 0
0
Q16
Q11
Q26
0
Q12 þ Q66
2Q16
0
Q66 9 8 0 > > > > > > > > > > > = < 2Q16 > ¼ Q11 > > > > > > > Q66 > > > > > ; : 0
0 Q12 þ Q66
Q16
0
Q26
0
0
y 2
38 9 a1 > > > > > > > 7> > Q16 a > 7> 3> 7< = 7 a4 0 7 > > > 7> a5 > > > Q12 þ Q66 5> > > > ; : > a 6 Q26
1
0
θ
x
ð7Þ
It is noted that, any other values of a2 could have been chosen leaving aside the homogeneous solutions. Different values of a2 will make the coefficients of the unknowns w different. The values of a2 which will lead to the diagonal dominance of the algebraic equations of unknown w at the nodal points are desirable in terms of increasing accuracy and shorter computational time. a2 = 1 leads to most simplified differential equations and hence more well-conditioned algebraic equations. The problem has now been reduced to the solution of a single function w from a single equation (Eq. (5b)). Let the equation be written as b1
39
o4 w o4 w o4 w o4 w o4 w þ b2 3 þ b3 2 2 þ b4 þ b5 4 ¼ 0 4 3 ox ox oy ox oy oxoy oy ð8Þ
then the coefficients bÕs are given by b1 ¼ a1 Q16 þ a4 Q66 b2 ¼ a1 Q12 þ Q66 þ a2 Q16 þ 2a4 Q26 þ a5 Q66
x,y,z - orthogonal loading axes
Z, 3
Fig. 1. Coordinate system for stress analysis of two-dimensional composite structure.
relations are thus obtained here as a special case of the general anisotropic material by substituting the value of h = 0. The stress–strain relations for orthotropic materials are expressed through the stiffness matrix, as follows [15]: 8 9 2 38 9 0 > Q11 Q12 > < rxx > < exx > = = 6 7 ryy ¼ 4 Q12 Q22 0 5 eyy ð9Þ > > > : : > ; ; cxy rxy 0 0 Q66 The two equilibrium equations, for the problems of orthotropic materials in terms of the displacement components, ux and uy, obtained from Eq. (1) are, Q11
o2 ux o2 uy o2 ux þ ðQ12 þ Q66 Þ þ Q66 2 ¼ 0 2 ox oxoy oy
ð10aÞ
Q66
o2 uy o2 ux o2 uy þ Q22 2 ¼ 0 þ ðQ12 þ Q66 Þ 2 ox oxoy oy
ð10bÞ
b3 ¼ a1 Q26 þ a2 Q12 þ Q66 þ a3 Q16 þ a4 Q22 þ 2a5 Q26 þ a6 Q66
b4 ¼ a2 Q26 þ a3 Q12 þ Q66 þ a5 Q22 þ 2a6 Q26 b5 ¼ a3 Q26 þ a6 Q22 Therefore, the governing differential equation (8) can be obtained explicitly once the values of aÕs of Eq. (7) are known. 3. Special cases of the general governing equation 3.1. Case-1: Orthotropic materials In this section, the governing equation of the w-formulation is deduced for the orthotropic materials where the principal material axes are aligned with the natural body axes, that is, the value of h is 0, in Fig. 1. All the
1,2,3 - principal material axes
The corresponding two equilibrium equations in terms of the function w, as obtained directly from Eq. (5) (or combining Eqs. (4) and (10)), are as follows: o4 w o4 w þ ða2 Q11 þ a4 Q12 þ a4 Q66 Þ 3 4 ox ox oy o4 w þ ða1 Q66 þ a3 Q11 þ a5 Q12 þ a5 Q66 Þ 2 2 ox oy o4 w o4 w þ ða2 Q66 þ a6 Q12 þ a6 Q66 Þ þ a3 Q66 4 ¼ 0 3 oxoy oy ð11aÞ
a1 Q11
40
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
o4 w o4 w þ ða1 Q12 þ a1 Q66 þ a5 Q66 Þ 3 4 ox ox oy o4 w þ ða2 Q12 þ a2 Q66 þ a4 Q22 þ a6 Q66 Þ 2 2 ox oy o4 w o4 w þ ða3 Q12 þ a3 Q66 þ a5 Q22 Þ þ a Q ¼0 6 22 oxoy 3 oy 4 ð11bÞ
a4 Q66
Likewise the case of general anisotropic material, assuming a2 = 1, the solutions for the remaining aÕs are obtained from the conditions required for the Eq. (11a) to be automatically satisfied. The values of aÕs thus obtained, are as follows: a1 ¼ a3 ¼ a5 ¼ 0 a2 ¼ 1 a4 ¼
a6 ¼
E211 Dn
ð12Þ
G12 ðE11 Dn
3.2. Case-2: Isotropic materials In addition to metallic materials, many of the particulate composites and also the randomly oriented discontinuous fiber composites can reasonably be treated as the isotropic materials. The corresponding equations for the problems of isotropic materials can readily be obtained from the general formulation of anisotropic materials by substituting the following conditions: 8 h ¼ 0o > > >
l212 E22 Þ
where, Dn ¼ l12 E11 E22 þ G12 ðE11 l212 E22 Þ. Note that the solutions for aÕs can also be obtained from the matrix Eq. (7) when the condition of orthotropic material is inserted. Knowing the values of aÕs, the coefficients of different derivatives of Eq. (11b), that is, bÕs, are obtained as follows: E2 G12 b1 ¼ a4 Q66 ¼ 11 Dn
lE=1 l2 E=1 l2 0
38 9 0 > < exx > = 7 e 0 5 yy > : > ; cxy E=2ð1 þ lÞ
ð15Þ
Similarly, the two equilibrium equations in terms of the displacement components, obtained from Eq. (1) are o2 ux 1 þ l o2 uy 1 l o2 ux þ ¼0 þ ox2 oxoy oy 2 2 2
ð16aÞ
o2 uy 1 þ l o2 ux 1 l o2 uy þ þ ¼0 2 2 oy 2 oxoy ox2
ð16bÞ
b2 ¼ b4 ¼ 0
The corresponding equilibrium equations in terms of the function w, obtained from Eqs. (5) and (14), are given by
b3 ¼ a2 ðQ12 þ Q66 Þ þ a4 Q22 þ a6 Q66
a1
¼
E11 E22 ð2l12 G12 E11 Þ Dn
b5 ¼ a6 Q22 ¼
E11 E22 G12 Dn
Substituting the values of bÕs in Eq. (8) the explicit expression of the governing differential equation for the case of orthotropic materials becomes
E11 G12
o4 w o4 w þ E22 ðE11 2l12 G12 Þ 2 2 4 ox ox oy
þ E22 G12
o4 w ¼0 oy 4
ð13Þ
E o4 w 1 l2 ox4 4 E lE E ow þ a2 þ a4 þ 2 2 1l 1l 2ð1 þ lÞ ox3 oy 4 E E lE E ow þ a5 þ a3 þ a1 þ 2 2 1l 2ð1 þ lÞ 1l 2ð1 þ lÞ ox2 oy 2 4 E lE E ow þ a6 þ a2 þ 2ð1 þ lÞ 1 l2 2ð1 þ lÞ oxoy 3 þ a3
E o4 w ¼0 2ð1 þ lÞ oy 4
ð17aÞ
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51 a4
E o4 w 2ð1 þ lÞ ox4 4 lE E E ow þ a5 þ a1 þ 2 1l 2ð1 þ lÞ 2ð1 þ lÞ ox3 oy 4 lE E E E ow þ a4 þ a2 þ þ a6 2 2 1l 2ð1 þ lÞ 1l 2ð1 þ lÞ ox2 oy 2 4 lE E E ow þ a þ a3 þ 5 1 l2 2ð1 þ lÞ 1 l2 oxoy 3 E o4 w þ a6 ¼0 1 l2 oy 4
n y φ
A
(90-φ)
ut
un O C
uy ux t
Boundary segment
ð17bÞ
As before, assuming a2 = 1, the corresponding values of the aÕs for the case of isotropic materials can be obtained from Eq. (12) or directly from Eq. (17a), and are as follows:
41
x Fig. 2. Components of displacement on the boundary segment, AOC.
a1 ¼ a3 ¼ a5 ¼ 0 a2 ¼ 1
a6 ¼
2 1þl
–Y ð18Þ
y
1l 1þl
φ
σxx
When the above values of aÕs and the conditions of Eq. (14) are substituted in the general expressions of bÕs in Eq. (8), we get b1 ¼
b5 ¼
σxy
–X a
φ C
E
a sinφ
t B
σxy σyy
ð1 þ lÞ2
(90-φ)
x
b2 ¼ b4 ¼ 0 b3 ¼
n
A
a cosφ
a4 ¼
Fig. 3. Components of stress on the boundary segment, AB.
2E ð1 þ lÞ2
4. Boundary conditions
E ð1 þ lÞ2
Substitution of the above values of bÕs in Eq. (8) gives the governing differential equation for the case of isotropic material as o4 w o4 w o4 w þ2 2 2þ 4 ¼0 4 ox ox oy oy
ð19Þ
Note that the above governing bi-harmonic equation for isotropic materials, obtained as a special case of the general formulation for anisotropic materials, is found to be the identical to that used in our displacement potential function formulation for mixed boundary-value problems of elasticity. Reliability as well as suitability of that formulation has already been verified in our previous researches [7–10,16,17].
The physical conditions specified along the edge or boundary of a body are either known restraints or forces or some kind of combination of the two. Both the restraints, that is, known displacements and the forces are defined by their respective components. They are (Figs. 2 and 3), (a) (b) (c) (d)
Normal displacement (un) Tangential displacement (ut) Normal stress (rn) Tangential stress (rt)
The solution governing equation requires two known boundary conditions at any point on the boundary. Therefore, the possible boundary conditions created out of the four components are (1) (rn, rt), (2) (rn, un), (3) (rn, ut), (4) (rt, un), (5) (rt, ut), and (6) (un, ut).
42
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
However, the combinations (rn, un) and (rt, ut) do not occur in practice. The remaining four possible boundary conditions are therefore (1) (rn, rt), (2) (rn, ut), (3) (rt, un), and (4) (un, ut).
Although the known conditions are usually specified in terms of the normal and tangential components of stress or displacement of a point on the boundary, the parameters within the body are known as x- and y-components of stress and displacement. Since our objective is to solve the problem in terms of a single function w, all the boundary conditions of interest are required to be expressed in terms of the function. The components of displacement in terms of w are ux ðx; y Þ ¼ a1
o2 w o2 w o2 w þ a þ a 2 3 ox2 oxoy oy 2
ð20aÞ
uy ðx; y Þ ¼ a4
o2 w o2 w o2 w þ a5 þ a6 2 2 ox oxoy oy
ð20bÞ
Expressions for the components of stress for the interior points of an anisotropic body are obtained when the displacement components in the stress–displacement relations (Eq. (3)) are replaced by Eq. (4). Thus, we get 3
o w rxx ðx; yÞ ¼ a1 Q11 þ a4 Q16 ox3 o3 w þ a1 Q16 þ a2 Q11 þ a4 Q12 þ a5 Q16 ox2 oy o3 w þ a2 Q16 þ a3 Q11 þ a5 Q11 þ a6 Q16 oxoy 2 o3 w þ a3 Q16 þ a6 Q12 oy 3 ð21aÞ 3
o w ryy ðx; yÞ ¼ a1 Q12 þ a4 Q26 ox3 o3 w þ a1 Q26 þ a2 Q12 þ a4 Q22 þ a5 Q26 ox2 oy o3 w þ a2 Q26 þ a3 Q12 þ a5 Q22 þ a6 Q26 oxoy 2 3 o w þ a3 Q26 þ a6 Q22 oy 3 ð21bÞ
o3 w rxy ðx; yÞ ¼ a1 Q16 þ a4 Q66 ox3 o3 w þ a1 Q66 þ a2 Q16 þ a4 Q26 þ a5 Q66 ox2 oy o3 w þ a2 Q66 þ a3 Q16 þ a5 Q26 þ a6 Q66 oxoy 2 3 o w þ a3 Q66 þ a6 Q26 oy 3 ð21cÞ 4.1. Normal and tangential components of displacement In Fig. 2, the positive normal direction on the boundary is outward, positive tangential direction is clockwise and the positive /, the angle drawn from positive normal to positive x-axis, is clockwise. With these conventions, the relations between interior and boundary components of displacement can be written as ux ¼ un cos / þ ut sin / uy ¼ un sin / ut cos /
ð22Þ
Solving Eqs. (22) for un and ut, we have, un ðx; yÞ ¼ ux cos / þ uy sin / ut ðx; yÞ ¼ ux sin / uy cos /
ð23Þ
Combining Eqs. (20) and (23), the expressions for the normal and tangential components of displacement in terms of w, on the boundary of an anisotropic material, are given by un ðx; yÞ ¼ R1
o2 w o2 w o2 w þ R3 2 þ R2 2 ox oxoy oy
ð24aÞ
ut ðx; yÞ ¼ R4
o2 w o2 w o2 w þ R þ R 5 6 ox2 oxoy oy 2
ð24bÞ
where R1 ¼ a1 cos / þ a4 sin / R2 ¼ a2 cos / þ a5 sin / R3 ¼ a3 cos / þ a6 sin / R4 ¼ a1 sin / a4 cos / R5 ¼ a2 sin / a5 cos / R6 ¼ a3 sin / a6 cos /
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
4.2. Normal and tangential components of stress
ð25Þ
Y ¼ ryy sin / þ rxy cos /
Again, the normal and tangential components of stress in terms of the components, X and Y , can be written as rn ¼ X cos / þ Y sin /
Q11 Q12 sin 2/ Q16 cos 2/ 2 Q Q22 T 5 ¼ 12 sin 2/ Q26 cos 2/ 2 Q Q26 T 6 ¼ 16 sin 2/ Q66 cos 2/ 2 T4 ¼
With reference to Fig. 3, if the condition of equilibrium is applied for the forces on a segment of boundary AB, the components of stress along the direction of coordinate axes are X ¼ rxx cos / þ rxy sin /
43
As far as numerical computation is concerned, it is thus evident from the above expressions that all the boundary conditions of interest can easily be discretized in terms of the single function w in the method of finite-difference.
ð26Þ
rt ¼ X sin / Y cos /
5. Special cases of the boundary expressions
or 5.1. Case-1: Orthotropic materials
rn ¼ rxx cos2 / þ ryy sin2 / þ rxy sin 2/ r r xx yy sin 2/ rxy cos 2/ rt ¼ 2
ð27Þ
Now, if the expressions of body parameters, namely rxx, ryy, and rxy of Eq. (21), are substituted in Eq. (27), the normal and tangential components of stress in terms of w are 3
3
3
ux ðx; yÞ ¼
3
ow ow ow ow rn ðx; yÞ ¼ P 1 3 þ P 2 2 þ P 3 þ P4 3 ox ox oy oxoy 2 oy
Explicit expressions for the body parameters in terms of the function w are obtained from the general equations (20) and (21) when the respective conditions of orthotropic materials are substituted in them. The conditions are,
ð28aÞ
o2 w oxoy
uy ðx; yÞ ¼
where,
o2 w 1 o2 w E211 2 þ G12 E11 l212 E22 Dn ox oy 2
P 1 ¼ a1 T 1 þ a4 T 3
ð29bÞ E11 G12 o3 w o3 w rxx ðx; yÞ ¼ E11 2 l12 E22 3 ox oy oy Dn
P 2 ¼ a1 T 3 þ a2 T 1 þ a4 T 2 þ a5 T 3 P 3 ¼ a2 T 3 þ a3 T 1 þ a5 T 2 þ a6 T 3 P 4 ¼ a3 T 3 þ a6 T 2
ryy ðx; yÞ ¼
ð29dÞ
T 2 ¼ Q12 cos2 / þ Q22 sin2 / þ Q26 sin 2/
rxy ðx; yÞ ¼
2
T 3 ¼ Q16 cos / þ Q26 sin / þ Q66 sin 2/ and rt ðx; yÞ ¼ P 5
o3 w o3 w o3 w o3 w þ P þ P þ P 6 7 8 ox3 ox2 oy oxoy 2 oy 3
where P 5 ¼ a1 T 4 þ a4 T 6 P 6 ¼ a1 T 6 þ a2 T 4 þ a4 T 5 þ a5 T 6 P 7 ¼ a2 T 6 þ a3 T 4 þ a5 T 5 þ a6 T 6 P 8 ¼ a3 T 6 þ a6 T 5
ð29cÞ
E11 E22 o3 w o3 w ðl12 G12 E11 Þ 2 G12 3 Dn ox oy oy
T 1 ¼ Q11 cos2 / þ Q12 sin2 / þ Q16 sin 2/
2
ð29aÞ
ð28bÞ
E11 G12 o3 w o3 w E11 3 l12 E22 ox oxoy 2 Dn
ð29eÞ
Similarly, evaluating the coefficients RÕs for the orthotropic materials and substituting in Eq. (24), the expressions for the normal and tangential components of displacement are E211 sin / o2 w o2 w þ cos / 2 Dn ox oxoy 2 G12 E11 l12 E22 sin / o2 w Dn oy 2
un ðx; yÞ ¼
ð30aÞ
44
ut ðx; yÞ ¼
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
E211 cos / o2 w o2 w þ sin / 2 Dn ox oxoy G12 E11 l212 E22 cos / o2 w þ Dn oy 2
rxy ðx; yÞ ¼ ð30bÞ
When the respective values of the coefficients PÕs and TÕs are substituted in the general equation (28), the expressions for the normal and tangential components of stress for points on the boundary of an orthotropic body are,
ð32eÞ
Evaluating the corresponding values of the coefficients RÕs for the case of isotropic material and substituting them in Eq. (24), the expressions for the normal and tangential components of displacement are 2 o2 w o2 w sin / 2 þ cos / 1þl ox oxoy
un ðx; yÞ ¼
E211 G12 sin 2/ o3 w Dn ox3 E11 E11 G12 cos2 / E11 E22 sin2 / þ l12 E22 G12 sin2 / þ Dn
rn ðx; yÞ ¼
1l o2 w sin / 2 1þl oy
ut ðx; yÞ ¼
o3 w l E11 E22 G12 sin 2/ o3 w 2 þ 12 ox oy Dn oxoy 2 E11 E22 G12 l12 cos2 / þ sin2 / o3 w oy 3 Dn
3 ow o3 w þ l ox3 oxoy 2 ð1 þ l Þ2 E
ð33aÞ
2 o2 w o2 w cos / 2 þ sin / 1þl ox oxoy þ
1l o2 w cos / 2 1þl oy
ð33bÞ
ð31aÞ E2 G12 cos 2/ o3 w rt ðx; yÞ ¼ 11 Dn ox3 E11 ðE11 G12 þ E11 E22 l12 E22 G12 Þ sin 2/ þ 2Dn
o3 w l E11 E22 G12 cos 2/ o3 w 12 2 ox oy Dn oxoy 2
Similarly, evaluating the coefficients PÕs and TÕs for the case of isotropic material, the expressions for the normal and tangential components of stress for points on the boundary of an isotropic body, are obtained from Eq. (28) as follows: rn ðx; yÞ ¼
E11 E22 G12 ð1 l12 Þ sin 2/ o3 w þ 2Dn oy 3
þ ð31bÞ þ
5.2. Case-2: Isotropic materials Here, the boundary conditions for the case of isotropic materials are derived as a special case of those for the general anisotropic material. Substituting the conditions of Eq. (14) along with the respective aÕs of Eq. (18), expressions for the body parameters in terms of the function w can be obtained from the general Eqs. (20) and (21) as follows: ux ðx; yÞ ¼
o2 w oxoy
uy ðx; yÞ ¼
rxx ðx; yÞ ¼
ð32aÞ
2 1 ow o2 w 2 2 þ ð1 lÞ 2 1 þ l ox oy
E ð1 þ lÞ2
ryy ðx; yÞ ¼
o3 w o3 w l 3 2 ox oy oy
o3 w o3 w þ ð 2 þ l Þ ox2 oy oy 3 ð1 þ lÞ2 E
ð32bÞ
ð32cÞ
ð32dÞ
E ð1 þ lÞ
sin 2/
2
E ð1 þ lÞ
2
El
E
cos2 / 2sin2 / lsin2 /
sin 2/
ð1 þ lÞ2 ð1 þ lÞ2
o3 w ox3 o3 w ox2 oy
o3 w oxoy 2
lcos2 / þ sin2 /
o3 w oy 3 ð34aÞ
rt ðx; yÞ ¼
E 2
cos 2/
o3 w Eð3 þ lÞ þ sin 2/ ox3 2ð1 þ lÞ2
ð1 þ l Þ o3 w El o3 w cos 2/ 2 ox oy ð1 þ lÞ2 oxoy 2 þ
Eð1 lÞ 2ð1 þ lÞ
2
sin 2/
o3 w oy 3
ð34bÞ
One can thus easily realize that all the boundary conditions for the problems of orthotropic as well as isotropic materials can be deduced from those of the general expressions of anisotropic material. It is noted that the differential equations of the boundary conditions for the case of isotropic material are found to be identical with those used in our w-formulation [10,16,17].
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51 Stencil for both orthotropic and isotropic material
(a)
45
(b)
i-forward & j-forward j
k h
k i j j
y
i
y i
h
x
Stencil for anisotropic material
x i-forward & j-forward
(c)
h k j
y i
x Fig. 4. (a) Finite-difference grid structure for the governing equilibrium equations; (b) Finite-difference stencil for the normal and tangential components of displacements on the boundary (Eq. (37)); (c) Finite-difference stencil for the normal and tangential components of stress on the boundary (Eq. (38)).
6. Finite-difference discretization of the governing equations The present computational scheme involves evaluation of a single function w at the mesh points of a uniform rectangular mesh network over a geometrically regular or irregular region. The governing differential equation (8), which is used to evaluate the function w at the internal mesh points of an anisotropic body, can be expressed in its corresponding difference form using central difference operators. The complete finite-difference expression of the governing equation is given by,
46
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
24b1 n4 þ 16b3 n2 þ 24b5 wði; jÞ 2
2
8n ð2n b1 þ b3 Þfwði þ 1; jÞ þ wði 1; jÞg þ b2 n3 fwði þ 2; j þ 1Þ wði þ 2; j 1Þ þwði 2; j 1Þ wði 2; j þ 1Þg þ b4 nfwði þ 1; j þ 2Þ wði 1; j þ 2Þþwði 1; j 2Þ wði þ 1; j 2Þg þ 4b1 n4 fwði þ 2; jÞ þ wði 2; jÞg þ 4b5 fwði; j þ 2Þ þwði; j 2Þg þ 2b2 n3 þ 4b3 n2 2b4 n fwði þ 1; j þ 1Þ wði 1; j 1Þg þ 2b2 n3 þ 4b3 n2 þ 2b4 n
expressions (for example, i- & j-forward, i- & j-backward, i-forward & j-backward, and i-backward & j-forward), for each of the boundary conditions, are developed for different regions of boundary. For example, the i- and j-forward difference formulae for the normal and tangential components of displacement on the boundary of an anisotropic body, corresponding to Eq. (24), are written as follows: R2 fwði þ 2; j þ 2Þ 4wði þ 2; j þ 1Þ 4hk þ 16wði þ 1; j þ 1Þ 4wði þ 1; j þ 2Þ R1 þ3wði; j þ 2Þ þ 3wði þ 2; jÞg þ 2 wði 1; jÞ h R1 3R2 R3 wði þ 1; jÞ þ 2 wði; j 1Þ þ 2 hk k h R3 3R2 wði; j þ 1Þ þ 2 hk k 9R2 2R1 2R3 2 2 wði; jÞ ð37aÞ þ 4hk h k
un ðx; yÞ ¼
fwði þ 1; j 1Þ þ wði 1; j þ 1Þg 8ðn2 b3 þ 2b5 Þ fwði; j þ 1Þ þ wði; j 1Þg ¼ 0 ð35Þ
where, n = k/h (Fig. 4). The corresponding grid structure of the governing equation (35) for an arbitrary internal mesh point (i, j) is shown in Fig. 4(a). It is thus clear that, application of the governing equation involves a total of 21 neighboring mesh points, and when this point of application (i, j) becomes an immediate neighbor to the physical boundary, the equation will involve mesh points not only interior but also exterior to the physical boundary. After substitution of the respective values of bÕs for the orthotropic or isotropic materials into Eq. (35), the corresponding central difference equations can easily be obtained. The finite-difference approximation of the governing equation thus obtained for isotropic materials is, 2n2 fwði þ 1; j þ 1Þ þ wði þ 1; j 1Þ þ wði 1; j þ 1Þ þwði 1; j 1Þg þ wði; j þ 2Þ 4n2 n2 þ 1 fwði þ 1; jÞ þ wði 1; jÞg 4 n2 þ 1 fwði; j þ 1Þ þ wði; j 1Þg þ n4 fwði þ 2; jÞ þwði 2; jÞg þ wði; j 2Þ þ 6n4 þ 8n2 þ 6 wði; jÞ ¼ 0 ð36Þ
The corresponding grid structure of Eq. (36), as shown in Fig. 4(a), involves a total of 13 neighboring mesh points. From the figures, it is clear that the use of the above central difference stencils of the governing equations, especially for the points in the immediate neighborhood of the physical boundary, demands the inclusion of an imaginary (false) boundary, exterior to the physical boundary of the domain concerned. As the differential equations associated with the boundary conditions contain second- and third-order partial derivatives of the function w, the use of centraldifference expressions eventually leads to the inclusion of points exterior to the false boundary. In order to avoid the occurrence, different versions of the finite-difference expressions like forward, backward and central have been adopted in a combined form to generate the difference equations for the boundary conditions. It is noted here that the order of local truncation error has been kept the same for all the expressions developed, that is, O(h2). Four different sets of finite-difference
R5 fwði þ 2; j þ 2Þ þ 16wði þ 1; j þ 1Þ 4hk 4wði þ 2; j þ 1Þ 4wði þ 1; j þ 2Þ R4 þ3wði; j þ 2Þ þ 3wði þ 2; jÞg þ 2 wði 1; jÞ h R4 3R5 R6 wði þ 1; jÞ þ 2 wði; j 1Þ þ 2 hk k h R6 3R5 wði; j þ 1Þ þ 2 hk k 9R5 2R4 2R6 2 2 wði; jÞ ð37bÞ þ 4hk h k
ut ðx; yÞ ¼
Similarly, from Eq. (28), the corresponding i- and j-forward difference formulae for the normal and tangential components of stress on the boundary are as follows: 5P 1 3P 2 3P 3 5P 4 rn ðx; yÞ ¼ wði; jÞ þ þ þ h3 h2 k hk 2 k3 2P 2 2P 3 wði þ 1; j þ 1Þ þ þ h2 k hk 2 6P 1 3P 2 4P 3 þ 2 þ 2 wði þ 1; jÞ h3 2h k hk 4P 2 3P 3 6P 4 wði; j þ 1Þ þ þ h2 k 2hk 2 k3 3P 1 3P 2 þ wði 1; jÞ 2h3 2h2 k 3P 3 3P 4 wði; j 1Þ þ 2hk 2 2k 3 3P 1 P 3 þ 2 wði þ 2; jÞ þ h3 hk
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
P 2 3P 4 P1 wði; j þ 2Þ 3 wði þ 3; jÞ þ 2 3 hk k 2h P2 þ 2 f4wði 1; j þ 1Þ wði 1; j þ 2Þ 2h k P4 wði þ 1; j þ 2Þg 3 wði; j þ 3Þ 2k P3 f4wði þ 1; j 1Þ wði þ 2; j þ 1Þ þ 2hk 2 wði þ 2; j 1Þg
½Am;m fX gm ¼ fBgm
þ
47
ð39Þ
for m unknowns; and
½Amþ1;mþ1 fX gmþ1 ¼ fBgmþ1
for m þ 1 unknowns
ð40Þ
To reduce Eq. (40) to (39), we have to perform (i) (m + 1)(m + 1) = (m + 1)2 divisions and (ii) (m + 1)m subtractions ð38aÞ
5P 5 3P 6 3P 7 5P 8 þ 2 þ 2 þ 3 wði; jÞ h3 h k hk k 2P 6 2P 7 þ wði þ 1; j þ 1Þ þ h2 k hk 2 6P 5 3P 6 4P 7 þ 2 þ 2 wði þ 1; jÞ h3 2h k hk 4P 6 3P 7 6P 8 wði; j þ 1Þ þ þ h2 k 2hk 2 k3 3P 5 3P 6 þ wði 1; jÞ 2h3 2h2 k 3P 7 3P 8 wði; j 1Þ þ 2hk 2 2k 3 3P 5 P 7 þ þ 2 wði þ 2; jÞ h3 hk P 6 3P 8 P5 þ 2 þ 3 wði; j þ 2Þ 3 wði þ 3; jÞ hk k 2h P6 þ 2 f4wði 1; j þ 1Þ wði 1; j þ 2Þ 2h k P8 wði þ 1; j þ 2Þg 3 wði; j þ 3Þ 2k P7 f4wði þ 1; j 1Þ wði þ 2; j þ 1Þ þ 2hk 2 wði þ 2; j 1Þg ð38bÞ
rt ðx; yÞ ¼
The corresponding finite-difference stencils for the normal and tangential components of displacement and stress for points on the boundary are illustrated in Fig. 4(b) and (c). It is interesting to note that the boundary stencils for the displacement and stress components take identical form for the general anisotropic, orthotropic and isotropic materials.
When m = 0, the increase in computation due to addition of one unknown is one division. That is, when we have only one unknown, computation involved is one division. When m = 1, the increase is 4 divisions and 2 subtractions. Let us consider the total number of nodes = N; Number of parameters at each node (like ux, uy) = P. Therefore, total number of unknowns, m, in Eq. (39) is given by m = NP. For a two-dimensional problem, in both existing finite-difference and finite-element techniques, P = 2. Here, in the present finite-difference scheme, we have to find only one parameter, w, i.e., P = 1. So the reduction of the unknowns between the two techniques are from 2N of finite-element to N of the present finite-difference scheme. The increase in computation due to the additional N unknowns is
(a)
Uniformly distributed load, w
D
x
0.1L
0.8L
0.1L
y
(b)
Physical boundary
Imaginary boundary
7. Scheme for reduction of computational work x, i
This section explains how the computational work is reduced due to reduction in the number of parameters to be evaluated at the nodal points. Both in finite-element and in finite-difference methods, 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
k h
y, j
Fig. 5. (a) Loading and geometry of a simply supported deep beam; (b) Finite-difference discretization of the beam domain.
48
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
(i) (N + 1)2 + (N + 2)2 + (N + 3)2 + + (2N)2 [division], and (ii) N(N + 1) + (N + 1)(N + 2) + (N + 2)(N + 3) + + (2N1)(2N) [subtraction] Let us now evaluate the percentage reduction in the computational efforts for the individual cases of division and subtraction. (i) Division:
ffi 87:5% when N is large
(ii) Subtraction: Saving in the number of subtraction 2N ð2N 1Þð4N 1Þ 2Nð2N 1Þ þ ¼ 6 2 N ðN 1Þð2N 1Þ N ðN 1Þ þ 6 2
Saving in the number of division ¼
Percentage saving in division alone ð2N Þð2N þ 1Þð4N þ 1Þ N ðN þ 1Þð2N þ 1Þ 100 ¼ 2N ð2N þ 1Þð4N þ 1Þ N þ1 100% ¼ 1 2ð4N þ 1Þ
ð2NÞð2N þ 1Þð4N þ 1Þ N ðN þ 1Þð2N þ 1Þ 6 6
1.5
1.0 σxx / (w/B)
σxx / (w/B)
1.0
1.5 FDM FEM Simple theory
0.5 0.0 -0.5 -1.0 -0.6
-0.4
-0.2
0.0 y/D
0.2
0.4
0.0
-1.0 -0.6
0.6
-0.4
-0.2
(b)
0.0 y/D
0.2
0.4
0.6
0.2
0.4
0.6
3
2 FDM FEM Simple theory
2
σxx / (w/B)
σxx / (w/B)
0.5
-0.5
(a)
1
FDM FEM Simple theory
0
FDM FEM Simple theory
1 0 -1
-1
-2 -2 -0.6
(c)
-0.4
-0.2
0.0 y/D
0.2
0.4
-3 -0.6
0.6
(d)
-0.4
-0.2
0.0 y/D
Fig. 6. Normalized bending stress distributions at the midspan section for the beam: (a) L/D = 0.5; (b) L/D = 1.0; (c) L/D = 1.5; (d) L/D = 2.0.
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
Percentage saving in subtraction alone " # 2 2N ð4N 2 1Þ N ðN3 1Þ 3 100 ¼ 2N ð4N 2 1Þ 3
ðN 2 1Þ 100% ¼ 1 2ð4N 2 1Þ ffi 87:5% when N is large It is thus verified that the present finite-difference method of solution can save approximately 87% of computational work for the solution of a two-dimensional problem, in comparison with that required for usual finite-element solution. 8. Example of application of the w-formulation Using the present w-formulation a classical problem of solid mechanics has been solved numerically and the results are compared with the available solutions in the literature. The present isotropic problem has been
49
solved here as a special case of the general anisotropic material. The problem is a deep beam of rectangular cross section, simply supported at the ends and carries a uniformly distributed load over its whole span. The geometry and loading are illustrated in Fig. 5(a). It is noted that the analysis of the present problem is being constantly looked into by several researchers in the field of solid mechanics [3,4,11,18,19]. However, each solution possesses certain limitations, and eventually none of the solutions are found to conform to all the physical characteristics of the problem appropriately. A relatively small number of mesh points (23 · 23) are considered here for the finite-difference solution of the problem, so that the solutions can be readily compared with the FEM results at identical sections of the beams. The finite-difference mesh-network used for the present problem is shown in Fig. 5(b), where an imaginary (false) boundary, exterior to the physical boundary, is included. It is noted that the error of the present computational approach is of the order of h2. It has been
0.6
0.5 FDM FEM Simple theory
0.4
0.4 σxy / (w/B)
σxy / (w/B)
0.3 0.2 0.1
0.2
0.0
0.0
FDM FEM Simple theory
-0.1 -0.6
-0.4
-0.2
(a)
0.0 y/D
0.2
0.4
-0.2 -0.6
0.6
-0.4
-0.2
(b)
0.0 y/D
0.2
0.4
0.6
0.4
0.6
1.0
0.6
0.4
σxy / (w/B)
σxy / (w/B)
0.8
0.2 FDM FEM Simple theory
0.0
0.6 0.4 0.2
FDM FEM Simple theory
0.0
-0.2 -0.6
(c)
-0.4
-0.2
0.0 y/D
0.2
0.4
0.6
-0.6
(d)
-0.4
-0.2
0.0 y/D
0.2
Fig. 7. Normalized shearing stress distributions at section: (a) x/L = 0.25 for the beam, L/D = 0.5; (b) x/L = 0.25 for the beam, L/D = 1.0; (c) x/L = 0.25 for the beam, L/D = 1.5; (d) x/L = 0.25 for the beam, L/D = 2.0.
50
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
Table 1 Comparison of maximum bending stress at section, x/L = 0.5 Beam size
Normalized maximum bending stress Simple theory
Tension side, at y/D = 0.5 L/D = 0.5 0.15 L/D = 1.0 0.60 L/D = 1.5 1.35 L/D = 2.0 2.40 Compression side, at y/D = 0.5 L/D = 0.5 0.15 L/D = 1.0 0.60 L/D = 1.5 1.35 L/D = 2.0 2.40
FEM
FDM
1.27 1.30 1.69 2.44
1.33 1.35 1.71 2.46
0.32 0.50 1.22 2.25
0.27 0.47 1.24 2.24
Table 2 (a) Comparison of maximum shearing stress at section, x/ L = 0.25; (b) Comparison of shearing stress at the section, y/ D = 0.5, x/L = 0.25 Beam size
(a) L/D = 0.5 L/D = 1.0 L/D = 1.5 L/D = 2.0
Normalized maximum shearing stress Simple theory
FEM
FDM
0.19 0.38 0.56 0.75
0.41 0.43 0.49 0.78
0.40 0.43 0.57 0.76
Normalized shearing stress (b) L/D = 0.5 L/D = 1.0 L/D = 1.5 L/D = 2.0
0.00 0.00 0.00 0.00
0.06 0.18 0.27 0.34
0.00 0.00 0.00 0.00
0.2 0.0
σyy / (w/B)
-0.2
L / D = 0.5 L / D = 1.0 L / D = 1.5 L / D = 2.0
-0.4 -0.6 -0.8 -1.0 -1.2 -0.6
-0.4
-0.2
0.0 y/D
0.2
0.4
0.6
Fig. 8. Finite-difference solutions for the normal stress component, ryy at section, x/L = 0.25, of the beams.
however confirmed by our FDM results that no significant change in accuracy is observed with decreasing h, when it is sufficiently small, as used in the present solution. The governing differential equation is used to evaluate the function w only at the internal mesh points. The discrete values of the displacement potential function w(x, y), at the mesh points of the domain concerned are obtained from a system of linear algebraic equations, resulting from the discretization of the governing equation and the associated boundary conditions, by the use of direct method of solution (CholeskeyÕs triangular-decomposition method). Finally, all the parameters of interest in the solution, namely, displacement and stress components, are obtained as soon as w is known at the mesh points of the domain, as all the components are expressed as summation of different derivatives of the function w. In the case of present problem, stresses are normalized with respect to the uniform loading w, and beam width B, consistent with that of Ref. [19]. In order to obtain the numerical results, values for YoungÕs modulus and PoissonÕs ratio of 209 GPa and 0.3, respectively, are used. The present solutions are compared with the solutions of simple beam theory and also with the finite-element solutions [19]. The finite-element solutions were obtained by using a regular mesh of 200, 8-noded, plane stress, isoparametric elements (20 · 10). Fig. 6 presents the bending stress distributions at the mid-sections of the beams of aspect ratio, L/D = 0.5, 1.0, 1.5, and 2.0. Along with the present finite-difference solutions, the FEM results [19] and also the results of simple bending theory of beam are included for the convenience of comparison. The mid-section has been chosen here for consideration since it is the section of maximum bending moment and also it is almost free from the local influence of the loading. It is seen that, for all the beams, the bending stress curves agree reasonably well with the finite-element solutions but differ significantly, as expected, from those of the simple flexure theory. For higher L/D ratios (L/D > 1.0) the finite-element predictions are found to be identical with those of our FDM results. As expected, the nonlinearity in the distributions as well as the deviations from the simple theory is found to decrease as the L/D ratio is increased. The distribution of shearing stress at section x/ L = 0.25 is shown in Fig. 7 for different beams. The parabolic distributions based on simple theory and the corresponding FEM solutions are included in the figures. For small L/D ratios, the shear stress distributions are significantly different from that obtained from the simple theory. The maximum shear stress predicted by our FDM approach is higher for all the beams inspected than those predicted by the simple theory. However, the deviation is found to decrease as L/D ratio increases. The FEM solutions, in general, agree well with the pre-
S.R. Ahmed et al. / Computers and Structures 83 (2005) 35–51
sent FDM results. However, it is interesting to note that the finite-element prediction of maximum shear stress for the beam of L/D = 1.5 is found to take a lower value even than that predicted by simple theory, which is in contrast to the solutions obtained for beam of other L/ D ratios. More importantly, the FEM predictions of shear stress in the neighborhood of both top and bottom surfaces are highly unreliable, even though the section concerned is free from singularity. This discrepancy is found to remain quite significant even in the case of the beam of L/D = 2.0. Therefore, the present investigation verifies the fact that the present finite-difference method of solution is free from the inability of FEM in predicting the actual state of stresses at the surfaces. Tables 1 and 2 list some of the numerical data for the bending and shear stresses at different sections of the beams; one can easily compare the accuracy of the present solutions with those of FEM and simple beam theory. Finally, Fig. 8 presents the distributions of normal stress component ryy at section x/L = 0.25. It should be noted here that the finite-difference analysis reproduces the normal stress exactly both at the top as well as at bottom boundaries of all the beams considered. 9. Conclusions The present paper demonstrates the extension of the displacement potential formulation of isotropic materials to the case of general anisotropic materials. It has been shown that the governing equations for the isotropic as well as orthotropic materials can readily be obtained as a special case of the proposed general formulation. The present approach of stress analysis for the plane problems of solid mechanics by finite-difference technique saves us approximately 87% of computational effort in comparison to that taken by exiting FE and FD methods. In addition, the quality of solutions will be improved significantly because of the fact that the solution of any algebraic system deteriorates with increasing number of unknowns and eventually fails when the number of unknowns become very large. The present approach not only increases our horizon of managing bigger problems in stress analysis but also in other fields of science and engineering. Summarizing, we would like to point out that the use of the present approach is expected to be the most promising one to save computational effort as well as to increase the horizon of capabilities in terms of solving bigger problems with higher accuracy.
51
References [1] Timoshenko D, Goodier JN. Theory of elasticity. 3rd ed. New York, NY: McGraw-Hill Book Company; 1979. [2] Uddin MW. Finite difference solution of two-dimensional elastic problems with mixed boundary conditions. M.Sc. thesis, Carleton University, Canada, 1966. [3] Conway HD, Chow L, Morgan GW. Analysis of deep beams. J Appl Mech Trans ASME 1951;18(2):163–72. [4] Chow L, Conway HD, Winter G. Stresses in deep beams. Trans ASCE, 1952, Paper No. 2557. [5] Chapel RE, Smith HW. Finite-difference solutions for plane stresses. AIAA J 1968;6(6):1156–7. [6] Conway HD, Ithaca NY. Some problems of orthotropic plane stress. J Appl Mech Trans ASME 1953:72–6. Paper No. 52-A-4. [7] Ahmed SR, Khan MR, Islam KMS, Uddin MW. Analysis of stresses in deep beams using displacement potential function. J Inst Engrs (India) 1996;77:141–7. [8] Idris ABM, Ahmed SR, Uddin MW. Analytical solution of a 2-D elastic problem with mixed boundary conditions. J Inst Engrs (Singapore) 1996;36(6):11–7. [9] Ahmed SR, Idris ABM, Uddin MW. Numerical solution of both ends fixed deep beams. Comp Struct 1996;61(1): 21–9. [10] Ahmed SR, Khan MR, Islam KMS, Uddin MW. Investigation of stresses at the fixed end of deep cantilever beams. Comp Struct 1998;69:329–38. [11] Durelli AJ, Ranganayakamma B. Parametric solution of stresses in beams. J Eng Mech 1989;115(2):401–15. [12] 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. [13] Smart J. On the determination of boundary stresses in finite elements. J Strain Anal 1987;22(2):87–96. [14] Dow JO, Jones MS, Harwood SA. A new approach to boundary modeling for finite difference applications in solid mechanics. Int J Numer Meth Eng 1990;30:99–113. [15] Jones RM. Mechanics of composite materials. McGrawHill Book Company; 1975. [16] Akanda MAS, Ahmed SR, Khan MR, Uddin MW. A finite-difference scheme for mixed boundary-value problems of arbitrary-shaped elastic bodies. Adv Eng Softw 2000;31(3):173–84. [17] Akanda MAS, Ahmed SR, Uddin MW. Stress analysis of gear teeth using displacement potential function and finitedifferences. Int J Numer Meth Eng 2001;53(7):1629–49. [18] Suzuki S. Stress analysis of short beams. AIAA J 1986;24:1396–8. [19] Hardy SJ, Pipelzadeh MK. Static analysis of short beams. J Strain Anal 1991;26(1):15–29.