Next Article in Journal
Threat-Oriented Collaborative Path Planning of Unmanned Reconnaissance Mission for the Target Group
Next Article in Special Issue
Critical Assessment of the Bonded Single Lap Joint Exposed to Cyclic Tensile Loading
Previous Article in Journal
Numerical Investigation on Aerodynamic Characteristics of an Active Jets-Matrix Serving as Pitch Control Surface
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Closed-Form Analysis of Thin-Walled Composite Beams Using Mixed Variational Approach

Department of Aerospace Information Engineering, Konkuk University, 120 Neungdong-ro, Gwangjin-gu, Seoul 05029, Korea
*
Author to whom correspondence should be addressed.
Aerospace 2022, 9(10), 576; https://doi.org/10.3390/aerospace9100576
Submission received: 1 September 2022 / Revised: 29 September 2022 / Accepted: 30 September 2022 / Published: 4 October 2022
(This article belongs to the Special Issue Composites in Aerospace Application)

Abstract

:
An analytical methodology is developed for thin-walled composite beams with arbitrary geometries and material distributions. The approach uses a mixed variational theorem which sets the shell displacements, and the shear stress resultants and hoop moments as the unknowns to obtain the stiffness constants at the level of Timoshenko–Vlasov beam. All the field equations and the continuity conditions that should be satisfied over the shell wall are derived in closed form as the part of the analysis. Numerical simulations are conducted to show the validity of the proposed analysis. The comparison of the predicted stresses and the stiffness constants indicates close correlation with those of the detailed finite element (FE)-based analysis and up-to-date beam approaches. Symbolically generated stiffness constants and the sectional warping deformation modes of coupled composite beams are presented explicitly to illustrate the strength of the proposed theory.

1. Introduction

Thin-walled beams made of composites exhibit complex mechanical behavior that cannot be observed for conventional isotropic counterpart structures. Particularly, anisotropic composites induce elastic couplings between various deformation modes such as extension, bending, torsion, and shear, which makes the analysis more involved. There is a need to model the local behavior of the shell wall properly, in response to the global deformation modes of the beam. Although three-dimensional (3D) solid or shell finite element (FE) approaches can lead to solutions of high resolution, these often demand heavy computational costs, which are unrealistic in many engineering applications (e.g., optimum structural design). During the past decades, development of accurate and efficient analysis methodologies for thin-walled composite beams has been a key research topic in the fields of the static and the dynamic analyses of helicopter and wind turbine blades [1,2,3].
The theoretical foundation for the modern composite beam analysis has been established by Berdichevskii [4] who justifies the decomposition of the 3D elasticity theory of the structure into 1D global beam problem and 2D local cross-sectional analysis, using the variational asymptotic arguments [3]. While the original 3D nature of the beam is essentially represented with the cross-sectional stiffness constants, a correct treatment of the sectional warpings is crucial for the successful development of a beam model. Many beam theories have been devised for generic cross-section analysis methods with detailed cross-sectional warpings and elastic coupling effects of composites. These range from finite element (FE)-based cross-section analyses [5,6,7,8,9], polynomial or series expansion methods [10,11,12], closed-form elasticity solutions [13,14], contour-based analytic formulations [15,16,17,18,19], and to semi-analytic approaches [20,21,22].
Each method has its own pros and cons depending upon kinematic foundations and numerical schemes adopted in the solution schemes. Among others, the analytical method has clear benefits in terms of computational efficiency and modeling simplicity. In addition, no serious discretization on 2D cross-sectional domain is needed which is typical in FE-based approaches. Instead, the thin-walled profile of a beam section is modeled as an assembly of 1D line segments which enables to compute the cross-sectional stiffness or flexibility constants with simple 1D line (analytic) integrals. Although degeneration of the cross-sectional information is inevitable, the beam response can be obtained in a set of closed-form solutions where the influences of beam sectional parameters can be clearly visible. Most of the analytical formulations have been made possible with simplifying (called ad hoc) assumptions such as the rigid cross-section in its own plane [1,2].
In the analytical beam approaches, there exist several attempts which are so general that none or the least level of beam assumptions is introduced where the nonclassical effects of composite beams are captured rigorously. Hodges and his coworkers [23,24,25] applied the variational-asymptotic method (VAM) to evaluate the stiffness constants and 3D warpings for anisotropic beams. In this method, the characteristic dimension of the cross-section with respect to the wavelength of the beam deformation mode is defined as a small parameter to find the desired representation of the beam analysis (e.g., Vlasov or Timoshenko level). The theory is generic and straightforward in mathematical terms, capturing all 3D warpings and elastic coupling effects precisely. However, an iterative asymptotic expansion is needed to systematically eliminate higher-order contributions from the strain energy. Furthermore, there appears to be a difficulty in reaching a beam representation equivalent to a Timoshenko-Vlasov level that takes into account the effects of the transverse shear and warping restraint simultaneously. Reissner and Tsai [26] presented a generalized solution procedure starting from the linear shell theory to obtain the sectional stiffness constants and 3D contour warpings for anisotropic beams under extension, torsion, and bending. Tsai [27] extended this concept to find the exact expressions for shell stress resultants and 3D warpings under the action of pure shear loads. Though the proposed method is generic and mathematically sound, the application is limited to simple cross-sections with homogeneity. Jung et al. [16] derived stiffness constants based on the Reissner’s mixed variational theorem [28] for extension, torsion, bending, shear, and restrained torsion. Despite the usage of simplifying assumptions of zero-hoop stress flows along the cross-sectional contour, the shell stress resultants obtained from the shell equilibrium equations resolved the limitations through the mixed formation to yield reasonable accuracy solutions as compared to those with the most sophistication. This feature has been recognized by the late Prof. Hodges [3] that the most accurate and powerful of the ad hoc methods to date appears to be Jung et al. (2002). The theory leads to [7 × 7] stiffness constants that include the effects of shell wall thickness, transverse shear, and restraint warping. However, only the torsion-related out-of-plane warpings are considered in the beam formulation with the lack of the stress recovery capability.
In the present study, the work of Jung et al. [16] is extended to combine the merits of Reissner and his associates [26,27,28] for variationally consistent, mixed beam analytical model that represents a Timoshenko-Vlasov beam. Given the shell theory as the starting point, all the stiffness coefficients as well as the field equilibrium equations and the boundary constraint conditions of the shell wall are derived in closed form. To facilitate the beam analysis, the same assumption adopted previously [16] is retained. The stress recovery part is included in the post-stage of the analysis to evaluate the layer-wise distribution of stresses (and strains) over the cross-sectional wall. The sectional warpings associated with different sets of loads are determined through the recovery process. The analysis is validated against the benchmark results available for composite beams in the literature. Symbolic expressions on the stiffness constants of an elastically coupled composite beam are presented to demonstrate the strength of the analytical beam theory developed based on the mixed variational formulation.

2. Theory

The theory is based on the mixed variational formulation developed by Reissner et al. [28,29] and Jung et al. [16]. In this method, either the displacements (active) or the stresses (reactive) are treated as the unknowns which will be determined simultaneously in the analysis. The details of the derivation leading to [7 × 7] flexibility or stiffness constants along with the recovery of stresses of thin-walled composite beams are described in this section.

2.1. Strain-Displacement Relations

Figure 1 shows the geometry and coordinate systems of a thin-walled beam section. Two systems of coordinate axes are employed: an orthogonal Cartesian system (x, y, z) for the global beam axis and a curvilinear coordinate system (x, s, n) for the local shell wall of the section. The displacements of an arbitrary material point on the beam cross-section, u 3 D = [ u v w ] T , are defined as the sum of the roto-translational displacements r of the beam reference line and the sectional warping function ω described along the midline contour as:
u 3 D = B d r + ω
where
B d = [ 1 0 0 0 z y 0 1 0 z 0 0 0 0 1 y 0 0 ] r = [ U V W ϕ β y β z ] T ω = [ ω x ω y ω z ] T
where ( U ,   V ,   W ) and ( ϕ ,   β y ,   β z ) are the translational displacements and rotation angles about x, y, z axes, respectively, and the superscript T denotes the transpose of an array. ( ω x ,   ω y ,   ω z ) are the contour warping functions along x, y, z axes, respectively. The shell midplane displacements u c 0 = [ u x 0 u s 0 u n 0 ] T are expressed via the transformation between the two reference frames as:
u c 0 = T s u 3 D
with the transformation matrix T s given by
T s = [ 1 0 0 0 y , s z , s 0 z , s y , s ]
where (   ) , s denotes the partial differentiation with respect to s. Assuming small strains and small local rotations, the strain- and curvature-displacement relations for the shell wall of an arbitrary geometry are written in symbolic forms as [16,30].
e = [ ϵ 0 κ γ ] T = s [ u c 0 Ψ ]
where ϵ 0 = [ ϵ x x 0 ϵ s s 0 γ x s 0 ] T are the in-plane strains at the midplane, κ = [ κ x x κ s s κ x s ] T the curvatures, γ = [ γ x n γ s n ] T the transverse shear strains across the shell wall, and Ψ = [ ψ x ψ s ] T the shell rotation angles about s and x axes, respectively. The operator matrix s above is given by:
s = [ (   ) , x 0 0 0 0 0 (   ) , s 1 / a 0 0 (   ) , s (   ) , x 0 0 0 0 0 0 (   ) , x 0 0 0 0 0 (   ) , s 0 0 0 (   ) , s (   ) , x 0 0 (   ) , x 1 0 0 1 / a (   ) , s 0 1 ]
where a is the local shell radius of the curvature. Inserting Equation (3) into (5), one can obtain the shell strain measures e in terms of the generalized beam displacement vector q and the derivatives of sectional warpings ω ˜ , x as [16]:
e = B e q + ω ˜ , x
with the arrays B e and q defined respectively as:
B e = [ 1 0 0 0 z y 0 0 0 0 0 0 0 0 0 y , s z , s r n 0 0 0 0 0 0 0 y , s z , s r t 0 0 0 0 0 0 0 0 0 0 2 0 0 0 0 z , s y , s 0 0 0 0 0 0 0 0 0 0 0 ]
q = [ U , x γ x y γ x z ϕ , x β y , x β z , x ϕ , x x ] T
where r n ,   r t are the offsets from the reference beam origin (see Figure 1), and γ x y ,   γ x z are the transverse shear strains of the cross-section, respectively. The sectional warping vector ω ˜ = [ ω 0 ] T in Equation (7) consists of the sectional warpings in its first three entries with the remaining ones being a null vector. The strain-displacement relations of Equation (7) form the basis of the displacement method for thin-walled beams.

2.2. Laminate Constitutive Relations

The classical lamination plate theory (CLPT) [31] is used to relate the stress resultants and couples with the strains and curvatures for the shell wall of the section as given by,
[ N M ] = [ A B B D ] [ ϵ 0 κ ] Q = k A γ
where N ,   Q ,   and   M are the stress resultants and couples of the laminate, and A ,   B ,   and   D are the laminate stiffness matrices for the extension, extension-bending coupling, and bending, and k  ( = 5 / 6 ) the shear correction coefficient, respectively [16]. Following Jung et al. [16], it is assumed that the hoop stress resultant N s s and the transverse shear resultant N s n are negligibly small. While shifting the reactive components of shell stress resultants and couples such as (Nxs, Mss, Nxn) into the right-hand side of the equation, one can express the shell constitutive relations in a semi-inverted form as:
[ p a ϵ r ] = [ C a a C a r C a r T C r r ] [ ϵ a p r ] γ x n = C x n N x n
with
p a = [ N x x M x x M x s ] T , p r = [ N x s M s s ] T ϵ a = [ ϵ x x 0 κ x x κ x s ] T ,   ϵ r = [ γ x s 0 κ s s ] T
where the subscripts a and r denote the active and reactive components, respectively [16,26,27,28]. The elements of the matrices C a a ,   C a r ,   C r r and C x n are found in Appendix A Equation (A1). Using the definition of shear strains, the active strains ε a are expressed as follows:
ε a = B q q b + ω , x
where q b = [ U , x ϕ , x β y , x β z , x ϕ , x x ] T are the generalized strains of the beam, and the matrix Bq is given by
B q = [ 1 0 z y ω x ϕ 0 0 y , s z , s r t 0 2 0 0 0 ]
where ω x ϕ is the out-of-plane warping function (unknown) under a unit twist of the cross-section. Equation (12) will be frequently recalled in the subsequent sections. It should be mentioned that only the sectional out-of-plane warping is survived due to the zero hoop stress resultant assumption ( N s s = 0 ), and thereby the in-plane warpings are dropped out from the warping deformation arrays in ω . It can be shown that the out-of-plane warping ω x is found by combining the definition of the shear strain γ x s 0 in Equations (5) and (1) as:
ω x = s 0 s γ x s 0 d s s 0 s r n d s   ϕ , x
The out-of-plane warping function ω x is not given a priori and obtained rigorously from the derivation. Inserting the constitutive relations of γ x s 0 into Equation (14), the out-of-plane warping deformation is obtained as:
ω x = ω q q b + ω q x q b , x
where ω q ,   ω q x are the row vectors that contain the corresponding out-of-plane warping function to each entry in q b ,   q b , x which will be determined later in the following section. It is noted that the sectional warpings are not only functions of the geometries but also the material distributions along the contour of the cross-section.

2.3. Governing Equations

The governing equations of the shell wall of the cross-section are found using a mixed variational theorem which is stated as [28]:
δ Π = δ U δ W = 0
where the virtual strain energy per unit length, δU, and the external virtual work done over the cross-sectional area A per unit length, δW, are respectively given by [26]
δ U = C [ δ ϵ a T p a ¯ + δ ( ϵ r k ) T p r + δ p r T ( ϵ r k ϵ r ) + δ Q T γ ] d s δ W = ( A δ u c T t d A ) , x
where C is the length of the sectional contour, ϵ r k is the kinematically defined reactive strains, and t = [ σ x n σ x s σ x x ] T is the surface traction distributed over the cross-section. It can be shown that the variational statement of Equation (16) leads to a set of Euler-Lagrange equations followed by
N x x , x + N x s , s = 0 N s s , s + N s n / a + N x s , x = 0 N s n , s N s s / a + N x n , x = 0 M x s , s N x n + M x x , x = 0 M s s , s N s n + M x s , x = 0
along with the constraint conditions for the reactive strain and curvature components. The compatibility conditions set for the reactive strains ϵ r are satisfied in a weak sense, while δ p r acting as the Lagrange multiplier to enforce the constraint conditions. While performing the integration by parts operation of Equation (18) leads to the following relations for the reactive shell stress resultants and couples p r T = [ N x s M s s ] as:
p r = B p p r 0 p r b
with
p r 0 = [ N x s 0 M s s 0 N y 0 N z 0 ] T ,         B p = [ 1 0 0 0 0 1 z y ] p r b = [ s 0 s N x x , x d s s 0 s ( 2 M x s , x + N x s , x r n ) d s ] = S r q b , x
where p r 0 is the unknown constant of integration defined for each cell of closed cross-sections, p r b is the vector that contains the shear flow resultant and hoop moment which vary along the cross-sectional contours, and S r is the corresponding branch flows vector, respectively. Equation (19) is inserted into δ p r T in Equation (17) to have the continuity conditions for the reactive strains:
B p T ϵ r d s = A s q b
with A s given by
A s = [ 0 r n d s 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
where the symbol denotes a closed-loop integral over the contour of the cross-section. It is remarked that there exist four sets of constraint equations for each cell of closed, multi-cell sections, similar to the cases of [16,23,26]. Inserting Equations (10) and (12) into (21), one obtains the constraint conditions expressed in a form as:
p r 0 = C 1 q b + C 2 q b , x
where C 1 ,   C 2 are given as closed-loop integrals that represent the material distributions and sectional warpings along the contour of the shell wall. Inserting Equations (23) into (19), one can find the reactive stress resultant vector p r as:
p r = B p C 1 q b + ( B p C 2 S r ) q b , x = f q b + g q b , x
where f indicates the circuit shear flows and hoop moments of a closed cell of the cross-section with respect to the generalized beam strains, and g is the resultant shear flows and hoop moments distributed along the sectional contour. The second part of Equation (24) is related with the transverse shear couplings in the Timoshenko-like beam which is neglected with the Euler–Bernoulli assumption. Having obtained the reactive stress resultant vector p r , the shell transverse shear resultant Nxn can be expressed using the fourth equilibrium equations of Equation (18) along with Equation (10) as
N x n = h q b + m q b , x
where h and m are defined similar to f and g in Equation (24), respectively. The effects of the transverse shear resultant Nxn should be negligible on the global beam behavior, and these have been neglected in most beam theories in the literature. Nevertheless, they are kept for more consistency and completeness in the present beam formulation.
Using Equations (7), (10), (24) and (25), the active part (singly underlined term) of the strain energy variation in Equation (16) leads to the cross-sectional stress resultants F b = [ F x M x M y M z M ω ] T expressed in terms of the generalized beam strains and their derivatives as:
F b = K b b q b + K b v q b , x
where K b b is the (5 × 5) cross-sectional stiffness matrix that represents the beam in the Euler–Bernoulli–Vlasov level, and K b v represents a coupling between F b and q b , x . The relationships between q b and the transverse shear forces F s on the cross-section are yet to be determined. The transverse shear forces affect F b , resulting in the variations of the shell stress resultants and generalized strains with respect to x. The relations are obtained under a pure shear condition [27]. When the higher-order terms ( q b , x x ) are neglected, the differentiation of Equation (26) by x gives
F b , x = K b b q b , x
Considering the equilibrium of forces and moments of the beam in the transverse directions, one obtains the following relations
F b , x = T b   F s
with
T b = [ 0 0 0 0 0 0 0 1 0 1 0 0 0 0 1 ] , F s = [ F y F z T ω ]
where T ω ( T ω = M ω , x ) is the warping torque [2,32]. It is remarked that only the transverse shear forces remain as the independent variables since the Vlasov torsion T ω is combined with the St. Venant torsion T S V to form the complete torsion equation: M x = T S V + T ω [32]. Equating Equations (27) and (28) yields the relations between q b , x and F s . Inserting this relation into Equations (26) yields the following equation:
F b = K b b q b + K b s F s
with the relations K b s = K b v S b b T b and S b b = K b b - 1 , respectively. The pure shear condition is introduced to relate q b and F s with a null vector for F b ,
q b = S b b K b s F s = S b s F s
A linear superposition for the bending and shear contributions from Equations (30) and (31) leads to the generalized beam strains q b expressed in terms of the generalized beam forces as:
q b = S b b F b + S b s F s
The beam cross-sectional stiffness constants taking into account the effects of the transverse shear of the section and warping inhibition can now be determined. Using Equations (10), (12), (15), (24), and (25), one obtains the variation of strain energy in the following form:
δ U = δ [ q b F s ] T C [ K b b K b s K b s T K s s ] d s [ q b F s ]
Using the relationship between q b and F s given in Equation (32), the strain energy variation in Equation (33) is transformed as:
δ U = δ [ F b F s ] T S ¯ 7 × 7 [ F b F s ]
where S ¯ 7 × 7 is a (7 × 7) flexibility matrix given by
S ¯ 7 × 7 = [ S b b S b s 0 I 2 × 2 ] T [ K b b K b v K b v T K v v ] [ S b b S b s 0 I 2 × 2 ]
where I 2 × 2 is a (2 × 2) identity matrix. An inversion of the flexibility matrix yields the (7 × 7) stiffness matrix: K ¯ 7 × 7 = S ¯ 7 × 7 - 1 . The elements in K ¯ 7 × 7 are properly rearranged so that the beam force–displacement relations of the cross-section are expressed in a convenient form:
F = K 7 × 7   q
where
F = [ F x F y F z M x M y M z M ω ] T
The beam generalized strains q are defined previously in Equation (8b). The (7 × 7) stiffness constants Kij describe the beam in a Timoshenko-Vlasov level while incorporating the effects of elastic couplings and shell wall thickness of composite beams. It should be noted that the transverse shear couplings essentially affect the elements of K b b so that the (7 × 7) stiffness constants associated with the extension, bending, torsion, and their coupled entities become substantially modified.

2.4. Recovery Relations

The stresses and strains over the shell wall thickness can be recovered through an inverse-procedure of obtaining the stiffness matrix K 7 × 7   . Once the (7 × 7) stiffness constants are determined, the cross-sectional stress resultants F are used to find q b via Equation (32). q b and F s are then used to evaluate the sectional warping (Equation (15)) and the reactive stress resultants p r (Equation (24)), which in turn gives the active strains ϵ a (Equation (12)). The reactive strains ϵ r are also recovered from Equation (10). With ϵ a and ϵ r are evaluated, the hoop strain ϵ s s 0 which has been omitted by the zero-hoop stress flow assumption ( N s s = 0 ) can be recovered using the laminate constitutive relations (Equation (9)) given by:
ϵ s s 0 = ( A 12 ϵ x x 0 + A 26 γ x s 0 + B 12 κ x x + B 22 κ s s + B 26 κ x s ) / A 22
where Aij and Bij (i, j = 1, 2, 6) are the laminate stiffness components defined in Equation (9). It is noted that the above ϵ s s 0 is constructed from the constitutive relations (not from the kinematics). The recovered strains and curvatures at the midplane are then used to evaluate the (layer-wise) off-axis strains ϵ off = [ ϵ x x ϵ s s γ x s ] T at any location of the shell wall with an offset n along the thickness direction as:
ϵ off = ϵ 0 + n κ
where ϵ 0 denotes the shell strain measures at the midplane. The off-axis stresses σ off = [ σ x x σ s s σ x s ] T are then evaluated using the stress–strain relations:
σ off = Q ¯ ϵ off
where Q ¯ is the off-axis lamina stiffness matrix [31]. Further recovery analysis is performed for on-axis stresses and strains using transformational relations given in [31]. It should be mentioned that the overall procedure to obtain the field variables, such as stresses and strains, are purely analytical, in the sense that all the cross-sectional stiffness constants (including sectional warpings) are obtained using a series of analytical integrations, which does not require interior domain discretizations, heavy numerical factorizations, and/or eigenvalue analyses through the entire derivation stage.

3. Numerical Examples

Numerical simulations are conducted to verify the current mixed beam formulation while demonstrating the strength and the accuracy of the proposed methodology. The examples include anisotropic rectangular strip section beam, four-layered laminated composite beam, and single-cell box section beam with circumferentially asymmetric (CAS) layup.

3.1. Anisotropic Rectangular Strip Beam

The first example is a composite thin rectangular strip section beam which has been studied by Yu and Hodges [24]. The strip is composed of a symmetric layup (bending-torsion coupling) with the dimension of the width b and the thickness h. In [24], the generalized expressions for (4 × 4) sectional stiffness constants are obtained using VAM based on the CLPT. Though their variational-asymptotic expansion has been limited by Euler–Bernoulli level, the symbolic expressions of the stiffness constants derived for an isotropic beam and the predicted values obtained for the composite strip beam are used as the reference for the cross-validation of the present results. For the composite strip beam, the non-zero stiffness constants in a Timoshenko-Vlasov level are obtained in symbolic form as:
K 11 = b { A 11 ( A 16 ) 2 A 66 + 5 ( A 16 ) 2 6 A 66 _ 5 36 ( ( A 16 ) 2 A 66 ) 2 A 11 5 6 ( A 16 ) 2 A 66 _ _ } K 12 = b { 5 6 A 16 _ 5 36 ( A 16 ) 3 A 66 A 11 5 6 ( A 16 ) 2 A 66 _ _ } ,   K 22 = b { 5 6 A 66 _ 5 36 ( A 16 ) 2 A 11 5 6 ( A 16 ) 2 A 66 _ _ } K 33 = 48 b [ ( D 16 D 12 D 26 D 22 ) 2 ( D 11 ( D 12 ) 2 D 22 ) ( D 66 ( D 26 ) 2 D 22 ) ] 2 b 2 ( D 11 ( D 12 ) 2 D 22 ) ( D 16 D 12 D 26 D 22 ) 2 + 12 C x n [ 3 ( D 16 D 12 D 26 D 22 ) 2 2 ( D 11 ( D 12 ) 2 D 22 ) ( D 66 ( D 26 ) 2 D 22 ) ] 2 K 44 = 4 b ( D 66 ( D 26 ) 2 D 22 ) ,   K 45 = 2 b ( D 16 D 12 D 26 D 22 ) ,   K 55 = b ( D 11 ( D 12 ) 2 D 22 ) K 66 = b 3 12 ( A 11 ( A 16 ) 2 A 66 ) ,   K 77 = b 3 12 ( D 11 ( D 12 ) 2 D 22 )
where K ij denotes the respective entities of the stiffness matrix of K 7 × 7 given in Equation (36). The off-diagonal elements K 12 and K 45 are the extension-shear and bending-torsion couplings, respectively, which are characterized for composite strip beams. The underlined terms are responsible for the shear deformation of the cross-section in x-y plane. The terms with single underline are due to the elastic couplings between the extension and shear while the terms with double underlines are due to the sectional warpings ω x related with the shear deformation. All the stiffness constants are derived in closed form which include the effects of the shell wall thickness, the transverse shear over both the shell wall and the beam, and the sectional warpings due to the couplings. It is noted that, as compared with the results of [24], the (4 × 4) stiffness coefficients expressed in symbolic form (isotropic beam) and numerical results (up to six digits) obtained for the composite beam remain identical with the present predictions (not shown explicitly), which confirms the soundness of the present mixed beam formulation.

3.2. Four-Layered Laminated Composite Beam

The next example is considered to illustrate the recovery analysis capability. A four-layered laminated composite beam employed in Liu and Yu [33] is considered for this purpose. Figure 2 shows the schematic of the beam with one end fixed and the other end subjected to a tip axial force (Fx = 10 kN). The beam is 1 m long with a width and height of 0.18 m and 0.04 m, respectively. The section consists of four layers having [90/0/90/0], which induces elastic coupling between extension and bending. The mechanical material properties are [33]: E1 = 110.5 GPa, E2 = E3 = 13.64 GPa, G12 = 3.92 GPa, G13 = G23 = 3.26 GPa, and ν12 = ν13 = 0.329.
In order to validate the current analysis, a detailed 3D, FE structural analysis is carried out using the commercial software MSC.NASTRAN (version: 2021.1, MSC software, Irvine, CA, United States, https://www.mscsoftware.com/product/msc-nastran (access date: 28 July 2022)). The beam is discretized using the 20-noded solid element (CHEXA) having 8 finite elements through the thickness, 36 elements along the width, and 200 elements along the beam length, leading to 57,600 in total. A RBE2 element is used to apply the axial load at the free tip. The predicted stresses through the thickness direction at the center of the beam, as indicated in Figure 2, are compared against the results by Liu and Yu [33] and the 3D solid FE analysis. Figure 3 shows the comparison between the three sets of results obtained respectively for the normal stresses ( σ x x ,   σ y y ) distributed through the thickness. It is observed that the present results are in excellent agreements with the detailed 3D FE results and the other beam prediction. It is emphasized that, as can be seen in Figure 3b, the in-plane normal stress σ y y along the shell wall can also be recovered using Equation (40), even though the deformation along the contour is kinematically suppressed due to the zero-hoop stress flow assumption.

3.3. Single-Cell Composite Box Beam with CAS Layup

The third example is a thin-walled, single-cell box-section beam with CAS layup, which has been exploited by many researchers for the validation works [34,35,36]. Figure 4 shows the geometry and the ply layup which induce couplings between bending-torsion and extension-shear motions. The outer dimensions of the box section are: b = 24.206 mm, h = 13.640 mm, with the thickness of an individual layer t = 0.127 mm. The effective length of the cantilevered beam is 762 mm. The positive directions of the local curvilinear coordinates at each wall are indicated in Figure 4. The mechanical material properties are given as [36]: E1 = 141.96 GPa, E2 = E3 = 9.79 Gpa, G12 = G13 = 6.00 Gpa, G23 = 4.80 Gpa, ν12 = ν13 = 0.42.
Table 1 shows the comparison of non-zero stiffness constants predicted between the present analysis and the other references [34,35,36] for the box-section beam when the ply orientation angle θ is set to 15 deg. The NABSA is a code name for detailed FE beam model constructed based on the work of Giavotto et al. [5], and the stiffness results are taken from Popescu and Hodges [34]. The work of Krenk and Couturier [35] uses 3D Fes to compute the stiffness parameters using a complementary form of the strain energy. In Dhadwal and Jung [36], 2D FE-based section analysis method is developed to compute (6 × 6) stiffness constants where 3D warpings and stress fields over the cross-sectional nodes are solved simultaneously. The comparison results indicate good agreements between the different sets of predictions. The maximum error of the present computation is less than 1.9%, as compared with the more involved, FE-based cross-sectional analysis predictions, which proves the validity of the present mixed beam approach.
The sectional warping deformation modes (continuous lines) associated with the six fundamental modes of extension, shear, torsion and bending are presented in Figure 5 along with their corresponding undeformed modes (dashed lines). As is indicated in the plots, all the warping deformations are influenced by the couplings between the extension-shear and bending-torsion motions. In Figure 6, the predicted stresses of the composite box beam subjected to a unit tip shear load (4.448 N) are compared against the ones of Dhadwal and Jung [36] and 3D FE analysis using the MSC.NASTRAN. The 3D FE model consists of 20-noded 1,152,000 CHEXA solid elements, along with a RBE2 element installed at the tip for the load application. Figure 6 presents the layer-wise distribution of stresses obtained along the A-A line of the left wall of the box section (see Figure 4). The in-plane normal stress σ s s is recovered as well as the other shell stress components ( σ x x ,   σ x s ). It is observed that the present predictions capture the general trends of the linear variation over each layer of the box section, with the amplitudes agreed well with those of the detailed FE-based method and 3D FE analysis. The comparison demonstrates again the validity of the present work for the stress recovery analysis capability of closed cross-sections. It should be noted that the present predictions are resulted from purely analytical treatments of the sectional integrals without resorting to any mesh discretization over the cross-section. No convergence check is needed as well, such as the case with the FE-based approaches, which greatly saves the computation cost.

4. Concluding Remarks

In this study, a variationally consistent, thin-walled composite beam model was developed combining the mixed beam theory and the generalized beam formulation. Starting from the shell theory, the Euler–Lagrange equations and the constraint conditions of the shell wall of the beam were obtained correctly as the part of the analysis. In order to facilitate the beam analysis, a zero-hoop stress flow assumption was retained for the sectional wall. The effects of the shell wall thickness and the out-of-plane warpings derived from the constitutive relations were incorporated in the theory, leading to the (7 × 7) stiffness constants equivalent to the level of Timoshenko-Vlasov beam. Three different, elastically coupled composite beams were tested to show the validity of the present analysis. First, symbolic expressions were derived explicitly for (7 × 7) stiffness constants of an anisotropic strip beam, identifying each of the contributing terms due to the couplings. Next, the cross-sectional warping deformation modes associated with the loads were identified while capturing the coupling behavior of the beam motions. All the stress fields including the in-plane stress were recovered and compared with 3D FE and other detailed cross-sectional analysis results, resulting in an excellent correlation between the predictions. Finally, the comparison of the stiffness constants for composite box section demonstrated a deviation less than 1.9%, against up-to-date FE-based methods. Given no mesh discretization nor heavy numerical factorization was involved in obtaining the warping and stiffness constants, the present work demonstrated excellent prediction capability which could be exploited further for an optimum structural design of composite helicopter blades and wind turbine blades.

Author Contributions

Conceptualization, J.S.B. and S.N.J.; software, J.S.B.; validation, J.S.B.; writing—original draft preparation, J.S.B. and S.N.J.; writing—review and editing, J.S.B. and S.N.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

The data presented in this study are available upon request from the corresponding author.

Acknowledgments

This paper was supported by Konkuk University in 2019.

Conflicts of Interest

The authors declare no conflict of interest.

Nomenclature

a Radius of curvature for the shell wall
F x Normal force acting on the cross-section along x axis
Fy, FzTransverse shear forces acting on the cross-section
k Shear correction coefficient
M x Beam torsional moment
Mxx, Mss, MxsShell moment resultants
My, MzBeam bending moments about y, z axes
M ω Beam torsional bi-moment
Nxx, Nss, NxsShell force resultants
Nxn, NsnTransverse shear resultants for the shell
ux, us, unShell displacements
U, V, WBeam displacements
ϵ x x 0 , ϵ s s 0 , γ x s 0 Shell membrane strain measures in x, s coordinates
ϕ ,   β y ,   β z Cross-sectional rotations of the beam
κ x x , κ s s , κ x s Shell curvature measures in x, s coordinates
γ x n , γ s n Transverse shear strain measures over the shell wall
ω x , ω y , ω z Sectional warping functions in x, y, z coordinates
ω x ϕ Torsion-related out-of-plane warping function
ψ x , ψ s Shell rotation angles about x, s axes
( ) T Transpose of an array
( ) , s Partial differentiation with s coordinate
( ) , x Partial differentiation with x coordinate

Appendix A

The submatrices C a a , C a r and C r r listed in Equation (9) are obtained as:
C a a = [ C n ϵ C n κ C n ϕ C n κ C m κ C m ϕ C n ϕ C m ϕ C ϕ ϕ ] ,         C a r = [ C n γ C n τ C m γ C m τ C ϕ γ C ϕ τ ] ,         C r r = [ C γ γ C γ τ C γ τ C τ τ ]
where the respective laminate stiffness constants are given by
C n ε = A 11 { A 16 ( A 16 D 22 B 12 B 26 ) B 12 ( A 16 B 26 A 66 B 12 } / Δ C n κ = B 11 { A 16 ( D 22 B 61 D 12 B 26 ) B 12 ( B 61 B 26 A 66 D 12 } / Δ C n φ = B 16 { A 16 ( D 22 B 66 D 26 B 26 ) B 12 ( B 26 B 66 A 66 D 26 } / Δ C n γ = ( A 16 D 22 B 12 B 26 ) / Δ   , C n τ = ( A 66 B 12 A 16 B 26 ) / Δ C m κ = D 11 { B 61 ( D 22 B 61 D 12 B 26 ) D 12 ( B 61 B 26 A 66 D 12 } / Δ C m φ = D 16 { B 61 ( D 22 B 66 D 26 B 26 ) D 12 ( B 26 B 66 A 66 D 26 } / Δ C m γ = ( B 61 D 22 D 12 B 26 ) / Δ   , C m τ = ( A 66 D 12 B 61 B 26 ) / Δ C φ φ = D 66 { B 66 ( D 22 B 66 D 26 B 26 ) D 26 ( B 26 B 66 A 66 D 26 } / Δ C φ γ = ( D 22 B 66 D 26 B 26 ) / Δ   , C φ τ = ( A 66 D 26 B 66 B 26 ) / Δ C γ γ = D 22 / Δ   , C γ τ = B 26 / Δ   , C τ τ = A 66 / Δ C x n = { k ( A 44 A 45 2 A 55 ) } 1
with Δ = A 66 D 22 ( B 26 ) 2 . The (   ) indicates that the values are obtained after the zero-hoop stress flow assumption ( N s s = 0 ) is applied to Equation (9).

References

  1. Hodges, D.H. Review of composite rotor blade modeling. AIAA J. 1990, 28, 561–565. [Google Scholar] [CrossRef]
  2. Jung, S.N.; Nagaraj, V.T.; Chopra, I. Assessment of composite rotor blade modeling techniques. J. Am. Helicopter Soc. 1999, 44, 188–205. [Google Scholar] [CrossRef]
  3. Hodges, D.H. Nonlinear Composite Beam Theory; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2006. [Google Scholar]
  4. Berdichevskii, V.L. On the energy of an elastic rod. J. Appl. Math. Mech. 1982, 45, 518–529. [Google Scholar] [CrossRef]
  5. Giavotto, V.; Borri, M.; Mantegazza, P.; Ghiringhelli, G.; Carmaschi, V.; Maffioli, G.; Mussi, F. Anisotropic beam theory and applications. Comput. Struct. 1983, 16, 403–413. [Google Scholar] [CrossRef]
  6. Morandini, M.; Chierichetti, M.; Mantegazza, P. Characteristic behavior of prismatic anisotropic beam via generalized eigenvectors. Int. J. Solids Struct. 2010, 47, 1327–1337. [Google Scholar] [CrossRef] [Green Version]
  7. Yu, W.; Hodges, D.H.; Ho, J.C. Variational asymptotic beam sectional analysis—An updated version. Int. J. Eng. Sci. 2012, 59, 40–64. [Google Scholar] [CrossRef]
  8. Han, S.; Bauchau, O.A. Nonlinear three-dimensional beam theory for flexible multibody dynamics. Multibody Syst. Dyn. 2015, 34, 211–242. [Google Scholar] [CrossRef]
  9. Dhadwal, M.K.; Jung, S.N. Multifield variational sectional analysis for accurate stress computation of multilayered composite beams. AIAA J. 2019, 57, 1702–1714. [Google Scholar] [CrossRef]
  10. Carrera, E.; Giunta, G. Refined beam theories based on a unified formulation. Int. J. Appl. Mech. 2010, 2, 117–143. [Google Scholar] [CrossRef]
  11. Carrera, E.; de Miguel, A.; Pagani, A. Hierarchical theories of structures based on Legendre polynomial expansions with finite element applications. Int. J. Mech. Sci. 2017, 120, 286–300. [Google Scholar] [CrossRef]
  12. Heyliger, P.R. Elasticity alternatives to generalized Vlasov and Timoshenko models for composite beams. Compos. Struct. 2016, 145, 80–88. [Google Scholar] [CrossRef]
  13. Yu, W.; Hodges, D.H. Elasticity solutions versus asymptotic sectional analysis of homogeneous, isotropic, prismatic beams. J. Appl. Mech. 2004, 71, 15–23. [Google Scholar] [CrossRef]
  14. Ho, J.C. Shear Stiffness of Homogeneous, Orthotropic, Prismatic Beams. AIAA J. 2017, 55, 4357–4363. [Google Scholar] [CrossRef]
  15. Smith, E.C.; Chopra, I. Formulation and evaluation of an analytical model for composite box beams. J. Am. Helicopter Soc. 1991, 36, 23–35. [Google Scholar] [CrossRef]
  16. Jung, S.N.; Nagaraj, V.T.; Chopra, I. Refined structural model for thin- and thick-walled composite rotor blades. AIAA J. 2002, 40, 105–116. [Google Scholar] [CrossRef]
  17. Librescu, L.; Song, O. Thin-Walled Composite Beams: Theory and Application; Springer: Berlin/Heidelberg, Germany, 2006. [Google Scholar]
  18. Zhang, C.; Wang, S. Structural mechanical modeling of thin-walled closed-section composite beams, Part 1: Single-cell cross section. Compos. Struct. 2014, 113, 12–22. [Google Scholar] [CrossRef]
  19. Vo, T.P.; Lee, J. Flexural–torsional behavior of thin-walled closed-section composite box beams. Eng. Struct. 2007, 29, 1774–1782. [Google Scholar] [CrossRef] [Green Version]
  20. Vieira, R.; Virtuoso, F.; Pereira, E. A higher order model for thin-walled structures with deformable cross-sections. Int. J. Solids Struct. 2014, 51, 575–598. [Google Scholar] [CrossRef] [Green Version]
  21. Han, S.; Bauchau, O.A. On the analysis of thin-walled beams based on Hamiltonian formalism. Comput. Struct. 2016, 170, 37–48. [Google Scholar] [CrossRef]
  22. Deo, A.; Yu, W. Thin-walled composite beam cross-sectional analysis using the mechanics of structure genome. Thin-Walled Struct. 2020, 152, 106663. [Google Scholar] [CrossRef]
  23. Volovoi, V.V.; Hodges, D.H. Theory of anisotropic thin-walled beams. J. Appl. Mech. 2000, 67, 453–459. [Google Scholar] [CrossRef] [Green Version]
  24. Yu, W.; Hodges, D.H. Best Strip-Beam Properties Derivable from Classical Lamination Theory. AIAA J. 2008, 46, 1719–1724. [Google Scholar] [CrossRef]
  25. Harursampath, D.; Harish, A.B.; Hodges, D.H. Model reduction in thin-walled open-section composite beams using variational asymptotic method. Part I: Theory. Thin-Walled Struct. 2017, 117, 356–366. [Google Scholar] [CrossRef]
  26. Reissner, E.; Tsai, W.T. Pure bending, stretching, and twisting of anisotropic cylindrical shells. J. Appl. Mech. 1972, 39, 148–154. [Google Scholar] [CrossRef]
  27. Tsai, W. On the problem of flexure of anisotropic cylindrical shells. Int. J. Solids Struct. 1973, 9, 1155–1175. [Google Scholar] [CrossRef]
  28. Reissner, E. On mixed variational formulations in finite elasticity. Acta Mech. 1985, 56, 117–125. [Google Scholar] [CrossRef]
  29. Murakami, H.; Reissner, E.; Yamakawa, J. Anisotropic beam theory with shear deformation. J. Appl. Mech. 1996, 63, 660–668. [Google Scholar] [CrossRef]
  30. Kraus, H. Thin Elastic Shells; Wiley: New York, NY, USA, 1967. [Google Scholar]
  31. Jones, R.M. Mechanics of Composite Materials, 2nd ed.; Taylor and Francis: Boca Raton, FL, USA, 1999. [Google Scholar]
  32. Jung, S.N.; Lee, J.Y. Closed-form analysis of thin-walled composite I-beams considering non-classical effects. Compos. Struct. 2003, 60, 9–17. [Google Scholar] [CrossRef]
  33. Liu, X.; Yu, W. A novel approach to analyze beam-like composite structures using mechanics of structure genome. Adv. Eng. Softw. 2016, 100, 238–251. [Google Scholar] [CrossRef]
  34. Popsecu, B.; Hodges, D.H. On asymptotically correct Timoshenko-like anisotropic beam theory. Int. J. Solids Struct. 2000, 37, 535–558. [Google Scholar] [CrossRef]
  35. Krenk, S.; Couturier, P.J. Equilibrium-based nonhomogeneous anisotropic beam element. AIAA J. 2017, 55, 2773–2782. [Google Scholar] [CrossRef] [Green Version]
  36. Dhadwal, M.K.; Jung, S.N. Generalized multifield variational formulation with interlaminar stress continuity for multilayered anisotropic beams. Compos. B Eng. 2019, 168, 476–487. [Google Scholar] [CrossRef]
Figure 1. Schematic of the reference coordinate frames.
Figure 1. Schematic of the reference coordinate frames.
Aerospace 09 00576 g001
Figure 2. Schematic of the four-layer laminated beam.
Figure 2. Schematic of the four-layer laminated beam.
Aerospace 09 00576 g002
Figure 3. Comparison of predicted normal stresses through the thickness between 3D FE, Liu and Yu (2016) [33] and the present: (a) axial stress; (b) transverse stress.
Figure 3. Comparison of predicted normal stresses through the thickness between 3D FE, Liu and Yu (2016) [33] and the present: (a) axial stress; (b) transverse stress.
Aerospace 09 00576 g003
Figure 4. Geometry and ply layup structure of the composite single-cell box section.
Figure 4. Geometry and ply layup structure of the composite single-cell box section.
Aerospace 09 00576 g004
Figure 5. Out-of-plane warping modes under unit loads: (a) Fx; (b) Fy; (c) Fz; (d) Mx; (e) My; (f) Mz.
Figure 5. Out-of-plane warping modes under unit loads: (a) Fx; (b) Fy; (c) Fz; (d) Mx; (e) My; (f) Mz.
Aerospace 09 00576 g005
Figure 6. Comparison of the layer-wise stress distributions at the left wall between the present, Dhadwal and Jung (2019) [36] and 3D FE: (a) axial stress; (b) hoop stress; (c) tangential stress.
Figure 6. Comparison of the layer-wise stress distributions at the left wall between the present, Dhadwal and Jung (2019) [36] and 3D FE: (a) axial stress; (b) hoop stress; (c) tangential stress.
Aerospace 09 00576 g006
Table 1. Comparison of non-zero stiffness constants of single-cell composite box section.
Table 1. Comparison of non-zero stiffness constants of single-cell composite box section.
StiffnessNABSA [34]Krenk et al. (2017) [35] (%) 1Dhadwal et al. (2019) [36] (%) 1Present (%) 1
K11 (N)6.0938 × 1066.1100 × 106 (0.27)6.1200 × 106 (0.43)6.1124 × 106 (0.31)
K12 (N)8.1843 × 1058.1900 × 105 (0.07)8.1600 × 105 (−0.30)8.1413 × 105 (−0.53)
K22 (N)3.9320 × 1053.9400 × 105 (0.20)3.9500 × 105 (0.46)3.8774 × 105 (−1.39)
K33 (N)1.7570 × 1051.7600 × 105 (0.17)1.7700 × 105 (0.74)1.7247 × 105 (−1.83)
K44 (N-m2)4.9645 × 1014.9800 × 101 (0.31)5.0000 × 101 (0.71)4.9093 × 101 (−1.11)
K45 (N-m2)−5.1654 × 101−5.1800 × 101 (0.28)−5.1500 × 101 (−0.30)−5.1503 × 101 (−0.29)
K55 (N-m2)1.7448 × 1021.7500 × 102 (0.30)1.7500 × 102 (0.30)1.7439 × 102 (−0.05)
K66 (N-m2)4.1036 × 1024.1000 × 102 (−0.09)4.1200 × 102 (0.40)4.0844 × 102 (−0.47)
1 Percentage differences are evaluated with respect to NABSA values as: (Analysis-NABSA)/NABSA × 100.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bae, J.S.; Jung, S.N. Closed-Form Analysis of Thin-Walled Composite Beams Using Mixed Variational Approach. Aerospace 2022, 9, 576. https://doi.org/10.3390/aerospace9100576

AMA Style

Bae JS, Jung SN. Closed-Form Analysis of Thin-Walled Composite Beams Using Mixed Variational Approach. Aerospace. 2022; 9(10):576. https://doi.org/10.3390/aerospace9100576

Chicago/Turabian Style

Bae, Jae Seong, and Sung Nam Jung. 2022. "Closed-Form Analysis of Thin-Walled Composite Beams Using Mixed Variational Approach" Aerospace 9, no. 10: 576. https://doi.org/10.3390/aerospace9100576

APA Style

Bae, J. S., & Jung, S. N. (2022). Closed-Form Analysis of Thin-Walled Composite Beams Using Mixed Variational Approach. Aerospace, 9(10), 576. https://doi.org/10.3390/aerospace9100576

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop