Bolted double-lap composite joints under mechanical and thermal loading

International Journal of Solids and Structures 38 (2001) 7801±7837

www.elsevier.com/locate/ijsolstr

Bolted double-lap composite joints under mechanical and thermal loading V. Kradinov a, A. Barut a, E. Madenci a,*, D.R. Ambur b a

Department of Aerospace and Mechanical Engineering, University of Arizona, P.O. Box 210119, Tucson, AZ 85721, USA b Mechanics and Durability Branch, NASA Langley Research Center, Hampton, VA 23665, USA Received 13 October 2000; in revised form 3 April 2001

Abstract This study concerns the determination of the contact stresses and contact region around bolt holes and the bolt load distribution in single- and double-lap joints of composite laminates with arbitrarily located bolts under general mechanical loading conditions and uniform temperature change. The unknown contact stress distribution and contact region between the bolt and laminates and the interaction among the bolts require the bolt load distribution, as well as the contact stresses, to be as part of the solution. The present method is based on the complex potential theory and the variational formulation in order to account for bolt sti€ness, bolt-hole clearance, and ®nite geometry of the composite laminates. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Bolt; Double lap; Joint; Composite; Mechanical; Thermal

1. Introduction Bolted joints provide the primary means for transferring load among composite components in the construction of aircraft and space structures. The stress state in a bolted joint is dependent primarily on the dimensions of the planar geometry, loading conditions, degree of material anisotropy, bolt-hole clearance, bolt ¯exibility, and friction between the laminates. Also, aircraft and space vehicles traveling at supersonic and hypersonic speeds can experience high temperature excursions. The in¯uence of thermal expansions can be signi®cant and may di€er signi®cantly among the materials for the bolts and laminates. As a result, high thermal stresses may develop as the temperature increases and may alter the bolt load distribution. Therefore, accurate determination of the stresses in bolted laminates under both mechanical and thermal loading is essential for reliable strength evaluation and failure prediction. A considerable amount of work on the behavior of composite joints with a single bolt exists in the literature. These studies investigated the stress distribution around a pin-loaded hole in laminated composites

*

Corresponding author. Tel.: +1-520-621-6113; fax: +1-520-621-8191. E-mail address: [email protected] (E. Madenci).

0020-7683/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 6 8 3 ( 0 1 ) 0 0 1 2 8 - 7

7802

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

based on either ®nite element analysis or analytical methods. Since the contact stress distribution and the contact region are not known a priori, a majority of the models did not directly impose the boundary conditions appropriate for modeling the contact and non-contact regions between the bolt and the boundary of the hole. These models usually assumed a cosinusoidal bearing stress distribution or zero radial displacements over the contact region of the hole boundary. In the case of multi-bolt joints, the commonly accepted approach is to ®rst determine the load distribution among the bolts in order to identify the critical (most highly loaded) bolt for a subsequent single-bolt analysis for local stress distribution. However, this type of analysis disregards the interaction among the bolts located in close proximity to each other. A representative study with extensive literature survey and experimental investigation was conducted by Ramkumar et al. (1986). They developed a special ®nite element for a ``loaded hole'' to predict the bolt load distribution in bolted lap joint and to determine the most critical bolt for a subsequent detailed stress analysis. The bolt load distribution analysis is essentially one dimensional. Therefore, the bolt load distribution is always equal among the bolts in a row perpendicular to the direction of loading. In order to account for this discrepancy, the load distribution among the various rows was predicted by invoking experimentally obtained ``joint sti€ness'' values. Also, this analysis is limited to particular bolt patterns in order to conform to the loaded hole element geometry. In order to eliminate these shortcomings, Madenci et al. (1998) developed a method for single-lap joints based on the boundary collocation technique. Their method determines the contact stresses and contact region, as well as the bolt load distribution, as part of the solution procedure. However, this method fails to provide converged solutions consistently depending on the number of bolts and their location in relation to each other or to the free boundaries. A detailed validation and demonstration of their approach, as well as an extensive review of previous analyses, were reported in detail by Madenci et al. (1997). In the literature, there are essentially no direct analyses of double-lap bolted joints for solid laminates under general loading conditions and appropriate boundary conditions arising from contact phenomenon. Madenci et al. (1999) extended their boundary collocation technique for single-lap joints to consider double-lap joints and thermal loading. This method provided converged results for particular con®gurations, but also su€ered from consistent convergence arising from the explicit partitioning of the domain. Xiong and Poon (1994) introduced an analytical approach utilizing a variational formulation in conjunction with the complex potential theory to single- and double-lap joints with many bolts. Their approach considers each laminate of the joint separately. The coupling of the laminates is achieved through bolt displacements, which are permitted only in the direction of loading. In their two-stage analysis, the ®rst stage provides the local deformation along the hole boundaries of one of the laminates subjected to the external boundary conditions and the prescribed cosinusoidal bearing stress representing the bolt load at each hole boundary. The local deformations and the bolt de¯ections are imposed as displacement constraints in the subsequent second stage to determine the contact stresses (bolt loads) and the contact region in the second laminate. Subsequently, these fastener loads are imposed as prescribed cosinusoidal bearing stress for the ®rst stage of the analysis, and the iterative process continues between the ®rst and second stage analyses until the constraint conditions are satis®ed. This study presents an analysis method for determining the bolt load distribution in single- and doublelap joints while accounting for the contact phenomenon and the interaction among the bolts explicitly under bearing and by-pass loading with or without thermal loading. Although an extension of the analysis introduced by Xiong and Poon (1994), it eliminates the requirement of a two-stage analysis and the associated iterative process. Furthermore, it considers each laminate of the joint simultaneously by coupling them through bolt displacements in both the horizontal and vertical directions. Also, it captures the e€ect of the clearance between the bolt and hole. Without resorting to a two-stage analysis, the resulting equations are solved in a coupled manner, leading to the contact stresses, contact region, and bolt load distribution.

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7803

2. Problem statement The geometry of a bolted single- and double-lap joint with composite laminates is described in Fig. 1. The joint can be subjected to a combination of bearing, by-pass and shear loads, and a uniform temperature change. Each laminate of the single- and double-lap joints, joined with L number of bolts, is subjected  …k† segment of the external boundary of each region deto traction components, ta…k† …a ˆ x; y† along the C …k† noted by C . The section of the external boundary subjected to displacement constraints, u~…k† q …q ˆ n; s†, ~ …k† . The subscripts ``n'' and ``s'' denote directions outward normal and tangent to the is denoted by C boundary, respectively. Each region with an area of A…k† can be under uniform temperature change, T …k† . The thickness of the laminates (regions) is denoted by h…k† . The contact region between the `th bolt and the ^ …k`† . The sub or superscripts ``…k†'' and ``…`†'' refer to the hole boundary in the kth region is denoted by C regions (laminates) and bolts, respectively. Their ranges are speci®ed by k ˆ 1; . . . ; K and ` ˆ 1; . . . ; L, with K and L being the total number of regions (laminates) and bolts, respectively. As illustrated in Fig. 2, the hole radius, a` , which is slightly larger than the bolt radius, R` , leads to a clearance of d` . The hole and bolt radii remain the same for each region. As shown in Fig. 2, the center of each hole, located at …X` ; Y` †, coincides with the origin of the Cartesian coordinates …x` ; y` †. As shown in Fig. 3, the free-body diagram of each component of a lap joint, the unknown boundary …k† ~…k† traction components, k~…k† u~…k† q ˆ …uq q † along a q , arise from the deformation of the boundary given by c …k† ~ portion of the external boundary, C . The unknown traction component in the outward normal direction, k^…k`† n , arises from the deformation of the contact zone between the `th bolt and the hole boundary in the kth …`† region laminate. This contact zone deformation is expressed by c^…k`† ˆ u…k† u^…k`† d…k`† along the n n …Dj † …k`† …k`† ^ . The extent of the contact region denoted by C ^ contact region denoted by C is dependent on the bolt …`† …k`† …k† displacement, u^n …Dj †, deformation of the hole boundary, un , and the gap, d…k`† . Because of the absence of friction between the bolt and the laminate, the tangential component of the bolt displacement, u^…k`† s , and ^…k`† the traction vector, k^…k`† ˆ 0 and k^…k`† ˆ 0. As shown in Fig. 4, at the point of initial contact s s , vanish, i.e., u s (prior to any deformation of hole boundary), the gap between the hole boundary and the bolt (distance PP 0 ) in the kth region is de®ned by PP 0 ˆ d…k`† ˆ d` ‰1

cos…h

h †Š

…1†

^ …k`† , is in which h speci®es the line of action, m, and d` is the clearance. The extent of the contact region, C de®ned by the angles hA and hB . The ¯exible bolts experience de¯ections given by n o T …`† …`† …`† …`† D…`† ˆ Dx1 ; Dx2 ; Dy1 ; Dy2 for a single-lap joint …2a† 

n o T …`† …`† …`† …`† …`† …`† for a double-lap joint D…`† ˆ Dx1 ; Dx2 ; Dx3 ; Dy1 ; Dy2 ; Dy3 …`† Dxi

…2b†

…`† Dyi

and …i ˆ 1; K† representing the bolt de¯ection components at the ith point along the length of with the `th bolt along the x- and y-directions, respectively. The material properties of each laminate are represented by the matrix A…k† relating the stress resultants, …k† …k† Nab , to strain resultants, eab , with a, b ˆ x, y, in the form 8 …k† 9 2 3…k† 8 …k† 9 A11 A12 A16 < Nxx = < exx = Nyy…k† ˆ 4 A12 A22 A26 5 e…k† or N…k† ˆ A…k† e…k† …3† : …k† ; : yy…k† ; A16 A26 A66 Nxy 2exy …k†

where Aij are the components of the in-plane sti€ness matrix A…k† of the kth region. The strain components …k† arising from temperature change,  eab , are expressed as      …k†  …k† …k† …k† …k† exx ; eyy ; 2 e…k† …4† ˆ a…k† xy xx ; ayy ; axy T

7804

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 1. Geometric description of single- and double-lap joints with many bolts. …k† …k† where the coecients a…k† xx , ayy , and axy represent the thermal expansion coecients of the kth region with respect to the global …X ; Y † coordinate system. The corresponding thermal stress resultants are de®ned …k† as  Nab ,

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7805

Fig. 2. Position of a bolt before and after the load is exerted.

8  …k† 9 2 A11 < Nxx =  …k† Nyy ˆ 4 A12 :  …k† ; A16 Nxy

A12 A22 A26

3…k† 8 …k† 9 A16 < axx = a…k† T …k† A26 5 : yy ; A66 a…k† xy

…5†

In matrix form, this relationship is expressed as 

N…k† ˆ A…k† e…k†

…6† T

…k† …k† in which  e…k† ˆ a…k† T …k† , with a…k† ˆ fa…k† xx ; ayy ; axy g. …`† The sti€ness matrix of the bolts, b , is given by " # bx…`† 0 …`† b ˆ 0 b…`† y

…7†

whose coecients are determined by modeling the bolt as a beam under concentrated forces. The explicit expressions for bolt sti€ness for a single- and double-lap joint, as well as the general lap con®gurations, are derived in Appendix A. These angles, the contact stresses, the components of bolt displacement, and the forces exerted by the bolts are the unknowns to be determined as part of the solution. Unless indicated otherwise, the subscripts a and b vary as a, b ˆ x, y, representing the …x; y† global coordinates. The subscript q varies as q ˆ n, s, representing the directions normal and tangent to the boundary, as shown in Fig. 1. Also, only repeated subscripts imply summation.

3. Solution method The solution method is based on the variational formulation in conjunction with the complex potential theory. The governing equations are derived by requiring the ®rst variation of the total potential energy, arising from thermal and mechanical loads, to vanish. The in-plane equilibrium equations in each region are satis®ed exactly by employing complex potential functions in the form suggested by Lekhnitskii (1968).

7806

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 3. Free-body diagram of each component in a bolted single-lap joint.

However, each of the bolt equilibrium equations and the boundary conditions are satis®ed by minimizing the total potential energy. 3.1. Governing equations The total potential energy for K regions connected with L number of bolts in the absence of friction between the laminates and the bolts along the contact region under mechanical and thermal loading can be expressed as pˆ

K X kˆ1

U …k† ‡

L X `ˆ1

B…`† ‡

K X L X kˆ1

`ˆ1

W^ …k`† ‡

K X

W~ …k† ‡

kˆ1

The strain energy of the kth laminate, U …k† , is given by

K X kˆ1

W …k†

…8†

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7807

Fig. 4. The gap between the bolt and hole boundary immediately before and after the load is exerted.

U …k† ˆ

1 2

Z

Z

…k† …k†

A…k†

Nab eab dA

…k†

A…k†

…k†

Nab  eab dA

…9†

with its ®rst variation derived in the form dU

…k†

Z ˆ

A…k† …k†

in which  ta…k† ˆ …Nab by

…k† Nab;b du…k† a dA …k†

Z ‡

 …k† C

 …k† ta du…k† a dC

Z ‡

~ …k† C

 …k† tq du…k† q dC

‡

L…k† Z X `ˆ1

^ …k`† C

 …k† tq du…k† q dC

…10†

2 Nab †nb with a, b ˆ x, y and q ˆ n, s. The strain energy of the `th bolt, B…`† , is given

7808

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837 …`† …`†

…`†

B…`† ˆ 12Di bij Dj

with i; j ˆ 1; 2K

…11†

with its ®rst variation …`†

…`†

…`†

dB…`† ˆ bij Dj dDi

…12†

The potential of the reaction forces, k^…k`† and k~…k† n q …q ˆ n; s†, arising from the contact between the bolt and the hole boundary and the applied displacement constraints along the external boundary are denoted by W^ …k`† and W~ …k† , respectively. They are expressed in the form Z

W^ …k`† ˆ

^ …k`† C

n J^…k`† k^…k`† u…k† n n

…`†

u^…k`† n …Dj †

o d…k`† dC

…13†

in which J^…k`† ˆ



1 contact between the kth plate and `th bolt 0 no contact between the kth plate and `th bolt

and W~ …k† ˆ

Z ~ …k† C

n k~…k† u…k† q q

o u~…k† dC q

with q ˆ n; s

…14†

The potential of the externally applied tractions, ta…k† , denoted by W …k† is expressed as W …k† ˆ

Z  …k† C

ta…k† u…k† a dC

with a ˆ x; y

…15†

in which ta…k† and u…k† a represent the applied traction and displacement components in the x- and y-directions, respectively. Their ®rst variations are obtained as Z Z n o …`†  …k`† …k`† …k`† …k† …k`† …k† ^ ^ ^ dW J J^…k`† k^…k`† ˆ un u^n …Dj † d…k`† dkn dC ‡ n dun dC ^ …k`† ^ …k`† C C Z o^ u…k`† …`† n J^…k`† k^…k`† dDi dC …16† n …`† …k`† ^ oDi C and dW~ …k† ˆ

Z ~ …k† C

n dk~…k† u…k† q q

Z o u~…k† dC ‡ q

~ …k† C

…k† k~…k† q duq dC

…17†

and dW …k† ˆ

Z  …k† C

ta…k† du…k† a dC

The ®rst variation of the total potential energy can be obtained as

…18†

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837 K Z X

dp ˆ

A…k†

kˆ1

‡

K Z X

 …k† ta

K X L Z X kˆ1

‡

n

 …k† C

kˆ1

‡

…k†

Nab;b du…k† a dA ‡

^ …k`† C

`ˆ1

K X L Z X kˆ1

^ …k`† C

`ˆ1

L X `ˆ1

( …`†

K Z X

…`†

bij Dj

kˆ1

K X L Z o X ta…k† du…k† dC ‡ a kˆ1

`ˆ1

K Z X

…k`†

o^ un …`† J^…k`† k^…k`† dC dDi n …`† …k`† ^ oDi C 

^ …k`† C

n

 …k† tn

 ‡ J^…k`† k^…k`† du…k† n n dC

K Z o X …k† dC ‡ ‡ k~…k† du q q

 …k† ts du…k† s dC

‡

n J^…k`† u…k† n

o …`†  ^…k`† …D † ‡ d u^…k`† j n …k`† dkn dC

kˆ1

~ …k† C

 …k† tq

7809

)

kˆ1

~ …k† C

n

u…k† q

o u~…k† dk~…k† q dC q …19†

…`† …k† ^…k`† are arbitrary independent quantities and requiring the u…k† u…k† Noting that du…k† a , duq , d^ n , d^ s , dDi , and dkn ®rst variation of the total potential energy to vanish lead to the equilibrium equations for the laminates and bolts as …k†

Nab;b ˆ 0

on A…k†

…20†

K Z X

…`† …`† bij Dj

kˆ1

…k`†

o^ un J^…k`† k^…k`† dC ˆ 0; n …`† …k`† ^ oDi C

i; j ˆ 1; 2K

…21†

Their associated boundary conditions are obtained as n o  …k†  …k† ta…k† ˆ 0 on C ta 

 …k† tn

 …k† ts

 ‡ J^…k`† k^…k`† ˆ0 n

ˆ0

n

…22a†

^ …k`† on C

…22b†

^ …k`† on C

…22c†

 …k† tq

o ‡ k~…k† ˆ0 q

~ …k† on C

…22d†

n u…k† q

o u~…k† ˆ0 q

~ …k† on C

…22e†

n J^…k`† u…k† n

o …`†  ^ …k`† u^…k`† ˆ 0 on C n …Dj † ‡ d…k`†

…22f†

where a, b ˆ x, y; q ˆ n, s; k ˆ 1; . . . ; K; and ` ˆ 1; . . . ; L. 3.2. Total potential energy The strain energy expression given in Eq. (9) can be rewritten in terms of displacement components as Z Z 1 …k† …k† …k† …k†  …k† …k† U ˆ N u dA Nab ua;b dA …23† 2 A…k† ab a;b …k† A Its integration by parts yields Z h  1  …k† 1 …k† …k† Nab;b 2 Nab;b u…k† U …k† ˆ dA ‡ Nab a 2 2 A…k†

 i …k† dA 2 Nab u…k† a ;b

…24†

7810

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Under uniform temperature distribution and applying Gauss' theorem, it reduces to Z Z   1 1 …k† …k† …k† …k† Nab;b u…k† dA ‡ Nab 2 Nab nb u…k† U …k† ˆ a a dC 2 A…k† 2 C…k†

…25†

This expression can be further simpli®ed by invoking the stress resultants, Nab;b , that satisfy the equilibrium equations, Nab;b ˆ 0, given by Eq. (20), as Z   1 …k† …k† …k† U …k† ˆ Nab …26† 2 Nab nb u…k† a dC 2 C…k† The explicit expressions for the stress resultants and displacement components satisfying both the equilibrium equations and the compatibility conditions are given by Lekhnitskii (1968). Finally, the total potential energy expression given in Eq. (19) is reduced to a form, free of area integrals, as Z  K L N Z  n X X 1 1X …k† …k† …k† …`† …`† …`†  …k† pˆ Nab bij Di Dj ‡ u…k† 2 Nab nb ua dC ‡ k~…k† q q 2 C…k† 2 `ˆ1 ~ …k† C kˆ1 kˆ1 K X L Z K Z n o X X …`†  …k`† ^…k`† …k† …k`† ta…k† u…k† J kn un ‡ u^n …Dj † ‡ d…k`† dC a dC kˆ1

`ˆ1

^ …k`† C

kˆ1

o u~…k† dC q

 …k† C

…27†

with a, b ˆ x, y; q ˆ n, s; k ˆ 1; . . . ; K; and ` ˆ 1; . . . ; L. It can be rewritten in matrix form as pˆ

Z K L K Z X X
T 1 1X ~k…k†T ~c…k† dC N…k† 2 N…k† n…k† u…k† dC ‡ b…`† D…`† ‡ 2 C…k† 2 `ˆ1 ~ …k† C kˆ1 kˆ1 Z Z K X L K X X T T t…k† u…k† dC J^…k`† ^ ‡ k…k`† ^c…k`† dC kˆ1

`ˆ1

^ …k`† C

kˆ1

 …k† C

…28†

where ~c…k† ˆ u…k†

~u…k† ;

n o T …k† u…k† ˆ u…k† ; x ; uy

^c…k`† ˆ u…k†

^ u…k`† …D…`† †

n o T ~ ~…k† u…k† ˆ u~…k† ; n ;u s

d…k`† n o T ^ u…k† ˆ u^…k† n ;0

…29a† …29b†

n o T t…k† ˆ tx…k† ; ty…k† ;

n o ~…k†T ˆ k~…k† ; k~…k† ; k n s

n o ^…k`† ˆ k^…k`† ; 0 k n

…29c†

n o T …k† ~c…k† ˆ c~…k† ~ ; c ; n s

n o T ^c…k`† ˆ c^…k`† ; 0 ; n

n o  dT ˆ d ; 0 …k`† …k`†

…29d†

and the matrix of unit normals, n…k† , is expressed as 2 …k† 3 0 nx 5 n…k† ˆ 4 0 n…k† y …k† …k† ny nx

…29e†

As given by Eq. (3), the unknown bolt displacements are contained in vector D…`† , and matrix b…`† represents the bolt sti€ness.

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7811

The constraint condition given by Eq. (22e), …k† c~…k† q ˆ uq

u~…k† q ˆ 0

with q ˆ n; s

…30†

is rewritten in vector form as ~c…k† ˆ u0

…k†

~ u0

…k†

ˆ0

…31†

in which u0

…k†

ˆ T…k† u…k†

…32†

…k†

with T representing the transformation matrix between the (x, y) and (n, s) coordinate systems. The …k† displacement vectors u0 and u…k† are de®ned in the (n, s) and (x, y) coordinate systems, respectively,  …k†T …k† …33† u0 ˆ u…k† n ; us n o T …k† u…k† ˆ u…k† x ; uy

…34†

Substituting for u…k† from Eq. (B.19) into Eq. (32) yields u0

…k†

T

ˆ T…k† U…k† a…k†

…35†

leading to T

~c…k† ˆ T…k† U…k† a…k†

~ u0

…k†

…36†

~…k† , de®ned in Eq. (29c) can be assumed as The unknown traction vector, k ~…k† ˆ k

J X jˆ0

~ …k† Pj K j

…37†

~ …k† where the matrix Pj and vector of unknown coecients K j are de®ned as   n o T p 0 ~ …k† ~ nj ; K ~ sj Pj ˆ j ˆ K and K j 0 pj

…38†

with pj being the jth-order Legendre polynomial. In matrix notation, this equation becomes ~ …k† ~…k† ˆ PK k

…39†

where P ˆ ‰P0 P1 P2


PJ Š and

n o ~ …k†T ˆ K ~ T; K ~ T; K ~ T; . . . ; K ~T K 0 1 2 J

…40†

~…k†T ~c…k† that appears in the expression for the From Eqs. (36) and (40), the boundary integral of the product k total potential energy is obtained as Z Z Z ~ …k†T C ~ …k†T ~f …k† …41† ~…k†T ~c…k† dC ˆ ~ …k†T PT T…k† U…k†T a…k† dC ~ …k†T PT T…k† ~u…k† dC ˆ K ~ …k† a…k† K k K K ~ …k† C

where ~ …k† ˆ C and

~ …k† C

Z ~ …k† C

T

PT T…k† U…k† dC

~ …k† C

…42†

7812

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

~f …k† ˆ

Z ~ …k† C

PT T…k† ~ u…k† dC

…43†

The constraint condition given by Eq. (22f), …`†

 u^…k`† n …Dj † ‡ d…k`† ˆ 0

c^…k`† ˆ u…k† n

…44†

is rewritten in vector form as ^c…k`† ˆ u0

…k†

…k`†

^u0 …D…`† † ‡ d…k`† ˆ 0

The bolt displacement vector, ^ u the form ^ u0

…k`†

0…k`†

, at the contact region can be expressed in terms of the bolt de¯ections in

^ …k`† D…`† ˆG

…46† ^ …k`†

in which the matrix G ^ …k`† ˆ G

with g…zi g…zi and

…45†



sin hkl 0

…`†

and the vector D are de®ned as  ... g…z1 zk † g…z2 zk †  cos hkl 0 0 ... 0 0 0 g…z1 zk † g…z2 zk †

zk † de®ned as  1 if i ˆ k zk † ˆ 0 if i ˆ 6 k

g…zK

zk †

0 ... ...



0 g…zK

…47†

zk †

…i ˆ 1; . . . ; K†

n o T …`† …`† …`† …`† …`† …`† D…`† ˆ Dx1 Dx2


DxK Dy1 Dy2


DyK

…48†

With substitutions from Eqs. (35) and (46), the expression for ^c…k`† becomes T

^c…k`† ˆ T…k† U…k† a…k†

^ …k`† D…`† ‡ d G …k`†

…49†

^…k`† , de®ned in Eq. (29c) can be assumed as The unknown traction vector, k ^…k`† ˆ k

I X iˆ0

^i Fi K

…k`†

…50†

^ …k`† are de®ned as where the matrix Fi and vector of unknown coecients K i   n o T F 0 ^ …k`† ; 0 ^ …k`† Fi ˆ i ˆ K and K i ni 0 Fi

…51†

in which Fi is an orthogonal trigonometric function satisfying the condition of zero stress at the endpoints, i.e., Fi …s ˆ s0 † ˆ Fi …s ˆ s1 † ˆ 0

…52†

These functions are formulated as 

ip…h…s† h…s0 †† Fi …s† ˆ sin …h…s1 † h…s0 ††

 with i ˆ 1; I

…53†

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7813

in which h…s0 † and h…s1 † represent the beginning and end angles of the contact region. These angles are measured with respect to the global X-axis, and any point in the contact region is identi®ed by the angle h…s†. The expression for the unknown traction vector at the contact region can be rewritten as ^ …k`† ^…k`† ˆ FK k

…54†

where F ˆ ‰F0 F1 F2


FI Š

and

n o ^ …k`†T ˆ K ^ …k`† ; K ^ …k`† ; . . . ; K ^ …k`† K I 0 1

…55† T

^…k† ^c…k† that appears in the expression for From Eqs. (49) and (54), the boundary integral of the product k the total potential energy is obtained as Z Z Z Z ^…k`†T ^c…k`† dC ^ …k`†T FT T…k† U…k†T a…k† dC ^ …k`†T FT G ^ …k`† D…`† dC ‡ ^ …k`†T FT d dC ^ˆ k K K K ^ …k`† C

^ …k`† C

^ …k`† C

^ …k`† C

…k`†

…56† or

Z ^ …k`† C

^ …k`† C ^…k`†T ^c…k`† dC ^ …k`† a…k† ^ ˆK k

where ^ …k`† ˆ C …k`† ^ g1

Z Z

…k`† ^ g0 ˆ

…57†

T

^ …k`† C

ˆ

…k`† ^ …k`†T ^ ^ …k`†T ^g…k`† g1 D…`† ‡ K K 0

^ …k`† C

Z ^ …k`† C

FT T…k† U…k† dC

…58†

^ …k`† dC FT G

…59†

FT d…k`† dC

…60†

 …k† is expressed as Also, the potential energy of the traction vector, t…k† , acting on C Z Z Z T T T T t…k† u…k† dC ˆ u…k† t…k† dC ˆ a…k† U…k†t…k† dC ˆ a…k† f  …k† C

where f ˆ

 …k† C

Z  …k† C

…61†

 …k† C

U…k†t…k† dC

…62†

Substituting from Eqs. (B.19), (41), (57) and (61) for the appropriate terms in Eq. (28), the expression for the total potential energy is obtained in matrix form as pˆ

K T 1X a…k† H…k† a…k† 2 kˆ1

‡

K X L X kˆ1

in which

`ˆ1

K X

 …k† …k†

h a

‡

kˆ1

^ …k`† C ^ …k`† a…k† K

K X L X kˆ1 `ˆ1

L K X T 1X ~ …k†T C ~ …k† a…k† D…`† b…`† D…`† ‡ K 2 `ˆ1 kˆ1

…k`† ^ …k`†T ^ g1 D…`† ‡ K

K X L X kˆ1

`ˆ1

^ …k`†T ^g…k`† K 0

K X kˆ1

K X

~ …k†T ~f …k† K

kˆ1

f T a…k†

…63†

7814

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

H…k† ˆ  …k†

h

Z

T

C…k†

Z ˆ

S…k† n…k† U…k† dC 

C…k†

T

…64†

T

N…k† n…k† U…k† dC

…65†

where S…k† and U…k† are given by Eqs. (B.19) and (B.20). The total potential energy in Eq. (63) is compacted to its ®nal form p ˆ 12aT Ha



~ T Ca ~ ha ‡ 12DT bD ‡ K

~ T~f ‡ K ^ Ca ^ K

^ T ^g ^ T ^g D ‡ K K 1 0

fa

in which the vectors and matrices are de®ned by n T o T T aT ˆ a…1† ; a…2† ; . . . ; a…K†

…67a†

n o T T T DT ˆ D…1† ; D…2† ; . . . ; D…L† 

~T ˆ K



^T ˆ K

~ …1†T ; K ~ …2†T ; . . . ; K ~ …K†T K ^ …1†T ; K ^ …2†T ; . . . ; K ^ …K†T K

…67b†

 …67c† 

^ …k†T ˆ with K



^ …k1†T ; K ^ …k2†T ; . . . ; K ^ …kL†T K

  f T ˆ f …1†T ; f …2†T ; . . . ; f …K†T  T

h ˆ

n

 …1†T  …2†T

h

; h

; . . . ;  h…K†

 ˆ

…1†T …2†T …K†T ^ g0 ; ^ g0 ; . . . ; ^ g0

 …67d† …67e†

T

o

…67f†

  ~f T ˆ ~f …1†T ; ~f …2†T ; . . . ; ~f …K†T ^ gT0

…66†

…67g†

 with

…k†T ^ g0

 ˆ

T …k2†T …kL†T ^g…k1† ; ^g0 ; . . . ; ^g0 0

 …67h†

and 2 6 6 Hˆ6 4 2 6 6 bˆ6 4

H…1†

b…1†

3 H…2†

..

.

7 7 7 5

…68a†

H…K† 3

b…2†

..

.

7 7 7 5 b…L†

…68b†

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

2 6 ~ ˆ6 C 6 4 2 6 ^ ˆ6 C 6 4

~ …1† C

^ …1† C

~ …2† C

..

.

~ …K† C

7 7 7 5

…68c†

3 ^ …2† C

..

.

^ …K† C

7 7 7 5

with

h ^ …k†T ˆ C ^ …k1†T C

and

2 ^ gT1

ˆ

h

…1†T ^ g1

7815

3

T ^g…2† 1






…K†T ^ g1

i

with

…k† ^ g1

6 6 ˆ6 6 4

^g…k1† 1

^ …k2†T C






^ …kL†T C

i

…68d†

3 ^g…k2† 1

..

.

^g…kL† 1

7 7 7 7 5

…68e†

The minimization of the total potential energy, i.e., dp ˆ 0, leads to the following equilibrium equations 9 2 38 9 8 ~T C ^ T > a > > f ‡  h > H 0 C > > > > = < = < 60 b 0 ^ gT1 7 0 6 7 D ˆ …69† ~ > > ~f > 4C ~ 0 K 0 0 5> > > > > ; ; : : ^ ^ ^ K ^ g0 0 C g1 0 Solving Eq. (69) for the unknown vectors permits the calculation of the stress and displacement components in each laminate, bolt de¯ections, and forces exerted by the bolts. Because the angles hA and hB de®ning the contact region are unknown, these equations are solved for their assumed values until convergence is achieved through an iterative scheme. The iterative scheme for solving the system of algebraic …0† …0† equations given in Eq. (69) begins with the initial estimates of hA and hB , shown in Fig. 5, de®ning the contact region. These angles are measured in the counterclockwise direction from the x…`† axis of the local coordinate system, …x…`† ; y …`† †. The initial estimate in most cases does not represent the true contact region, ^ …k`† , and the fact that the start and the end for which the radial stresses are all compressive r…k† rr < 0 on C …k† points of the contact region have zero radial stresses, rrr …hA † ˆ 0 and r…k† rr …hB † ˆ 0. As shown in Fig. 5, three distinct cases exist, depending on whether the radial stress near the start angle, …0† hA , is larger or smaller than the true value of the start angle or equal to the true value, hA , of the contact ^ …k`† . The initial guess of the starting angle, h…0† region, C A , is smaller than its true value, hA , if the compressive …0† radial stresses change sign and become tensile near the start point of the contact region, hA . The initial …0† guess of the starting angle, hA , is larger than its true value, hA , if the compressive radial stresses do not change sign and remain compressive near the start point of the contact region. The initial estimate of the …0† end angle, hB , is larger than its true value, hB , if the compressive stresses change sign and become tensile …0† near the end angle. If the compressive stresses do not change sign, then the initial guess, hB , is smaller than its true value. …0† …0† During the iteration process, an initially guessed contact region de®ned by hA and hB converges to the true contact region de®ned by hA and hB through incrementally changing the values of the start and end angles. The increment is forced to decrease each time the true value of the start or end point is passed, and the direction of angle change is altered. The convergence of the iterative process is achieved when the incremental value of the angle reaches a pre-de®ned value.

7816

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 5. The behavior of radial stress near the start point of a contact region.

In order to avoid the case of radial stresses having zero values at the start or end points of the contact region but with tensile and compressive stresses along the contact region, two auxiliary points are considered as shown in Fig. 6. These points are located inside the contact region near the start and end points of the contact region. In order to achieve convergence for the contact region, the radial stresses at these two auxiliary points must be compressive.

4. Numerical results Three di€erent types of load transfer through bolted joints are considered in order to validate the present analysis. The ®rst con®guration is a pin-loaded square plate considered by Ireman et al. (1993). The second con®guration is a single-lap joint with four bolts investigated by Xiong and Poon (1994). The third con®guration is a joint of dissimilar bars with a single bolt under di€erent uniform temperature changes analyzed previously by Gatewood (1957). Then, the capability of the present analysis is demonstrated by considering two di€erent double-lap joint con®gurations: one with three bolts under mechanical loading only, uniform temperature change only, and their combination, and the other with seven bolts under mechanical loading only.

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7817

Fig. 6. Auxiliary points along a contact region.

Fig. 7. Pin-loaded plate con®guration.

A pin-loaded square plate was investigated through detailed ®nite element analysis by Ireman et al. (1993). As shown in Fig. 7, the geometry of the plate is de®ned by the parameters W ˆ 24 mm, D ˆ 6 mm, d ˆ 0:021 mm, and thickness t ˆ 3:046 mm. Three di€erent sets of laminate material properties considered in the analysis are given in Table 1. The bolt was assumed to be rigid and the applied load, P, was taken as 5483 N. In the present analysis, the number of terms, N, retained in the series representation of the complex …k`† …k`† …k`† potential functions, U…k`† r …zr † and ur …zr †, is taken as 80. The number of terms retained in the series ^…k`† , denoted by J and I, respectively, is taken as 24. Figs. ~…k† and k representation of the reaction tractions, k 8±10 shows the favorable comparison of the normalized radial and tangential stresses obtained through the ®nite element analysis reported by Ireman et al. (1993) with those of the present calculations. The slight Table 1 Material properties of laminates Case

Ex (GPa)

Ey (GPa)

Gxy (GPa)

mxy

A B C

99.2 35.5 51.5

35.5 99.2 51.5

8.5 8.5 19.3

0.24 0.08 0.33

7818

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 8. Comparison of normalized stress around hole for pin-loaded plate for laminate type A.

oscillations appear due to the trigonometric series representation of unknown tractions along the contact region. The radial and tangential stress components are normalized with respect to the applied bearing stress, P =…dt†, and the applied stress, P =…Wt†, respectively. The geometry and loading conditions of a single-lap joint of aluminum and composite plates with four bolts are shown in Fig. 11. As shown in this ®gure, the geometry of the joint is de®ned by the dimensions of h…1† ˆ 0:31 in., h…2† ˆ 0:117 in., D ˆ 0:3125 in., W ˆ 3:125 in., s ˆ 1:25 in., e ˆ 0:9375 in., ` ˆ 2:75 in. There exists no bolt-hole clearance, d ˆ 0. The composite lay-up is ‰…45°=0°= 45°=0°†2 =0°=90°ŠS with lamina

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7819

Fig. 9. Comparison of normalized stress around hole for pin-loaded plate for laminate type B.

properties of EL ˆ 18:5  106 psi, ET ˆ 1:9  106 psi, GLT ˆ 0:85  106 psi, and mLT ˆ 0:3. In obtaining the present predictions, the series representations of the functions are truncated at N ˆ 30 and J ˆ I ˆ 5. Although the entire geometry of the lap joint is considered in the present analysis, only the predictions for the bolt loads exerted by Bolt 1 and Bolt 3 are presented due to the presence of symmetry. As compared in Table 2, the present analysis predictions are in remarkable agreement with those calculated by Xiong and Poon (1994).

7820

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 10. Comparison of normalized stress around hole for pin-loaded plate for laminate type C.

The geometry of the steel and aluminum plates connected with a single bolt is shown in Fig. 12. Each plate is subjected to a di€erent uniform temperature change. The bolt is assumed to be rigid, with no clearance between the bolt and the bolt holes. The material properties and the temperature change for each bar is presented in Table 3. In the present analysis, the bars are assumed as narrow plates whose geometry is speci®ed by W ˆ 0:5 in., e ˆ 3:5 in., s ˆ 0:5 in., D ˆ 0:2 in., and h…1† ˆ h…2† ˆ 0:06 in. The present analysis predictions are obtained by considering N ˆ 30 and J ˆ I ˆ 5 in the series representations of the functions. The present analysis predicts the force exerted by the bolt to be 442 lb. The bolt force obtained by Gatewood (1957) according to the strength of materials approach is 461 lb. The di€erence of 4.1% in bolt

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7821

Fig. 11. Single-lap joint with four bolts.

Table 2 Bolt load distribution P1 =P P3 =P

Present analysis

Xiong and Poon (1994)

0.23 0.27

0.23 0.27

load prediction is possibly due to the fact that the strength of materials approach does not account for the presence of the bolt hole. The geometry of a double-lap joint with three bolts is shown in Fig. 13. The geometry of the plates all having the same dimensions are described by W ˆ 50 mm, s ˆ 25 mm, e ˆ 12:5 mm, D ˆ 6 mm, h ˆ 5 mm, h…1† ˆ h…2† ˆ 4 mm, and d ˆ 0:005 mm. The material properties for the steel and aluminum plates, respectively, are: ES ˆ 200 GPa, mS ˆ 0:3, aS ˆ 11:7  10 6 …°F† 1 , EA ˆ 70 GPa, mA ˆ 0:3, and aA ˆ 23:0  10 6 …°F† 1 . All of the bolts in this joint are assumed to be rigid. For the ®rst case of mechanical loading only, the applied tensile stress, r , is 200 N/mm. In the second loading case, the plates are under uniform temperature change only. Accordingly, the steel plates are free from temperature change, DTS ˆ 0°C, while the aluminum plate is subject to DTA ˆ 125°C. The third case represents a combination of the mechanical and thermal loads, with r ˆ 200 N/mm, DTS ˆ 0°C, and DTA ˆ 125°C. The number of terms retained in the series representation of the functions in the present analysis is speci®ed by N ˆ 30 and J ˆ I ˆ 5. Due to the presence of symmetry, only the results concerning Bolt 2 and Bolt 3 are presented in Figs. 14± 16 and Table 4. The tangential and radial stresses around the bolt holes, shown in Figs. 14±16, correspond to mechanical loading, uniform temperature change, and their combination, respectively. As observed in these ®gures, the location and extent of the contact region, as well as the magnitudes of the bearing stresses, change signi®cantly. The bolt load distributions for each of these loading cases are given in Table 4. As expected, the force exerted by Bolt 2 alters its direction when the nature of the loading changes from mechanical to thermal. This behavior is caused by the thermal expansion of the aluminum plate constrained

7822

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 12. Steel and aluminum plates connected with a single bolt under thermal loading.

Table 3 Material properties and temperature changes Material Plate 1 Plate 2

Steel Aluminum

E (psi) 7

3  10 107

m

a ……°F † 1 †

DT (°F)

0.3 0.3

6:5  10 6 12  10 6

500 100

with three bolts. Consequently, under thermal loading, the bolt forces with their directions toward the middle of the plate resist the thermal expansion of the plate. Under combined mechanical and thermal loading, the forces exerted by Bolt 2 and Bolt 3 are approximately equal to the sum of the forces acting on these bolts in the cases of mechanical and thermal loads. This shows the e€ect of the non-linearity arising from the contact analysis. A complex double-lap joint with seven bolts, shown in Fig. 17, is subjected to a tensile loading. The outer plates are steel and the inner plate is a composite. The bolt-hole clearance is speci®ed as 1% of the hole diameter. The composite plate is made of graphite and ®berglass with material properties EL ˆ 4:7  106 psi, ET ˆ 4:75  106 psi, GLT ˆ 1:2  106 psi, and mLT ˆ 0:24. Steel properties for the plates and bolts are taken as E ˆ 30  106 psi and m ˆ 0:3. The applied load is 70,000 lbs, corresponding to r ˆ 10,000 psi. In the calculation of the results, the number of terms retained in the series representation of the functions in the present analysis is speci®ed by N ˆ 30 and J ˆ I ˆ 7. Because of the presence of symmetry, only the stress distributions around Bolts 1, 2, 4, and 5 are shown in Fig. 18. The forces exerted by these bolts are given in Table 5. As expected, the bolts in the ®rst row share more of the load than the bolts in the second row. Also, the outer bolts share more of the load than the inner bolts. Accordingly, Bolt 1 and Bolt 3 are the most highly loaded bolts for this con®guration.

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7823

Fig. 13. Steel plates bolted to an aluminum plate with three bolts.

5. Conclusions In this study, a new approach based on a complex potential theory in conjunction with a variational formulation has been introduced for the thermo-elastic contact analysis of a general bolted-joint con®guration containing multiple laminates joined by multiple bolts. The total potential energy of the joint is formulated by using a solution in the form of a complex potential series that automatically satis®es the stress equilibrium equations and compatibility conditions, thus avoiding the necessity to perform area integrations and resulting in boundary integral expressions for the strain energy of the laminates. The total potential energy also includes the strain energy of the bolts based on a shear deformable beam theory. In order to capture high gradient variations of stresses near the free or bolted holes, the stress ®eld is de®ned as the superposition of complex potential series originating from each hole. Hence, not only are continuous stress and displacement ®elds obtained but the modeling of the entire joint is also simpli®ed considerably. By only entering boundary information, hole size and locations, and the number of terms to be used in complex and other series, the solution provides all the stress, displacement, and contact force distributions at any point in the joint. Contact between the bolts and the laminates is established by enforcing displacement continuity along the contact region between the bolts and the plates. This is established by incorporating the work done by the unknown contact forces over the contact displacements into the total potential energy expression. The contact displacements are de®ned by constraint equations that take into account the gap between the bolts and plates. The contact forces are assumed in the form of trigonometric series that satisfy stress-free conditions at the ends of the contact regions. Since the contact regions are unknown a priori, an iterative scheme is adopted in order to determine the beginning and end angles of the contact regions. Starting with an initial guess, the system matrix is generated to solve for unknown plate and bolt displacements and contact forces simultaneously. The simultaneous determination of bolt displacements and contact forces, along with the plate displacements, is

7824

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 14. Stresses around bolt holes in aluminum plate due to mechanical loading only.

a unique feature of the present formulation. A new guess is then obtained by monitoring the stress distribution along the contact regions. The iterative scheme is continued until a con®guration for contact regions is reached where all the contact forces become compressive. The validation problems show excellent agreement of the present formulation against those reported by other investigators. The pin-loaded panel in the ®rst validation problem provides a comparison of contact angles and contact force distribution. For all laminate con®gurations, remarkable agreement is obtained between the present analysis and the re®ned ®nite element solution. In the case of a single-lap con®guration subjected to thermal loading, the strength of material solution is available. As expected, the present analysis achieves the right bolt load and compares well with the strength of material solution if full contact is assumed around the bolt. In the case of a single- and a double-lap joint containing multiple holes, the present analysis captures the correct load distribution shared by each bolt. The versatility of the present formulation has been demonstrated by solving a double-lap joint con®guration containing three bolts and two laminates with di€erent material properties, and the joint is subjected to thermal, mechanical, and thermo-mechanical loadings. All of these cases are solved assuming variable contact regions. Therefore, the rule of superposition is invalid for these problems since the contact regions are changing. This can be clearly observed in the tabulated results where the summation of the bolt load distribution arising from thermal and mechanical loads is signi®cantly di€erent from those arising from the combined loading.

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7825

Fig. 15. Stresses around bolt holes in aluminum plate due to thermal loading only.

Appendix A The bolt sti€ness matrix is derived based on the Timoshenko's zeroth-order shear deformable beam theory. The cross-section of a bolt connecting K laminates (regions) is shown in Fig. 19. The bolt number is denoted by ` and the regions are numbered from bottom to top in sequential order. The bolt has a circular uniform cross-section, A` , moment of inertia, I` , and Young's and shear moduli, E` and G` , respectively. The nuts at the ends of the bolt are assumed to represent clamped boundary conditions, thus preventing rotations but creating reaction moments at the ends of the bolt. The bolt is subjected to forces arising from the contact between the bolt and the laminates. Because of the variation in laminate thickness and sti€ness, these forces exerted by the laminates vary through the length of the bolt. Because of the variable contact forces along the length of the bolt, the large in-plane bolt sti€ness compared to those of the laminates and the small ratio of the bolt diameter to its length, the most suitable and accurate representation of the bolt can be achieved by discretizing the bolt into small Timoshenko beam elements connected at nodal points, as shown in Fig. 19. The bolt discretization is based on the most e€ective locations of the contact forces. Thus, two of the nodes are selected at the top and bottom ends of the bolt in order to obtain the largest in-plane bolt de¯ection. The contact forces exerted on the bolt by the top and bottom laminates are assigned to these end nodal points. The remaining nodal points are chosen at the intersections of mid-planes of the inner laminates and the bolt longitudinal axis, as shown in Fig. 19. Hence, the contact forces exerted on the bolt by the inner laminates are assigned to these intermediate nodal points. The nodal de¯ections and rotations

7826

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 16. Stresses around bolt holes in aluminum plate due to combined thermal-mechanical loading.

Table 4 Bolt load distribution in a double-lap joint of steel plates and an aluminum plate with three bolts Bolt load (N) Mechanical loading Bolt # 2

Fx Fy

Bolt # 3

Fx Fy

2313:99 54:74 0

5374:18

Thermal loading

Combined loading

6906.94 6597.65

4383.20 6372.89

0

13 815:46

0

18 769:19

along the length of the bolt permit the determination of the bolt de¯ection at any point along the bolt by utilizing interpolation functions. Because the bolt material is homogeneous and isotropic and its cross-section is circular, the moment of inertia of the bolt is homogeneous on the …x; y† plane. Hence, there exists no coupling between the deformations of the bolt on the …x; z† and …y; z† planes, leading to the uncoupled sti€ness matrices of the bolt associated with the …x; z† and …y; z† planes. However, their forms are identical since the sti€ness properties of the bolt are homogeneous on the …x; y† plane. Therefore, the derivation of the sti€ness matrix associated with the …x; z† plane applies to the sti€ness matrix associated with the …y; z† plane.

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7827

Fig. 17. Steel plates bolted to a composite laminate with seven bolts.

As shown in Fig. 19, the bolt connecting K laminates is modeled by …K 1† number of beam elements …`† …`† …`† and K nodes. Each node is assigned a de¯ection, Dxi ˆ D…`† x …zi †, and a rotation, /xi ˆ /x …zi †, with the subscript i ˆ 1, K representing the node numbers. The positive directions of the de¯ections and rotations are shown in Fig. 19. Also, the length of each beam element is denoted by Li , with i ˆ 1, K. Based on the geometry and material properties of the `th bolt, the strain energy arising from the bolt deformation associated with the …x; z† plane can be written as Ux…`† ˆ

K 1 X iˆ1

Ux…i`†

…A:1†

The strain energy of each element is expressed as Z Z
2
1 1 …i`† 2 Ux…i`† ˆ E` I` j…i`† dz ‡ G` Af ` czx dz zz 2 Li 2 Li

…A:2†

…i`† in which the strains j…i`† zz and czx are based on Timoshenko's zeroth-order shear deformation theory, and they are de®ned as

j…i`† zz ˆ

d2 D…i`† x dz2

and

c…i`† zx ˆ

dD…i`† x dz

/…i`† x

…A:3†

In Eq. (A.2), Af ` denotes the corrected area of the bolt and is de®ned as Af ` ˆ c 2 A`

…A:4†

7828

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

Fig. 18. Stresses around bolt holes in composite plate with seven bolts.

Table 5 Bolt load distribution in a double-lap joint of steel plates and a composite laminate with seven bolts Bolt

Fx (lb)

1 2 4 5

1181:5

0 2.5 238.3

Fy (lb)

Fy =Fapplied

11983.7 11340.9 9947.8 7427.6

0.1712 0.1620 0.1421 0.1061

in which c2 represents the shear correction factor, which is a correction to the shear strain energy due to uniform transverse shear deformations. The displacement and rotation ®eld of the bolt is represented by piecewise continuous interpolation functions. These functions are de®ned individually over each element as …`†

…`†

D…i`† x …z† ˆ H1 …s†Dxi ‡ H2 …s†Dx…i‡1† ‡ H3 …s† …`†

…`†

dDx…i‡1† dDxi ‡ H4 …s† dz dz

…`† …`† /…i`† x …z† ˆ N1 …s†/xi ‡ N2 …s†/x…i‡1† ‡ N3 …s†/xm

…`†

…A:5a† …A:5b†

where the subscript m denotes the mid-point of the element. The variable s ˆ z zi is a local coordinate system for the ith element. The interpolation functions Hj …s† …j ˆ 1; . . . ; 4† represent the cubic Hermitian polynomials de®ned as

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7829

Fig. 19. A close view of the section in the vicinity of a bolt, and the discretization of the bolt into Timoshenko beam elements.

3s2 2s3 ‡ 3 L2i Li

H1 …s† ˆ 1 H2 …s† ˆ

3s2 L2i

H3 …s† ˆ s H4 …s† ˆ

…A:6a†

2s3 L3i

…A:6b†

2s2 s3 ‡ 2 Li Li

…A:6c†

s2 s3 ‡ Li L2i

…A:6d†

Also, the functions Nj …s† …j ˆ 1; . . . ; 3† represent the Lagrangian shape functions de®ned as   N1 …s† ˆ …s Li =2†…s Li †= L2i =2 N2 …s† ˆ s…s

  Li †= L2i =2

N3 …s† ˆ s…Li

  s†= L2i =4 D…i`† x

…A:7a† …A:7b† …A:7c†

/x…i`†

and in terms of the nodal unknowns de®ned at the end points of the element, In order to express two successive steps are performed. In the ®rst step, constraint of uniform shearing strain along the beam element is enforced into the kinematic ®eld, thus leading to

7830

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

dc…i`† d2 Dx…i`† zx ˆ dz dz2

d/x…i`† ˆ0 dz

…A:8†

Substituting from Eqs. (A.3) and (A.5) into Eq. (A.8) and grouping the terms as coecients of z0 and z, this constraint equation produces two algebraic equations of the form …`†

…`† 4 dDxi 2 dDxi‡1 ‡ ˆ Li dz Li dz …`†

…`† 6 dDxi 6 dDxi‡1 ‡ 2 ˆ L2i dz Li dz

6 …`† 6 …`† 3 1 …`† D ‡ 2 Dx…i‡1† ‡ /…`† / xi ‡ 2 xi Li Li Li Li x…i‡1†

4 …`† / Li xm

12 …`† 12 …`† 4 4 …`† D ‡ 3 Dx…i‡1† ‡ 2 /xi…`† ‡ 2 /x…i‡1† Li Li L3i xi Li …`†

8 …`† / L2i xm

…A:9a†

…A:9b†

…`†

These two equations are then solved for dDxi =dz and dDx…i‡1† =dz as …`†

dDxi ˆ dz

1 …`† 1 …`† 5 Dxi ‡ Dx…i‡1† ‡ /…`† Li Li 6 xi

…`†

dDx…i‡1† dz

1 …`† 1 …`† D ‡ Dx…i‡1† Li xi Li

ˆ

1 …`† / 6 x…i‡1†

1 …`† 5 …`† / ‡ /x…i‡1† 6 xi 6

2 …`† / 3 xm

…A:10a†

2 …`† / 3 xm

…A:10b†

Substituting from Eqs. (A.10a) and (A.10b) into Eq. (A.5a) and rearranging the terms, the in-plane de¯ection component, D…i`† x , is obtained as       s s …`† 5s 3s2 2s3 s s2 2s3 …`† …`† …`† D…i`† ˆ 1 ‡ D ‡ ‡ ‡ ‡ D / /x…i‡1† xi x xi Li Li x…i‡1† 6 2Li 3L2i 6 2Li 3L2i   2s 2s2 4s3 ‡ ‡ …A:11† /…`† xm 6 Li 3L2i Thus, the total number of unknowns in the expressions for the de¯ection and rotation is reduced from 7 to 5. Substituting from Eqs. (A.11) and (A.5) for Dx…i`† and /…i`† x in the strain de®nitions, Eq. (A.3), and carrying out the integrations in Eq. (A.2), analytically result in the following strain energy expression for the ith element of the `th bolt in matrix form ( )T " #( ) …i`† …i`† …i`† 1 q…i`† b11 b12 qx1 …i`† x1 Ux ˆ …A:12† …i`† …i`†T …i`† 2 q…i`† qx2 b12 b22 x2 where …i`†T

qx1

ˆ

n

…`†

Dxi

…`†

Dx…i‡1†

/xi…`†

…`†

/x…i‡1†

o

…A:13a†

…i`†

…`† qx2 ˆ /xm

2 6 6 …i`† b11 ˆ 6 6 4

…i`†T b12

 ˆ

G` Af ` Li G` Af ` Li G` Af ` Li G` Af ` Li

2 G` Af ` 3

…A:13b† G` Af ` Li G` Af ` Li G` Af ` Li G` Af ` Li

G` Af ` 6 G` Af ` 6 G` Af ` Li 7E` I` ‡ 3Li 36 G` Af ` Li E` I` ‡ 3Li 36

2 G` Af ` 3

G` Af ` 6 G` Af ` 6 G` Af ` Li E` I` ‡ 3Li 36 G` Af ` Li 7E` I` ‡ 3Li 36

8E` I` G` Af ` Li ‡ 3Li 9

3 7 7 7 7 5

8E` I` G` Af ` Li ‡ 3Li 9

…A:13c†

 …A:13d†

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837 …i`†

b22 ˆ

16E` I` 4G` Af ` Li ‡ 3Li 9

7831

…A:13e†

A further reduction of the total number of unknowns from 5 to 4 is achieved in the next step. The slope de®ned at the mid-point of the beam, /…`† xm , can be condensed out of the strain energy expression, Eq. (A.12), by employing the static condensation procedure. In the absence of nodal forces at the mid-point of the …i`† element, the ®rst variation of the strain energy with respect to qx2 yields the equilibrium equation at the mid-point as …i`†T …i`†

…i`† …i`†

b12 qx1 ‡ b22 qx2 ˆ 0 …i`† qx2

(i.e., Solving for pression reduces to

/…`† xm )

…A:14†

in the above equation and substituting into Eq. (A.12), the strain energy ex-

T

b…i`† q…i`† Ux…i`† ˆ 12q…i`† x x in which q…i`† is identical x 2 …i`† …i`† b11 b12 6 …i`† …i`† b22 6b b…i`† ˆ 6 12 …i`† 4 b13 b…i`† 23 …i`† …i`† b14 b24 or

2 6 6 b…i`† ˆ 6 4

…i`†

12b1 …i`† 12b1 …i`† 6b1 Li …i`† 6b1 Li

…A:15† …i`†

to qx1 , and b…i`† is de®ned as 3 …i`† …i`† b13 b14 …i`† …i`† 7 b23 b24 7 …i`† …i`† 7 b33 b34 5 …i`† …i`† b34 b44 …i`†

12b1 …i`† 12b1 …i`† 6b1 Li …i`† 6b1 Li

…i`†

6b1 Li …i`† 6b1 Li …i`† 4b2 …i`† 2b3

…A:16†

3 …i`† 6b1 Li 7 …i`† 6b1 Li 7 …i`† 7 2b3 5 …i`†

4b2

where Af ` E` G` I`
Li 12E` I` ‡ Af ` G` L2i

…A:17a†


E` I` 3E` I` ‡ Af ` G` L2i
ˆ Li 12E` I` ‡ Af ` G` L2i

…A:17b†


E` I` 6E` I` Af ` G` L2i
ˆ Li 12E` I` ‡ Af ` G` L2i

…A:17c†

…i`†

b1 ˆ …i`† b2

…i`† b3

Substituting from Eq. (15) into Eq. (A.1), the strain energy of the `th bolt is expressed as Ux…`† ˆ

K X 1 …i`†T …i`† …i`† q b qx 2 x iˆ1

…A:18†

As mentioned previously, the presence of nuts prevents the bolt from rotating at the end points; thus, it …`† …`† resembles clamped-type boundary conditions, requiring that /x1 ˆ 0 and /xK ˆ 0. Invoking this condition into Eq. (A.18), the strain energy expression is modi®ed as K 2 T T 1 1X 1 0…1`† 0…1`† Ux…`† ˆ q0…1`† b q ‡ q…i`† b…i`† qx…i`† ‡ q0……K x x x 2 2 iˆ2 2 x

1†`†T 0……K 1†`† 0……K 1†`† b qx

…A:19†

7832

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

where the vectors and matrices with a prime are de®ned as n o T …`† …`† …`† q0…1`† ˆ D ; D ; / x1 x2 x2 x q0……K x

1†`†T

n …`† ˆ Dx…K

…`† …`† 1† ; DxK ; /x…K 1†

…A:20a†

o

…A:20b†

…A:20c† and

…A:20d†

The terms on the right-hand side of Eq. (A.19) can be rearranged such that all the unknown de¯ection and rotation components and the coecients of the element sti€ness matrices are assembled in large vectors and matrices. Hence, the strain energy expression is compacted to ( )T " …`† #( ) …`† …`† …`† b b 1 q q DD D/ xD xD Ux…`† ˆ …A:21† …`† …`† …`†T …`† 2 qx/ qx/ bD/ b// where

n o …`†T …`† …`† …`† qxD ˆ Dx1 Dx2


DxK t;

n …`†T …`† …`† …`† qx/ ˆ /x2 /x3


/x…K

o 1†

…A:22a†

…A:22b†

…A:22c†

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7833

…A:22d†

The bolt is subjected to contact forces that are assigned to the nodal points and the reaction moments due to clamped support at the end points of the bolt. None of the internal nodal points, where the nodal rotations are active, is subjected to external moments. Therefore, it is appropriate to reduce the total number of unknowns by statically condensing out the rotational components of the bolt. Because no external moments are acting on the internal nodal points, the …`† ®rst variation of the strain energy with respect to the vector qx/ must vanish, thereby yielding the moment equilibrium equations in matrix form as …`†T …`†

…`†

…`†

bD/ qxD ‡ b// qx/ ˆ 0 Solving

…`† qx/

…A:23†

in the above equation and substitution into Eq. (A.21), after rearranging the terms, lead to T

…`† …`† Ux…`† ˆ 12D…`† x bx D x

…A:24† …`†

where the vector Dx…`† is identical to qxD and the matrix b…`† x is de®ned as 1 …`†T …`† 1 …`† b b b 2 D/ // D/

…`†

bx…`† ˆ bDD

…A:25†

As mentioned previously, the sti€ness properties of the bolt in the …x; z† and …y; z† planes are identical because the bolt material is isotropic and its cross-section is circular. Thus, the strain energy of the bolt in the …y; z† plane can be expressed as T

…`† …`† Uy…`† ˆ 12D…`† y by D y

…A:26†

in which the vector D…`† y contains n o T …`† …`† …`† ˆ D D


D D…`† yK y1 y2 y

…A:27†

…`† and the matrix b…`† y is identical to bx , i.e., T 1 …`† …`† b b 2 D/ //

…`†

by…`† ˆ bDD

1

…`†

bD/

…A:28†

The total strain energy of the bolt becomes T

T

…`† …`† 1 …`† …`† …`† U …`† ˆ 12D…`† x bx Dx ‡ 2Dy by Dy

or, in a more compact form, ( )T " #( ) …`† …`† …`† T 1 D D b 0 x x x U …`† ˆ ˆ 12D…`† b…`† D…`† …`† …`† 0 b D 2 D…`† y y y

…A:29†

…A:30†

The analytical derivation of Eq. (A.30) depends on the number of plates connected with the `th bolt. For example, the analytic derivation of matrix bx…`† (ˆ by…`† ) for bolts used in single- and double-lap joints is obtained as

7834

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

 …`† b…`† x ˆ by ˆ

b b

b b

 …A:31†

with b ˆ

12E` I` G` Af ` L3i G` Af ` ‡ 12E` I` Li

…A:32†

for a single-lap joint, and

…A:33†

…1`†

…2`†

with bij and bij de®ned by Eq. (A.16) for a double-lap joint. If a bolt connects more than three laminates, the analytic derivation would be too lengthy to present herein. For this case, numerical calculation of the matrices becomes more appropriate. Appendix B …k† For the kth region (laminate), the displacement components, u…k† x and uy , can be rewritten in terms of real vector quantities as





…k† u…k† x ; uy



L…k† X 2 N   X X …k`†T …k`†T Ux…rn† ; Uy…rn† a…k`† ˆ rn `ˆ1

…B:1†

rˆ1 nˆ N

L…k† X 2 N   X  X …k`†T …k`†T …k`†T Sxx…rn† ; Syy…rn† ; Sxy…rn† a…k`† Nx…k† ; Ny…k† ; Nxy…k† ˆ rn `ˆ1

…B:2†

rˆ1 nˆ N

where L…k† is the number of bolt holes in the kth laminate and the real vectors are given by      …k`†T Ux…rn† ˆ 2 Re prk …k`† Urn ; 2 Im prk …k`† Urn

…B:3a†

   …k`†T Uy…rn† ˆ 2 Re qrk …k`† Urn ;

  2 Im qrk …k`† Urn

…B:3b†

   …k`†T Sxx…rn† ˆ 2 Re l2rk …k`† urn ;

  2 Im l2rk …k`† urn

…B:4a†

  …k`†T Syy…rn† ˆ 2 Re

2 Im

  …k`†T Sxy…rn† ˆ 2 Re lrk

…k`†

 urn ;

…k`†

   T ˆ Re …k`† arn ; a…k`† rn

 urn ; Im



…k`†

 2 Im lrk

…k`†

arn



urn

…k`†



urn



…B:4b† …B:4c† …B:5†

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

The functions U…k`† and u…k`† are de®ned as rn rn    n d  ˆ n…k`† n…k`† and u…k`† ˆ …k`† U…k`† U…k`† r r rn rn rn dzr

7835

…B:6†

based on the explicit form of the complex potential functions given by Lekhnitskii (1968) N  n
X U…k`† z…k`† n…k`† a…k`† ˆ r r r rn

…B:7†

nˆ N n6ˆ0

…k`† map the `th circular hole in the kth laminate to in which a…k`† rn are the complex unknown coecients and nr a unit circle in the mapped plane, thus permitting Laurent series representations to approximate the ®eld …k`† variables. The prime denotes di€erentiation with respect to n…k`† introduced by r . The mapping functions nr Bowie (1956) are in the form r  

n…k`† r

ˆ

z…k`†  r

…k`†

zr

ak` …1

2

a2k` …1

l2rk †

…B:8†

ilrk † 1=2

in which ak` is the of the `th hole in the kth region and i ˆ … 1† . The sign of the square-root term is radius P 1. The complex parameters z…k† and z…k`† are de®ned as chosen so that n…k`† r r r z…k† r ˆ X ‡ lrk Y

and

z…k`† ˆ x…`† ‡ lrk y …`† r

in which X and Y are the global coordinates and x…`† and y …`† are the coordinates associated with the bolt hole. The complex parameters l1k and l2k are the roots of the characteristic equation derived by Lekhnitskii (1968). …k`† …k`† Terms arising from the expansion of Ua…rn† , Sab…rn† , and a…k`† rn for r ranging from 1 to 2 can be contained in the following vectors: n o …k`†T …k`†T …k`†T Ua…n† ˆ Ua…1n† ; Ua…2n† …B:9† n o …k`†T …k`†T …k`†T Sab…n† ˆ Sab…1n† ; Sab…2n†

…B:10†

n o T …k`†T …k`†T ˆ a1n ; a2n a…k`† n

…B:11†

The terms arising from the expansion of these vectors for n ranging from vectors de®ned as n o T …k`†T …k`†T …k`†T …k`†T …k`†T ˆ Ua… N † ; Ua… N ‡1† ; . . . ; Ua… 1† ; Ua…1† ; . . . ; Ua…N † U…k`† a …k`†T

Sab

n o …k`†T …k`†T …k`†T …k`†T …k`†T ˆ Sab… N † ; Sab… N ‡1† ; . . . ; Sab… 1† ; Sab…1† ; . . . ; Sab…N †

n o T …k`†T …k`†T …k`†T …k`†T …k`†T a…k`† ˆ a… N † ; a… N ‡1† ; . . . ; a… 1† ; a…1† ; . . . ; a…N †

N to N are contained in the …B:12† …B:13† …B:14†

with …a; b ˆ x; y†. Thus, using Eqs. (B.12) and (B.13), the series expansions for displacements and resultant stresses in Eqs. (B.1) and (B.2), respectively, can be rewritten as

7836

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

u…k† a ˆ

…`†

Nab ˆ

L…k† X `ˆ1

T

U…k`† a…k`† a

L…k† X `ˆ1

…k`†T

Sab

a…k`†

…B:15†

…B:16†

with …a ˆ x; y†. These expressions can be recast as u

…k†

ˆ

L…k† X

T

U…k`† a…k`†

…B:17†

`ˆ1

S…k† ˆ

L…i† X

T

S…k`† a…k`†

…B:18†

`ˆ1

by de®ning the following vectors n o T …k† u…k† ˆ u…k† ; u x y n o T N…k† ˆ Nxx…k† ; Nyy…k† ; Nxy…k† h i …k`† U U…k`† ˆ U…k`† x y h i …k`† …k`† S…k`† ˆ S…k`† xx Syy Sxy Finally, these equations can be expressed in a more compact form as T

u…k† ˆ U…k† a…k† T

N…k† ˆ S…k† a…k†

…B:19† …B:20†

where

h i T T T …k† T U…k† ˆ U…k1† U…k2†


U…kL †

h i T T T …k† T S…k† ˆ S…k1† S…k2†


S…kL † n o T T T …k† T a…k† ˆ a…k1† ; a…k2† ; . . . ; a…kL †

References Bowie, O.L., 1956. Analysis of an in®nite plate containing radial cracks originating at the boundary of an internal circular hole. Journal of Mathematics and Physics 35, 60±71. Gatewood, B.E., 1957. Thermal Stresses. McGraw-Hill, New York. Ireman, T., Nyman, T., Hellbom, K., 1993. On design methods for bolted joints in composite aircraft structures. Composites and Structures 25, 567±578.

V. Kradinov et al. / International Journal of Solids and Structures 38 (2001) 7801±7837

7837

Lekhnitskii, S.G., 1968. Anisotropic Plates. Gordon and Breach, New York. Madenci, E., Shkarayev, S., Sergeev, B., 1997. Analysis of composite laminates with multiple fasteners. Report to the Federal Aviation Administration. University of Arizona, Tucson, AZ. Madenci, E., Shkarayev, S., Sergeev, B., Oplinger, D.W., Shyprykevich, P., 1998. Analysis of composite laminates with multiple fasteners. International Journal of Solids and Structures 35, 1793±1811. Madenci, E., Sergeev, B., Shipilov, Y., Scotti, S., 1999. Analysis of double-lap composite joints with multiple-row bolts under combined mechanical and thermal loading. International Conference on Joining and Repair of Plastics and Composites, The Institution of Mechanical Engineers, London, England, 1999, pp. 55±64. Ramkumar, R.L., Saether, E.S., Appa, K., 1986. Strength analysis of laminated and metallic plates bolted together by many fasteners. AFWAL-TR-85-3034. Xiong, Y., Poon, C., 1994. A design model for composite joints with multiple fasteners. Aeronautical Note IAR-AN-80, NRC no. 3216, National Research Council, Canada.