Estimating of the Static and Dynamic Behaviors of Orthogrid-Stiffened FRP Panel Using Reduced-Order Plate Model

The orthogrid-stiffened FRP panel (OSFP) is a generic structural element in weight-sensitive structure applications. Based on the variational asymptotic method, a 2D reduced-order plate model (2D-RPM) of OSFP was constructed through matching the strain energy of the original panel for static and dynamic analyses. The local field distributions were recovered using the recovery relationship and global response. The relative influences of select parameters on the effective performance of the OSFP were revealed by parametric studies. The comparative results showed that the effective performance of the OSFP predicted by the 2D-RPM were consistent with those predicted by the 3D finite element model, but the computational efficiency was greatly improved. The stiffener height had the greatest influence on the natural frequency of the panel. The layup configurations of laminates had significant influences on the equivalent stiffness and buckling load of the OSFP but had little effect on the vibration modes, which could be varied by adjusting the stiffening forms.


Introduction
Fiber-reinforced polymer (FRP) has the advantages of being lightweight and having high strength, corrosion resistance, and tailorability. It is a new type of material that exhibits excellent performances, and has been extensively used in civil engineering, shipbuilding, aerospace and so on. Under the same load capacity, the weight of an FRP bridge is only 30% of that of a steel bridge and 5% of that of a reinforced concrete bridge. Due to the excellent properties of FRP, it has gradually become a substitute for traditional building materials (e.g., steel and concrete). Compared with traditional steel and concrete bridges, FRP bridges are not only lightweight but are also more convenient to manufacture, transport, and install. Furthermore, they have long service lives of more than 100 years [1].
The orthogrid-stiffened FRP panel (OSFP) is characterized by a lattice of rigid, interconnected stiffeners, which is a generic structural element in weight-sensitive structure applications [2]. Its stiffness, buckling, and vibration characteristics are not only related to the stiffening forms but are also closely related to the structural and material parameters, which increases the difficulty of analysis. At present, most analysis methods for stiffened panels can be summarized as numerical methods, analytical methods, and a combination of the two methods. The calculation models for the effective analysis of stiffened panels can be divided into finite element models (FEMs) [3][4][5][6][7][8], discrete stiffener models (DSMs) [9], and smeared stiffener models (SSMs) [10]. The stiffeners are modeled as members with axial bending/torsion stiffnesses on the attached skin in DSM, and they can only be effectively applied for the analysis of stiffened panels with thin skins and rigid stiffeners. The geometric features of the stiffeners and skin are contained in FEM, which is more flexible and accurate than DSM. However, it not only requires a long simulation time due to the large computational workload, but also the inability to set different configurations and materials during preliminary design and optimization [11]. The basic idea of SSM is to smear the stiffnesses of the stiffeners into the panel and calculate the effective properties [12][13][14]. This may be applicable to the overall analysis of stiffened panels with intensive interconnected stiffeners but not to the stress and strain analysis of stiffeners [15].
The global buckling is considered as the main failure mode of stiffened panels under axial compression or/and external pressure according to the failure mode map [16,17]. Among the methods used to study the global buckling of stiffened panels, extensive studies have concentrated on the FEM and SSM [18,19]. The extension-bending coupling interaction caused by the eccentricity of the stiffeners was neglected in the work by Chen et al. [20], which resulted in the imprecise buckling prediction results [21]. Abhijit et al. [22] used the FEM to study the free vibration characteristics of stiffened plates with symmetric stiffeners. Later, the isoparametric stiffened-plate element was introduced to analyze the free vibration of eccentrically stiffened plates [23,24]. Lam et al. [25] explored the vibrations of stiffened plates by dividing the plate domain into appropriate rectangular segments.
On the basis of literature review, there were some issues in the existing methods for predicting the effective performance of stiffened panels. First, the shear stiffness matrix cannot be predicted using the homogenized elastic constants and plate thickness together with classic plate theory. Second, due to the assumptions for defining the kinematics, the models were usually applicable to some specific form of stiffeners. Third, few plate models can accurately predict the local stress and strain distributions of stiffeners, which were of great significance to the failure analysis of stiffened panels.
Recently, Yu and Zhong [26][27][28][29] put forward a new multiscale modeling technique to deal with dimensionality-reducing structures based on the variational asymptotic method (VAM) [30]. With this method, the small structural parameters (such as the thicknesswidth ratio) were used to asymptotically expand the energy functional, and the higherorder terms were removed to obtain the approximate energy of different orders and the corresponding dimension reduction model, which achieved a good tradeoff between accuracy and effectiveness. In this work, a VAM-based reduced-order plate model of stiffened FRP panel was established to solve the three issues mentioned above. The influences of geometric parameters (such as height, thickness, length-width ratio, and period length) on the effective performance of the OSFP were investigated. Finally, the displacements and free-vibration modes of stiffened FRP panels with different stiffening forms, such as the orthogrid-, T-and blade-stiffened panels, were compared. To the author's knowledge, this method has never been used for this purpose. Figure 1, if the sizes of the whole panel (denoted by the macro-coordinates x i ) are much larger than those of a unit cell (denoted by the micro-coordinates y i ), then y i = x i /ξ (ξ is a small parameter), and the derivative of the function f ξ (x α ) with respect to

As shown in
where i, j = 1, 2, 3; α, β = 1, 2. To construct a reduced-order plate model of the OSFP using VAM, the 3D displacement field of the original OSFP u i need to be represented by using 2D plate variables v i such as where w i is the fluctuating function to be solved, and the underlined terms should meet the following constraints where · represents the volume integration over the unit cell.
The non-underlined terms should satisfy the following conditions The 3D strain field can be expressed as Plugging Equation (2) into Equation (5) gives the 3D strain field as Γ 11 = ε 11 + ξy 3 κ 11 + w 1,1 2Γ 12 = 2ε 22 + 2ξy 3 κ 12 + w 1,2 + w 2,1 where ε αβ and κ αβ can be expressed as The 3D strain field can be obtained as where Γ e , Γ s , Γ t are strain matrices of 3D-FEM; As shown in Figure 2, the unit cell within the OSFP can be divided into three parts to facilitate the integral solution. Then we obtain the strain energy of the panel as where with the subscripts A, B, and C representing the skin, longitudinal stiffener, and transverse stiffener, respectively. Equation (10) can be rewritten as where D e , D es , D et , D s , D st , and D t are the corresponding sub-matrices of three-dimensional 6 × 6 material matrix. The virtual work done by the applied loads is where f i is the body force, α i is the traction force applied on the lateral surfaces, β i and τ i denote the traction forces on the bottom and top surface, respectively. Plugging Equation (2) into Equation (13) gives where Φ 1 = v 3,2 , Φ 2 = −v 3,1 , and The values of p i , q α , P i , and Q α in Equation (14) can be calculated as According to VAM, E * can be ignored, and the total potential energy is

VAM-Based Reduction Analysis of OSFP 2.2.1. Zeroth-Order Approximation
Plugging Equation (8) into Equation (17) gives the total potential energy density as where the underlined items and the double-underlined item can be ignored according to VAM.
To impose the constraints on the fluctuating function, we introduce the Lagrange multipliers λ i , such as The zeroth-order approximate variational expression is The corresponding Euler-Lagrange equations are where The boundary conditions of the top and bottom of the panel can be defined as where the superscript "+ / −" indicates the items at the top and bottom of the panel. From these conditions, we can solve w || and w 3 as where Plugging Equation (23) into Equation (18) gives the zeroth-order strain energy as where A, D and B are tensile, bending, and coupling stiffness sub-matrix, respectively, and can be expressed as

Transforming into Reissner-Mindlin Model
There are two additional transverse shear strains γ = 2γ 13 2γ 23 T in the Reissner-Mindlin model. To transform Equation (25) into the Reissner-Mindlin model, we must eliminate the coupled stiffness terms between and γ as follows: where R is Reissner-Mindlin generalized strains. The final form of the total energy can be expressed as where F is a load-related term and The resultant stress of the panel can be expressed as Due to the symmetry of the axis and plane, some stiffness components disappear, and the constitutive relation of the OSFP can be obtained as The original 3D geometric nonlinear problem in Equation (17) is mathematically decomposed into constitutive modeling over the unit cell in Equation (31) and geometric nonlinear plate analysis. That is to say, as an alternative to the direct numerical simulation using 3D nonlinear finite element analysis, the global analysis of the OSFP can be reduced to 2D plate analysis using the linear solver in ABAQUS, with the constitutive relation obtained from the constitutive modeling of the unit cell.

Equivalent Density of 2D-RPM
For heterogeneous materials and structures, the equivalent density of the unit cell is usually used to characterize the weight of the reduced-order plate model. The total mass of the unit cell within the OSFP can be expressed as where ρ 1 , ρ 2 , and ρ 3 are the density of the skin, the longitudinal stiffener and transverse stiffener, respectively. The total volume of the unit cell is According to the equivalent density formula, we obtain The equivalent density of the 2D-RPM can be expressed as

Validation Example
To verify the accuracy and effectiveness of the present reduced model, the static and dynamic behaviors of OSFP predicted by the present model were compared with those of 3D-FEM. The 3D-FEM had 15 unit cells in the x 1 and x 2 direction as shown in Figure 3. The equivalent stiffness of OSFP was obtained by variational asymptotic analyzing over the unit cell shown in Figure 4b and inputted into the 2D-RPM (300 mm × 300 mm) using shell elements, as shown in Figure 4c

Static Displacement Analysis
Six typical boundary conditions shown in Figure 6, including CCCC, CCSS and CSCS, CSSS, SSSS and FFCC, were used for static displacement analysis. The naming convention of boundary conditions is four letters, where S denotes simply supported constraint, C for fixed constraint, and F for free constraint.
To verify the effectiveness of the 2D-RPM, a uniform load of 5 kPa was applied to the top surface of the OSFP, and the displacement distributions along Path 1 of the 3D-FEM and 2D-RPM were compared in Figure 7. The comparative results show that the displacement distributions predicted by the 2D-RPM were basically in agreement with those of the 3D-FEM. The differences were due to the different meshing methods used in the two models. The maximum displacement error under the CCSS boundary condition was the largest, but it was still within 5%. It is worth noting that the differences between the 2D-RPM and 3D-FEM in Figure 7c-e were much greater than other cases, which may be due to the gradual enhancement of boundary constraints from SSSS to CCCC. It was concluded that the 2D-RPM can predict the static displacement of the stiffened FRP panel with high accuracy and effectiveness, and the equivalent stiffness obtained from the VAM had sufficient accuracy.

Local Displacement and Stress Field
According to Equation (2), the local displacement distributions within the unit cell at the geometric center of the 2D-RPM can be recovered as shown in Figure 8. It can be observed that the maximum and minimum values of U 1 and U 2 were located on opposite sides of the stiffener, and the overall displacement presented a centrosymmetric distribution trend. The smallest value of U 3 was located at the center of the stiffener, while the maximum value of U 3 was on the skin. The displacement distribution was centrally symmetric about the intersection point of the stiffeners. The maximum value of U 3 within the unit cell was 10.27 mm, and the error was about 1% compared with that of 3D-FEM presented in Section 3.1 (U 3 = 10.17 mm), indicating that the recovered displacement distribution is accurate and can be used to evaluate the location of maximum local displacement.  Figure 9 shows the local stress fields within the unit cell at the geometric center of the plate recovered from Equation (8) and the 3D Hooke's law. It can be observed that the stiffeners played an important role in the process of load transfer, and there was a large stress concentration at the intersection of the stiffeners and the skin, showing a significant skin-stiffener effect. The stress distribution on the skin was relatively uniform, and there was no evident mutation. It was concluded that the stiffeners improved the bearing capacity of the panel, and the OSFP had a low weight and high strength compared to the ordinary panel. Figure 10 shows the local von Mises stress and displacement distribution along Path 1 of the skin within the unit cell (as shown in Figures 8c and 9a) predicted by 2D-RPM and 3D-FEM. It can be observed that the local stress and displacement curves predicted by 2D-RPM and 3D-FEM agreed well, and the maximum error was less than 5%. The local stress at the junctions between the skin and the stiffeners decreased significantly, indicating that these regions were very incidental to be damaged.  Table 1 shows the the first four vibration modes and natural frequencies of the OSFP under the CCCC boundary condition predicted by 2D-RPM and 3D-FEM. The vibration modes of the 3D-FEM and 2D-RPM were in good agreement. For example, there were one and two half-waves along the x 1 direction, two half-waves along the x 2 direction, and two half-waves along the x 1 and x 2 directions for the first, second, third, and forth mode shapes, respectively, for both the 2D-RPM and 3D-FEM results. The natural frequencies of the 3D-FEM and 2D-RPM were also highly consistent, and the maximum error of the natural frequency was less than 6.83%. It was worth noting that the first-order vibration frequency showed relative big error compared with the third and fourth vibration frequencies, which may be because the first-order frequency was more sensitive to different meshing between 2D-RPM and 3D-FEM. The 2D-RPM was more time-efficient than the 3D-FEM in analyzing the vibration modes: the 2D-RPM required 30 s with one CPU as opposed to the nearly 20 min required for the 3D-FEM with four CPUs. It can be concluded that the 2D-RPM had high accuracy in free-vibration analysis of OSFP under CCCC boundary condition. To further verify the effectiveness of 2D-FEM in analyzing vibration modes, the vibration modes of OSFP under different boundary conditions are analyzed as shown in Table 2. The 3D-FEM and VAM-based 2D-RPM had good consistency in the prediction of natural frequencies and vibration modes under various boundary conditions, and the maximum error of natural frequency was less than 6% under the CCSS boundary condition. The stronger the boundary condition is, the higher the natural frequency is. The natural frequency under CCCC boundary condition was about twice that under SSSS boundary condition. The natural frequencies under CCSS and CSCS boundary conditions were almost the same, but the asymmetry of CSCS boundary condition led to the asymmetry of vibration mode.

Global Buckling Analysis
To verify the accuracy and effectiveness of 2D-RPM, the buckling modes and critical loads of OSFP under different boundary and load conditions illustrated in Figure 11 are listed in Table 3 Table 3 shows that the buckling load in Case 3 (471.83 N) was about 9 times that in Case 2 (57.96 N) and 1.577 times that in Case 4 (300.67 N). The buckling loads in Case 1 and Case 2 were basically the same. The error of the critical buckling load under various boundary conditions was less than 5%, indicating that the VAM-based 2D-RPM and 3D-FEM predictions of the global buckling of OSFP agreed closely.

Parameter Study
The 2D-RPM was selected to conduct parametric study. The material properties and layup configurations of the skin and the stiffeners were the same as those in Section 3, except where explicitly indicated. The boundary conditions for uniaxial buckling analysis were SFFF in Case 1, while the boundary conditions for the free-vibration analysis were CCCC. Figure 12a shows the effect of the stiffener thickness on the equivalent stiffness of the OSFP when the other parameters remained unchanged. The equivalent stiffness A ij and D ij increased with increasing stiffener thickness, and in particular, the bending stiffness components D 11 and D 22 increased significantly. This may be because A ij was directly proportional to the cross-sectional area, which increased with increasing stiffener thickness. In contrast, D ij was proportional to the moment of inertia, which was linearly related to the stiffener thickness. Figure 12b shows the effect of stiffener height on the equivalent stiffness when other parameters remain unchanged. It can be seen that A ij increased linearly and D ij increased nonlinearly with the increasing stiffener height. The reason was that A ij was proportional to the sectional area, which was linear with the stiffener height, while D ij was proportional to the moment of inertia, resulting in a parabolic growth trend. Figure 12c shows the effect of length-width ratio on the equivalent stiffness. It can be observed that A 11 and D 11 had no obvious change with the increasing length-width ratio, while D 22 decreased significantly. The reason was that the extension area and moment of inertia in x 1 direction remained unchanged, while the extension area and moment of inertia in the x 2 direction gradually decreased with the increasing length-width ratio, which led to the nonlinear decrease in equivalent bending stiffness. Figure 12d shows the effect of the periodic length on the equivalent stiffness of the stiffened FRP panel. It can be observed that A ij and D ij decreased nonlinearly with the increase in periodic length. This was because A ij and D ij were, respectively, proportional to the extension area and the moment inertia, and there was a negative nonlinear relationship between the extension area/moment inertia and the periodic length, resulting in a parabolic downward trend with the increasing periodic length.

Influence of Structural Parameters on Buckling Loads and Natural Frequencies
To further investigate the influence of structural parameters on the effective performance of OSFP, the first four buckling loads and natural frequencies of OSFP with different stiffener height, thickness, periodic length, and length-width ratio were calculated by using 2D-RPM, as shown in Figure 13.
The first four natural frequencies of the OSFP increased with increasing stiffener thickness and height and decreased with increasing length-width ratio and periodic length. The effect of the stiffener height h on the natural frequency was much greater than that of other structural parameters. The reason is that the variation trends of the equivalent stiffness and equivalent mass were consistent with those of structural parameters, and their influences on the natural frequency might counteract each other. However, the effect of the stiffener height h on the equivalent stiffness was much greater than that on the equivalent mass. The buckling load of the OSFP increased with the increase in stiffener thickness and height but decreased with increasing length-width ratio and periodic length, which was the same as the change trend of equivalent stiffness.

Influence of Layup Configuration on the Effective Performance of OSFP
The layup configurations of the laminates would affect the effective performance of the OSFP due to the anisotropy and heterogeneity. In this section, the influences of the layup configuration on the equivalent stiffness, free vibrations, and buckling mode of the OSFP are analyzed. The layup configuration was set to [0/θ/0/θ/0] s , where θ increased from 0 • to 90 • at 15 • intervals. The boundary condition was fixed on one side and simply supported on three sides (CSSS). Figure 14a shows the effect of the layup configuration on the equivalent stiffness of the OSFP. With the gradual increase in the ply angle, the stiffness components A 11 and D 11 showed nonlinear downward trends, while A 22 and D 22 showed significant nonlinear increases when the ply angle was greater than 45 • . Figure 14b shows the effect of the layup configuration on the first four natural frequencies and buckling loads of the OSFP. It can be observed that the layup configuration had little effect on the natural frequency, and the first natural frequency first decreased and then increased with the increasing ply angle, reaching the minimum value in the range of 30-60 • ply angle. The buckling load first increased and then decreased with the increasing ply angle, and the buckling load of each order reached the maximum value at 30-45 • ply angle.  Figure 16a shows that the influences of the 0 • -ply ratio on A 11 and D 11 were much greater than those on the other stiffness components because the stiffness along the fiber direction (0 • direction) was stronger. Figure 16b shows that the first natural frequency increased with increasing 0 • -ply ratio and reached a maximum value when the 0 • ply ratio was 100%. The third to fourth natural frequencies increased first and then decreased and reached the maximum value when the 0 • -ply ratio was between 0.2 and 0.6. With the increase in the 0 • -ply ratio, the first to fourth buckling loads of the OSFP increased gradually and reached a maximum value when the 0 • ply ratio was 100%. In engineering applications, the effective performance of the OSFP in the corresponding direction can be improved by adjusting the 0 • -ply ratio.

Comparison with Other Stiffened FRP Panels with Different Stiffening Forms
To compare the effects of different stiffening forms on the effective performance of the OSFP, the 3D FE models and 2D reduced-order plate models of orthogrid-, T-, and blade-stiffened FRP panel were established. The 3D FE models were obtained by repeating the unit cell 15 times in the x 1 and x 2 directions as shown in Figure 17 The static displacements along the center line of the stiffened FRP panels under the CCCC boundary condition and 5 kPa of uniform load were analyzed. The comparative results in Figure 18 show that the displacement of the blade-stiffened FRP panel was the largest, followed by the orthogrid-and T-stiffened FRP panels due to the fact that the equivalent bending stiffness of T-stiffened FRP panel was greater than the other two stiffened FRP panels.  Table 4 shows the first four natural frequencies of the stiffened FRP panel with different stiffening forms under the CCCC boundary condition. The first natural frequency of the T-stiffened FRP panel was the largest. From the second order, the natural frequency of the orthogrid-stiffened FRP panel was the largest, followed by the T-and blade-stiffened FRP panels. The comparative results showed that the natural frequencies of the OSFP increased faster with increasing modal order, while those of the T-and blade-stiffened FRP panels showed little change. The vibration modes of the T-and blade-stiffened FRP panels were basically the same, but the vibration modes of the orthogrid-stiffened FRP panel were very different (there were one and two half-waves along the x 1 direction, two half-waves along the x 2 direction, and two half-waves along the x 1 , and x 2 directions for the first, second, third, and forth mode shapes, respectively). It was concluded that the vibration modes of the stiffened FRP panel could be changed by adjusting the stiffening forms.

Conclusions
In this work, a VAM-based reduced-order plate model was established to analyze the effective performance of the orthogrid-stiffened FRP panel (OSFP). The influences of the material and structural parameters of the OSFP were investigated by parametric studies, and the following conclusions were drawn.
(1) The results of static displacement, local field distributions, natural frequencies and buckling loads predicted by 2D-RPM were consistent with those of 3D FE model, but the computational efficiency was greatly improved, which verifies the accuracy and effectiveness of the VAM-based reduced-order plate model.
(2) The equivalent stiffness increased gradually with increasing stiffener thickness or height and decreasing length-width ratio or periodic length. The influences of the structural parameters on the buckling load were closely related to the equivalent stiffness. The effect of the stiffener height on the natural frequency of the OSFP was much greater than those of the other structural parameters.
(3) Different layup configurations had significant influences on the equivalent stiffness and buckling load of the OSFP, while they had little effect on the vibration modes. With the increase in the 0 • ply ratio, the equivalent stiffness in the fiber direction increased significantly, and the first buckling load and natural frequency increased gradually. The static displacement and vibration modes of the orthogrid-stiffened FRP panel were different from those of the blade-and T-stiffened FRP panels, indicating that the vibration modes of the stiffened FRP panel could be varied by adjusting the stiffening forms. Data Availability Statement: Data available on request due to restrictions, e.g., privacy or ethical. The data presented in this study are available on request from the corresponding author. The data are not publicly available due to subsequent analyzes and publications.

Conflicts of Interest:
The authors declare no conflict of interest.