Set of field equations for thick shell of revolution made of functionally graded materials in curvilinear coordinate system/ Lauko lygciu taikymas aprasant storus sukinio formos kevalus, pagamintus is funkciskai pagerintu medziagu kreivineje koordinaciu sistemoje.

Nejad, M. Zamani ; Rahimi, G.H. ; Ghannad, M. 等

1. Introduction

Nonhomogeneous material systems with gradual variation in properties are collectively referred to as functionally graded materials or FGMs. (FGMs) used initially as thermal barrier materials for aerospace structural applications and fusion reactors are now developed for the general use as structural components in high temperature environments and being strongly considered as a potential structural material candidate for the design of high speed aerospace vehicles. Further, these materials are inhomogeneous, in the sense that the material properties vary smoothly and continuously in one or more directions, and obtained by changing the volume fraction of the constituent materials [1]. In the past few years, some researchers began to pay attention to contact problem of functionally graded materials.

Thick shells of revolution of functionally graded materials have been of great interest in many engineering applications. A body of revolution is obtained by rotating a plane area 360[degrees] about an axis in its plane. These bodies may be either solid or hollow. If they are hollow, they may be called shells of revolution. Circular cylinders, spheres, cones, ellipsoids, paraboloids, and hyperboloids are just some of the shapes which bodies of revolution may take and are frequently used in various industries because of their great structural efficiency.

There are a number of publications on field equations, equations of motion, for thin or moderately thick shells of revolution. Flugge [2], Timoshenko and Woinowsky-Krieger [3] have provided an adequate description of the geometry and the coordinate system for thin shells of revolution. In these derivations, simple kinematic assumptions are made about the variation of the displacements through the thickness, the result being that a three-dimensional theory is reduced to a two-dimensional one by the middle surface displacements. The twodimensional models which already exist have a number of shortcomings. This is especially the case when we deal with thick and very thick shells. In order to analyze static or dynamic displacements and stresses, and free and forced vibrations for such shells, there must be a threedimensional theory. The equations of equilibrium or motion derived could be solved by the exact method or the approximate method. Using tensor analysis, which is extremely useful especially for three-dimensional problems in any curvilinear coordinate system, a number of researchers, including Fung [4], Green and Zerna [5] have developed the general theory of elasticity.

A lot of research has been carried out into general thick-walled shells [6, 7], thick-walled cylindrical shells [8-16], thick-walled spherical shells [17, 18], thick-walled conical shells [19-21], thick-walled elliptical [22], helicoidal [23] and toroidal shells [24]. The field equations derived for thick shells of revolution were considered for the homogenous condition. A distinguishing feature of the present paper is that for a thick shell of revolution with arbitrary curvature and variable thickness made of functionally graded materials, a set of field equations has been developed by tensor analysis in the curvilinear coordinate system. However, to the best of the researchers' knowledge, no extensive study has yet been carried out on threedimensional displacement equations for thick shell of revolution with of functionally graded materials.

2. Analysis

A surface of revolution is the result of the rotation of a plane curve around an axis in its plane. The resulting curve is technically called a meridian. The meridian of a thick shell of revolution with an arbitrary curvature and variable thickness z is explained by Eq. (1), or alternatively by Eq. (2)

R = H (x) (1)

R = R([phi]), x = x([phi]) (2)

According to Fig. 1, R represents the distance of one of its points from the axis of rotation x - axis , and x represents the distance of its points from the R - axis [phi] is the angle between the normal to the meridian curve and the axis of revolution. [R.sub.1] and [R.sub.2], shown in Fig. 1 are the two principal radii of curvature.

If Eq. (1) is used, these can be obtained using the following equations

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3)

It is also possible to calculate them using Eq. (2), as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (4)

Using Fig. 1, we can define the following equations

R = [R.sub.2] sin[phi] (5)

ds = [R.sub.1]d[phi] (6)

dR = dscos[phi] (7)

dx = dssin[phi] (8)

where ds is the line element of meridian. Based on Eqs. (5) to (8), we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

[FIGURE 1 OMITTED]

In Fig. 2, the cross-section of an arbitrary shell element of variable thickness z is shown. The coordinate used is curve linear coordinate system ([phi], [OMEGA], [THETA]), where [THETA] is the meridional coordinate; [OMEGA] the normal distance from the midsurface to an arbitrary point P; and [THETA] the circumferential angle. The two faces on the top and bottom, shown in Fig. 2, are flat and normal to the midsurface. The arbitrary point of P within the shell element could be determined by [phi], [OMEGA] and [THETA].

In Fig. 2, r shows the radial distance of the arbitrary point P , which depends on meridional angle and the normal distance [OMEGA], as

r([phi], [OMEGA]) = R([phi]) + [OMEGA]sin[phi] (10)

By substituting R([phi]) defined in Eq. (5) in Eq. (10), we have

r([phi], [OMEGA] = [[R.sub.2]([phi] + [OMEGA]]sin[phi] (11)

The curvilinear coordinates ([phi], [OMEGA], [THETA]), using tensor analysis, could be transformed to Cartesian coordinate ([x.sup.1], [x.sup.2], [x.sup.3])

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (12)

In this equation, [OMEGA] and [THETA] must be within the following ranges

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (13)

[FIGURE 2 OMITTED]

3. Field equations of elasticity

In order to develop field equations, base vectors, metric coefficients and Christoffel symbols are used in the curvilinear coordinates. If r is the position vector of an arbitrary point P , then prior to deformation, the covariant base vectors [[??].sub.i] are defined as follows

[[??].sub.i] = [[??].sub.j] (14)

where comma (,) represents partial differentiation with respect to [phi], [OMEGA], and [THETA].

In the Cartesian coordinate systems, the covariant base vectors [[??].sub.i] are given as follows.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (15)

where

[R.sub.[OMEGA]] = [R.sub.1]([phi]) + [OMEGA] (16)

For the unstrained body, the contravariant base vectors [[??].sub.i] are defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (17)

where the Kronecker delta [[delta].sub.ij] is defined as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (18)

and i and j could take any positive integer value.

The contravariant base vectors [[??].sup.1], [[??].sup.2] and [[??].sup.3] could be calculated using the following equations

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (19)

where

J = [[??].sub.1]([[??].sub.2] x [[??].sub.3]) (20)

By using the Eqs. (15), (19) and (20) we can obtain the contravariant base vectors.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (21)

Dot products of the base vector give a set of symmetric numbers which are called metric coefficients. The covariant and contravariant metric coefficients ([g.sub.ij] and [g.sup.ij]) are calculated using the following equations

[g.sub.ij] = [[??].sub.i][[??].sub.j] (22)

[g.sub.ij] = [[??].sup.i][[??].sup.j] (23)

If the coordinate system is orthogonal, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (24)

The nonvanishing covariant and contravariant metric coefficients, which take the values other than zero, are given as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (25)

By definition, the Christoffel symbols of the second kind, which have symmetric components, [[SIGMA].sup.k.sub.ij] = [[SIGMA].sup.k.sub.ji], are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (26)

where [g.sup.j] represents spatial curvilinear coordinates in which j = 1, 2, 3 correspond to [phi], [OMEGA], and [theta]. The non-vanishing Christoffel symbols of the second kind, are defined by Eqs. (15), (21) and (26) as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (27)

All other components of the Christoffel symbols of the second kind become zero.

Physical components related to stress tensor ([[sigma].sup.*.sub.if]) in general coordinate system are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (28)

where ^ = aJi are the components of stress tensor. In addition

cj = gkck (29)

where [[sigma].sup.k.sub.j] [not equal to] [[sigma].sup.j.sub.k] are the mixed components of stress tensor.

Using the Eqs. (25), (28) and (29) mixed components of stress tensor are obtained as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (30)

Components of displacement [u.sub.i] in the general coordinate system are

[u.sub.i] = [square root of ([g.sub.ii][u.sup.*.sub.i])] (31)

where [u.sup.*.sub.i] are physical components of displacement. Therefore

[u.sub.1] = [R.sub.[OMEGA]][u.sub.[psi]]

[u.sub.2] = [u.sub.[OMEGA]]

[u.sub.3] = r[u.sub.[theta]] (32)

Physical components of strain tensor [[epsilon].sup.*.sub.if] in the general coordinate system are

[[epsilon].sup.*.sub.if] = [square root of ([g.sup.ii])][square root of ([g.sup.jj])] [[epsilon].sub.ij] (33)

where [[epsilon].sub.ij] are the components of strain tensor and are defined from the following equations

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (34)

Given the Eqs. (25), (27), (32) and (34), physical components of strain tensor are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (35)

The physical stress-strain ([[sigma].sup.*.sub.ij] - [[epsilon].sup.*.sub.ij]) relations in tensor form for a linearly elastic material are given by

[[sigma].sup.*sub.ij] = [lambda][[delta].sub.ij][[epsilon].sup.*.sub.kk] + 2 [mu][[epsilon].sup.*.sub.ij] (36)

where [lambda] and [mu] are the Lame constants. The constants are related to Young's modulus E and Poisson's ratio u for an isotropic solid by

[lambda] = E[upsilon]/(1+[upsilon])(1 - 2[upsilon]) (37)

[mu] = E/2(1 + [upsilon]) (38)

By using Eq. (36), the stress-strain relations in terms of the physical components are obtained as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (39)

The equations of translational motion could be written compactly in tensor form

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (40)

where the vertical bar (|) represents covariant differentiation. Also, [f.sub.i] are components of the body force vector [??], per unit volume; [rho] is the mass density per unit volume; [[??].sub.j] are the covariant components of the acceleration of the volume in the deformed body.

The covariant derivatives of a tensor of order two are also tensors and are represented as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (41)

where [A.sup.i.sub.j], are the mixed components of a typical tensor of order two.

By substituting Eqs. (27) and (30) into Eq. (40), the equations of motion in terms of the physical components could be derived as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (42)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (43)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (44)

In deriving the Eqs. (42) to (44), was used the following

[r.sub.,[phi] = [R.sub.[OMEGA]] cos[phi] (45)

Combining Eqs. (35), (39) and (42) to (44), the displacement equations of motion are obtained

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (46)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (47)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (48)

where the Lame constant derivatives are related to Young's modulus E and Poisson's ratio [upsilon], as shown below.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (49)

4. Solution for axisymmetric circular cylindrical state

In this section, a special case of thick shell of revolution, i.e. thick-walled FGM cylinder ([partial derivative]/[partial derivative][OMEGA] = [partial derivative]/[partial derivative]r, [phi] = [pi]/2, [R.sub.[OMEGA]] [right arrow] [infinity]) with an inner radius [r.sub.i], and thickness t, subjected to an internal pressures P, that is axisymmetric, (Fig. 3), will be considered.

[FIGURE 3 OMITTED]

For axisymmetric case, in plane elasticity theory (PET), we have: [u.sub.[theta]] = 0, [partial derivative]/[partial derivative][theta] = 0, [u.sub.[phi],i] = 0 and the Eqs. (43) yield the following equation:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (50)

where [[sigma].sub.rr] and [[sigma].sub.[theta][theta]] are the radial and circumferential stress components, respectively. Equation (50) is equilibrium equation in the radial direction, neglecting the body force components.

By using Eqs. (35), the radial strain [[epsilon].sub.rr] and circumferential strain [[epsilon].sub.[theta][theta]] are related to the radial displacement [u.sub.r] by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (51)

By using Eqs. (39) and (49), the stress-strain relations for nonhomogenous and isotropic materials are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (52)

where A and B are related to Poisson's ratio u as Plane strain condition:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (53)

Plane stress condition:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (54)

Using Eqs. (50) to (54), the Navier equation in terms of the radial displacement is

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (55)

The material properties are assumed to be radially dependent. Given that the radial coordinate r is normalized as [bar.r] = r/[r.sub.i], the module of elasticity through the wall thickness is assumed to vary as follows.

E = [E.sub.i] [[bar.r].sup.[beta]] (56)

here [E.sub.i] is the module of elasticity at the inner surface r = [r.sub.i], and [beta] is the inhomogeneity constants which are determined empirically.

Since the analysis was carried out for thick wall cylindrical pressure vessel of isotropic FGM, and given that the variation of Poisson's ratio, [upsilon], for engineering materials is small, the Poisson's ratio is assumed as constant.

By substituting Eq. (56) into Eq. (55), the Navier equation would be

[r.sup.2][u.sub.r,rr] + ([beta] + 1)r[u.sub.r,r] + ([[upsilon].sup.*][beta] - 1) [u.sub.r] = 0 (57)

where [[upsilon].sup.*] = B/A

Equation (57) is the nonhomogeneous Euler-Caushy equation whose complete solution is

[u.sub.r] = [C.sub.1][r.sup.m1] + [C.sub.2][r.sup.m2] (58)

where [m.sub.1] and [m.sub.2] are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (59)

By substituting Eq. (58) into Eqs. (51) and (52), the radial stress is obtained as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (60)

To determine the constants [C.sub.1] and [C.sub.2], consider the boundary conditions for stresses given by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (61)

Substituting the boundary conditions (61) into Eq. (60) and solving for [C.sub.1] and [C.sub.2], we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (62)

where k = 1 + t/[r.sub.i].

Hence, the radial stress, circumferential stress and radial displacement are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (63)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (64)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (65)

The analytical solution obtained in this section may be checked for one example.

Example. Consider a hollow functionally graded cylinder of inner radius [r.sub.i] = 40 mm , and the thickness t = 20 mm. The modulus of elasticity at inner radius has the value of [E.sub.i] = 200 GPa . It is also assumed that the Poisson's ratio, u , has a constant value of 0.3 . The applied internal pressure is 80MPa . In addition, [beta] ranges from -2 to 2 . The range -2 [less than or equal to] [beta] [less than or equal to] 2 to be used in the present study covers all the values of coordinate exponent encountered in the references cited earlier.

Displacement and stress distributions depending on an inhomogeneity constant are compared with the solutions of the finite element method (FEM) and are presented in the form of graphs.

The radial displacement along the radius for the conditions of plane strain, and plane stress, is plotted in Figs. 4 and 5.

There is a decrease in the value of the radial displacement as [beta] increases. Besides, for similar values of [beta], the value of radial displacement is the highest for the plane stress condition and for the plane strain the lowest. Figs. 6 and 7 show the distribution radial and circumferential stresses in the radial direction. As [beta] increases, so does the magnitude of the radial stress. For [beta] > 1, the circumferential stress increases as the radius increases whereas for [beta] < 1 the circumferential stress along the radius decreases. Given that [beta] = 1, the circumferential stress remains nearly constant along the radius.

[FIGURE 4 OMITTED]

[FIGURE 5 OMITTED]

[FIGURE 6 OMITTED]

[FIGURE 7 OMITTED]

5. Conclusions

Assuming small strains, for thick shells of revolution, with arbitrary curvature, variable thickness, and of functionally graded materials, a set of partial differential equations in terms of displacement components was developed, which could be useful for analyzing the static and dynamic behavior. These consist of six stress-strain relations, six strain-displacement relations, and three equations of motion in terms of physical components. The equations are expressed in terms of coordinates tangent and normal to the shell middle surface. The relationships are combined to yield equations of motion in terms of orthogonal displacement components taken in the meridional, normal and circumferential directions. For various thick shells of revolution, approximate solutions could be obtained by using proper numerical methods. In addition, exact solutions of displacement equations of motion are possible for some shell configurations, such as constant thick-walled cylindrical and spherical shells. It is apparent that exact solutions are highly significant in simplified versions of real engineering problems.

Using the field equations derived in the present study, exact solutions for stresses in a functionally graded materials pressurized thick-walled hollow circular cylinder are obtained under generalized plane strain and plane stress assumptions, respectively. The material properties are assumed to vary nonlinearly in the radial direction, and the Poisson's ratio is assumed constant. Depending on an inhomogeneity constant, the displacements and stress distributions are compared with the solutions of the finite element method and good agreement are found. To show the effect of inhomogeneity on the stress distributions, different values of [beta] were considered. Results show that the inhomogeneous constant [beta] presented in the current study is a useful parameter from a design point of view in that it can be tailored to specific applications to control the stress distributions. Thus, by selecting a proper value of [beta], it is possible for engineers to design FGM pressurized thick hollow cylinder that can meet some special requirements. It is also possible to find an optimum value for the power law index [beta] such that the variation of stresses along the radial direction is minimized, yielding optimum use of material.

Received April 16, 2009

Accepted May 29, 2009

References

[1.] Suresh, S., Mortensen, A. Functionally graded metals and metal-ceramic composites-Part 2, thermomechanical behavior. -Int. Materials Reviews, 1997, v.42, No3, p.85-116.

[2.] Flugge, W. Stresses in Shells. -Berlin: Springer-Verlag, 1973.-18p.

[3.] Timoshenko, S., Woinowsky-Krieger, S. Theory of Plates and Shells. -New York: McGraw-Hill 1959. -445p.

[4.] Fung, Y.C. Foundations of Solid Mechanics. -Prentice-Hall, 1965.-87p.

[5.] Green, A.E, Zerna, W. Theoretical Elasticity. -Oxford, 1968.-148p.

[6.] Friedrich, C.M. An elastic stress analysis for thick shells of revolution. -Nuclear Engineering and Design, 1967, v.6, p.367-375.

[7.] Kang, J.H., Leissa, A.W. Three-dimensional field equations of motion, and energy functionals for thick shells of revolution with arbitrary curvature and variable thickness. -J. of Applied Mechanics-Transactions of the ASME, 2001, v.68, No6, p.953-954.

[8.] Dorosevas, V., Volkovas, V. Analysis the interaction of two cylindrical surfaces under shock impact load. -Mechanika. -Kaunas: Technologija, 2007, Nr.6(68), p.49-52.

[9.] Ghannad, M., Rahimi, G.H., Khadem, S.E. General shear deformation solution of axisymmetric functionally graded thick cylindrical shells. -J. of Modares Technology and Engineering, 2009, Accepted. (in Iran).

[10.] Horgan, C.O., Chan, A.M. The pressurized hollow cylinder or disk problem for functionally graded isotropic linearly elastic materials. -J. of Elasticity, 1999, v.55, No1, p.43-59.

[11.] Partaukas, N., Bareisis, J. The stress state in two-layer hollow cylindrical bars. -Mechanika. -Kaunas: Technologija, 2009, Nr.1(75), p.5-12.

[12.] Rahimi, G.H., Zamani Nejad, M. Exact solutions for thermal stresses in a rotating thick-walled cylinder of functionally graded materials. -J. of Applied Sciences, 2008, v.8, No18, p.3267-3272.

[13.] Shao, Z.S., Ma, G.W. Thermo-mechanical stresses in functionally graded circular hollow cylinder with linearly increasing boundary temperature. -Composite Structures, 2008, v.83, No3, p.259-265.

[14.] Shi, Z.F., Zhang, T.T., Xiang, H.J. Exact solutions of heterogeneous elastic hollow cylinders. -Composite Structures, 2007, v.79, No1, p.140-147.

[15.] Zamani Nejad, M., Rahimi, G.H. Deformations and stresses in rotating FGM pressurized thick hollow cylinder under thermal load. -Scientific Research and Essays, 2009, v.4, No3, p.131-140.

[16.] Zamani Nejad, M., Rahimi, G.H. Elastic analysis of FGM rotating cylindrical pressure vessels. - Journal of the Chinese Institute of Engineers, Accepted.

[17.] Kang, J.H., Leissa, A.W. Three-dimensional vibrations of thick spherical shell segments with variable thickness. -Int. J. of Solids and Structures, 2006, v.43, No 9, p.2848-2851.

[18.] Poultangari, R., Jabbari, M., Eslami, M.R. Functionally graded hollow spheres under non-axisymmetric thermo-mechanical loads. -Int. J. of Pressure Vessels and Piping, 2008, v.85, No5, p.295-305.

[19.] Eipakchi, H.R., Khadem, S.E., Rahimi, G.H. Axisymmetric stress analysis of a thick conical shell with varying thickness under nonuniform internal pressure. -J. of Engineering Mechanics-ASCE, 2008, v.134, No8, p.601-610.

[20.] Panferov, I.V. Stresses in a transversely isotropic conical elastic pipe of constant thickness under a thermal load. -PMM J. of Applied Mathematics and Mechanics, 1992, v.56, No3, p.410-415.

[21.] Srinivasan, R.S., Hosur, V. Axisymmetric Vibration of thick conical shells. -J. of Sound and Vibration, 1989, v.135, No1, p.171-176.

[22.] Patel, B.R., Shukla, K.K., Nath, Y. Thermo-elastic stability behavior of laminated cross-ply elliptical shells. -Structural Engineering and Mechanics, 2005, v.19, No6, p.749-755.

[23.] McGee, O.G. The three-dimensional vibration analysis of a cantilevered skewed helicoidal thick shell. -J. of the Acoustical Society of America, 1993, v.93, No3, p.1431-1444.

[24.] Buchanan, G.R., Liu, Y.J. An analysis of the free vibration of thick-walled isotropic toroidal shells. -Int. J. of Mechanical Sciences, 2005, v.47, No2, p.277-292.

M. Zamani Nejad *, G.H. Rahimi **, M. Ghannad ***

* Mechanical Engineering Department, Tarbiat Modares University, Tehran, Iran, E-mail: mzamani@modares.ac.ir

** Mechanical Engineering Department, Tarbiat Modares University, Tehran, Iran, E-mail: gh.rahimi.s@gmail.com

*** Mechanical Engineering Faculty, Shahrood University of Technology, Shahrood, Iran, E-mail: ghannad.mehdi.@gmail.com