Analytical solution for the fully coupled static response of variable stiffness composite beams

Analytical Solution for the Fully Coupled Static Response of Variable Stiffness Composite Beams

Journal Pre-proof

Analytical Solution for the Fully Coupled Static Response of Variable Stiffness Composite Beams Pedram Khaneh Masjedi, Paul M. Weaver PII: DOI: Reference:

S0307-904X(19)30749-8 https://doi.org/10.1016/j.apm.2019.12.010 APM 13196

To appear in:

Applied Mathematical Modelling

Received date: Revised date: Accepted date:

25 July 2019 27 November 2019 5 December 2019

Please cite this article as: Pedram Khaneh Masjedi, Paul M. Weaver, Analytical Solution for the Fully Coupled Static Response of Variable Stiffness Composite Beams, Applied Mathematical Modelling (2019), doi: https://doi.org/10.1016/j.apm.2019.12.010

This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Highlights • An analytical solution is proposed for the fully coupled variable stiffness composite beams. • Direct integration technique and series expansion representation are employed to obtain the closedform solution. • Chebyshev collocation method is also applied as an alternative numerical approach.

1

Analytical Solution for the Fully Coupled Static Response of Variable Stiffness Composite Beams Pedram Khaneh Masjedia,˚, Paul M. Weavera a Bernal

Institute, School of Engineering, University of Limerick, Limerick, Ireland

Abstract An analytical solution is presented for the 3D static response of variable stiffness non-uniform composite beams. Based on Euler-Bernoulli theory, a set of governing differential equations are obtained, in which four degrees of freedom are fully coupled. For the variable stiffness beam, the governing field equations have variable coefficients reflecting the stiffness variation along the beam. Using the direct integration technique, the general analytical solution is derived in the integral form and the closed-form expressions of the obtained solutions are presented employing a series expansion approximation. The series expansion representation enables the proposed approach to be applicable for variable stiffness composite beams with arbitrary spanwise variation of properties. As an alternative solution, the Chebyshev collocation method is applied to the proposed formulation to verify the results obtained from the analytical solution. A number of variable stiffness composite beams made by fibre steering with various boundary conditions and stacking sequences are considered as the test cases. The static response are presented based on the analytical solution and Chebyshev collocation method and excellent agreement is observed for all test cases. The proposed model presents a reliable and efficient approach for capturing the complicated behaviour of variable stiffness nonuniform composite beams. Keywords: Variable Stiffness, Composite Beam, Analytical Solution, Static Response

1. Introduction Variable stiffness non-uniform structures are used as real-life engineering designs in which multiple design criteria due to variable and state-dependent characteristics are required. Applications include modern high aspect ratio aircraft wings, helicopter rotor blades and wind turbine rotor blades which are considered as slender beam-like structures. Therefore, the development of methodologies capable of capturing the structural behaviour of variable stiffness non-uniform beams accurately and efficiently is of great importance for ˚ Corresponding

author: Phone:+353 834603273 Email addresses: [email protected]; [email protected] (Pedram Khaneh Masjedi), [email protected] (Paul M. Weaver)

Preprint submitted to Applied Mathematical Modelling

December 13, 2019

the engineering community. In this context, analytical solutions play an important role for the development of efficient preliminary design tools and for the benchmarking purposes. In the case of beam-like structures, the variable stiffness properties can be achieved through the spanwise distribution of material properties (e.g. by gradation of material or by using fibre steering), the change of geometry (e.g. by tapering the cross-section) or even by using both strategies. Based on this premise, tapered beams, axially or bidirectional graded beams and beams made by fibre steering are considered as variable stiffness beams. Early attempts to address the analytical solution of variable stiffness non-uniform beams can be traced back to the classical works of Hetenyi [1], Jones et al. [2], Jones [3], Gundersen and Iversen [4], Boley [5]. In the context of isotropic beams, Medwadowski [6] considered a shear deformable beam and obtained closed-form expressions for the plane deflection of a non-uniform beam with linear variation of the properties along the beam span. Rao and Spyrakos [7] proposed a Fourier-polynomial series for the analytical solution of linear partial or ordinary differential equations with variable coefficients which was applied for the static analysis of a tapered beam of variable rectangular cross section. Lee et al. [8] assuming that the bending stiffness is second-order differentiable, obtained the exact solution for the deflection of plane Euler-Bernoulli beams. Later, the same approach was used for plane Timoshenko beams by Lee and Kuo [9] to obtain the static deflection. Eisenberger [10] presented the exact stiffness matrix for non-uniform beams including axial, torsional and bending stiffness. Friedman and Kosmatka [11] derived the exact axial, torsional and stiffness matrices for arbitrary non-uniform beam finite elements. By transforming the fourth-order differential equations with variable coefficients into fourth-order differential equations with constant coefficients, Romano [12] presented closed-form solutions for the bending of Timoshenko beams with linearly varying depth and width. Baker [13] proposed closed-form series solutions for the small inplane deflection of non-prismatic beams. Wang et al. [14] obtained the relationships concerning the bending solution between non-uniform Timoshenko and Euler-Bernoulli beams. Al-Gahtani and Khan [15] used the boundary integral method to obtain an analytical solution for the bending of non-prismatic Euler-Bernoulli beams. Closed-form exact expressions were provided for the two cases of linear and parabolic span-wise variation of the beam depth. Employing the theory of generalised functions, Yavari et al. [16] studied the bending of Euler-Bernoulli and Timoshenko beams with so called jump discontinuities in the deflection and slope and also with sudden changes in the beam stiffness properties. Hodges et al. [17], using the variationalasymptotic method, obtained the exact expressions for the strain energy and the stiffness properties of linearly tapered strip-beams and then they presented the exact solutions for the classical problems of pure extension, pure bending and flexure. Attarnejad et al. [18] derived the exact shape functions for arbitrarily tapered Timoshenko beams, based on the concept of basic displacement functions. Balduzzi et al. [19] proposed a model for non-prismatic planar beams based on Timoshenko theory and considered the effects of non-prismatic geometry. They also presented the analytical solution for the simple case of a tapered beam. This approach was further extended for multilayer Timoshenko beams by Balduzzi et al. [20]. 3

In the context of non-uniform variable stiffness anisotropic beams, Shahba et al. [21] presented exact shape functions for axially functionally graded plane beams for the static, stability and free vibration analysis. Nguyen et al. [22] obtained closed-form analytical solution for the static response of tapered axially functionally graded beams using the power series methodology. Pydah and Sabale [23] presented analytical solution for the in-plane bending of circular bidirectional functionally graded beams based on Euler-Bernoulli theory. However, the analytical results are limited to the case of statically determinate beams. Pydah and Batra [24] analytically solved the static deflection of circular bidirectional graded shear deformable beams including a logarithmic function. More recently, based on the variational-asymptotic method, Sachdeva and Padhee [25] obtained closed-form analytical solutions to predict the nonlinear mechanical behaviour of bidirectional functionally graded cylindrical beams. Soltani and Asgarian [26] evaluated the static and buckling stiffnesses exactly for axially functionally graded Timoshenko beams. One can also find various numerical solutions for the variable stiffness beam problem. To name just a few, Karamanlı [27] considered the in-plane static deflection of bidirectional functionally graded beams using Symmetric Smoothed Particle Hydrodynamics method (SSPH). Using a generalized differential quadrature method (DQM), Li et al. [28] solved the bending, buckling and vibration problems of plane axially FG beams. Xie et al. [29] presented an integrated collocation method for the static and free vibration analysis of plane Euler-Bernoulli beams with axially variable cross section, modulus of elasticity, and mass density. G¨ unay and Timarci [30] studied the static behaviour of variable stiffness thin-walled laminated composite beams with closed cross-sections made by fibre steering. A finite element method is used as the numerical solution of the displacements and rotations. Using the Ritz technique, Ai and Weaver [31] investigated the effects of a combination of geometric taper and variable stiffness of the core on the static response of sandwich beams. Rajasekaran and Bakhshi Khaniki [32] studied the static and dynamic behaviour of nonuniform size-dependent axially functionally graded beams using finite element method. Zappino et al. [33] considered the static analysis of tapered thin-walled composite beams based on a refined one-dimensional unified formulation using finite element method. Macquart et al. [34] presented enhanced beam elements for variable stiffness beams to capture the static behaviour more accurately. While the exact analytical solutions for the static deflection of fully coupled beams with constant stiffness are presented by Masjedi et al. [35] and Doeva et al. [36], according to the author’s best knowledge no closed-form analytical solution exists in the literature for the 3D static response of fully coupled variable stiffness composite beams. Thus, the main purpose of this paper is to present a formulation for fully coupled variable stiffness composite beams and to obtain an analytical closed-form solution for the 3D static response. Alternatively, the Chebyshev collocation method, which has been shown to be efficient and accurate in beam problems (Masjedi and Ovesy [37, 38], Masjedi and Maheri [39], Masjedi et al. [40]), is also applied in order to validate the analytical results. Various elastic coupling mechanisms introduced by the anisotropy of composite structures are utilised in 4

real-life engineering designs for performance enhancement. For example, axial(extension)-twist is proposed for passive blade twist control by using the high level centrifugal forces to achieve optimum blade twist distribution in order to improve aeroelastic performance of rotorcrafts [41] or using out-of-plane bend-twist coupling is well-known for aeroelastic tailoring of fixed wings [42]. Thus, to be able to capture all possible coupled behaviour of variable stiffness composite beams, in the mathematical formulation presented in this work all four degrees of freedom, including in-plane, out-of-plane and axial displacements, and twist are considered and are fully coupled. While in the present work all coupling terms are considered, simpler solutions already exist in which coupling between bend and twist occurs. Possibly the most thorough of which is due to Kollar and Springer [43](Section 6.2.2) and also Stodieck et al. [44]. The former considers effective decoupled bending from twist but assumes transverse bending moment is non-zero, whilst the latter considers that transverse curvature is zero and gives a solution for decoupled twist deformation. Both works highlight coupling effects of bending moments and torque and decouple their contributions effectively. Indeed the latter authors show that bending moment can augment or indeed diminish effective torque depending on its direction and the sign of bend twist coupling stiffness. In light of these works, our solutions can be considered as generalised versions of them. It is also noted that, in order to keep the formulation as general as possible the stiffness matrix entries are expressed in terms of the engineering constants. Consequently, the proposed method can be used for the analysis of various types of variable stiffness non-uniform composite beams. For example, the proposed approach can be employed to retrieve 3D deflections of simpler cases of axially or bidirectional functionally graded beams, for which, due to constitutive relations, some coupling terms are not present (and as a result the number of coupled differential equations reduces). While the proposed solution is 3D and closed-form, the available solutions for functionally graded beams in the literature, such in references [21], [22], [32] and [26] are neither closed-form nor 3D. However, the focus of current work is on the variable stiffness composite beams made by fibre steering. The content of this paper is outlined as follows: In Section 2, the governing equations of a fully coupled variable stiffness composite beam are presented. In Section 3, a general exact solution is derived in an integral form using a direct integration technique. In Section 4, in order to obtain the closed-form expressions, series expansion representation is employed. In Section 5, the Chebyshev collocation method is introduced as an alternative numerical solution. In Section 6, several benchmark test problems are considered and the results obtained from the analytical approach are verified by those obtained from the Chebyshev collocation method. Some conclusions and remarks are drawn in the last section.

2. Governing equations In this section, using the principle of virtual work, a linear mathematical formulation is obtained for a fully coupled variable stiffness composite beam based on Euler-Bernoulli theory. It is noted that in the 5

derivation of the beam model it is assumed that the material is linearly elastic, displacements, rotations and strain are small and the effects of cross-sectional warping are neglected. Assume a straight beam of length ` for which the x coordinate is along the beam longitudinal axis and the coordinates y and z define the cross-sectional planes. The principle of virtual work is expressed as: ż`´ ¯ δWint ´ δWext dx “ 0,

(2.1)

ż`

(2.2)

0

where δWint and δWext are the variations of internal and external works respectively. Variation of internal work of the beam can be expressed as:

0

δWint dx “

ż`

δεT Sεdx,

0

where vector of strains ε and stiffness matrix S are written as:

ε “ ru1 »

ϕ1

´ w2

v 2 sT ,

EA pxq SET pxq SEF pxq SEL pxq — — — SET pxq GJ pxq SF T pxq SLT pxq S“— — — SEF pxq SF T pxq EIy pxq SF L pxq – SEL pxq SLT pxq SF L pxq EIz pxq

fi

ffi ffi ffi ffi , ffi ffi fl

(2.3a)

(2.3b)

where u, ϕ, w and v are the axial elongation, twist, out-of-plane displacement and in-plane displacement, respectively, pq1 denotes the derivative with respect to x, EA pxq is the distribution of extensional stiffness,

GJ pxq is the distribution of twist stiffness, EIy pxq is the distribution of out-of-plane bending stiffness,

EIz pxq is the distribution of in-plane bending stiffness, SET pxq is the distribution of coupling between

axial elongation and twist, SEF pxq is the distribution of coupling between out-of-plane bending and axial

elongation, SEL pxq is the distribution of coupling between in-plane bending and axial elongation, SF T is

the distribution of coupling between out-of-plane bending and twist, SLT pxq is the distribution of coupling

between in-plane bending and twist, and SF L pxq is the distribution of coupling between out-of-plane and in-plane bending. Variation of external work can be expressed as: ż` 0

δWext dx “

ż`

δU QT dx,

(2.4)

0

where displacement vector U and vector of external loads Q are written as:

U “ ru Q “ rqx

ϕ

vsT ,

w

qϕ 6

qz

qy sT ,

(2.5a) (2.5b)

where qx , qy , qz are the distributed loads in the x, y, z directions respectively, and qϕ is the distributed torque. Now Eqs. (2.2) and (2.4) can be written as: ż` 0

δWint dx “

ż` 0

` ˘ tδu1 EA pxq u1 ` SET pxq ϕ1 ` SEF pxq w2 ` SEL pxq v 2

` ˘ `δϕ1 SET pxq u1 ` GJ pxq ϕ1 ` SF T pxq w2 ` SLT pxq v 2 ˘ ` ´δw2 SEF pxq u1 ` SF T pxq ϕ1 ` EIy pxq w2 ` SF L pxq v 2 ` ˘ `δv 2 SEL pxq u1 ` SLT pxq ϕ1 ` SF L pxq w2 ` EIz pxq v 2 udx,

and ż` 0

δWext dx “

ż` 0

pδuqx ` δϕqϕ ` δwqz ` δvqy q dx.

(2.6)

(2.7)

Substituting Eqs. (2.6) and (2.7) into Eq. (2.1), the set of governing differential equations with variable coefficients are obtained as:

` ˘1 δu : ´EA pxq u1 ´ SET pxq ϕ1 ` SEF pxq w2 ´ SEL pxq v 2 “ qx ` ˘1 δϕ : ´SET pxq u1 ´ GJ pxq ϕ1 ` SF T pxq w2 ´ SLT pxq v,x “ qϕ ` ˘2 δw : ´SEF pxq u1 ´ SF T pxq ϕ1 ` EIy pxq w2 ´ SF L pxq v 2 “ qz ` ˘2 δv : SEL pxq u1 ` SLT pxq ϕ1 ´ SF L pxq w2 ` EIz pxq v 2 “ qy .

(2.8a) (2.8b) (2.8c) (2.8d)

At x “ 0 and x “ `, boundary conditions can be expressed as follows:

u“0

or

EA pxq u1 ` SET pxq ϕ1 ´ SEF pxq w2 ` SEL pxq v 2 “ Fx

(2.9a)

ϕ“0

or

SET pxq u1 ` GJ pxq ϕ1 ´ SF T pxq w2 ` SLT pxq v 2 “ Mx

(2.9b)

w1 “ 0

or

v1 “ 0

or

w“0

or

v“0

or

´ SEF pxq u1 ´ SF T pxq ϕ1 ` EIy pxq w2 ´ SF L pxq v 2 “ ´My SEL pxq u1 ` SLT pxq ϕ1 ´ SF L pxq w2 ` EIz pxq v 2 “ Mz ` ˘1 SEF pxq u1 ` SF T pxq ϕ1 ´ EIy pxq w2 ` SF L pxq v 2 “ Fz ` ˘1 ´SEL pxq u1 ´ SLT pxq ϕ1 ` SF L pxq w2 ´ EIz pxq v 2 “ Fy ,

(2.9c) (2.9d) (2.9e) (2.9f)

where Fx , Fy , Fz are the tip loads in the x, y and z direction respectively, Mx is the tip torque, and My , Mz are the tip moments about the y and z axes respectively. 7

x pxq be a block matrix such In order to express Eqs. (2.8) and (2.9) in a compact matrix form, let K

that:

»

x pxq “ – K

where

»

fl , B T pxq D pxq

SEF pxq

´SEL pxq

»

EIy pxq

D pxq “ –

SF T pxq

´SF L pxq

K pxq “ –

KA pxq

GJ pxq

(2.10)

fi

»

SET pxq

»

´1

fi

SET pxq

B pxq “ –

x such that K pxq “ K

B pxq

EA pxq

A pxq “ –

Now consider the matrix

A pxq

fl ,

fi

fl ,

´SLT pxq

´SF L pxq EIz pxq

KB pxq

(2.11a)

fi

fl .

fi

fl , KB T pxq KD pxq

(2.11b)

(2.11c)

(2.12)

x pxq. Entries of K pxq can be determined pxq, i.e. K pxq is an inverse of matrix K

using the blockwise inversion formula (Bernstein [45]):

” ı´1 B T pxq A´1 pxq , KA pxq “ A´1 pxq ` A´1 pxq B pxq D pxq ´ B T pxq A´1 pxq B pxq ” ı´1 KB pxq “ ´A´1 pxq B pxq D pxq ´ B T pxq A´1 pxq B pxq , ” ı´1 KB T pxq “ ´ D pxq ´ B T pxq A´1 pxq B pxq B T pxq A´1 pxq , ı´1 ” . KD pxq “ D pxq ´ B T pxq A´1 pxq B pxq

(2.13a) (2.13b) (2.13c) (2.13d)

Using matrix (2.10), the governing equations (2.8) and boundary conditions (2.9) are written in a compact matrix form: » fi2 ˛1 » fi » fi1 u w q ˝´A pxq – fl ` B pxq – fl ‚ “ – x fl , ϕ v qϕ ¨ » fi1 » fi2 ˛2 » fi q u w ˝´B T pxq – fl ` D pxq – fl ‚ “ – z fl . ϕ v qy ¨

8

(2.14a)

(2.14b)

» fi » fi u 0 – fl “ – fl , ϕ 0 » fi » fi w 0 – fl “ – fl , v 0 » fi1 » fi w 0 – fl “ – fl . v 0

» fi1 » fi2 » fi u w Fx A pxq – fl ´ B pxq – fl “ – fl , ϕ v Mx » fi1 » fi2 » fi u w ´M y fl , ´B T pxq – fl ` D pxq – fl “ – ϕ v Mz ¨ » fi1 » fi2 ˛1 » fi u w F ˝B T pxq – fl ´ D pxq – fl ‚ “ – z fl . ϕ v Fy

(2.15a)

(2.15b)

(2.15c)

(2.16a)

(2.16b)

(2.16c)

In the following section, the procedure of obtaining the closed-form analytical solution of the system of differential Eqs. (2.14) with variable coefficients is presented. 3. General solution To obtain the exact general solution of Eqs. (2.14), Eq. (2.14a) is integrated and rearranged. It is assumed that the loads at the right hand side of Eqs. (2.14) are uniformly distributed, thus the first derivative of the vector ru

ϕsT can be obtained as follows:

» fi1 » fi2 » fi u w q – fl “ A´1 pxq B pxq – fl ´ xA´1 pxq – x fl ´ A´1 pxq C 1 . ϕ v qϕ

(3.1)

Substituting Eq. (3.1) into Eq. (2.14b) and rearranging it, we obtain:

» fi2 » fi ˛2 » fi ¨ ´ ¯ w q q ˝ D pxq ´ B T pxq A´1 pxq B pxq – fl ` xB T pxq A´1 pxq – x fl ` B T pxq A´1 pxq C 1 ‚ “ – z fl . qϕ qy v (3.2) By successive integration of Eq. (3.2), the following expressions are derived: ¨ » fi2 » fi ˛1 » fi ´ ¯ w q q x ˝ D pxq ´ B T pxq A´1 pxq B pxq – fl ` xB T pxq A´1 pxq – fl ` B T pxq A´1 pxq C 1 ‚ “ x – z fl ` C 2 , v qϕ qy (3.3) 9

» fi » fi » fi2 ´ ¯ w 2 q q x x – z fl`xC 2 `C 3 . D pxq ´ B T pxq A´1 pxq B pxq – fl `xB T pxq A´1 pxq – fl`B T pxq A´1 pxq C 1 “ 2 qϕ qy v (3.4) Rearranging Eq. (3.4) and using equalities introduced in Eq. (2.13), the second derivative of the vector rw

vsT is obtained: » fi2 ¨ » fi ˛ ¨ » fi ˛ 2 w qz q x x – fl “ KD pxq ˝ – fl ` xC 2 ` C 3 ‚` KB T pxq ˝x – fl ` C 1 ‚. 2 qy v qϕ

(3.5)

Integrating Eq. (3.5) twice, the solution for the out-of-plane displacement w and in-plane displacement v can be obtained: » fi » fi1 ż 2 ż ż qz w x – – fl “ fl KD pxq dx ` xKD pxq dx C 2 ` KD pxq dx C 3 2 qy v » fi ż ż qx T – fl ` xKB pxq dx ` KB T pxq dx C 1 ` C 4 , (3.6) qϕ » fi » fi ij ij ij 2 qz w x – fl – fl “ KD pxq dx dx ` xKD pxq dx dx C 2 ` KD pxq dx dx C 3 2 qy v » fi ij ij qx T – fl ` xKB pxq dx dx ` KB T pxq dx dx C 1 ` x C 4 ` C 5 , (3.7) qϕ

where C i , i “ 1, 2, ..., 5, are the vectors of unknown integrating constants to be determined.

Substituting Eq. (3.5) into Eq. (3.1) and using Eq. (2.13), the the first derivative of the vector ru

ϕsT

is obtained: » fi ˛ » fi1 ¨ » fi 2 qz u qx x – fl “ ´ ˝xKA pxq – fl ` KA pxq C 1 ` KB pxq – fl ` xKB pxq C 2 ` KB pxq C 3 ‚. 2 ϕ qϕ qy

(3.8)

Integrating Eq. (3.8), the exact solution for the displacement u and twist ϕ can be derived:

» fi ¨ » fi ż ż q u – fl “ ´ ˝ xKA pxq dx – x fl ` KA pxq dx C 1 ϕ qϕ » fi ˛ ż 2 ż ż qz x ` KB pxq dx – fl ` xKB pxq dx C 2 ` KB pxq dx C 3 ‚` C 6 , (3.9) 2 qy 10

where C i , i “ 1, 2, 3, 6, are the vectors of unknown coefficients to be determined. Substituting Eqs. (3.8) and (3.5) into Eqs. (2.16a), (2.16b) and (2.16c), the natural boundary conditions can be simplified as: »

´x –

qx qϕ

fi

»

fl ´ C 1 “ –

Fx Mx

fi

fl ,

» fi » fi ´My x2 –qz fl fl , ` xC 2 ` C 3 “ – 2 qy Mz » fi » fi qz Fz ´x – fl ´ C 2 “ – fl . qy Fy

(3.10a)

(3.10b)

(3.10c)

Eqs. (3.9) and (3.7) are the general solutions for the 3D static deflection of a fully coupled variable stiffness composite beam subject to uniformly distributed and tip loads. These general solutions have six vectors of unknown constants to be determined. Once the boundary conditions are given, these unknown vectors can be obtained and the particular solutions can be found. Eqs. (3.9) and (3.7) are expressed in integral form, however if the functions inside the integrals are integrable, exact closed-form expressions can be retrieved. In the following section a series expansion representation is proposed to obtain the closed-form solutions in the case of arbitrary span-wise distribution of beam stiffness properties.

4. Series expansion representation For arbitrarily span-wise distributed stiffness properties, it is impossible to guarantee that the exact closed-form expressions can be derived from Eqs. (3.9) and (3.7). However, it is shown in this section that a series expansion representation can be employed to obtain the analytical closed-form expressions for arbitrary smooth variation of beam stiffness properties. Toward this goal, using a series expansion representation the functions inside the integrals in Eqs. (3.9) and (3.7) are expressed as:

0 KA pxq “ WA φ pxq ,

1 xKA pxq “ WA φ pxq ,

0 KB pxq “ WB φ pxq ,

1 xKB pxq “ WB φ pxq ,

0 φ pxq , KD pxq “ WD

1 xKD pxq “ WD φ pxq ,

(4.1a) x2 2 KB pxq “ WB φ pxq , 2 2 x 2 KD pxq “ WD φ pxq , 2

(4.1b) (4.1c)

where, Wαβ , α “ A, B, D and β “ 0, 1, 2 are block matrices containing the vectors of unknown weight

coefficients and φ pxq is a known function of x. Substituting Eqs. (4.1) into Eqs. (3.7) and (3.9), we obtain:

11

» fi » fi ij ij ij w qz 1 0 – fl “ W 2 – fl φ pxq dx dx ` WD φ pxq dx dx C 2 ` WD φ pxq dx dx C 3 D v qy » fi ij ij qx T 0T 1 – fl φ pxq dx dx C 1 ` x C 4 ` C 5 , (4.2) φ pxq dx dx ` WB ` WB qϕ » fi ¨ » fi ż ż u q – fl “ ´ ˝W 1 φ pxq dx – x fl ` W 0 φ pxq dx C 1 A A ϕ qϕ » fi ˛ ż ż ż q z 2 1 0 φ pxq dx – fl ` WB φ pxq dx C 2 ` WB φ pxq dx C 3 ‚` C 6 , (4.3) `WB qy

where the matrices of weight coefficients are expressed as follows: »

0 WA “–

fi

0 wA 11

0 wA 12

0 wA 21

0 wA 22

0 wB 11

0 wB 12

0 wB 21

0 wB 22

0 wD 11

0 wD 12

0 wD 21

0 wD 22

»

0 WB “–

»

0 WD “–

»

1 wA 11

1 wA 12

1 wA 21

1 wA 22

1 wB 11

1 wB 12

1 wB 21

1 wB 22

1 wD 11

1 wD 12

1 wD 21

1 wD 22

fl ,

1 WA “–

fl ,

1 WB “–

fl ,

1 WD “–

fi

fi

»

»

fi

fl ,

(4.4a)

fi

»

2 wB 12

2 wB 21

2 wB 22

2 wD 11

2 wD 12

2 wD 21

2 wD 22

fl ,

2 WB “–

fl ,

2 WD “–

fi

fi

2 wB 11

»

fl

(4.4b)

fl ,

(4.4c)

fi

β where wα , α “ A, B, D, β “ 0, 1, 2 and i, j “ 1, 2 are the vectors of unknown weighting coefficients to be ij

calculated following a standard curve fitting algorithm. By substituting the expression of φ pxq, it is possible to obtain closed-form expressions from Eqs. (4.2) and “ ‰T (4.3). Using a power series representation, i.e. φ pxq “ 1, x, x2 , ¨ ¨ ¨ , xN where N is the highest degree of polynomials, we obtain:

» fi » fi w q – fl “ W 2 Φii pxq – z fl ` W 1 Φii pxq C 2 ` W 0 Φii pxq C 3 D D D v qy » `

1 T ii WB Φ

12

pxq –

qx qϕ

fi

fl ` W 0 T Φii pxq C 1 ` x C 4 ` C 5 , (4.5) B

» fi ¨ » fi u q – fl “ ´ ˝W 1 Φi pxq – x fl ` W 0 Φi pxq C 1 A A ϕ qϕ

˛ » fi q z 0 i 1 i 2 i Φ pxq C 3 ‚` C 6 , (4.6) Φ pxq C 2 ` WB Φ pxq – fl ` WB `WB qy

where,

i

Φii pxq “

ij

Φ pxq “

ż

φ pxq dx dx “

„ φ pxq dx “ x,

„

x2 , 2

x3 , 6

x2 , 2

¨¨¨ ,

¨¨¨ ,

xN `1 N `1

xN `2 pN ` 1q pN ` 2q

T

T

,

(4.7a)

.

(4.7b)

Eqs. (4.5) and (4.6) are the closed-form analytical solutions for the 3D deflection of variable stiffness composite beams. Applying boundary conditions, it is possible to determine constants of integration Ci , i “ 1, 2, ¨ ¨ ¨ , 6. Sample calculations are given in Appendix A for the case of cantilever and simply supported beams. It is noted that 20 terms in the power series were sufficient to obtain convergence. It is also worth mentioning that the proposed approach is not limited to power series and any type of appropriate series can be used in this approach.

5. Chebyshev collocation method The Chebyshev collocation method (CCM) is applied as an alternative numerical solution for Eq. (2.8). For this purpose, unknown variables u, v, w and ϕ are discretised using the Chebyshev polynomials, and the Chebyshev points are employed as the collocation points. The Chebyshev polynomial of the first kind in x and of degree n is defined as (Mason and Handscomb [46]): Tn pxq “ cospnθq when

x “ cos θ,

´1 ď x ď `1,

n “ 0, 1, . . .

(5.1)

It is noted that for any arbitrary range a ď x ď b the Chebyshev polynomials can be shifted by replacing the independent variable x in Eq.(5.1) by: x“

2 b`a x´ . b´a b´a

For an arbitrary interval a ď x ď b the transformed Chebyshev points are: ˆ ˙ 1 1 2i ´ 1 xi “ pa ` bq ´ pb ´ aq cos π , i “ 1, 2, . . . , N ` 1, 2 2 2N where N is the highest degree of the Chebyshev polynomials. 13

(5.2)

(5.3)

Using the Chebyshev polynomials, unknown variables are discretised as: u“

N ÿ

i“0

ai Ti pxq,

ϕ“

N ÿ

i“0

bi Ti pxq,

w“

N ÿ

i“0

ci Ti pxq,

v“

N ÿ

i“0

di Ti pxq.

(5.4)

Eqs. (5.4) are substituted into Eqs. (2.8) and the residuals at the Chebyshev points are set equal to zero. It is worth mentioning that in order to have a well-posed system of 4 ˆ pN ` 1q equations in conjunction with the boundary conditions (2.9), the same approach proposed by Masjedi et al. [35] is followed herein: in the case of u and ϕ, the equations associated with the residuals at the first and the last Chebyshev points are eliminated, whereas in the case of w and v, those equations associated with the first two and the last two Chebyshev points from Eq. (5.3) are discarded systematically. Finally, unknown coefficients ai , bi , ci and di where i “ 0, 1, 2, ..., N are determined solving the system of linear equations. Using a discretisation of N “ 20 gives converged results and all numerical results from the Chebyshev collocation method in the following sections are obtained on this premise.

6. Numerical results In this section for verifying the current approach a constant stiffness composite beam is considered first and the numerical results are compared to exact results available in the literature. Then a number of variable stiffness composite beam examples under the action of tip load or distributed load are presented for different boundary conditions. A slender composite beam made by fibre steering with a rectangular cross-section is considered for all numerical samples. The material and geometric properties used in the samples are listed in Table 1. It is noted that the distribution of stiffness properties is calculated as functions of beam span-wise coordinate ”x” using the closed-form expressions presented by Yu and Hodges [47](see Appendix B). To investigate the effects of different coupling terms on the deflection of variable stiffness composite beams, the following stacking sequences are considered: 1. Symmetric rθ3 ss with bend-twist coupling, 2. Antisymmetric rθ3 { ´ θ3 s with extension-twist coupling, 3. Unsymmetric rθ3 {03 s with bend-twist, extension-twist and extension-bend coupling. where the fibre angle θ varies linearly along the beam span as: θ pxq “

θt ´θr ` x

` θr , where θr is the fibre

angle at x “ 0 and θt is the fibre angle at x “ `. In this work it is assumed that θr “ 0˝ and θt “ 45˝ . The span-wise distribution of stiffness properties for each stacking sequence is depicted in Figs.(1), (2) and (3). In the following sections numerical results are presented for the static deformations of variable stiffness composite beams with various boundary conditions and stacking sequences. The numerical results are presented using 10 decimal places for the comparison and benchmarking purposes, even though they are not considered to be meaningful in real-life problems. 14

Table 1: Beam properties

Length

1

m

Width

0.01

m

Number of Layers

6

Ply Thickness

0.000125

m

E11

135.64

GP a

E22

10.14

GP a

G12

5.86

GP a

ν12

0.29

50

0

45

-0.02

40

-0.04

35

-0.06

30

-0.08

25

-0.1

20

-0.12

15

-0.14

10

-0.16

5

-0.18

0

0

0.2

0.4

0.6

0.8

-0.2

1

(a)

0

0.2

0.4

0.6

(b)

Figure 1: Distribution of Stiffness Properties, Symmetric rθ3 ss

15

0.8

1

0

30

-10

25

-20 20

-30 -40

15

-50

10

-60 5

0

-70

0

0.2

0.4

0.6

0.8

-80

1

0

0.2

0.4

(a)

0.6

0.8

1

0.6

0.8

1

(b)

Figure 2: Distribution of Stiffness Properties, Asymmetric rθ3 { ´ θ3 s

16

0 -10

14

-20 12

-30

10

-40

8

-50 -60

6

-70 4 2

-80 0

0.2

0.4

0.6

0.8

-90

1

(a)

0

0.2

0.4

(b)

Figure 3: Distribution of Stiffness Properties, Unsymmetric rθ3 {03 s

16

6.1. Cantilevered constant stiffness composite beam under tip load A constant stiffness cantilevered composite beam made from unidirectional straight fibres under the action of a tip load is considered in this section. Numerical results obtained from the proposed approach is compared to the exact solution presented by Hodges [48] which is based on an intrinsic formulation (IF) for a cantilever beam subject to tip shear forces and moments. Four different stacking sequences are considered herein: 1. Symmetric r453 ss with bend-twist coupling, 2. Antisymmetric r453 { ´ 453 s with extension-twist coupling, 3. Cross-ply r03 {903 s with extension-bend coupling, 4. Unsymmetric r603 {303 s with bend-twist, extension-twist and extension-bend coupling. The stiffness matrices obtained from closed-form expressions presented in Appendix B, are given for each stacking sequence in Table (6). Numerical results for the tip deflections are presented in Tables (2) through (5). In all cases, excellent agreement is observed between different sets of results. Table 2: Tip deformation of bend-twist coupled cantilevered beam (r453 ss ) under tip shear load (Fz “ 0.001N )

upmq

ϕpradq

wpmq

vpmq

Exact (IF)

0

-0.0324433326

0.0645604090

0

Analytical

0

-0.0324433326

0.0645604090

0

CCM

0

-0.0324433326

0.0645604090

0

Table 3: Tip deformation of extension-twist coupled cantilevered beam (r453 { ´ 453 s) under tip twist (Mx “ 0.001N.m)

upmq

ϕpradq

wpmq

vpmq

Exact (IF)

0.0000070724

0.0455822430

0

0

Analytical

0.0000070724

0.0455822430

0

0

CCM

0.0000070724

0.0455822430

0

0

17

Table 4: Tip deformation of extension-bend coupled cantilevered beam (r03 {903 s) under tip shear load (Fz “ 0.001N )

upmq

ϕpradq

wpmq

vpmq

Exact (IF)

0.0000070724

0

0.0292097534

0

Analytical

0.0000070724

0

0.0292097534

0

CCM

0.0000070724

0

0.0292097534

0

Table 5: Tip deformation of the unsymmetric cantilevered beam (r603 {303 s) under tip shear load (Fz “ 0.001N )

upmq

ϕpradq

wpmq

vpmq

Exact (IF)

-0.0000048161

-0.0292043302

0.0520136990

0

Analytical

-0.0000048161

-0.0292043302

0.0520136990

0

CCM

-0.0000048161

-0.0292043302

0.0520136990

0

18

Table 6: Stiffness matrix for different stacking sequences

Stacking Sequence

r453 ss

r453 { ´ 453 s

r03 {903 s

r603 {303 s

» — — — — — — –

»

S 110155.0

0

0

0

0.0176

0

´0.0059

´0.0059

0

0

152399.0

— — —´23.6451 — — — 0 – 0 »

548117.8

— — — 0 — — — 88.4751 – 0

»

163771.4

— — — 15.7916 — — —´16.0203 – 0

0

0.0071 0

´23.6451

0

0.0256

0

0

0.0061

0

0

ffi ffi ffi ffi ffi 0 ffi fl 0.9179 0

ffi ffi ffi ffi ffi 0 ffi fl 1.2700 0

fi

88.4751

0

0.0082

0

0

0

0.0257

0

0

0

4.5676

0.0167 ´0.0072 0

fi

0

0

15.7916

fi

´16.0203

ffi ffi ffi ffi ffi ffi fl 0

´0.0072

0

0.0101

0

0

1.3648

fi

ffi ffi ffi ffi ffi ffi fl

6.2. Cantilevered variable stiffness composite beam under tip and distributed loads A cantilevered composite beam subject to concentrated tip and distributed load is considered in this section. Numerical results for the tip deflections are presented from the analytical solution and the Chebyshev collocation method (CCM). Tables (7) through (9) for the case of tip load and Tables (10) through (12) for the case of distributed load, show the tip deflections for different stacking sequences. Figs. (4) and (5) show the composite beam deformations for different stacking sequences based on the analytical solution and also CCM. Deformations with zero value are not shown. A very good agreement is observed between the two sets of results. The effects of different coupling terms is clearly observed on the deformation behaviour of the variable stiffness composite beam. 19

Table 7: Tip deformation of bend-twist coupled cantilevered beam (rθ3 ss ) under tip shear load (Fz “ 0.005N )

upmq

ϕpradq

wpmq

vpmq

Analytical

0

-0.1092212560

0.0741482245

0

CCM

0

-0.1092212560

0.0741482245

0

Table 8: Tip deformation of extension-twist coupled cantilevered beam (rθ3 { ´ θ3 s) under tip twist (Mx “ 0.001N.m)

upmq

ϕpradq

wpmq

vpmq

Analytical

0.0000069634

0.0753015244

0

0

CCM

0.0000069634

0.0753015100

0

0

Table 9: Tip deformation of the unsymmetric cantilevered beam (rθ3 {03 s) under tip shear load (Fz “ 0.005N )

upmq

ϕpradq

wpmq

vpmq

Analytical

-0.0000072537

-0.0473716972

0.0500256940

0

CCM

-0.0000072537

-0.0473717014

0.0500256984

0

20

0.08

0.09

0.06

0.08

0.04

0.07

0.02

0.06

0

0.05

-0.02

0.04

-0.04

0.03

-0.06

0.02

-0.08

0.01

-0.1

0

-0.12

0

0.2

0.4

0.6

0.8

-0.01

1

0

(a) rθ3 ss , Fz “ 0.005N

0.2

0.4

0.6

0.8

(b) rθ3 { ´ θ3 s, Mx “ 0.001N.m

0.08 0.06 0.04 0.02 0 -0.02 -0.04 -0.06 -0.08

0

0.2

0.4

0.6

0.8

1

(c) rθ3 {03 s, Fz “ 0.005N Figure 4: Deformation of a cantilevered beam under the action of tip load for different stacking sequences

21

1

Table 10: Tip deformation of bend-twist coupled cantilevered beam (rθ3 ss ) under uniformly distributed load (qz “ 0.005N {m)

upmq

ϕpradq

wpmq

vpmq

Analytical

0

-0.0297228069

0.0232070638

0

CCM

0

-0.0297228069

0.0232070638

0

Table 11: Tip deformation of extension-twist coupled cantilevered beam (rθ3 { ´ θ3 s) under uniformly distributed load (qϕ “ 0.001N )

upmq

ϕpradq

wpmq

vpmq

Analytical

0.0000030192

0.0449962310

0

0

CCM

0.0000030192

0.0449962380

0

0

Table 12: Tip deformation of the unsymmetric cantilevered beam (rθ3 {03 s) under uniformly distributed load (qz “ 0.005N {m)

upmq

ϕpradq

wpmq

vpmq

Analytical

-0.0000015907

-0.0134062849

0.0171096123

0

CCM

-0.0000015907

-0.0134062868

0.0171096143

0

22

0.03 0.05 0.02

0.045 0.04

0.01

0.035 0.03

0

0.025 -0.01

0.02 0.015

-0.02

0.01 0.005

-0.03 0

0.2

0.4

0.6

0.8

0

1

0

0.2

(a) rθ3 ss , qz “ 0.01N {m

0.4

0.6

0.8

1

(b) rθ3 { ´ θ3 s, qϕ “ 0.01N

0.025 0.02 0.015 0.01 0.005 0 -0.005 -0.01 -0.015 -0.02

0

0.2

0.4

0.6

0.8

1

(c) rθ3 {03 s, qz “ 0.01N {m Figure 5: Deformation of a cantilevered beam under the action of uniformly distributed load for different stacking sequences

6.3. Simply supported variable stiffness composite beam under distributed load In this section a variable stiffness composite beam simply supported at both ends subject to uniformly distributed load is considered. The values of the maximum deformations obtained from the analytical 23

solution and the CCM are given in Tables (13) through (15). Figs. (6) depict the deformations of the variable stiffness composite beam for different stacking sequences based on the two methods. In all cases an excellent agreement is observed between the analytical results and those of the CCM. For this problem, the loading and boundary conditions are symmetric, however, due to anisotropy and non-uniform distribution of the beam stiffness properties the symmetry of the deformations is broken. Table 13: Maximum deformation of bend-twist coupled simply supported beam (rθ3 ss ) under uniformly distributed load (qz “ 0.1N {m)

upmq

ϕpradq

wpmq

vpmq

Analytical

0

0.1214569346

0.0735328464

0

CCM

0

0.1214569346

0.0735328464

0

Table 14: Maximum deformation of extension-twist coupled simply supported beam (rθ3 { ´ θ3 s) under uniformly distributed load (qϕ “ 0.01N )

upmq

ϕpradq

wpmq

vpmq

Analytical

0.0000049513

0.0824035794

0

0

CCM

0.0000049513

0.0824035798

0

0

Table 15: Maximum deformation of unsymmetric simply supported beam (rθ3 {03 s) under uniformly distributed load (qz “

0.1N {m)

upmq

ϕpradq

wpmq

vpmq

Analytical

0.0000051278

0.0319742931

0.0428916957

0

CCM

0.0000051278

0.0319742931

0.0428916958

0

24

0.09 0.1

0.08 0.07

0.05

0.06 0.05

0

0.04 0.03

-0.05

0.02 0.01

-0.1

0 -0.15

0

0.2

0.4

0.6

0.8

-0.01

1

0

0.2

(a) rθ3 ss , qz “ 0.1N {m

0.4

0.6

0.8

1

(b) rθ3 { ´ θ3 s, qϕ “ 0.01N

0.06 0.04 0.02 0 -0.02 -0.04 -0.06

0

0.2

0.4

0.6

0.8

1

(c) rθ3 {03 s, qz “ 0.1N {m Figure 6: Deformation of a simply supported beam at both ends under the action of uniformly distributed load for different stacking sequences

25

6.4. Composite beams with different boundary conditions under combined distributed load In order to show the capabilities of the proposed beam model at capturing 3D static deformations under the action of different loading cases with various boundary conditions, a variable stiffness composite beam subject to a combined distributed load is considered. Uniformly distributed shear loads qy and qz and a distributed torque qϕ are applied to the variable stiffness composite beam with an unsymmetric stacking sequence (rθ3 {03 s). For the cases of simply supported - simply supported, clamped - simply supported and clamped - clamped boundary conditions, the numerical results for the maximum deformations of the beam are given, while for the case of the cantilevered beam the values for the tip deformations are presented. Table (16) shows the results obtained from the analytical solution and the CCM, and the deformations of the variable stiffness composite beams are depicted in Figs. (7) for different boundary conditions. Excellent agreement is observed between the two sets of results for all cases. In the case of symmetric simply supported - simply supported and clamped - clamped boundary conditions, due to the combined loading, non-uniform distribution of stiff properties and anisotropy the symmetry of the deformations is broken. Table 16: Deformation of the unsymmetric beam (rθ3 {03 s) under combined uniformly distributed loads (qy “ 0.1N {m, qz “

0.02N {m, qϕ “ 0.0025N )

B.C. C-F

S-S

C-S

C-C

upmq

ϕpradq

wpmq

vpmq

Analytical

-0.0000018797

0.0804209075

0.0550321640

0.0016945527

CCM

-0.0000018797

0.0804208996

0.0550321725

0.0016945493

Analytical

0.0000135773

0.0280105004

0.0082840413

0.0002243471

CCM

0.0000135773

0.0280105005

0.0082840414

0.0002243471

Analytical

0.0000006974

0.0243149416

0.0028549117

0.0000870663

CCM

0.0000006974

0.0243149416

0.0028549117

0.0000870663

Analytical

0.0000013336

0.0271275318

0.0020362427

0.0000428771

CCM

0.0000013336

0.0271275318

0.0020362427

0.0000428771

26

0.12

0.04

0.1

0.035 0.03

0.08

0.025 0.06

0.02

0.04

0.015 0.01

0.02

0.005 0 0 -0.02

-0.005 0

0.2

0.4

0.6

0.8

1

(a) Cantilevered

0.04

0.035

0.035

0.03

0.03

0.025

0.025

0.02

0.02

0.015

0.015

0.01

0.01

0.005

0.005

0

0 0

0.2

0.4

0.6

0.2

0.4

0.6

0.8

1

(b) Simply Supported - Simply Supported

0.04

-0.005

0

0.8

-0.005

1

(c) Clamped - Simply Supported

0

0.2

0.4

0.6

0.8

1

(d) Clamped - Clamped

Figure 7: Deformation of composite beam with different boundary conditions under the action of combined distributed load (qy “ 0.1N {m, qz “ 0.02N {m, qϕ “ 0.0025N )

27

6.5. Effects of fibre path In order to show the effects of fibre path on the static behaviour of variable stiffness beams, various test cases are considered in this section. A cantilevered beam subject to uniformly distributed loads; qz “ 0.02N {m and qϕ “ 0.0025N is considered. Four different values are assumed for θt ; 15˝ , 30˝ , 45˝ and 60˝

which in turn introduce different distribution of stiffness properties along the beam reference axis. The distribution of the beam stiffness properties is depicted in Figs.(8) and (9) for an unsymmetric lay-up. Figs.(10) show the deformations of beam for various values of θt . By increasing θt from 15˝ to 60˝ it is observed that the twist deflection ϕ decreases which is due to the increase of torsional stiffness. However, while the out-of-plane stiffness decreases continuously by increasing θt from 15˝ to 60˝ , the out-of plane deflection w first decreases from the case of 15˝ to the case of 30˝ and then increases by increasing θt to 45˝ and 60˝ . This behaviour can be attributed to the complex effects of anisotropy. It is also observed from the in-plane deflections that while for the cases of θt “ 15˝ , 30˝ , 45˝ , the beam is under tensile stress this can change to a compressive stress for the case of θt “ 60˝ . This example clearly shows the highly complicated

influence of fibre path on the static behaviour of variable stiffness beams on one hand and the design space offered by fibre steering on the other hand. It is worth mentioning that the results shown in this section are highly influenced by anisotropy, distribution of stiffness properties and boundary conditions and the study of the complicated behaviour of variable stiffness composite beams are almost impossible intuitively. However, the proposed model presents a reliable and efficient method for the future engineering optimisation and design studies for obtaining deeper insights.

28

10.5

14

10 9.5

13

9 12

8.5 8

11

7.5 7

10

6.5 9

6 5.5

0

0.2

0.4

0.6

0.8

1

0

0.2

(a) Axial stiffness

0.4

0.6

0.8

1

0.8

1

(b) Torsional stiffness

46

8

44

7.5

42 40

7

38 6.5

36 34

6

32

5.5

30

5

28 26

0

0.2

0.4

0.6

0.8

1

0

(c) Out-of-plane bending stiffness

0.2

0.4

0.6

(d) In-plane bending stiffness

Figure 8: Distribution of different stiffness properties for different values of θt , Unsymmetric rθ3 {03 s

29

0

0 -2

-10

-4 -20

-6

-30

-8 -10

-40

-12 -50

-14

-60

-16 -18

-70

-20 -80

-22 0

0.2

0.4

0.6

0.8

1

0

0.2

(a) Extension-bend coupling

0.4

0.6

0.8

(b) Extension-twist coupling

0 -10 -20 -30 -40 -50 -60 0

0.2

0.4

0.6

0.8

1

(c) Bend-twist coupling Figure 9: Distribution of different coupling properties for different values of θt , Unsymmetric rθ3 {03 s

30

1

0.03 0.12 0.02 0.1 0.01 0.08 0

0.06

-0.01

0.04

-0.02 -0.03

0.02

0

0.2

0.4

0.6

0.8

0

1

0

0.2

0.4

(a) Axial elongation

0.6

0.8

1

(b) Twist

0.5 0.45 0.4 0.35 0.3 0.25 0.2 0.15 0.1 0.05 0

0

0.2

0.4

0.6

0.8

1

(c) Out-of-plane bending Figure 10: Deformation of variable stiffness cantilevered beam for different values of θt , Unsymmetric rθ3 {03 s

7. Conclusions The closed-form analytical solution for the fully coupled 3D static response of variable stiffness composite beams is presented. The fully coupled governing differential equations with variable coefficients are derived 31

based on Euler-Bernoulli theory and are expressed in a compact matrix form. The entries of the stiffness matrix are denoted by the engineering constants which allows this model to be used for laminated, fibrereinforced or functionally graded beams with arbitrary cross-sectional shapes. The general and closed-form analytical solutions are derived by the method of direct integration and series expansion representation. In addition, in order to verify the obtained results from the analytical solution, the Chebyshev collocation method, which is well-known to be accurate and computationally efficient, is applied to the problem. Numerical results are obtained for a slender beam with rectangular cross-section made by fibre steering, subject to tip and uniformly distributed loads with different boundary conditions. The fibre angles are assumed to vary linearly along the beam span. Three types of stacking sequences showing various combinations of coupling terms are considered to show the effects of anisotropy on the static response. It is shown that the obtained results from the analytical solution and the Chebyshev collocation method are in a very good agreement. The proposed variable stiffness non-uniform composite beam model is fully coupled and not limited by the shape of cross-section and boundary conditions, and it can capture complicated full 3D static deformations, thus it is a promising tool for engineering designs and applications.

Acknowledgements The authors would like to thank Science Foundation Ireland (SFI) for funding Spatially and Temporally VARIable COMPosite Structures (VARICOMP) Grant No. (15/RP/2773) under its Research Professor programme.

References References [1] Hetenyi M. Deflection of beams of varying cross section. ASME Trans J Appl Mech 1937;4:A49–52. [2] Jones E, et al. The flexure of a non-uniform beam. Pacific Journal of Mathematics 1955;5(5):799–806. [3] Jones E.

Finite fourier transform analysis of the flexure of a non-uniform beam.

The Aeronautical Journal

1956;60(552):805–6. [4] Gundersen RM, Iversen LE. Small deflections of non-uniform beams. Journal of the Franklin Institute 1962;273(5):365–74. [5] Boley BA. On the accuracy of the bernoulli-euler theory for beams of variable section. Journal of applied mechanics 1963;30(3):373–8. [6] Medwadowski SJ. Nonprismatic shear beams. Journal of Structural Engineering 1984;110(5):1067–82. [7] Rao HVG, Spyrakos C. Closed form series solutions of boundary value problems with variable properties. Computers & Structures 1986;23(2):211–5. [8] Lee SY, Ke HY, Kuo YH. Exact static deflection of a non-uniform bernoulli-euler beam with general elastic end restraints. Computers & structures 1990;36(1):91–7. [9] Lee S, Kuo Y. Static analysis of nonuniform timoshenko beams. Computers & structures 1993;46(5):813–20. [10] Eisenberger M. Exact solution for general variable cross-section members. Computers & Structures 1991;41(4):765–72.

32

[11] Friedman Z, Kosmatka J. Exact stiffness matrix of a nonuniform beami. extension, torsion, and bending of a bernoulli-euler beam. Computers & structures 1992;42(5):671–82. [12] Romano F. Deflections of timoshenko beam with varying cross-section. International journal of mechanical sciences 1996;38(8-9):1017–35. [13] Baker G. Exact deflections in nonprismatic members. Computers & structures 1996;61(3):515–28. [14] Wang CM, Chen C, Kitipornchai S. Shear deformable bending solutions for nonuniform beams and plates with elastic end restraints from classical solutions. Archive of Applied mechanics 1998;68(5):323–33. [15] Al-Gahtani H, Khan M. Exact analysis of nonprismatic beams. Journal of engineering mechanics 1998;124(11):1290–3. [16] Yavari A, Sarkani S, Reddy J. On nonuniform euler–bernoulli and timoshenko beams with jump discontinuities: application of distribution theory. International Journal of Solids and Structures 2001;38(46-47):8389–406. [17] Hodges D, Ho J, Yu W. The effect of taper on section constants for in-plane deformation of an isotropic strip. Journal of Mechanics of Materials and Structures 2008;3(3):425–40. [18] Attarnejad R, Shahba A, Semnani SJ. Analysis of non-prismatic timoshenko beams using basic displacement functions. Advances in Structural Engineering 2011;14(2):319–32. [19] Balduzzi G, Aminbaghai M, Sacco E, F¨ ussl J, Eberhardsteiner J, Auricchio F. Non-prismatic beams: a simple and effective timoshenko-like model. International Journal of Solids and Structures 2016;90:236–50. [20] Balduzzi G, Aminbaghai M, Auricchio F, F¨ ussl J. Planar timoshenko-like model for multilayer non-prismatic beams. International Journal of Mechanics and Materials in Design 2018;14(1):51–70. [21] Shahba A, Attarnejad R, Hajilar S. A mechanical-based solution for axially functionally graded tapered euler-bernoulli beams. Mechanics of Advanced Materials and Structures 2013;20(8):696–707. [22] Nguyen N, Kim N, Cho I, Phung Q, Lee J. Static analysis of transversely or axially functionally graded tapered beams. Materials Research Innovations 2014;18(sup2):S2–260. [23] Pydah A, Sabale A. Static analysis of bi-directional functionally graded curved beams. Composite structures 2017;160:867– 76. [24] Pydah A, Batra R. Shear deformation theory using logarithmic function for thick circular beams and analytical solution for bi-directional functionally graded circular beams. Composite Structures 2017;172:45–60. [25] Sachdeva C, Padhee SS. Analysis of bidirectionally graded cylindrical beams using variational asymptotic method. AIAA Journal 2019;:1–13. [26] Soltani M, Asgarian B. Finite element formulation for linear stability analysis of axially functionally graded nonprismatic timoshenko beam. International Journal of Structural Stability and Dynamics 2019;19(02):1950002. [27] Karamanlı A. Elastostatic analysis of two-directional functionally graded beams using various beam theories and symmetric smoothed particle hydrodynamics method. Composite Structures 2017;160:653–69. [28] Li X, Li L, Hu Y, Ding Z, Deng W. Bending, buckling and vibration of axially functionally graded beams based on nonlocal strain gradient theory. Composite Structures 2017;165:250–65. [29] Xie X, Zheng H, Zou X. An integrated spectral collocation approach for the static and free vibration analyses of axially functionally graded nonuniform beams. Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science 2017;231(13):2459–71. [30] G¨ unay MG, Timarci T. Static analysis of thin-walled laminated composite closed-section beams with variable stiffness. Composite Structures 2017;182:67–78. [31] Ai Q, Weaver PM. Simplified analytical model for tapered sandwich beams using variable stiffness materials. Journal of Sandwich Structures & Materials 2017;19(1):3–25. [32] Rajasekaran S, Bakhshi Khaniki H. Finite element static and dynamic analysis of axially functionally graded nonuniform small-scale beams based on nonlocal strain gradient theory. Mechanics of Advanced Materials and Structures 2018;:1–15.

33

[33] Zappino E, Viglietti A, Carrera E. Analysis of tapered composite structures using a refined beam theory. Composite Structures 2018;183:42–52. [34] Macquart T, Pirrera A, Weaver PM. Finite beam elements for variable stiffness structures. AIAA Journal 2018;56(8):3362– 8. [35] Masjedi PK, Doeva O, Weaver PM. Exact analytical solution for the static deflection of fully coupled composite beams. International Journal of Solids and Structures 2019;Under Consideration. [36] Doeva O, Masjedi PK, Weaver PM. Exact solution for the deflection of composite beams under non-uniformly distributed loads. In: AIAA Scitech Forum. 2020,. [37] Masjedi PK, Ovesy HR. Chebyshev collocation method for static intrinsic equations of geometrically exact beams. International Journal of Solids and Structures 2015;54:183–91. [38] Masjedi PK, Ovesy HR. Large deflection analysis of geometrically exact spatial beams under conservative and nonconservative loads using intrinsic equations. Acta Mechanica 2015;226(6):1689–706. [39] Masjedi PK, Maheri A. Chebyshev collocation method for the free vibration analysis of geometrically exact beams with fully intrinsic formulation. European Journal of Mechanics-A/Solids 2017;66:329–40. [40] Masjedi PK, Maheri A, Weaver PM. Large deflection of functionally graded porous beams based on a geometrically exact theory with a fully intrinsic formulation. Applied Mathematical Modelling 2019;76:938–57. [41] Amoozgar M, Shaw A, Zhang J, Friswell M. Composite blade twist modification by using a moving mass and stiffness tailoring. AIAA Journal 2018;:1–8. [42] Shirk MH, Hertz TJ, Weisshaar TA. Aeroelastic tailoring-theory, practice, and promise. Journal of Aircraft 1986;23(1):6– 18. [43] Kollar LP, Springer GS. Mechanics of composite structures. Cambridge university press; 2003. [44] Stodieck O, Cooper JE, Weaver PM. Interpretation of bending/torsion coupling for swept, nonhomogenous wings. Journal of Aircraft 2016;53(4):892–9. [45] Bernstein DS. Scalar, Vector, and Matrix Mathematics: Theory, Facts, and Formulas-Revised and Expanded Edition. Princeton university press; 2018. [46] Mason JC, Handscomb DC. Chebyshev polynomials. Chapman and Hall/CRC; 2002. [47] Yu W, Hodges DH. Best strip-beam properties derivable from classical lamination theory. AIAA journal 2008;46(7):1719– 24. [48] Hodges DH. Nonlinear composite beam theory. American Institute of Aeronautics and Astronautics; 2006.

Appendix A. Sample derivation of constants of integration Appendix A.1. Cantilevered beam subject to tip loads Consider a cantilevered beam under the action of free end tip loads. The expressions for the vectors C i , i “ 1, 2, . . . , 6, are obtained by applying boundary conditions (2.15) at x “ 0 and boundary conditions (3.10) at x “ `. qu , qϕ , qz and qy are also set to be zero. The vectors C i , i “ 1, 2, . . . , 6, are obtained as:

34

»

C1 “ ´ –

Fx

fi

fl Mx » fi Fz C 2 “ ´ – fl Fy » fi » fi ´My F fl ` ` – z fl C3 “ – Mz Fy C 4 “ C 5 “ C 6 “ r0

T

0s

(A.1a)

(A.1b)

(A.1c) (A.1d)

Appendix A.2. Cantilevered beam subject to distributed loads Consider a cantilevered beam subject to uniformly distributed loads. Applying boundary conditions (2.15) at the clamped end and boundary conditions (3.10) at the free end and setting tip loads Fx , Fy , F z, Mx , My and Mz to be zero we obtain: »

C 1 “ ´` –

qx

fi

fl qϕ » fi qz C 2 “ ´` – fl qy » fi `2 –qz fl C3 “ 2 qy

(A.2a)

(A.2b)

(A.2c)

C 4 “ C 5 “ C 6 “ r0

T

0s

(A.2d)

Appendix A.3. Simply supported beam subject to distributed loads For a composite beam simply supported at both ends boundary conditions (2.15a), (2.15b) and (2.16b) should be applied at both x “ 0 and x “ `. The following expressions for the unknown vectors C i , i “ 1, 2, . . . , 6, are derived:

» fi˛ q 0 ˝W 1 Φi` – fl ` W 2 Φi` – z fl‚ C 1 “ ´ WA Φ` B A qy qϕ ¨ » fi » fi˛ qx qz 1 C 4 “ ˝P – fl ` R – fl‚ ` qϕ qy `

˘ i ´1

¨

»

C 2 “ C 3 “ C 5 “ C 6 “ r0

35

qx

T

0s

fi

(A.3a)

(A.3b) (A.3c)

where, ` 0 i ˘´1 1 i 0 T ii 1 ii Φ` WA Φ` Φ` P “ WB WA Φ` ´ WB ` 0 i ˘´1 2 i 0 T ii 2 ii R “ WB Φ` WA Φ` Φ` WB Φ` ´ WD

and

ii Φii ` “ Φ p`q “ r

`3 , 6

`2 , 2

`3 , 3

`2 , 2

Φi` “ Φi p`q “ r`,

`4 , 12

`N `1 T s N `1

¨¨¨ ,

¨¨¨ ,

`N `2 sT pN ` 1q pN ` 2q

Appendix B. Stiffness properties of a composite strip [47]

˜

2

B EA “ b A11 ´ 12 D22

SET

ˆ

¸

,

˜

2

D GJ “ 4b D66 ´ 26 D22

B 12 D26 “ ´2b B 16 ´ D22

˙

,

SEF

¸

,

˜

2

D EIy “ b D11 ´ 12 D22

ˆ

B 12 D12 “ b B 11 ´ D22

˙

,

SF T

¸

,

ˆ

A216 A22 ´ 2A12 A16 A26 ` A212 A66 A226 ´ A22 A66

B 11 “ B11 `

A12 A66 B12 ` A16 A22 B16 ´ A26 pA16 B12 ` A12 B16 q A226 ´ A22 A66

B 12 “ B12 `

A12 A66 B22 ` A16 A22 B26 ´ A26 pA16 B22 ` A12 B26 q A226 ´ A22 A66

B 16 “ B16 `

A12 A66 B26 ` A16 A22 B66 ´ A26 pA16 B26 ` A12 B66 q A226 ´ A22 A66

D11 “ D11 `

2 2 A66 B12 ´ 2A26 B12 B16 ` A22 B16 2 A26 ´ A22 A66

36

b2 EA 12

D12 D26 “ ´2b D16 ´ D22

SEL “ SLT “ SF L “ 0 A11 “ A11 `

EIz “

˙

D12 “ D12 `

A66 B12 B22 ` A22 B16 B26 ´ A26 pB16 B22 ` B12 B26 q A226 ´ A22 A66

D22 “ D22 ` D16 “ D16 ` D26

2 2 A66 B22 ´ 2A26 B22 B26 ` A22 B26 2 A26 ´ A22 A66

A66 B12 B26 ` A22 B16 B66 ´ A26 pB16 B26 ` B12 B66 q A226 ´ A22 A66

` 2 ˘ A66 B22 B26 ` A22 B26 B66 ´ A26 B26 ` B22 B66 “ D26 ` A226 ´ A22 A66 D66 “ D66 `

2 2 A66 B26 ´ 2A26 B26 B66 ` A22 B66 2 A26 ´ A22 A66

where b is the width of the strip and Aij , Bij and Dij , i, j “ 1, 2, 6 are the entries of the well-know A, B and D matrices of classical lamination theory.

37