An accurate Ritz approach for analysis of cracked stiffened plates

Applied Mathematical Modelling 73 (2019) 598–614

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

An accurate Ritz approach for analysis of cracked stiffened plates B.H.S. Oliveira a, E. Lucena Neto b,∗, F.A.C. Monteiro b a b

EMBRAER S.A., 12227-901 São José dos Campos, São Paulo, Brazil Instituto Tecnológico de Aeronáutica, 12228-900 São José dos Campos, São Paulo, Brazil

a r t i c l e

i n f o

Article history: Received 10 September 2018 Revised 11 March 2019 Accepted 4 April 2019 Available online 16 April 2019 Keywords: Stiffened plate Crack Ritz method Hierarchic set Polynomial set Stress intensity factor

a b s t r a c t The present paper promotes the capabilities of the Ritz method in accurately predicting the stress intensity factor of cracked plates with either continuously or discretely attached stiffeners. The original plate domain is initially treated as an assembly of a small number of disjoint subdomains, chosen according to the crack and stiffener locations, in order to properly construct Ritz bases accounting for discontinuities. A complete set of hierarchic polynomials is adopted to locally approximate the displacement field within each subdomain, and then the displacement continuity is exactly imposed between subdomains whenever required. The accuracy of the present solutions are confirmed through comparison with published data, when available, and with converged results obtained from developed finite element models. The proposed procedure is further employed to investigate the effect that a propagating crack has on the stress intensity factor as the crack approaches and crosses a stiffener. © 2019 Elsevier Inc. All rights reserved.

1. Introduction Stiffened plates have a wide range of applications in aerospace, civil and marine structures [1,2]. They are used, for instance, in wing structures, box girders and ship hulls. The main feature of stiffened-plate constructions lies in their high resistance-to-weight ratio. Such thin-walled structures have cracks as one of the most common defects, whose presence gives rise to singularities in the stress field around the crack tip. From the standpoint of fracture mechanics [3], the knowledge of this local stress field is the quantity of primary interest. It is closely related to the structure residual strength and rate of crack propagation, and has intensity characterized by the stress intensity factor. There are specialized handbooks [4,5] devoted to collect as many stress intensity factors as possible in order to provide their values for various geometry and loading. Unfortunately, stress intensity factors can only be evaluated analytically for some restricted cases [6,7]. Even a small change on the stiffener spacing or fastener pitch may be enough to considerably alter the stress distribution pattern and, thus, the magnitude of the stress intensity factor. A general purpose evaluation of such a factor must be necessarily carried out by sophisticated numerical procedures [8–10]. The p-version of the finite element method, in which convergence is sought by progressively increasing the number of basis functions used in the element formulation, has been acclaimed as one of the most efficient methods to deal with this problem [11,12].

∗

Corresponding author. E-mail address: [email protected] (E. Lucena Neto).

https://doi.org/10.1016/j.apm.2019.04.014 0307-904X/© 2019 Elsevier Inc. All rights reserved.

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

599

Although the superiority of the finite element method is well-established for complex geometry, the use of this method for structural analysis of rectangular cracked stiffened plates may be computationally expensive compared with a dedicated Ritz scheme. However, special care should be exercised in a Ritz solution of this problem when setting up the Ritz basis because all the discontinuities introduced by stiffeners and cracks into the displacement field must be properly handled. Another difficulty to be faced lies in establishing integrals that can be promptly evaluated. To investigate the vibrations of a cracked unstiffened rectangular Kirchhoff plate, Yuan and Dickinson [13] simplify the construction of the Ritz basis by decomposing the original plate domain into several subdomains. A set of shape functions is chosen to locally approximate the displacement field within each subdomain, and then artificial springs are used to join each subdomain to the adjacent ones along their interconnecting boundaries. A similar approach is followed by Liew et al. [14], but now the degree of continuity required by the Ritz method for the displacement field along the interconnecting boundaries is imposed in an integral sense. Taking also advantage of the Ritz method with a domain decomposition, Barrette et al. [15] investigate the vibration of stiffened plates without cracks, while Kumar and Paik [16] predict the buckling of cracked unstiffened plates. The inclusion of stiffeners in the buckling analysis of cracked plates has been handled by Dang and Kapania [17]. To locally approximate the displacement field within each subdomain, the references [15–17] adopt the hierarchic trigonometric set proposed by Beslin and Nicolas [18]. More recently, the above Ritz solution scheme has also been used in the analysis of composite omega stiffened panels by Vescovini and Bisagni [19] and in the analysis of buckling and vibration of stiffened composite panels with debonding defect by Castro and Donadon [20]. The degree of displacement continuity between subdomains is exactly imposed in references [15–17], while this constraint is enforced in references [19,20] by means of a penalty-based approach. In this paper, the Ritz method with a domain decomposition is extended to create a flexible approach to accurately predict the stress intensity factor of cracked stiffened plates. Each skin subdomain is supposed to behave as a membrane and the stiffeners are taken as stretching bars either continuously or discretely attached to the skin. A complete set of hierarchic polynomials [21] is adopted to locally approximate the displacement field within each subdomain. The first and second elements of the set are the linear Lagrange polynomials, and the remaining ones stem from a single integration of Legendre polynomials. A polynomial set is chosen herein because of its superior convergence properties in capturing localized features [22]. The Ritz method is then used to discretize the potential energy of a subdomain. The displacement continuity is exactly imposed between subdomains, whenever required, and used to identify the discrete potential energy of the whole structure. The satisfaction of this continuity in an exact manner is an important issue to preserve some properties of the Ritz method, such as monotonic convergence of the strain energy from below even in the presence of singularities [23]. The accuracy of the present solutions are confirmed through comparison with published data, when available, and with converged results obtained from developed finite element models. In addition, the proposed procedure is employed to investigate the effect that a propagating crack has on the stress intensity factor as the crack approaches and crosses a stiffener. 2. Fundamentals Fig. 1a depicts a cracked stiffened plate subjected to distributed loads qx , qy and qxy per unit length applied to the edges, and concentrated loads Qx and Qy applied to the stiffener ends. The plate skin is supposed to behave as a membrane and the stiffeners are taken as stretching bars either continuously or discretely attached to the skin. 2.1. Potential energy In order to simplify the construction of Ritz bases that properly handle the discontinuities introduced by stiffeners and cracks into the displacement field, the original plate domain is decomposed into a small number of disjoint subdomains [13–17,19,20]. The decomposition illustrated in Fig. 1a consists of 20 subdomains, each of which is delimited by a stiffener, crack, plate edge, or a skin section. It corresponds to the decomposition with the minimum number of subdomains for that stiffened plate. A typical rectangular subdomain of size ai × bi , shown in Fig. 1b, has one stiffener attached to the skin at x = ai and another at y = bi . For the typical subdomain i, the strain energy of the skin is given by

Uskin =

1 2

 0

bi



ai 0

mT Am dxdy

(1)

where the strain and the membrane rigidity are

m = u,x

v,y

T

u,y + v,x 

A=

⎡

1

Eh ⎢ν ⎣ 1 − ν2 0

ν 1 0

⎤

0 0 ⎥. ⎦ 1−ν 2

(2)

The quantities u and v are the skin midsurface displacements in the x and y directions, a comma followed by x (or y) indicates differentiation with respect to x (or y), and E, ν and h are the Young’s modulus, the Poisson’s ratio and the skin

600

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

(a)

(b) Fig. 1. Cracked stiffened plate: (a) original plate domain decomposed into a minimum number of 20 subdomains; (b) a typical subdomain i.

thickness. The contribution of the stiffeners to the strain energy is

Ustiffener =

1 2

 0

ai

Ex Ax u¯ 2,x dx +

1 2

 0

bi

Ey Ay v¯ 2,y dy.

(3)

In this expression, Ex and Ey are the Youngs´ modulus of the horizontal and vertical stiffeners, Ax and Ay are their crosssection areas, u¯ (x ) and v¯ (y ) are the stiffener displacement fields.

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

601

At the subdomain level, a general expression for the potential energy of the externally applied loads takes the form

 Vi = −



ai 0



(qy v + qxy u )y=bi − (qy v + qxy u )y=0 dx −

bi 0



(qx u + qxy v)x=ai − (qx u + qxy v)x=0 dy

− Qx u¯ |x=ai + Qx u¯ |x=0 − Qy v¯ |y=bi + Qy v¯ |y=0

(4)

from which loads not actually applied must be removed. The load applied to a skin edge (or end of a stiffener) is not necessarily equal to the load applied to the opposite edge (or end), but for simplicity they are denoted here by the same symbol. Letting Ui = Uskin + Ustiffener , the potential energy  of the entire plate decomposed into subdomains reads

=



(Ui + Vi )

(5)

provided that u(x, y) and v(x, y) are continuous between subdomains, u¯ (x ) and v¯ (y ) are continuous between stiffeners [11]. A C 0 continuity is thus needed in the assembly of the subdomains, except at the crack locations. 2.2. Discretization of the potential energy For convenience, the nondimensional coordinates

ξ=

2 x−1 ai

η=

2 y−1 bi

(6)

are introduced so that − 1 ≤ ξ ≤ 1 and − 1 ≤ η ≤ 1 in a subdomain, as shown in Fig. 1b. The displacement fields are then locally approximated by

u= u¯ =

 

α j G j ( ξ , η ) = αT G α¯ j G¯ j (ξ ) = α¯ T G¯

v= v¯ =





β j H j (ξ , η ) = β T H

β¯ j H¯ j (η ) = β¯ T H¯

(7)

¯ , H and H ¯ collect the Ritz basis functions, and the vectors α, α¯ , β and β¯ collect unknown constant where the vectors G, G coefficients. In view of (6) and (7), expressions (1) and (3) read

1 T c k i ci 2 i

Uskin = where

1 ki = ai bi

Ustiffener =

1 1 −1 −1

⎡

¯ ,ξ G ¯ T dξ G ,ξ

−1

ai

⎤

0

⎡

(8)

1 1

NTu ANu dξ dη

−1 −1

NTu ANv dξ dη

−1 −1

NTv ANv dξ dη

1 1

sym.

Ex Ax 1

k¯ i = 2

1 T¯ c¯ ki c¯ i 2 i



0 Ey Ay bi

⎤

1 −1



¯ ,η H ¯ T,η dη H

    α¯ α T ⎦ T ⎣ ⎦ ⎣ a H Nu = Nv = ci = c¯ i = 0 i ,η ¯ β β T T bi GT,ξ

0T

(9)

bi H,ξ

ai G,η

assuming that the rigidities Ex Ax and Ey Ay are constant along the lengths of the stiffeners. Similarly, expression (4) reduces to

Vi = −cTi fi − c¯ Ti f¯i

(10)

with

ai fi = 2





1

−1

qxy G qy H

η=1 η=−1

bi dξ + 2



1 −1



Substitution of (8) and (10) into (5) yields

=

1 2

cTi ki ci +

ξ =1

qx G qxy H



 dη

f¯i =

ξ =−1

ξ =1 

¯ Qx G

ξη==1−1 .

¯ Qy H

(11)

η=−1





1 T¯ 1 KC − F . c¯ ki c¯ i − cTi fi − c¯ Ti f¯i = CT 2 i 2

(12)

602

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

The stiffness matrix K, the applied force vector F and the vector of unknown constant coefficients C of the whole plate can be written as

⎡

k1 ⎢0 ⎢. ⎢ .. ⎢

K=⎢

⎤

··· ··· .. .

0 k2 .. .

k¯ 1 0 .. .

⎢ ⎢ ⎣

⎧ ⎫ f1 ⎪ ⎪ ⎪ ⎪ ⎪ f2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ... ⎪ ⎬

⎧ ⎫ c1 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ c2 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ .. ⎪ ⎬

.

.

⎥ ⎥ ⎥ . ⎥ C= . ⎥ F = f¯ · · ·⎥ c¯ 1 ⎪ 1⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎥ ⎪ ⎪ ⎪ ⎪ ¯ c¯ 2 ⎪ ⎪ ⎪ f2 ⎪ · · ·⎦ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ .. ⎩ .. ⎪ ⎭ ⎩ .. ⎪ ⎭

0 k¯ 2 .. .

.

(13)

2.3. Continuity constraints The Ritz basis G (and H) is constructed in the separable form with p (and r) functions φ i (ξ ) and q (and s) functions φ i (η) in the ξ and η directions, respectively, as

G=

⎧ ⎫ φ ξ  ⎪ ⎪ 1 ( ) gη ⎪ ⎪ ⎨ φ2 (ξ )gη ⎬

H=

.

.. ⎪ ⎪ ⎪ ⎪ ⎩ ⎭ φ p (ξ )gη pq×1

⎧ ⎫ φ ξ  ⎪ ⎪ 1 ( ) hη ⎪ ⎪ ⎨ φ2 (ξ )hη ⎬ ⎪ ⎪ ⎩

.. .

φr (ξ )hη

(14)

⎪ ⎪ ⎭ rs×1

where

gη = φ1 (η )

φ2 ( η )

φq (η )T hη = φ1 (η )

···

φ2 ( η )

···

φs (η )T .

(15)

The functions φ i (ξ ) are selected from the following complete set of hierarchic polynomials:

φ1 ( ξ ) = φi ( ξ ) =

1 (1 − ξ ) 2 (i−1 )/2 n=0

1 2

φ2 ( ξ ) = ( 1 + ξ )

(−1 )n (2i − 2n − 5)!! i−2n−1 ξ 2n n! (i − 2n − 1 )!

i = 3, 4, . . .

(16)

Here i!! = i(i − 2)(2 or 1) and 0!! = ( − 1)!! = 1. Moreover, only the integer part of the summation upper bound (i − 1)/2 is retained and the functions φ i (η) follow directly from φ i (ξ ) by replacing ξ with η. While the first and second elements in (16) are the linear Lagrange polynomials, the remaining ones stem from a single integration of Legendre polynomials [21]. The polynomial set (16) is chosen herein because of its superior convergence properties in capturing localized features [22]. The C 0 continuity between subdomains and all the geometric boundary conditions are controlled by the first two functions φ i (ξ ) at ξ = ±1, and by the first two functions φ i (η) at η = ±1. Higher order functions (i > 2) possess zero displacement at each end, being significant in the solution refinement. The expressions (14) may be rewritten in a more compact form by defining the following operator, called Kronecker product [24,25], between two matrices P and Q:

⎡

p11 Q p P  Q = ⎣ 21 Q ...

⎤

··· · · ·⎦ .. .

p12 Q p22 Q .. .

(17)

where pij are the elements of P. Thus,

G = gξ  gη

H = hξ  hη

gξ = φ1 (ξ )

φ2 ( ξ )

(18)

with

···

φ p (ξ )T hξ = φ1 (ξ )

φ2 ( ξ )

···

φr (ξ )T .

(19)

For the purpose of illustration, suppose that the C 0 continuity of the displacement u(ξ , η) is enforced between the subdomain i and the subdomain i + 1 on its right:

ui (1, η ) = ui+1 (−1, η ).

(20)

Substituting (7) and (18) into (20) gives

 i T  gξ α

 gη

ξ =1

=



αi+1

T 

gξ  gη

ξ =−1

.

(21)

Since the set of functions (16) reduces gξ (1) and gξ ( − 1) to the unit vectors

gξ (1 ) = 0

1

0

···

0T

gξ (−1 ) = 1

0

0

···

0T ,

(22)

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

603

Eq. (21) may be rewritten as q 

αqi + j φ j (η ) =

j=1

q 

α ij+1 φ j (η )

(23)

j = 1, 2, . . . , q.

(24)

j=1

which holds for any η if

αqi + j = α ij+1

Three other groups of similar constraints must be generated by meeting the continuity of u(ξ , η) between the subdomain i and adjacent subdomains on the left, above and below. Moreover, continuity requirements on v(ξ , η) between subdomains, u¯ (ξ ) between horizontal stiffeners and v¯ (η ) between vertical stiffeners generate additional constraints. While the continuity of u(ξ , η) creates linear dependency between the coefficients α j , the continuities of v(ξ , η), u¯ (ξ ) and v¯ (η ) create linear dependency between β j , α¯ j and β¯ j respectively. Attachment of the stiffeners to the skin also generates constraints involving linear dependency between α j and α¯ j for horizontal stiffeners, and β j and β¯ j for vertical stiffeners. Collecting all the independent coefficients in the vector C¯ , one may state the relation

C = TC¯

(25)

where matrix T contains only binary entries 1 and 0. Substitution of (25) into (12) yields

 = C¯ T

1 2



¯ C¯ − F¯ K

(26)

with

¯ = TT KT F¯ = TT F K

(27)

being the stiffness matrix and force vector associated with the independent coefficients. Finally, the discrete form of the equations that govern the response of the cracked stiffened plate is given by the stationary condition of  with respect to C¯ :

∂ ¯ C¯ = F¯ . =0 ⇒ K ∂ C¯

(28)

2.4. Integral evaluation The integration to be performed in the force vector fi depends on the type of load. For instance, suppose that the evaluation of

bi 2



1

ξ =1

−1

qx G|ξ =−1 dη

(29)

is to be carried out for qx constant. From (18), one writes

bi qx 2



1



−1

gξ  gη

bi qx ξ =1 dη = ξ =−1 2

!

1 −1

" gξ dη  gη |ξξ =1 . =−1

(30)

1 Thus, the integral to be computed is simply −1 gξ dη. The force vector fi may be evaluated in advance for a sufficiently large number of functions and stored. Once low order vectors are embedded in high order vectors, lower approximations are accomplished by simply removing exceeding lines from the previously stored vector. The above precomputation approach can also be extended to the evaluation of the stiffness matrices ki and k¯ i . For instance, suppose the evaluation of the submatrix



1



−1

1 −1

NTu A Nv dξ dη

(31)

of ki is to be carried out. Matrices Nu and Nv may be decomposed into the products

⎡

⎤⎡

GT,ξ

bi

0

0

Nu = ⎣ 0

1

0 ⎦⎣ 0T

⎢

0

0

ai

⎥⎢

⎤

⎡

1

⎥ ⎢ ⎦ = Ju N¯ u Nv = ⎣0

GT,η

0

⎤⎡

0T

⎤

0

0

ai

¯v 0 ⎦⎣HT,η ⎦ = Jv N

0

bi

⎥⎢

⎥

(32)

HT,ξ

¯ u and N ¯ v . Eq. (31) then with the subdomain geometry isolated in Ju and Jv , and the coordinates ξ and η isolated in N becomes



1 −1



1 −1

¯ Tu JTu A Jv N ¯ v dξ dη. N

(33)

604

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

Fig. 2. Centrally cracked plate with the solution domain decomposed into two subdomains.

In order to perform the integration, it is convenient to introduce the “vec” operator which converts a matrix into a vector by stacking its columns one underneath the other. For a given matrix P, for instance,

vec P =  p11

p21

···

p1n

· · · T .

p2n

(34)

The definitions of “vec” operator and Kronecker product stated in (17) lead to the following useful property



vec(PQR ) = RT  P vec Q

(35)

whenever the product PQR of these three matrices is defined [24,25]. Applying the “vec” operator to (33), followed by the use of the above property, gives

! 1 

1

vec −1

−1

"

¯ Tu JTu A Jv N ¯ v dξ dη N

=

! 1  −1

1

−1

"



¯ Tv  N ¯ Tu dξ dη vec JTu A Jv . N

(36)

Now the integrand is clearly independent of both material and geometry, and then can be promptly precomputed for a generic rectangular subdomain. Appendix A shows that all the integrals present in ki and k¯ i can be completely evaluated by means of the three smaller matrices



L00 cd =

1 −1

c Td dx L10 cd =



1 −1

c,x Td dx L11 cd =



1

−1

c,x Td,x dx,

(37)

computed in advance for a sufficiently large number of functions and stored. Note that the variable x stands for either ξ or

1 η, indices c and d indicate the number of functions in c and d , and L01 = −1 d Tc,x dx = (L10 )T . As explained for fi , lower dc cd approximations are accomplished by simply removing exceeding lines and columns from the previously stored matrices L. Such a procedure reduces computational cost during the solution refinement. 3. Numerical tests A code in MATLAB [26] language has been written to carry out the numerical tests. Matrices L in (37) have been evaluated by symbolic computing to circumvent round-off errors, and stored in advance for the number of functions c = d = 100. Stiffened plates with single or multiple cracks of length 2a have been chosen to illustrate the performance of the proposed Ritz scheme. Stiffeners may be either continuously or discretely attached. All the plates have length 2L = 60, width 2b = 20, thickness h = 1, and skin material given by E = 1 and ν = 0.3 (data in consistent but unspecified units). The displacement fields (7) are locally approximated within all subdomains with the same number p of functions (16) in the ξ and η directions. 3.1. Centrally cracked plate In this first example, the centrally cracked plate to be treated is unstiffened and subjected to a uniformly distributed load q as shown in Fig. 2. Due to symmetry, only the highlighted quarter of the plate is defined as the solution domain. The geometric boundary conditions

u=0

at x = L, b + a ≤ y ≤ 2b

v = 0 at y = b

(38)

are exactly satisfied by (7). As the plate is supposed to be in a state of plane stress, the stress intensity factor KI is evaluated by means of [3]

KI2 = −

E ∂ . h ∂a

(39)

The derivative ∂ /∂ a of the potential energy is approximated by the central finite difference

∂ (a + a ) − (a − a ) ≈ ∂a 2 a

(40)

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

605

Table 1 Values of β /β ref for a centrally cracked plate. a/b

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 a

β ref

p

[5]

10 (351)a

20 (1501)

30 (3451)

40 (6201)

50 (9751)

1.006 1.024 1.058 1.109 1.186 1.303 1.487 1.814

0.932 0.970 0.977 0.979 0.980 0.978 0.974 0.964

0.985 0.992 0.994 0.995 0.995 0.995 0.994 0.992

0.993 0.997 0.998 0.998 0.998 0.998 0.998 0.997

0.996 0.998 0.999 0.999 0.999 0.999 0.999 0.999

0.998 0.999 0.999 0.999 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0

Number of degrees of freedom NDF.

Fig. 3. Plate with a central crack a/b = 0.5: error ε β = |β /β ref − 1| versus number of degrees of freedom NDF for h- and p-refinements. The five points on Nastran curve correspond to 1 × 2, 2 × 4, 4 × 8, 8 × 16 and 16 × 32 meshes.

with a = 0.01a. Table 1 shows the results for several crack lengths in terms of the parameter

β=

KI , K0

(41)

√

where K0 = q π a/h corresponds to the stress intensity factor for an infinite plate (L = b = ∞). The number of degrees of freedom NDF (dimension of C¯ after inclusion of all the geometric boundary conditions) is also shown in the table. Normalized by the empirical formula

Kref = βref

q√ πa h

#  a 2  a 4 $% π a βref = 1 − 0.025 + 0.06 sec b

b

2b

(42)

found in [5], the results are obtained with the quarter of the plate decomposed into the minimum number of two subdomains, as sketched in Fig. 2, under p-refinement (i.e., when convergence is sought by progressively increasing the number p of functions φ i ). The results converge toward the reference values monotonically and from below. The error

   β   εβ =  − 1 βref

(43)

for a crack length a/b = 0.5 is plotted against the number of degrees of freedom NDF on a log-log scale in Fig. 3 using different refinements. The p-refinement is carried out with the quarter of the plate decomposed into two subdomains. In the h-refinement (i.e., when convergence is sought by progressively increasing the number of subdomains), it is adopted p = 2 or p = 5 and the number of subdomains varies from two to a higher value by successively splitting each subdomain into four smaller ones. It is also included in the figure results from Nastran models of CQUAD4 elements (1 × 2, 2 × 4, 4 × 8, 8 × 16 and 16 × 32 meshes) [27]. As the slope of each curve is proportional to the convergence rate, it can be seen that this uniform h-refinement, including that related to Nastran results, have similar rates which are smaller than the rate achieved by the p- refinement. In this sense, the advantage of using the proposed methodology with the minimum number of subdomains under p-refinement is clear.

606

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614 Table 2 Normalized strain energy U/U0 for a centrally cracked plate. a /b

p

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 a

10 (351)a

20 (1501)

30 (3451)

40 (6201)

50 (9751)

1.004 1.019 1.046 1.087 1.146 1.232 1.357 1.555

1.005 1.021 1.049 1.091 1.153 1.241 1.371 1.579

1.005 1.021 1.049 1.092 1.154 1.242 1.374 1.584

1.005 1.021 1.050 1.092 1.154 1.243 1.374 1.585

1.005 1.021 1.050 1.093 1.154 1.243 1.375 1.586

Number of degrees of freedom NDF.

Fig. 4. Plate with two coplanar cracks and solution domain decomposed into three subdomains.

Table 2 shows the strain energy obtained with the quarter of the plate decomposed into two subdomains under prefinement and normalized by the strain energy U0 = 2q2 bL/Eh of the corresponding uncracked plate. The expected monotonic convergence from below associated with the conventional Ritz solution is preserved even in the presence of singularities. The remarkable approximation capability exhibited by the proposed methodology when using the minimum number of subdomains under p-refinement motivates the adoption of similar approach in all the remaining examples.

3.2. Plate with two coplanar cracks The plate of Fig. 2 is now taken with two identical coplanar cracks at distance s apart, as sketched in Fig. 4. Due to symmetry, only the highlighted quarter of the plate is defined as the solution domain which will be decomposed into the minimum number of three subdomains. The geometric boundary conditions, this time given by



u = 0 at x = L,

b ≤ y ≤ b + 2s b + 2s + 2a ≤ y ≤ 2b

v = 0 at y = b,

(44)

are exactly satisfied by (7). For a fixed value of the crack length 2a/b = 0.5, Table 3 shows how the parameter β evaluated at the crack tips A and B changes with the distance s. As the unstiffened plate with a single crack, the convergence toward the reference values is monotonic from below. For the purpose of conveniently observing the interaction effect between the cracks on β , the numerical results are also plotted in Fig. 5. The β solution √ at tip B increases asymptotically as s → 0. At tip A, on the other hand, the solution is bounded by the value 1.677 (= 2 × 1.186) as s → 0 because the two cracks become a single crack with twice the original length of each crack (see Table 1 for a/b = 0.5). 3.3. Centrally cracked plate with edge stiffeners continuously attached Stiffeners with rigidity Ex Ax = 1 are now continuously attached to the edges of the plate of Fig. 2, as sketched in Fig. 6, and submitted to the end load Q = Ex Ax q/Eh to satisfy strain compatibility between stiffeners and skin [6]. The geometric boundary conditions to be satisfied are given in (38). Table 4 shows the results for several crack lengths in terms of the parameter β obtained with the quarter of the plate decomposed into two subdomains, as sketched in Fig. 6, under p-refinement. The results are normalized by the empirical formula

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

607

Table 3 Values of β /β ref at the crack tips A and B shown in Fig. 4 (2a/b = 0.5). s/2a

β ref

p

Nastran (961,600) 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 a

a

1.203∗ 1.401∗∗ 1.169 1.249 1.156 1.193 1.154 1.167 1.161 1.155 1.179 1.152 1.216 1.158 1.293 1.176 1.489 1.218

10 (522)

20 (2242)

30 (5162)

40 (9282)

50 (14,602)

0.957 0.914 0.964 0.948 0.965 0.960 0.964 0.964 0.963 0.965 0.960 0.965 0.956 0.964 0.947 0.962 0.930 0.958

0.990 0.984 0.992 0.989 0.992 0.991 0.992 0.992 0.992 0.992 0.991 0.992 0.990 0.992 0.988 0.992 0.979 0.990

0.997 0.994 0.997 0.996 0.997 0.997 0.997 0.997 0.997 0.997 0.997 0.998 0.997 0.997 0.996 0.997 0.994 0.997

0.999 0.998 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.999 0.998 0.999

1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0

Number of degrees of freedom NDF

∗

∗∗

crack tip A

crack tip B.

Fig. 5. Plate with two coplanar cracks of length 2a/b = 0.5: β versus s/2a at the crack tips A and B evaluated by the Ritz model for p = 50.

Kref =

βref

q√ πa h

 a 2  a 4  a 6  a 8 βref = 1 + 0.4970 + 0.3366 + 0.2306 + 0.1611 b b b b  a 10  a 12  a 14 + 0.1156

b

+ 0.0843

b

+ 0.0627

b

(45)

608

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

Fig. 6. Centrally cracked plate with edge stiffeners continuously attached. The solution domain is decomposed into two subdomains. Table 4 Values of β /β ref for a centrally cracked plate with edge stiffeners continuously attached. a /b

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 a

β ref

p

[6]

10 (351)a

20 (1501)

30 (3451)

40 (6201)

50 (9751)

1.005 1.020 1.048 1.089 1.150 1.237 1.366 1.564

0.928 0.965 0.972 0.975 0.976 0.976 0.974 0.969

0.980 0.987 0.989 0.990 0.990 0.990 0.990 0.990

0.989 0.992 0.993 0.993 0.993 0.993 0.993 0.993

0.991 0.993 0.994 0.994 0.994 0.994 0.994 0.995

0.993 0.994 0.994 0.994 0.994 0.994 0.994 0.995

Number of degrees of freedom NDF.

Fig. 7. Centrally cracked plate with edge stiffeners continuously attached: β versus a/b. The Ritz models have NDF = 351 and 9751 for p = 10 and 50, respectively, and the Nastran model has NDF = 961,600.

found in Isida et al. [6]. Plots of β against a/b for p = 10 (NDF = 351) and 50 (NDF = 9751) are given in Fig. 7, which includes the results obtained from a Nastran model (NDF = 961,600) made up of 1200 × 400 CQUAD4 elements for the skin and 1200 CROD elements for the stiffener. It is observed that the convergence toward the reference values is monotonic from below, as before. 3.4. Centrally cracked plate with internal stiffeners continuously attached The stiffeners of the plate of Fig. 6 are now supposed to be internally attached, as depicted in Fig. 8. The minimum number of subdomains into which the quarter of the plate must be decomposed changes, in this case, from two to three, as shown in Fig. 8. Results are compared in Table 5 for several crack lengths with those obtained from a Nastran model (NDF = 240,800) made up of 600 × 200 CQUAD4 elements for the skin and 600 CROD elements for the stiffener. The stress intensity factor decreases (β decreases) as the crack approaches a stiffener, indicating that the stiffener aids to restrain the crack or slow

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

609

Fig. 8. Centrally cracked plate with internal stiffeners continuously attached and the solution domain decomposed into three subdomains. Table 5 Values of β /β ref for a centrally cracked plate with internal stiffeners continuously attached. a/b

β ref Nastran (240,800)

0.1 0.2 0.3 0.4 0.42 0.44 0.46 0.48 0.5 0.52 0.54 0.56 0.58 0.6 0.7 0.8 a

0.999 1.012 1.030 1.038 1.034 1.023 0.997 0.930 – 1.653 1.532 1.501 1.497 1.508 1.666 2.007

p a

10 (522)

20 (2242)

30 (5162)

40 (9282)

50 (14,602)

0.937 0.973 0.981 0.988 0.991 0.995 1.004 1.039 – 1.066 1.012 0.997 0.991 0.988 0.978 0.967

0.990 0.995 0.996 0.998 0.998 0.999 1.001 1.010 – 1.016 1.004 1.0 0 0 0.998 0.998 0.996 0.994

0.999 0.999 0.999 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.001 – 1.004 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 0.999 0.999

1.001 1.001 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 0.999 – 0.999 1.0 0 0 1.0 0 0 1.0 0 0 1.0 0 0 1.001 1.001

1.003 1.001 1.001 1.001 1.0 0 0 1.0 0 0 1.0 0 0 0.997 – 0.997 1.0 0 0 1.0 0 0 1.001 1.001 1.001 1.002

Number of degrees of freedom NDF.

Fig. 9. Centrally cracked plate with internal stiffeners continuously attached: β versus a/b. The Ritz models have NDF = 522 and 14,602 for p = 10 and 50, respectively, and the Nastran model has NDF = 240,800.

down the propagation. Once the crack reaches the stiffener, it is assumed that the stiffener is completely and suddenly severed. The load previously carried by the stiffener is then shed to the remaining net section, increasing immediately the stress intensity factor. For a better visualization of the interaction effect between the crack and stiffeners on β , Fig. 9 depicts the solution for p = 10 (NDF = 522) and 50 (NDF = 14,602), highlighting the jump in the β value as a result of the stiffener being completely and suddenly severed. It can be seen that the convergence is monotonic from below provided that the crack is not very close to the stiffener (severed or not), otherwise (a/b = 0.46, 0.48, 0.52, 0.54) the convergence is still monotonic but from above. Although not

610

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

Fig. 10. Centrally cracked plate with internal stiffeners discretely attached. The solution domain is decomposed into three subdomains in the y direction and a number of subdomains in the x direction that depends on the pitch d.

Fig. 11. Centrally cracked plate with discretely attached stiffeners: β versus a/b. The Ritz and Nastran models have NDF = 145,182 and 962,790 for a/b = 0.3, NDF = 435,262 and 962,770 for a/b = 0.1, and NDF = 870,382 and 962,740 for a/b = 0.05. Table 6 Values of β /β ref for a centrally cracked plate with internal stiffeners discretely attached (d/b = 0.3). a /b

β ref Nastran (962,790)

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 a

1.002 1.015 1.034 1.053 1.063 1.075 1.123 1.239

p a

10 (5057)

20 (22,092)

30 (51,122)

40 (92,152)

50 (145,182)

0.996 0.997 0.996 0.994 1.036 0.988 0.980 0.973

1.001 1.0 0 0 1.0 0 0 0.999 1.011 0.999 0.998 0.995

1.002 1.001 1.001 1.001 1.008 1.005 1.006 1.006

1.002 1.001 1.001 1.002 1.008 1.009 1.011 1.012

1.002 1.001 1.001 1.003 1.008 1.011 1.015 1.017

Number of degrees of freedom NDF.

shown here for the sake of brevity, the strain energy converges monotonically and from below regardless of the crack location. 3.5. Centrally cracked plate with internal stiffeners discretely attached This last example takes the stiffeners of the plate of Fig. 8 as being discretely attached. The quarter of the plate is decomposed into three subdomains in the y direction, as before. Since the fasteners must be located at the corners of the subdomains, the minimum number of subdomains in the x direction depends on the pitch d, as illustrated in Fig. 10. Thus, the smaller the pitch d, the greater this minimum number. Results for d/b = 0.3 are compared in Table 6 for several crack lengths with those obtained from a Nastran model (NDF = 962,790) made up of 1200 × 400 CQUAD4 elements for the skin and 1200 CROD elements for the stiffener. The chosen value of d implies ten subdomains in the x direction.

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

611

The rate of reduction at which the stress intensity factor decreases as the crack approaches a stiffener is smaller than that for continuously attached stiffeners, as can be observed from Figs. 9 and 11. Once the crack reaches the stiffener (a/b = 0.5), it is assumed that the stiffener detaches from the skin at that point but remains intact. No jump is observed in the β value, indicating that the stiffener continues to limit the crack growth. The smaller the pitch d, the greater the action of the stiffener as a crack stopper. As in the previous example, the convergence is monotonic from below provided that the crack is not very close to the stiffener, otherwise (a/b = 0.5) the convergence is still monotonic but from above. 4. Conclusions The methodology proposed in this paper to evaluate stress intensity factors of cracked stiffened plates by the Ritz method is simple and remarkably accurate. This has been confirmed by comparing its results with those either available in the literature or obtained from converged finite element solutions. The methodology is also versatile in the sense that can accommodate both continuously and discretely attached stiffeners. Some points deserve to be highlighted: (a) the decomposition of the solution domain into the minimum number of subdomains, uniquely defined according to the crack and stiffener locations, is recommended; (b) the subdomain stiffness matrix should be derived using only the three matrices L for the sake of computational efficiency; (c) the convergence of the strain energy monotonically and from below, even in the presence of cracks, indicates that the proposed methodology preserves this important property of the conventional Ritz method; (d) the convergence of the stress intensity factor is monotonic from below provided that the crack is not very close to the stiffener, otherwise the convergence is still monotonic but from above; (e) the stress intensity factor decreases as the crack approaches a stiffener, as expected [28], indicating that the stiffener aids to restrain the crack or slow down its propagation. The rate of reduction before the crack reaches the stiffener is greater for continuously attached stiffeners. However, discretely attached stiffeners are better suited as crack stoppers since they remain intact after the crack passage; (f) the adoption of triangular subdomains, whose integrals could still be evaluated in closed form, should provide the methodology proposed here with the necessary tools to model irregularly shaped plates with inclined cracks. It is reasonable to assume that the numerical behavior observed for rectangular subdomains will remain valid for triangular subdomains, but this demands further detailed development and examination. Appendix A The purpose of this appendix is to demonstrate that the matrices ki and k¯ i given in (9) can be completely evaluated by means of the three smaller matrices



L00 cd =

1

−1

c Td dx L10 cd =



1 −1

c,x Td dx L11 cd =



1

−1

c,x Td,x dx,

(A.1)

where the variable x stands for either ξ or η, indices c and d indicate the number of functions in c and d , and L01 = dc

1 T dx = (L10 )T .   d c,x −1 cd The Kronecker product (17) satisfies the properties

(A  B )T = AT  BT (A  B )(C  D ) = AC  BD  T

a  b = vec ba ,

(A.2)

where A, B, C, D are matrices and a, b are vectors [24,25]. Of course, the product AC and BD must be defined for the second of the above properties to make sense. In general, ab = ba. However, there exists a matrix

Knm =

n  

ei  Im  eTi ,

(A.3)

i=1

where the unit vector ei denotes the i-th column of the identity matrix In of order n, such that

Knm (a  b ) = b  a

(A.4)

for the two vectors a and b of dimensions m and n, respectively. It is now convenient to introduce the matrix

Pmnpq = Im  Knp  Iq

(A.5)

which has the useful property

Pmnpq (a  c  b  d ) = a  b  c  d.

(A.6)

612

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

The vectors a, b, c and d have dimensions m, n, p and q, respectively. It is a simple matter to verify (A.6) using the second of properties (A.2) and (A.4):

Pmnpq (a  c  b  d ) = (Im  Knp  Iq )(a  c  b  d ) = Im a  (Knp  Iq )(c  b  d ) = Im a  Knp (c  b )  Iq d = a  b  c  d.

(A.7)

The matrix to be integrated in (36) is

⎡



⎛⎡

¯v N ¯u N

⎤

⎡

⎤⎞T

0T GT,ξ T = ⎝⎣H,η ⎦  ⎣ 0T ⎦⎠ HT,ξ GT,η

T

⎤T

0T  GT,ξ ⎢ 0T  0T ⎥ ⎢ T ⎥ ⎢ 0  GT,η ⎥ ⎢ T ⎥ ⎢H,η  GT,ξ ⎥ ⎢ T ⎥ T ⎥ =⎢ ⎢ HT,η  0T ⎥ . ⎢H,η  G,η ⎥ ⎢ ⎥ ⎢HT  GT ⎥ ,ξ ⎥ ⎢ ,ξ ⎣ HT  0T ⎦ ,ξ HT,ξ  GT,η

(A.8)

The Ritz bases (18) may be substituted into the last column of (A.8) to give



HT,ξ  GT,η

T

= H,ξ  G,η



= hξ  hη

,ξ





gξ  gη

,η

= hξ ,ξ  hη  gξ  gη,η



= Prspq hξ ,ξ  gξ  hη  gη,η







= Prspq vec gξ Thξ ,ξ  vec gη,η Thη



(A.9)

and



1



−1

1



−1

HT,ξ



T GT,η d

#

!

ξ dη = Prspq vec

1 −1

" gξ 

T hξ ,ξ d

ξ

!  vec

1

−1

"$ gη,η 

T hη d

η

.

(A.10)

Recalling that the vectors gξ , gη , hξ and hη have dimensions p, q, r and s, as defined in (15) and (19), the expression (A.10) takes the form



1

−1



1 −1



HT,ξ  GT,η

T



10 dξ dη = Prspq vec L01 pr  vec Lqs

(A.11)

in accordance with the notation (A.1). Repeating the above process for the remaining columns of (A.8), one gets

⎡ ⎢ ⎢ ⎢  ⎢ ⎢ Prspq vec  1 1 ⎢ T ¯v N ¯ u dξ dη = ⎢ N  ⎢ −1 −1 ⎢ Prspq vec ⎢  ⎢ Prspq vec ⎢ ⎣  

0T 0T 0T 01 L10 pr  vec Lqs 0T 11 L00 pr  vec Lqs 00 L11 pr  vec Lqs T 0

10 Prspq vec L01 pr  vec Lqs

Similarly, it follows that

⎤T ⎥ ⎥ T ⎥ ⎥ ⎥ ⎥ T ⎥ ⎥ . ⎥ T ⎥ ⎥ ⎥ ⎦ T

(A.12)

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

⎡



00 P pqpq vec L11 pp  vec Lqq 0T

⎢ ⎢ ⎢ P vec ⎢ pqpq ⎢  1 1 ⎢ T ¯u N ¯ u dξ dη = ⎢ N ⎢ −1 −1 ⎢ ⎢  ⎢ P pqpq vec ⎢ ⎣  

10 L01 pp  vec Lqq 0T 0T 0T 01 L10 pp  vec Lqq T 0

11 P pqpq vec L00 pp  vec Lqq

⎡ ⎢ ⎢ ⎢ ⎢ ⎢  1 1 ⎢ Prsrs vec T ¯ ¯ Nv  Nv d ξ d η = ⎢  ⎢ −1 −1 ⎢ Prsrs vec ⎢ ⎢  ⎢ ⎣ Prsrs vec  

0T 0T 0T 0T 11 L00 rr  vec Lss 01 L10 rr  vec Lss T 0 10 L01 rr  vec Lss

00 Prsrs vec L11 rr  vec Lss

613

 T ⎤T ⎥ T ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ T ⎥ ⎥ ⎥ ⎦ T ⎤T

⎥ ⎥ ⎥ ⎥ T ⎥ ⎥ T ⎥ ⎥ . ⎥ ⎥ T ⎥ ⎥ ⎦ T

(A.13)

Once the matrices (A.12) and (A.13) are known, the matrix ki may be directly evaluated. ¯ and H ¯ in (7) be vectors of dimensions m and The matrices L can also be employed in the evaluation of matrix k¯ i . Let G n, respectively. Thus,

!

1

vec −1

"

¯ ,ξ G ¯ T dξ G ,ξ

 =  =

1 −1 1 −1



¯ ,ξ G ¯ T dξ vec G ,ξ ¯ ,ξ  G ¯ ,ξ dξ G

= vec L11 mm

(A.14)

= vec L11 nn .

(A.15)

and, in a similar manner,

!

1

vec −1

"

¯ ,η H ¯ T,η dξ H

The evaluation of k¯ i follows directly from (A.14) and (A.15). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20]

O. Bedair, Analysis and limit state design of stiffened plates and shells: a world view, Appl. Mech. Rev. 62 (2) (2009) 020801-1–020801-16. O. Bedair, Recent developments in modeling and design procedures of stiffened plates and shells, Recent Pat. Eng. 7 (3) (2013) 196–208. T.L. Anderson, Fracture Mechanics: Fundamentals and Applications, 4th ed., CRC Press, Boca Raton, 2017. D.P. Rooke, D.J. Cartwright, Compendium of Stress Intensity Factors, HMSO, London, 1976. H. Tada, P.C. Paris, G.R. Irwin, The Stress Analysis of Cracks Handbook, 3rd ed., ASME Press, New York, 20 0 0. M. Isida, S. Tagami, Y. Itagaki, Stress concentration due to a central transverse crack in a strip reinforced on either sides, J. Jpn. Soc. Aeronaut. Eng. 10 (100) (1962) 141–149 (in Japanese). C.C. Poe Jr, The effect of broken stringers on the stress intensity factor for a uniformly stiffened sheet containing a crack, NASA TM X-71947, 1973. B. Andersson, U. Falk, I. Babuška, T. von Petersdorff, Reliable stress and fracture mechanics analysis of complex components using a h-p version of FEM, Int. J. Numer. Meth. Eng. 38 (13) (1995) 2135–2163. S.A. Fawaz, B. Andersson, Accurate stress intensity factor solutions for corner cracks at a hole, Eng. Fract. Mech. 71 (9–10) (2004) 1235–1254. A.K. Srivastava, P.K. Arora, H. Kumar, Numerical and experiment fracture modeling for multiple cracks of a finite aluminum plate, Int. J. Mech. Sci. 110 (2016) 1–13. B. Szabó, I. Babuška, Introduction to Finite Element Analysis: Formulation, Verification and Validation, John Wiley, Chichester, 2011. F.L.S. Bussamra, E. Lucena Neto, W.M. Ponciano, Simulation of stress concentration problems by hexahedral hybrid-Trefftz finite element models, CMES 99 (3) (2014) 255–272. J. Yuan, S.M. Dickinson, The flexural vibration of rectangular plate systems approached by using artificial springs in the Rayleigh-Ritz method, J. Sound Vib. 159 (1) (1992) 39–55. K.M. Liew, K.C. Hung, M.K. Lim, A solution method for analysis of cracked plates under vibration, Eng. Fract. Mech. 48 (3) (1994) 393–404. M. Barrette, A. Berry, O. Beslin, Vibration of stiffened plates using hierarchical trigonometric functions, J. Sound Vib. 235 (5) (20 0 0) 727–747. Y.V.S. Kumar, J.K. Paik, Buckling analysis of cracked plates using hierarchical trigonometric functions, Thin-Walled Struct. 42 (5) (2004) 687–700. T.D. Dang, R.K. Kapania, Ritz approach for buckling prediction of cracked-stiffened structures, J. Aircraft 50 (3) (2013) 965–974. O. Beslin, J. Nicolas, A hierarchical functions set for predicting very high order plate bending modes with any boundary conditions, J. Sound Vib. 202 (5) (1997) 633–655. R. Vescovini, C. Bisagni, Semi-analytical buckling analysis of omega stiffened panels under multi-axial loads, Compos. Struct. 120 (2015) 285–299. S.G.P. Castro, M.V. Donadon, Assembly of semi-analytical models to address linear buckling and vibration of stiffened composite panels with debonding defect, Compos. Struct. 160 (2017) 232–247.

614

B.H.S. Oliveira, E. Lucena Neto and F.A.C. Monteiro / Applied Mathematical Modelling 73 (2019) 598–614

[21] D.C. Zhu, Development of hierarchical finite element methods at BIAA, in: Proceedings of the International Conference on Computational Mechanics, Vol. 1, Tokyo, 1986, pp. I-123–I-128. [22] L.N. Yshii, E. Lucena Neto, F.A.C. Monteiro, R.C. Santana, Accuracy of the buckling predictions of anisotropic plates, J. Eng. Mech. 144 (8) (2018) 04018061-1–04018061-9. [23] B.A. Szabó, Some recent developments in finite element analysis, Comp. Maths. Appl. 5 (2) (1979) 99–115. [24] H.V. Henderson, S.R. Searle, The vec-permutation matrix, the vec operator and Kronecker product, Lin. Multilin. Alg. 9 (4) (1981) 271–288. [25] K.M. Abadir, J.R. Magnus, Matrix Algebra, Cambridge University Press, Cambridge, 2005. [26] MATLAB, version R 2016a, MathWorks, Natick, 2016. [27] NX Nastran, Element Library Reference, Siemens PLM Software, Plano, 2014. [28] R.J. Dexter, H.N. Mahmoud, Predicting Stable Fatigue Crack Propagation in Stiffened Panels, Report no. SSC-435, Ship Structure Committee, Washington, D.C., March 2004.