An Improved C0 FE Model for the Sandwich Lattice Composite Panel

Combining the improved C0 plate element using high-order zigzag theories and the beam element degenerated from the plate element, a type of analysis model for the sandwich lattice composite panel was developed. Compared with the actual test results including the mid-span deflections and the surface sheet normal stresses, the outstanding of that method was presented through numeric calculation. The results showed that the model has great potential to become an excellent and highly efficient analysis and design tool for sandwich lattice composite panel to avoid the conventional three-dimension hybrid element model, which usually may lead to the complex program establishment, and the coupling degrees of freedom among the different types of elements.


Introduction
Composite materials and structures based on fiber-reinforced polymer (FRP) are widely used in many areas, such as aerospace, transportation, and building due to their excellent performance qualities of high strength, lightweight, anti-corrosion, easy design, and multi-forms.
One of the abovementioned is the laminated sandwich plate (LSP), which is composed of lamina face sheets stacked with each other and thick core materials, shows superior mechanical properties, and is applied in engineering practices [1][2][3][4]. Face sheets are mainly made of fiber materials such as glass fiber-reinforced polymer (GFRP), carbon fiber-reinforced polymer (CFRP), and core materials that include polyvinyl chloride (PVC) foam, polyurethane (PU) foam, balsa wood, and bamboo. Plenty of investigations have been conducted on the mechanical properties of LSP such as flexure, pressure, fatigue, delamination, etc. [5][6][7][8][9][10][11]. The research showed that the delamination usually occurred at the bonding interface between face sheets and core materials due to poor constraint of the foam core. Weiqing Liu, Hai Fang, et al. [12][13][14][15] presented a kind of sandwich-latticed composite panel (SLCP) manufactured by vacuum-infusing resin process, as shown in Figure 1. The lattice webs were added into the foam core to improve the peel resistance between face sheet and core surface. However, the crisscross reinforced webs destroyed the continuity of core and led to a complicated analysis model which cannot be modeled as LSP simply. Therefore, this paper presented a type of analysis model to calculate SLCP accurately and efficiently.
LSP can be calculated by several methods, including: (i) Equivalent modulus analysis [16] based on the classical laminate theory, which is hard for the complicated SLCP due to existing webs, which neglect many details of cores and webs and omit the transverse shear stress. (ii) Three-dimension finite element method (FEM), which is usually considered to be omnipotent for complex model construction. However, several issues may block its convenient application, including huge programming workload, calculation efficiency, Kirchhoff's classic laminate plate theory (CLT) [17] was the first comprehensive theory developed for plates but it did not take the transverse shear strains into consideration. The single-layer theory is also known as first-order shear deformation theory (FSDT) and the disadvantage is lack of a shear correction factor to predict the actual parabolic variation of shear stress and shear locking phenomenon for thin plates, which is still being investigated now [18,19]. To overcome the shortcomings of FSDT, higher order shear deformation theory (HSDT) [20,21] was proposed by Reddy, which is accurate and accounts for the transverse shear deformation and the transverse shear traction-free conditions on the top and the bottom surfaces of the plate without shear correction factor. However, the C 1 continuity issue along the interelement boundary promotes the theory development to some degree. The layer-wise theories were presented to surmount that issue. Discrete layer theories were concerned and proposed by A.Toledano, H.Murakami [22], and Reddy [23], which take unknown displacement components through all the layer interfaces. This plate theory possesses good performance, but the numbers of unknowns increase rapidly with the increase of layers that, in fact, lead to huge unaccepted calculations. So, zigzag theories (ZZT, known as refined plate theories) were developed by H. Murakami [24], S.P. Lee, et al. [25] for solving the aforementioned question though the unknowns at different interfaces linking to those at the reference plan unknowns. Improved versions about these theories have been suggested continuously. However, the zigzag theory still faced the C 1 continuity issue of the transverse displacement at the nodes by FEM. Combining the benefits of the discrete layer wise and higher order zigzag theories (HOZT), sub-laminate models and penalty stiffness multiplier were propounded to overcome C 1 continuity successively for LSP [26][27][28][29][30].
In fact, the improved C 0 FE model could receive accurate results efficiently [31][32][33][34]. Chalak et al. [35][36][37][38][39] proposed an improved C 0 finite element model for LSP analysis with a soft, compressible core using HOZT. In this model, the in-plane displacement fields as a combination of a linear zigzag function with different slopes at each layer and cubically varying function over the entire thickness are assumed. The out-of-transverse displacement within the core and the surface of face sheets is presumed to be quadratic and constant. This model satisfies the transverse shear stress continuity conditions at the layer interfaces and the conditions of zero transverse shear stress at the top and bottom of the plate, and has made significant contributions: (i) overcoming the C 1 continuity problem associated with HOZT to implement a C 0 formulation, (ii) possessing great effects on the core compressibility in the formulation, and (iii) eliminating the requirement of using a penalty multiplier in the stiffness matrix formation. The model can also be applied for the ordinary and medium-thick LSP.
In this paper, the improved C 0 plate finite element with soft compressible core made use of high-order zigzag theories and the beam element, which is the degeneration form of the above plate element (IC 0 FEM-HOZT). This finite element was combined to calculate the sandwich lattice composite panel under the out-of-transverse quasi-static loading. Therefore, a type of high efficient analysis method was developed rather than the ordinary three-dimension FEM composed of conventional shell or plate, beam, and solid elements. Nine-nodes-isoparametric-quadratic plate element was used for the face sheets and core. Three-nodes-isoparametric-quadratic beam element, regarded as a type of plate model degeneration, was adopted for the lattice webs cutting apart the core. The spatial occupied core areas by the web volumes can be subtracted by the same virtualized volume web with the core material, or even omitted since the web is very thin and the elastic modulus is far higher than that of the core material. The outstanding performance of that method was presented through numeric calculation compared with the actual test results, including the mid-span deflections and the surface feet normal stresses.

Model Displacements
The in-plane displacement fields [28] for the plate, shown in Figure 2, were chosen as follows: where u 0 and v 0 are on behalf of the in-plane displacements of any point at the midsurface, θ x and θ y are the rotations of normal to the middle plane about the yand x-axis respectively; n u and n l are the number of upper and lower layers respectively; β x , β y , η x , η y are the higher unknowns; α i xu , α i yu , α j xl , α j yl are the slopes of i th and j th layer corresponding to upper and lower layers respectively, and H z − z u i and H −z + z l j are the unit step functions. The out-of-transverse displacement was assumed to vary quadratically through the core thickness and be constant over the face sheets, which is expressed as Equation (3) and shown in Figure 3: The out-of-transverse displacement was assumed to vary quadratically through the core thickness and be constant over the face sheets, which is expressed as Equation (3) and shown in Figure 3: where w u , w 0 , and w l are the values of the out-of-transverse displacement at the top layer, middle layer, and bottom layer of the core respectively, and l 1 , l 2 , l 3 , are Lagrangian interpolation functions in the thickness co-ordinate as defined in [35,36]. The out-of-transverse displacement was assumed to vary quadratically through the core thickness and be constant over the face sheets, which is expressed as Equation (3) where u w , 0 w , and l w are the values of the out-of-transverse displacement at the top layer, middle layer, and bottom layer of the core respectively, and l1, l2, l3, are Lagrangian interpolation functions in the thickness co-ordinate as defined in [35,36].

Model of Constitutive Relations
The constitutive relationship of any kth orthotropic layer for plate or beam having any fiber orientation with the respect to structural axes system (x-y-z or x-z) is depicted as

Model of Constitutive Relations
The constitutive relationship of any kth orthotropic layer for plate or beam having any fiber orientation with the respect to structural axes system (x-y-z or x-z) is depicted as where {σ} and {ε} are the stress vector and the strain vector.
[Q] k is the transformed rigidity matrix of kth layer that can be evaluated with the material properties (E, elastic modulus; v, poisson ratio; G, shear modulus) and fiber orientation of kth layer on mechanics of composite structure. The detailed Equation (4) for plate is: The detailed Equation (4) for beam is: On the conditions of zero transverse shear stress at the top and bottom surfaces of the plate or beam, and the transverse shear stress continuity at the interfaces between the layers with the condition, u = u u and v = v u at the top, and u = u l and v = v l at the bottom of the plate, β x , η x , β y , η y , α i xu , α i xl , α i yu , α i yl , ∂w u /∂x, ∂w l /∂x, ∂w u /∂y, and ∂w l /∂y can be expressed by the displacements u 0 , v 0 , θ x , θ y , u u , u l , v u , and v l as where for plate as

And elements of [A]
are dependent on material properties. The last derivatives entries of the vector {B} help overcome the problem of C 1 continuity as mentioned before, including last four items of the plate and last two items of the beam.
According to the above equations, the in-plane displacement fields for plate as given in Equations (1) and (2) may be expressed as For plate For beam where the coefficients b i and c i are functions of thickness coordinates, unit step functions, and material properties as defined in [35,36]. Equations (6)-(8) do not contain any firstorder derivative terms of out-of-transverse displacements and avoid the requirements of C 1 continuity efficiently without new field variables [28] and penalty method [40,41]. The generalized displacement vector {δ} for the plate and beam model can be presented as For plate {δ} plate = u 0 v 0 w 0 θ x θ y u u v u w u u l v l w l For beam {δ} beam = {u 0 w 0 θ x u u w u u l w l } With the linear constitutive relation and Equations (1) where Equation (9) where Equation (10) And the elements of [H] plate and [H] beam are functions of z and unit-step functions, as given in [35,36] Polymers 2021, 13, 4200 6 of 17 The potential energy of the system can be expressed as where U s is the strain energy and W est is the work due to the elemental out-of-transverse static load. Equations (4), (9), and (10), U s can be presented by where for plate For beam The work due to the elemental transverse static load P can be calculated by To solve this problem, a nine-node quadratic element with 11 field variables (u 0 , v 0 , w 0 , θ x , θ y , u u , v u , w u , u l , v l , w l ) per node was employed for the plate. A three-node quadratic element with seven field variables (u 0 , w 0 , θ x , u u , w u , u l , w l ) per node was employed for the beam, which coordinates the plate element conveniently. The generalized displacement vector at any point for any plate or beam can be expressed as where {δ} plate = u 0 , v 0 , w 0 , θ x , θ y , u u , v u , w u , u l , v l , w l T for plate, and {δ} beam = {u 0 , w 0 , θ x , u u , w u , u l , w l } T for beam. {δ} i is the displacement vector corresponding to node i of plate or beam element; N i is the shape function associated with node i and n is the number of nodes per element, that is, nine for the plate or three for the beam.
With the help of Equation (16) where [B] plate or [B] beam are the strain-displacement matrices in the Cartesian coordinate system. The elemental potential energy as given in Equation (11) can be rewritten as where [N w ] T is the shape function like matrix with non-zero terms associated only with the corresponding out-of-transverse nodal displacements. Though the beam element node lacks partial field variables, the existing ones are all included in the plate-element-node field variables and are easy to extend to the plate's.
In accordance with Principle of Minimum Potential Energy, minimizing ∏ e as given in Equation (19) with respect to {δ} plate , the equilibrium equation is where [K e ] is the element stiffness matrix, and {P e } is the nodal load vector. The global stiffness matrix was formed by taking the contribution of all the plate elements and beam elements. The formation of the global load vector for the whole SLCP was just formed in consideration of the plate elements that contain the all-beam nodes. Then, the global linear simultaneous equations were formed and solved for the SLCP incorporating appropriate boundary conditions. In order to improve the displacementscalculation efficiency by the FE model, the sparse-matrix technique was utilized to store the global stiffness matrix. The stresses were calculated with the constitutive relationship by using the condition of stress continuity as in Equation (5). Meanwhile, this model, combined with the improved C 0 Zigzag plate model and its degeneration beam model, naturally circumvented the DOF coupling, made the beam and plate deformation compatible, and simplified the programming process.

Comparison Cases
To verify the effectiveness of the aforementioned improved C 0 finite element method on the ground of HOZT (IC 0 FEM-HOZT), the corresponding test results were compared here. Two cases about SLCP flexural experiments in the literature are presented as the referenced examples in the following sections [14,15].

. Specimen and Experiment Introduction
Five two-way reinforced SLCPs, composed of GFRP face sheets, GFRP webs, and rigid polyurethane foam cores, by vacuum-assisted resin infusion process, were tested under the concentrated load, which was loaded by the hydraulic actuator. The face sheet was formed by 0 • /90 • GFRP clothes, and the web sheet consisted of −45 • /45 • GFRP clothes. That length and width of all SLCP specimens were 1000 mm, while the effective support spans were both 910 mm. The side length of the core and the lattice web spacing varied from 75 mm, 125 mm, and 175 mm. The geometric details of SLCP specimens are described in Table 1, and the material properties of SLCPs are listed in Table 2. The experimental scheme is shown in Figure 4, where LVDT for the out-of-transverse deflection was set up under the bottom of the SLCP, as shown in Figure 5.   rally circumvented the DOF coupling, made the beam and plate deformation compatible, and simplified the programming process.

Comparison Cases
To verify the effectiveness of the aforementioned improved C0 finite element method on the ground of HOZT (IC0FEM-HOZT), the corresponding test results were compared here. Two cases about SLCP flexural experiments in the literature are presented as the referenced examples in the following sections [14,15].

Specimen and Experiment Introduction
Five two-way reinforced SLCPs, composed of GFRP face sheets, GFRP webs, and rigid polyurethane foam cores, by vacuum-assisted resin infusion process, were tested under the concentrated load, which was loaded by the hydraulic actuator. The face sheet was formed by 0°/90°GFRP clothes, and the web sheet consisted of −45°/45°GFRP clothes. That length and width of all SLCP specimens were 1000 mm, while the effective support spans were both 910 mm. The side length of the core and the lattice web spacing varied from 75 mm, 125 mm, and 175 mm. The geometric details of SLCP specimens are described in Table 1, and the material properties of SLCPs are listed in Table 2. The experimental scheme is shown in Figure 4, where LVDT for the out-of-transverse deflection was set up under the bottom of the SLCP, as shown in Figure 5.

Results and Discussion
The out-of-transverse displacements and face-plate normal stresses of the flexural experiment of SLCPs are discussed here. Figure 6 shows the test load-deflection curves compared to the top, mid, and bottom surface out-of-plane vertical deflections by FEM. Figure 7 shows the load-stress curves on the mid bottom surface. Table 3 lists the deflections and surface normal stresses corresponding to the ultimate bearing capacity.
Polymers 2021, 13, x FOR PEER REVIEW 10 of 18 Figure 7 shows the load-stress curves on the mid bottom surface. Table 3 lists the deflections and surface normal stresses corresponding to the ultimate bearing capacity.   In accordance with Table 3, the bottom deflection errors between the test and FEM for two-way simply supported SLCPs were within 52-61%. The bottom surface norm stress errors were within 8-22%. However, Figure 6 presents that the test load-deflectio curves were evidently elastoplastic. During the elastic stage, the FEM load-deflectio curves including the top, mid surface, and bottom were very close to the tests'. Figure  presents that the test load-stress curves were weak nonlinear until the ultimate stage. Du ing the whole loading history, the FEM bottom stress curves became close to the tests'.   In accordance with Table 3, the bottom deflection errors between the test and FEM for two-way simply supported SLCPs were within 52-61%. The bottom surface normal stress errors were within 8-22%. However, Figure 6 presents that the test load-deflection curves were evidently elastoplastic. During the elastic stage, the FEM load-deflection curves including the top, mid surface, and bottom were very close to the tests'. Figure 7 presents that the test load-stress curves were weak nonlinear until the ultimate stage. During the whole loading history, the FEM bottom stress curves became close to the tests'.

Specimen and Experiment Introduction
Two unidirectional-web and six bidirectional-web reinforced SLCPs, manufactured by vacuum-assisted resin infusion process as in Case 1, were tested under four-point load action. The face sheet was formed by 0 • /90 • GFRP clothes, and the web sheet was composed of −45 • /45 • GFRP clothes. That length of all SLCP specimens were 1000 mm, while the effective support spans were both 800 mm, and the width was 300 mm. The side length of the core and the lattice web spacing were 75 mm. The geometric details of SLCP specimens are described in Table 4, and the material properties of SLCPs are listed in Table 5. The experimental scheme is shown in Figure 8, where three LVDTs for the out-of-transverse displacements were set up under the bottom of the SLCP and the top specimen surface at support; the longitudinal and in-plane transverse normal stress on the top and bottom surface center at the midspan.

Results and Discussion
To verify the results with the aforementioned IC 0 FEM-HOZT, the corresponding test results were compared. The out-of-transverse displacements and bilateral normal stresses of the top and bottom midspan face sheets by the flexural experiment for the single-way simply supported SLCPs are discussed here. Figure 9 shows the load-midspan out-of-transverse deflection curves. Figure 10 shows the load-bilateral stress curves on the midspan top and bottom surfaces. Table 6 lists the deflections corresponding to the ultimate bearing capacity. Table 6 lists the deflection corresponding to the ultimate bearing capacity of the test. Table 7 lists that the bilateral midspan normal stresses on the top and bottom surfaces under the ultimate bearing capacity.
Two unidirectional-web and six bidirectional-web reinforced SLCPs, manufactured by vacuum-assisted resin infusion process as in Case 1, were tested under four-point load action. The face sheet was formed by 0°/90°GFRP clothes, and the web sheet was composed of −45°/45°GFRP clothes. That length of all SLCP specimens were 1000 mm, while the effective support spans were both 800 mm, and the width was 300 mm. The side length of the core and the lattice web spacing were 75 mm. The geometric details of SLCP specimens are described in Table 4, and the material properties of SLCPs are listed in Table 5. The experimental scheme is shown in Figure 8, where three LVDTs for the out-of-transverse displacements were set up under the bottom of the SLCP and the top specimen surface at support; the longitudinal and in-plane transverse normal stress on the top and bottom surface center at the midspan.      In accordance with Table 6, the bottom deflection errors between the test and FEM for single-way simply supported SLCPs were within 1-31%. Based on Table 7, the longitudinal normal stress errors on the top surface were within 11-35%; the transverse errors were within 2-35% except for SXG-2-2-75, which cannot be estimated by the current error method; the longitudinal normal stress errors at the bottom surface were within 14-50%; the transverse errors were within 3-140%. Both the test load-midspan deflection curves and load-stress curves took on elastic, generally, according to Figures 9 and 10. During the whole loading stage, the FEM results including on the top, at the midsurface, and at the bottom were close to the tests'. However, the sudden change of test data at the ultimate loading end showed the nonlinear effect evidently, which can lead to large errors on the relative deflection and normal stress.

Conclusions
In this study, nine-nodes-isoparametric-quadratic plate element and three-nodes-isoparametric-quadratic beam element were combined to simulate the sandwich lattice composite panel with the improved C0 finite element method with the soft compressible core using high-order zigzag theories. The deflections and normal results of SLCPs under the out-of-plane quasi-static loading with FEM were compared to those by the test.
As a whole, the combination method by the improved C0 finite plate element and the beam element degenerated from the improved C0 finite plate element, based on high-order zigzag theories, is a suitable method for the analysis of the sandwich lattice composite In accordance with Table 6, the bottom deflection errors between the test and FEM for single-way simply supported SLCPs were within 1-31%. Based on Table 7, the longitudinal normal stress errors on the top surface were within 11-35%; the transverse errors were within 2-35% except for SXG-2-2-75, which cannot be estimated by the current error method; the longitudinal normal stress errors at the bottom surface were within 14-50%; the transverse errors were within 3-140%. Both the test load-midspan deflection curves and load-stress curves took on elastic, generally, according to Figures 9 and 10. During the whole loading stage, the FEM results including on the top, at the midsurface, and at the bottom were close to the tests'. However, the sudden change of test data at the ultimate loading end showed the nonlinear effect evidently, which can lead to large errors on the relative deflection and normal stress.

Conclusions
In this study, nine-nodes-isoparametric-quadratic plate element and three-nodesisoparametric-quadratic beam element were combined to simulate the sandwich lattice composite panel with the improved C 0 finite element method with the soft compressible core using high-order zigzag theories. The deflections and normal results of SLCPs under the out-of-plane quasi-static loading with FEM were compared to those by the test.
As a whole, the combination method by the improved C 0 finite plate element and the beam element degenerated from the improved C 0 finite plate element, based on high-order zigzag theories, is a suitable method for the analysis of the sandwich lattice composite panel. IC 0 FEM-HOZT can avoid the conventional three-dimension hybrid element model composed by cube element, shell element, and beam, which usually may lead to a complicated building program, and the coupling degrees of freedom among the different types of elements. Though some deviation still exists between the calculation results by IC 0 FEM-HOZT and those by test due to the multiple causes such as nonlinear, IC 0 FEM-HOZT has great potential to become an excellent and highly efficient analysis and design tool for the sandwich lattice composite panel, if appropriate modifications are adopted according to the actual work.