JOURNAL OF THEORETICAL
AND APPLIED MECHANICS
42, 4, pp. 739-753, Warsaw 2004
ON THE STRESS INTENSITY FACTORS FOR TRANSIENT
THERMAL LOADING IN AN ORTHOTROPIC THIN PLATE
WITH A CRACK
Bogdan Rogowski
Chair of Mechanics of Materials, Technical University of Łódź
e-mail: brogowsk@p.lodz.pl
This paper is concernedwith anorthotropic thinplate containinga crack
perpendicular to its surfaces. It is assumed that the transient thermal
stress is set up by the application of a heat flux as a function of time
and position along the crack edge and the heat flow by convection from
the plate surfaces. The exact analytical solutions for the stress intensity
factor and crack-opening displacement are derived. Numerical examples
show, among others, a dependence of the stress intensity factor on the
thermal and elastic constants of the orthotropic material.
Key words: orthotropic plate, crack problem, transient thermal loading,
stress intensity factor
1. Introduction
The study of thermoelastic problemshas always been an important branch
in solidmechanics (seeNowacki, 1986;Nowiński, 1978). Inparticular, the ther-
moelastic fracture problems subjected to various types of thermal boundary
conditions have been discussed extensively in the literature. Most of rese-
arch works discuss the steady-state crack problems and axisymmetric cases
for which the Hankel transform technique and the theory of dual integral
equations were usually employed (Sneddon, 1966). Recently, a report on a
penny-shaped or external crack subjected to temperature and heat flux, arbi-
trarily acting in a transversely isotropicmedium,was presented by the author
(Rogowski, 2003). The corresponding fundamental solution canplay an impor-
tant role in the boundary element method of thermoelastic fracture analysis.
Some metallic materials, such as zinc, magnesium, cadmium are transversely
740 B.Rogowski
isotropic (Hearmon, 1961). Many fibrous composites may also be modeled as
transversely isotropicmaterials (Christensen, 1979). Therehave beenmany re-
ports on crack analysis in transversely isotropic and orthotropic thermoelastic
materials. Among the studies, Tsai (1983a,b) calculated the stress intensity
factors of a penny shaped crack in a transversely isotropic material due to a
thermal loading, while Rogowski (2001a,b) presented analysis of a crack sys-
tem in transversely isotropic materials. Many of research works discuss the
two-dimensional thermal crack problem in the literature. Sumi (1981, 1982),
Aköz andTauchert (1972), Atkinson andClement (1977), Ghosh (1977), Cle-
ments and Tauchert (1979), Clements (1983), Tsai (1983a,b) and Rogowski
(1982) solved various problems in anisotropic thermoelastic solids. Gladwell et
al. (1983) considered the radiation boundary conditions. But, perhapsbecause
of mathematical complexity, the three-dimensional crack problem of an ani-
sotropic medium under transient thermal loading have not yet received much
attention. Among the studies, Koizumi and Niwa (1977) performed the ana-
lysis of an edge crack in a semi-infinite plate under transient thermal loading.
Noda and Matsunaga (1986) investigated the transient crack problem in an
infinite medium, while Ishida (1987) calculated the stress intensity factor for
a transient thermal loading in a transversely isotropic material. Ting and Ja-
cobs (1979) solved the problem for transient thermal stress in a cracked solid.
Many problems of thermoelasticity were solved in a book by Podstrigach and
Kolyano (1972).
This paper considers the transient thermal problem of a crack in an or-
thotropic thin plate. The method of solution involves the use of Fourier and
Laplace’s transforms and displacement potentials to reduce the mixed boun-
daryvalue problemto apair of dual integral equations.The solution is given in
an exact analytical form. The stress intensity factor of mode I and the crack-
opening displacement for a heat flux arbitrarily acting on the crack surface,
are determined. The numerical results are shown graphically to demonstrate
the influence of thermal andmechanical anisotropic parameters.
2. Analysis
2.1. Temperature field
Consider an orthotropic thin plate of thickness 2h containing a crack. Fi-
gure 1a shows the geometry of the problem where the position of the point
is defined by Cartesian co-ordiantes (x,y,z). In this co-ordinate system, the
On the stress intensity factors... 741
crack occupies the region y = 0, |x| ¬ a, |z| ¬ h. We shall suppose that the
crack is opened out by the heat flux depending on time and position applied
to its surfaces. Referring to the semi-infinite region y 0, the boundary con-
ditions in the problem can be assumed as shown in Fig.1b, since the thermal
and mechanical conditions on y = 0+ are identical with those on y = 0−.
Additionally, for a thin plate the unknown temperature distribution T(x,y,t)
is assumed to be constant over the thickness, giving the heat exchange by
convection on both surfaces of the plate, which equals −2γT , where T is the
temperature change and γ is theheat transfer coefficient on theplane surfaces.
Fig. 1. Geometry and co-ordinate system (a) and boundary conditions (b)
The equation heat conduction governing an unsteady-state temperature
field in an orthotropic thin plate with heat dissipation at both plane surfaces
is (Nowacki, 1986)
λ11
∂2T
∂x2
+λ22
∂2T
∂y2
−
γ
h
T = cρ
∂T
∂t
(2.1)
where c is the specific heat, ρ is the mass density and λ11 and λ22 are the
thermal conductivities in the x- and y-directions, respectively.
The initial and boundary conditions for the temperature field are
T =0 at t =0
λ22
∂T
∂y
= q0g(t)f(x)H(a−|x|) on y =0
(2.2)
where q0 is the heat flux per unit area and unit time, and H(·) denotes
Heaviside’s step function. The problem is symmetric with respect to the plane
y = 0, the temperature T(x,y,t) is an even function of y and differentiable
with respect to y at y = 0; in consequence, the heat flux is equal to zero
742 B.Rogowski
for y = 0, |x| > a. Applying Laplace’s transform to time and Fourier cosine
transform to the variable x, and using the convolution theorem for inverse
Laplace’s transform, the solution to (2.1) which satisfies (2.2) and (2.3) may
be expressed by
T=
t
∫
0
g(t− τ)
[
−
4
π2
q0χλ
2
λ22
∞
∫
0
f(s)cos(sx)ds
∞
∫
0
cos(py)e−χ(m
2+s2+λ2p2)τdp
]
dτ =
(2.3)
=
∞
∫
0
[
∞
∫
0
θ(s,p,t)cos(py) dp
]
cos(sx) ds
where
θ(s,p,t)=−
4
π2
q0χλ
2
λ22
f(s)
t
∫
0
g(t− τ)e−χ(m2+s2+λ2p2)τ dτ
f(s)=
a
∫
0
f(x)cos(sx) dx (2.4)
m2 =
γ
λ11h
χ =
λ11
cρ
λ2 =
λ22
λ11
From (2.4)2 it follows that only the symmetric problemwith respect to the y
axis is considered, since it is assumed that f(x) is an even function. For the
general case of the function f(x), its oddpartwill be associatedwithFourier’s
sine transform and the solution can be obtained in a similarmanner; formally
by replacement of cos(sx) with sin(sx) functions in Eq. (2.4)2.
2.2. Thermal stress and displacement
We consider the stress and displacement field. The stress-strain equations
for an orthotropic medium under a plane stress state are
σxx = c11exx+ c12eyy −β1T
σyy = c12exx+ c22eyy −β2T (2.5)
σxy =2Gexy
where eij are the strain components, σij are the stress components, cij are the
moduli of elasticity of thematerial, G is the shearmodulus, β1 = c11α1+c12α2,
On the stress intensity factors... 743
β2 = c12α1 + c22α2 and α1, α2 are the thermal expansion coefficients along
the x- and y-directions, respectively. The strain components are
exx =
∂ux
∂x
eyy =
∂uy
∂y
exy =
1
2
(∂ux
∂y
+
∂uy
∂x
)
(2.6)
where ux and uy are the displacement components along the axis. The equ-
ations of equilibrium for the plane stress in the absence of the body forces
are
∂σxx
∂x
+
∂σxy
∂y
=0
∂σxy
∂x
+
∂σyy
∂y
=0 (2.7)
From (2.5), (2.6) and (2.7), it follows
c11
∂2ux
∂x2
+G
∂2ux
∂y2
+(c12+G)
∂2uy
∂x∂y
=β1
∂T
∂x
(2.8)
G
∂2uy
∂x2
+ c22
∂2uy
∂y2
+(c12+G)
∂2ux
∂x∂y
= β2
∂T
∂y
The general solution to equilibriumequations (2.8)maybe obtained as the
superposition of two fields. The first corresponds to the solution to homoge-
neous equation (2.8), for which (Rogowski, 1975)
ux =
∂
∂x
(kϕ1+ϕ2) uy =
∂
∂y
(ϕ1+kϕ2)
σxx =−G(k+1)
∂2
∂y2
(ϕ1+ϕ2) σyy =−G(k+1)
∂2
∂x2
(ϕ1+ϕ2)
σxy = G(k+1)
∂2
∂x∂y
(ϕ1+ϕ2)
∂2ϕi
∂x2
+
1
s2i
∂2ϕi
∂y2
=0 (i =1,2)
(2.9)
where (i =1,2)
Gc22s
4
i − (c11c22− c212−2c12G)s2i +Gc11 =0 k =
c22s
2
1−G
c12+G
(2.10)
The second may be obtained in terms of the thermoelastic displacement po-
tential function ψ(x,y,t), defined as follows
ux =
∂ψ
∂x
uy = l
∂ψ
∂y
(2.11)
744 B.Rogowski
Equations (2.8) are satisfied if
c11
∂2ψ
∂x2
+G
∂2ψ
∂y2
+ l(c12+G)
∂2ψ
∂y2
= β1T
(2.12)
Gl
∂2ψ
∂x2
+ c22l
∂2ψ
∂y2
+(c12+G)
∂2ψ
∂x2
= β2T
A suitable expression for ψ defined by (2.12) for temperature distribution in
(2.3) is in the form
ψ =
∞
∫
0
[
∞
∫
0
C(s,p,t)cos(py) dp
]
cos(sx) ds (2.13)
This satisfies both equations (2.12) providing
C(s,p,t)[c11s
2+Gp2+ lp2(c12+G)] =−β1θ(s,p,t)
(2.14)
C(s,p,t)[l(c22p
2+Gs2)+s2(c12+G)] =−β2θ(s,p,t)
i.e.
l(s,p)=
β1s
2(c12+G)−β2(c11s2+p2G)
β2p
2(c12+G)−β1(c22p2+s2G)
(2.15)
C(s,p,t)= θ(s,p,t)
β2p
2(c12+G)−β1(c22p2+s2G)
(c11s2+Gp2)(c22p2+Gs2)− (c12+G)2p2s2
Appropiate solutions to Eqs (2.9)6 are
ϕ1(x,y)=−
s2
G(k+1)(s1−s2)
∞
∫
0
s−1A(s)e−s1sy cos(sx) ds
(2.16)
ϕ2(x,y)=
s1
G(k +1)(s1−s2)
∞
∫
0
s−1B(s)e−s2sy cos(sx) ds
On the stress intensity factors... 745
Using the above obtained potentials, we find
ux(x,y,t) =
1
G(k+1)(s1−s2)
∞
∫
0
[ks2A(s)e
−s1sy −s1B(s)e−s2sy]sin(sx)ds−
−
∞
∫
0
[
∞
∫
0
sC(s,p,t)cos(py) dp
]
sin(sx) ds
(2.17)
uy(x,y,t) =
s1s2
G(k+1)(s1−s2)
∞
∫
0
[A(s)e−s1sy −kB(s)e−s2sy]cos(sx) ds−
−
∞
∫
0
[
∞
∫
0
pC(s,p,t)l(s,p)sin(py) dp
]
cos(sx) ds
σxx(x,y,t) =
s1s2
s1−s2
∞
∫
0
s[s1A(s)e
−s1sy−s2B(s)e−s2sy]cos(sx) ds+
+ G
∞
∫
0
{
∞
∫
0
p2C(s,p,t)[l(p,s)+1]cos(py) dp
}
cos(sx) ds
σyy(x,y,t) = −
1
s1−s2
∞
∫
0
s[s2A(s)e
−s1sy−s1B(s)e−s2sy]cos(sx) ds+
(2.18)
+ G
∞
∫
0
[
∞
∫
0
s2C(s,p,t)[l(p,s)+1]cos(py) dp
]
cos(sx) ds
σxy(x,y,t) = −
s1s2
s1−s2
∞
∫
0
s[A(s)e−s1sy−B(s)e−s2sy]sin(sx) ds+
+ G
∞
∫
0
[
∞
∫
0
psC(s,p,t)[l(p,s)+1]sin(py) dp
]
sin(sx) ds
Themechanical boundary conditions on the plane y =0 are
σxy =0 (2.19)
σyy =0 on |x| < a uy =0 on |x| a (2.20)
Applying (2.18)3 to boundary condition (2.19), we obtain A = B.
746 B.Rogowski
Substituting (2.17)2 and (2.18)2 into boundary conditions (2.20) andusing
A = B, we obtain the following dual integral equations for A(s)
∞
∫
0
sA(s)cos(sx) ds =
∞
∫
0
s2F(s,t)cos(sx) ds on |x| < a
−
1
GC
∞
∫
0
A(s)cos(sx) ds =0 on |x| a
(2.21)
where
F(s,t)=−G
∞
∫
0
C(s,p,t)[l(p,s)+1] dp
(2.22)
C =(k+1)(k−1)−1(s−12 −s
−1
1 )
Equation (2.21)1 may be replaced by the following equation
∞
∫
0
A(s)sin(sx) ds =
∞
∫
0
sF(s,t)sin(sx) ds (2.23)
We introduce an integral representation of the function A(s)
A(s)=
a
∫
0
x′h(x′)J0(sx
′) dx′ (2.24)
where J0(sx
′) is theBessel function of the first kind and zero order, and h(x′)
is a new unknown function. This representation satisfies equation (2.21)2 and
converts equation (2.23) to the Abel integral equation for h(x′)
x
∫
0
x′h(x′)√
x2−x′2
dx′ =
∞
∫
0
sF(s,t)sin(sx) ds |x| < a (2.25)
The solution to this equation is
h(x′)=
∞
∫
0
s2F(s,t)J0(sx
′) ds (2.26)
where the following integral were employed
2
π
1
x′
d
dx′
x′
∫
0
xsin(sx)√
x′2−x2
dx = sJ0(sx
′) (2.27)
On the stress intensity factors... 747
Thus, the solution to the dual integral equations of form (2.21) is
A(s)=
a
∫
0
x′J0(sx
′)
[
∞
∫
0
q2J0(qx
′)F(q,t) dq
]
dx′ (2.28)
The above formula is exactly the same as that obtained by Sneddon (1966,
p.98) for dual integral equations of type (2.21). Therefore, we obtain the com-
plete solution to the problem by substituting (2.28) and (2.22)1 into (2.17)
and (2.18).
The singular stress σyy(x,0) is obtained as follows
σyy(x,0)=
xh(a)√
x2−a2
=
x√
x2−a2
∞
∫
0
s2F(s,t)J0(sa) ds as x → a+
(2.29)
Therefore, the stress intensity factor KI of mode I is defined as
KI = lim
x→a+
√
2π(x−a)(σyy)y=0 =
√
πa
∞
∫
0
s2F(s,t)J0(sa) ds (2.30)
Note that
GC(s,p,t)[l(p,s)+1]=−E1
α1p
2+α2p
2
δ2s4+2µs2p2+p4
θ(s,p,t)=
(2.31)
=−E1
( α1+ c0
s21s
2+p2
− c0
s22s
2+p2
)
θ(s,p,t)
where s21 and s
2
2 are the roots of algebraic equation (2.10)1, which may be
written in an equivalent form
s4i −2µs2i + δ2 =0 µ =
E1
2G
−ν21 δ2 =
E1
E2
(2.32)
and
c0 =
α1s
2
2−α2
s21−s22
=
α1(µ−
√
µ2− δ2)−α2
2
√
µ2−δ2
(2.33)
Here E1 and E2 are Young’s moduli in the x- and y-directions, respectively,
and ν21 is Poisson’s ratio.
748 B.Rogowski
The stress intensity factor KI is calculated from (2.30), and we get the
formula
KI =−
4
√
a
π
√
π
q0χλ
2E1
λ22
∞
∫
0
s2f(s)J0(as)
∞
∫
0
( α1+c0
s21s
2+p2
−
c0
s22s
2+p2
)
·
(2.34)
·
t
∫
0
g(t− τ)e−χ(m2+s2+λ2p2)τ dτ ds dp
For an isotropic material we have s1 = s2 = 1, α1 = α2 = α, λ = 1, and
stress intensity factor (2.34) assume the form
(KI)iso =−
4
√
a
π
√
π
q0χEα
λiso
∞
∫
0
s2f(s)J0(as)
∞
∫
0
1
s2+p2
·
(2.35)
·
t
∫
0
g(t− τ)e−χ(m2+s2+p2)τ dτ ds dp
The displacement uy(x,0) is obtained in the form
uy(x,0)=
4
√
2
π2
q0χλ
2δ
√
µ+ δ
λ22
a
∫
x
x′ dx′√
x′2−x2
∞
∫
0
s2f(s)J0(sx
′) ·
(2.36)
·
∞
∫
0
( α1+ c0
s21s
2+p2
−
c0
s22s
2+p2
)
t
∫
0
g(t− τ)e−χ(m2+s2+λ2p2)τ dτ ds dp
The displacement ux(x,0) is given by the formula
ux(x,0)=−
4
π2
q0χλ
2(δ−ν21)
λ22(s
2
2+ν21)
∞
∫
0
sf(s)sin(sx) ·
(2.37)
·
∞
∫
0
( c1+ c2
s21s
2+p2
−
c2
s22s
2+p2
)
t
∫
0
g(t− τ)e−χ(m2+s2+λ2p2)τ dτ ds dp
On the stress intensity factors... 749
where
a1 =2−
G
E1
ν12 b1 =
G
E1
c1 = a1α1− b1α2
a2 =1−
G
E2
ν12 b2 =
G
E2
c2 =
c1s
2
2− b2α1−a2α2
s21−s22
ν12
E2
=
ν21
E1
(2.38)
3. Numerical results
In calculating the temperature and the stress intensity factors, the follo-
wing dimensionless quantities are introduced
ξ =
x
a
η =
y
a
t′ = χ
t
a2
M2 = a2m2 =
γa2
λ11h
λ2 =
λ22
λ11
α =
α2
α1
δ2 = E =
E1
E2
µ =
E1
2G
−ν21 T =
Tλ22
q0a
KI =
KIλ22
α1E1aq0
√
a
σyy =
σyyλ22
α1E1aq0
Numerical calculations were carried out for two types of the heat supply q
Case 1: q = q0g(t)f(x) = q0
Case 2: q = q0g(t)f(x) = q0e
−t′
where t′ is the Fourier number.
Figure 2 shows the temperature at η = 0 for various values of λ2. The
temperature at η =0 increases with the ratio of thermal conductivity. Figures
3a,b show the effects of λ2 and α on the normal stress σyy at η =0 for case 1
of thermal loading.
Figures 4a-d show the effects of anisotropies of the material constants on
the stress intensity factor for case 1 and case 2. It is assumed that only one
of thematerial constants λ2, α, E, µ indicates various anisotropies, while the
other constants are kept equal to those of isotropic conditions.
Figures 2-4 show that the anisotropy effects of the material constants λ2,
α, E and µ on the stress intensity factor are large. In the figureswe can notice
750 B.Rogowski
Fig. 2. Variation of temperature at η =0with ξ for various λ2
Fig. 3. Variation of normal stress σyy at η =0with ξ for various values of λ
2 (a)
and α (b)
that the stress intensity factor increaseswith the thermal conductivity,Young’s
modulus and thermal expansion coefficient in the direction perpendicular to
the crack plane.
On the stress intensity factors... 751
Fig. 4. Variation of stress intensity factor with t′ for various values of λ2 (a), α (b),
E (c) and µ (d)
References
1. AközA.Y., TauchertT.R., 1972,Thermal stresses in an orthotropic elastic
semi-space, J. Appl. Mech. ASME, 39, 87-90
2. Atkinson C., Clements D.L., 1977, On some crack problem in anisotropic
thermoelasticity, Int. J. Solids Struct., 13, 855-864
3. Christensen R.M., 1979, Mechanics of Composite Materials, Wiley, New
York
4. ClementsD.L., 1983,A thermoelastic problem for a crack between dissimilar
anisotropic media, Int. J. Solids Struct., 19, 121-130
5. Clements D.L., Tauchert T.R., 1979, A thermoelastic crack problem for
an anisotropic slab, J. Aust. Math. Soc., Ser. B, 21, 243-255
6.
6. Ghosh S., 1977, A note on two-dimensional thermoelastic crack problem for
an aelotropic solid, Indian J. Mech. Math., 15, 6-15
752 B.Rogowski
7. GladwellG.M.L.,Barber J.R.,OlesiakZ., 1983,Thermalproblemswith
radiation boundary conditions,Q. Jl. Mech. Appl. Math., 36, 387-401
8. Hearmon R.F.S., 1961, An introduction to applied anisotropic elasticity,
Oxford Univ. Press
9. Ishida R., 1987, Crack problem in a transversely isotropic medium with a
penny-shapedcrackunder transient thermal loading,Zeit.Angew.Math.Mech.,
ZAMM, 67, 2, 93-99
10. KoizumiT.,NiwaH., 1977,Transient thermal stressproblem ina semi-infinite
platewith anedge crackperpendicular to the edge surface,Trans. of the Society
of Japan Mech. Eng., 43, 422-447
11. Noda N., Matsunaga Y., 1986, Transient thermoelastic problem in an infi-
nite body containing a penny-shaped crackdue to time and position dependent
temperature condition, Zeit. Angew. Math. Mech., ZAMM, 66, 6, 233-239
12. Nowacki W., 1986,Thermoelasticity, PWN,Warszawa
13. Nowiński J.L., 1978, Theory of Thermoelasticity with Applications, Sijthoff
and Noordhoff Int. Publ., The Netherlands
14. Podstrigach Ya.L., Kolyano Y.M., 1972,Transient Temperature Stresses
in Thin Plates, (in Russian), NaukowaDumka, Kiev
15. Rogowski B., 1975, Displacement functions for transversely isotropic media
(in Polish),Mechanika Teoretyczna i Stosowana, 13, 1, 69-83
16. Rogowski B., 1982, Thermal stresses in a transversely isotropic layer conta-
ining an annular crack. Tensile – and shear type crack, J. Theor. Appl. Mech.,
22, 3-4, 473-492
17. RogowskiB., 2001a,Heat conductionproblem inatransversely isotropicbody
with a circular crack system, Mathematical Methods and Physicomechanical
Fields, 44, 3, 128-134
18. Rogowski B., 2001b,Thermal stresses problem in a transversely isotropic bo-
dywith a circular crack system,Mathematical Methods and Physicomechanical
Fields, 44, 4, 85-92
19. Rogowski B., 2003, Fundamental solutions related to thermal stress intensity
factors of modes I and II – the axially symmetric problem, J. Theor. Appl.
Mech., 41, 2, 241-269
20. Sneddon I.W., 1966, Mixed Boundary Value Problems in Potential Theory,
Amsterdam: North-Holland
21. Sumi N., 1981, Thermal stresses in an orthotropic rectangular plate with a
rigid ribbonlike inclusion,Nucl. Eng. Des., 66, 405-411
22. SumiN., 1982,The thermoelastic problem for an orthotropic rectangular plate
with an oblique crack,Theor. Appl. Mech., 31, 135-148
On the stress intensity factors... 753
23. Ting V.C., Jacobs H.R., 1979, Stress intensity factors for transient thermal
loadings of a semi-infinite medium, J. Thermal Stresses, 2, 1, 1-13
24. Tsai Y.M., 1983a, Thermal stress in a transversely isotropic medium conta-
ining a penny-shaped crack,Trans. ASME J. Appl. Mech., 50, 24-28
25. Tsai Y.M., 1983b, Transversely isotropic thermoelastic problem of uniform
heat flow disturbed by a penny-shaped crack. J. Thermal Stresses, 6, 379-389
O współczynnikach intensywności naprężenia dla nieustalonego
termicznego obciążenia w ortotropowej cienkiej płycie ze szczeliną
Streszczenie
Rozpatrzono zagadnienie ortotropowej cienkiej płyty zawierającej szczelinę pro-
stopadłą do jej brzegów. Założono, że nieustalone naprężenia termiczne powstają
wwyniku przepływuprzez powierzchnie szczeliny strumienia ciepła będącego funkcją
czasu,miejsca i konwekcyjnegoprzepływuciepłaprzez powierzchniepłyty. Znaleziono
ścisłe, analityczne rozwiązanieokreślającewspółczynnik intensywnościnaprężeń i roz-
warcie szczeliny. Przykłady numeryczne pokazują zależności temperatury, naprężeń
i współczynnika intensywności naprężenia od parametrów geometrycznych i stałych
określającychwłasności termiczne i sprężyste materiału.
Manuscript received February 10, 2004; accepted for print April 24, 2004