Buckling of Rectangular Composite Pipes under Torsion

Featured Application: Thin-walled rectangular composite pipes, because of their light weight, high stiffness, and strength, can be applied to ultra-lightweight airplanes, satellites, spacecraft, automobiles, and civil structures. Abstract: The numerical buckling load of rectangular composite pipes under torsional load was derived by using the energy method. The authors found no available simple design method or chart for the buckling loads of rectangular composite pipes, which are often used airplanes, spacecraft, and other lightweight structures, through their involvement in a Mars exploration airplane project. Thus, numerical results were obtained for length-to-width ratios ( l / b ) from 1 to 20, width-to-height ratios ( h / b ) from 1 to 6, and [0/90] layer ratios ( r ) from 0 to 1, which means [(0/90) r ,( ± 45) 1-r ] s . The layups were assumed to be symmetric, and tension-bending, torsion-bending, and tension-shear coupling stiffnesses were ignored. To establish a simple design method, a closed-form polynomial equation for the buckling load factor was derived by minimizing the weighted residuals of the safe and non-safe side errors, which were obtained by comparing the derived numerical results with the polynomial equations. As a result, the errors of the polynomial equation for the buckling load factor were 4.95% for the non-safe side and 12.4% for the safe side. The errors are sufﬁciently good for preliminary design use and for parametric design studies and optimization.


Introduction
The Mars exploration airplane project [1,2] is a joint project of JAXA (Japan Aerospace Exploration Agency) and Japanese universities. It is intended to be the first extraterrestrial atmospheric flight exploration. Aerial exploration can avoid complex terrain where a rover cannot explore and can also explore where a satellite cannot observe, like the sides of a valley. An image of the Mars exploration airplane in flight is shown in Figure 1. To build the airplane, the authors designed a wing to be tested at high altitude, where the density of the air is close to the density of the Martian atmosphere. A rectangular pipe, which can easily align the attack angle of the wing, was used. In addition, because the density of Martian atmosphere is 1/100th that of the Earth, the airplane must be ultra-lightweight; thus, the rectangular pipe was designed using thin carbon fiber reinforced plastics (CFRP). Thin structures usually fail not for lack of material strength or delamination but through buckling load because the buckling load of thin structure lower than material failure load. The buckling load can be increased by careful design of the layup sequence of the composite layers; however, no available simple design equation or chart were found in spite of there being many studies on composite structures [3][4][5][6][7][8][9][10]. The buckling load under shear force has been studied; however, those studies were not for composite rectangular tubes, but rather for long anisotropic plates [11,12]. An approximate expression for the torsional buckling load was derived by Omidvari and Hematiyan [13]; however, the solution is for an isotropic tube. Banks and Rhodes investigated postbuckling behavior of composite box sections [14]. Loughlan analyzed buckling of composite stiffened box sections under compression and bending by finite strip method [15]. Vo and Lee conducted numerical analysis for flexural-torsional buckling of thin-walled composite box beams by finite element method [16]. Moreover, while a comprehensive design handbook for the buckling of composite structures [17] has been published, it does not contain a method for the torsional buckling load of composite rectangular tubes. Given this state of affairs, it was decided that a numerical analysis of the torsional buckling load of cantilever rectangular orthotropic pipes should be conducted by using the energy method including a consideration of layup anisotropy. In addition, to easily use this analysis for design, a simple closed-form polynomial equation for the buckling load factor was derived by minimizing the weighted residuals of the safe and non-safe side errors of the derived numerical results. Although a numerical analysis like a finite element analysis (FEA) can determine the buckling load of anisotropic structures, a simple closed-form design equation is especially useful for the initial design and optimum design of the whole system that require a huge number of iterations for analysis of not only the structure but also its aerodynamics, guidance, and control. however, the solution is for an isotropic tube. Banks and Rhodes investigated postbuckling behavior of composite box sections [14]. Loughlan analyzed buckling of composite stiffened box sections under compression and bending by finite strip method [15]. Vo and Lee conducted numerical analysis for flexural-torsional buckling of thin-walled composite box beams by finite element method [16]. Moreover, while a comprehensive design handbook for the buckling of composite structures [17] has been published, it does not contain a method for the torsional buckling load of composite rectangular tubes. Given this state of affairs, it was decided that a numerical analysis of the torsional buckling load of cantilever rectangular orthotropic pipes should be conducted by using the energy method including a consideration of layup anisotropy. In addition, to easily use this analysis for design, a simple closed-form polynomial equation for the buckling load factor was derived by minimizing the weighted residuals of the safe and non-safe side errors of the derived numerical results. Although a numerical analysis like a finite element analysis (FEA) can determine the buckling load of anisotropic structures, a simple closed-form design equation is especially useful for the initial design and optimum design of the whole system that require a huge number of iterations for analysis of not only the structure but also its aerodynamics, guidance, and control.

Energy Method
To use a rectangular thin-walled pipe as a wing spar, a buckling analysis using the energy method was conducted on the cantilevered rectangular pipe under torque T shown in Figure 2. The rectangular pipe consisted of two vertical plates and two horizontal plates. The vertical and horizontal plates were assumed to be simply supported. The outplane displacements v and w are described as follows: Note that the study discusses about local wall buckling of rectangular pipe, not global buckling of long column. Thus, warping is not included. As in reference [18], the plates' edge angles were assumed to be constant, as shown in Figure 3.

Energy Method
To use a rectangular thin-walled pipe as a wing spar, a buckling analysis using the energy method was conducted on the cantilevered rectangular pipe under torque T shown in Figure 2. The rectangular pipe consisted of two vertical plates and two horizontal plates. The vertical and horizontal plates were assumed to be simply supported. The out-plane displacements v and w are described as follows: Note that the study discusses about local wall buckling of rectangular pipe, not global buckling of long column. Thus, warping is not included. As in reference [18], the plates' edge angles were assumed to be constant, as shown in Figure 3.  In this case, the out-plane displacements v and w are described using mn as follows: where n = 1, 2, …; m = 1, 2, …; and mn mn mn mn is a common coefficient to determine amplitudes of w and v and Equation (5) is the condition of constant the plates' edge angles. This assumption will be validated later. The total potential energy  can be expressed in terms of the strain energy U and work of the external force W as The contributions of a single horizontal plate (in x-y plane) to U and W are written as  In this case, the out-plane displacements v and w are described using mn as follows: where n = 1, 2, …; m = 1, 2, …; and mn mn mn mn is a common coefficient to determine amplitudes of w and v and Equation (5) is the condition of constant the plates' edge angles. This assumption will be validated later. The total potential energy  can be expressed in terms of the strain energy U and work of the external force W as The contributions of a single horizontal plate (in x-y plane) to U and W are written as In this case, the out-plane displacements v and w are described using Φ mn as follows: where n = 1, 2, . . . ; m = 1, 2, . . . ; and Φ mn is a common coefficient to determine amplitudes of w and v and Equation (5) is the condition of constant the plates' edge angles. This assumption will be validated later. The total potential energy Π can be expressed in terms of the strain energy U and work of the external force W as The contributions of a single horizontal plate (in x-y plane) to U and W are written as For vertical plates, y replaced z in Equations (7) and (8) and related equations appearing hereafter. For horizontal plates (in x-y plane), by considering the bending strain due to the out-plane displacement w, the strains can be described as follows: Considering Hooke's Law for anisotropic layers and substituting Equation (10) into Equation (7), using Equation (9), and integrating Equation (7), U p becomes For a symmetric layup, D xx , D xy , D yy , and D ss are Generally, D xs and D ss are non-zero even in a symmetric layup; when the number of layups is large enough, the following assumption can be used.
For vertical plates, "y" is changed to "z" in Equations (7)- (13). The total strain energy U and the work W of the four plates is written as Substituting Equations (3) and (4) into Equation (14) yields the total strain energy U and the work W of the external torsional stress τ cr : where τ cr = T/(2bht), i = 1, 2, . . . ; and j = 1, 2, . . . . The buckling torsional stress is obtained by minimizing the total potential energy.
where summation indicates only those values of i and j for which m + i and n + j are odd. Equation (18) can be rewritten as following matrix form for m + n is even is and for m + n is odd is To obtain a non-zero solution of Φ mn , the determinant of the coefficient matrix for Φ mn should be zero. The coefficient matrix separates into one for which m + n is odd and one for which m + n is even. It is well known that the thin isotropic plate under shear load, in which m + n is even has a lower buckling load than the one in which m + n is odd except for long plate. Numerical calculations, however, both the case that m + n is even and m + n is odd, were conducted because orthotropic stiffness may lead to lower buckling load for the case that m + n is odd.
The normalized torsional buckling load factor k s is defined as Appl. Sci. 2021, 11, 1342 6 of 12

Case of a Constant Edge Angle
To solve Equations (19) and (20), the terms of m and n should be limited to a finite number, m max and n max . To check the accuracy relative to the number of terms, a numerical calculation was conducted on a quasi-isotropic layup of T300 CFRP. The material properties of the quasi-isotropic layup are listed in Table 1. m max and n max = 10, 20, . . . , 100, and b = h = 100 mm. The results are shown in Figure 4.  To check the assumption of a constant edge angle, an FE analysis was conducted using NASTRAN (version2018.1). The dimensions of the pipe were b = h = 100 mm, t = 1 mm, and l = 400 mm. A square mesh with a size of 10 was used. A rigid body element (RBE2) connected the end of the plates to the center of the section and torsional load (torque) was applied to the center of the section. Figure 5 shows the buckling mode and the changes in the angles of the pipe from 90°.

Case of a Variable Edge Angle
Supposing that the edge angle is variable, substituting Equations (1) and (2) into Equation (14) yields the total strain energy U and the work W of the external torsional stress cr as follows: The error of m max = n max = 20 is within 1% of m max = n max = 100 for l/b ≤ 6, and the buckling stress τ cr is constant for l/b > 6 because of effect of boundary constraints at x = 0 and x = l vanishing; thus, m max = n max = 20 for l/b ≤ 6 and buckling stress τ cr (l/b = 6) for l/b > 6 were used in the subsequent numerical calculations.
To check the assumption of a constant edge angle, an FE analysis was conducted using NASTRAN (version2018.1). The dimensions of the pipe were b = h = 100 mm, t = 1 mm, and l = 400 mm. A square mesh with a size of 10 was used. A rigid body element (RBE2) connected the end of the plates to the center of the section and torsional load (torque) was applied to the center of the section. Figure 5 shows the buckling mode and the changes in the angles of the pipe from 90 • .  To check the assumption of a constant edge angle, an FE analysis was conducted using NASTRAN (version2018.1). The dimensions of the pipe were b = h = 100 mm, t = 1 mm, and l = 400 mm. A square mesh with a size of 10 was used. A rigid body element (RBE2) connected the end of the plates to the center of the section and torsional load (torque) was applied to the center of the section. Figure 5 shows the buckling mode and the changes in the angles of the pipe from 90°.

Case of a Variable Edge Angle
Supposing that the edge angle is variable, substituting Equations (1) and (2) into Equation (14) yields the total strain energy U and the work W of the external torsional stress cr as follows:

Case of a Variable Edge Angle
Supposing that the edge angle is variable, substituting Equations (1) and (2) into Equation (14) yields the total strain energy U and the work W of the external torsional stress τ cr as follows: The buckling torsional stress is obtained by minimizing the total potential energy.
The following equation is obtained from Equation (17).
The normalized torsional buckling load factor k s was obtained by setting the determinant of the coefficient matrix for B mn and H mn to zero. To check the accuracy relative to the number of terms of n and m, a similar numerical calculation as shown in Figure 4 was conducted, and the same conclusion was obtained for the terms of n and m.
The numerical results for Equations (18) and (26), and FEM for h/b = 0.7 are shown in Figure 6. The results for Equations (18) and (26) for h/b = 1.0 are also plotted for comparison. In the case of l/b = 1, the FEM result of h/b = 0.7 is close to the constant-angle case; however, in the region l/b ≥ 2, the FEM result of h/b = 0.7 is close to the variable-angle case. There are no apparent differences between the results of Equation (18) for h/b = 1.0 and those of Equation (26) for h/b = 1.0 because of the section is rectangular and thus the edge angle is constant. These results show that the assumption of a constant edge angle is a stiff constraint for the region h/b < 1.0. Thus, the assumption of a variable edge angle, Equation (26), was employed in the subsequent analysis because of its more accurate and safe-side estimation. Note that both FE analysis and the method of minimizing the total potential energy are essentially the same because they are based on energy principle. If the difference appears, the reason is convergence due to mesh size of FE model or number of terms of displacement function of Equations (1) and (2). The convergence is discussed in Section 2.2 and Figures 4 and 6 shows the present results gives conservative result comparing with FE results, thus the accuracy of the solution is reliable. Appl. Sci. 2021, 11, x FOR PEER REVIEW 9 of 14

Closed-Form Polynomial Equation
To make quick and convenient estimates for torsional buckling of a composite tube, a simple closed form polynomial equation was derived by fitting the numerical solution of the orthotropic [(0/90)r,(±45)1-r]s layups. Here, r is the ratio of the 0/90 layer number to the total layer number. Because the [(0/90)r,(±45)1-r]s layup is often used, it is thus chosen. The numerical analysis described above is hard to use for design purposes and would entail a huge calculation cost if it were to be used for optimum design of the whole system. Thus, it is considered that a simple approximate design equation would still be useful. Table 2 shows the material properties of unidirectional pre-preg tape, and Table 3 shows the layup sequence used in the numerical analysis. The total ply number was 32, and the binding stiffness was systematically controlled by varying r. Assuming that the number of layups was large enough, the bending stiffness was calculated as Eij t 3 /12 = Dij (i, j = x, y, s). Advanced stiffness properties like through thickness

Closed-Form Polynomial Equation
To make quick and convenient estimates for torsional buckling of a composite tube, a simple closed form polynomial equation was derived by fitting the numerical solution of the orthotropic [(0/90) r ,(±45) 1-r ] s layups. Here, r is the ratio of the 0/90 layer number to the total layer number. Because the [(0/90) r ,(±45) 1-r ] s layup is often used, it is thus chosen. The numerical analysis described above is hard to use for design purposes and would entail a huge calculation cost if it were to be used for optimum design of the whole system. Thus, it is considered that a simple approximate design equation would still be useful. Table 2 shows the material properties of unidirectional pre-preg tape, and Table 3 shows the layup sequence used in the numerical analysis. The total ply number was 32, and the binding stiffness was systematically controlled by varying r. Assuming that the number of layups was large enough, the bending stiffness was calculated as E ij t 3 /12 = D ij (i, j = x, y, s). Advanced stiffness properties like through thickness properties discussed in Refs. [19][20][21] are not used because of to make quick and convenient estimates. Figures 7-9 show the buckling load factor k s for l/b, h/b, and r. The ranges of l/b, h/b, and r are 1 ≤ l/b ≤ 6, 0.7 ≤ h/b ≤ 1, and 0 ≤ r ≤ 1. properties discussed in Refs. [19][20][21] are not used because of to make quick and convenient estimates. Figures 7-9 show the buckling load factor ks for l/b, h/b, and r. The ranges of l/b, h/b, and r are 1 ≤ l/b ≤ 6, 0.7 ≤ h/b ≤ 1, and 0 ≤ r ≤ 1.    properties discussed in Refs. [19][20][21] are not used because of to make quick and convenient estimates. Figures 7-9 show the buckling load factor ks for l/b, h/b, and r. The ranges of l/b, h/b, and r are 1 ≤ l/b ≤ 6, 0.7 ≤ h/b ≤ 1, and 0 ≤ r ≤ 1.    Figure 9 can be fitted by a first-order polynomial; thus, the following polynomial was employed.