Flexural behavior of composite structural insulated panels with magnesium oxide board facings

The current report is devoted to the flexural analysis of a composite structural insulated panel (CSIP) with magnesium oxide board facings and expanded polystyrene (EPS) core, that was recently introduced to the building industry. An advanced nonlinear FE model was created in the ABAQUS environment, able to simulate the CSIP’s flexural behavior in great detail. An original custom code procedure was developed, which allowed to include material bimodularity to significantly improve the accuracy of computational results and failure mode predictions. Material model parameters describing the nonlinear range were identified in a joint analysis of laboratory tests and their numerical simulations performed on CSIP beams of three different lengths subjected to three- and four-point bending. The model was validated by confronting computational results with experimental results for natural scale panels; a good correlation between the two results proved that the proposed model could effectively support the CSIP design process.


Introduction
Sandwich structures are becoming increasingly popular in civil engineering applications, as their use allows for a reduction of dead weight, improvement of sustainability, and overall cost-efficiency [1][2][3]. The low weight makes even large modular elements easy to handle, significantly improving the speed of transport and assembly [2,4]. Some notable examples of sandwich structure applications in civil engineering can be found in the housing industry [2,4] and footbridges [5,6].
The general idea of a sandwich structure is to combine two thin high-strength facings and a light structural core, to create a panel of considerable strength and stiffness [7,8]. An essential advantage of the sandwich concept is its versatility. It can be tailored, through the use of different geometric proportions and various combinations of facing and core materials, to acquire properties required for a particular purpose [9,10].
Structural insulated panel (SIP) is a type of sandwich structure designed for use in single-and two-story buildings. SIPs are typically composed of a thick foam core and woodbased facings (e.g., oriented strand board), and distinguish themselves by the threefold role they play in a structure: (1) enveloping functionality, (2) thermal insulation, (3) transfer of structural loads. They are applied as structural parts of walls, floors, and roofs for a quick assembly of energyefficient, "green" buildings resistant to deterioration even in extreme climate [11][12][13]. However, the use of wood-based facings makes them prone to biological corrosion and fire. A developed version of SIP was proposed to overcome these disadvantages. Composite structural insulated panel (CSIP) employs composite material facings, which additionally results in higher strength and stiffness, making it useful in a broader range of civil structures [14][15][16].
The main subject of the current research is a CSIP with magnesium oxide board (MgO board) facings and an expanded polystyrene (EPS) core (Fig. 1). Employing the MgO board as a cladding material, makes the CSIP retain the functionality of a traditional SIP and eliminates its main disadvantages. The MgO board is an emerging new type of cladding that is not only fire-resistant and immune to mold, fungus, and insects, but it also has many "green features," making it an environmentally sustainable material. It also has high compressive and flexural strength, which is an important quality for a load-bearing composite [17,18]. Combined with exceptional thermal insulation properties of EPS [19], it creates an attractive new construction material that requires research attention.
The intrinsically heterogeneous nature of CSIP, and the complexity of interactions between its layers, makes it challenging to foresee its load-bearing capacity and failure mode. Reaching the ultimate stress in one of the facings is the most desirable type of failure; however, it might be preceded by global buckling, local facing buckling (wrinkling), delamination, core shear, or core crushing [7,8,20].
To meet market expectations, CSIP manufacturers modify their products regularly, and extensive cycles of full-scale laboratory tests are required after each adjustment. To reduce the costs of experimental procedures, it is crucial to discard in advance any changes that may lead to undesirable failure mechanisms. One of the most common changes is the increase in the core thickness. It is motivated by the need to improve thermal insulation, strength, and stiffness properties of CSIPs; however, it also makes them more prone to local forms of failure when subjected to transverse loads. For this reason, the flexural behavior of the panel under consideration is a fundamental matter to investigate.
Examples of different computational approaches to this class of problems can be found in the literature. Classical sandwich theory [7,8] and other equivalent single-layer models [21,22] were used when detailed core deformation analysis was not necessary [5,9,23,24]. A fully threedimensional approach was utilized to obtain more precise core strain state information, by discretizing all layers with solid elements [14,24,25], or using a combination of the shell (facings) and solid elements (core) to reduce the computational cost [16,[26][27][28][29][30][31]. Two-dimensional models were formulated within assumptions of plane stress conditions, allowing for enhanced meshing and more detailed results [32]. In delamination failure analysis, contact at the interface had to be included [27]; however, when this behavior was not observed before failure initiation, perfect bonding between layers was assumed instead [25][26][27][28][29][30][31][32]. Nonlinear behavior was included in the material models of facings [29], core [16,25], or both [26][27][28]30] depending on failure progression. It was shown, that the consideration of the postyield core behavior has a significant role in the modeling of soft foam panels [16,25] and that the complex phenomena occurring at the supports before failure depend strongly on the combination of normal and shear stress distribution [28]. The geometrical nonlinearity was included, when significant deformation affected the load distribution [14,[25][26][27][28][29].
Earlier studies allowed to determine fundamental material properties of the constituent materials, as well as to observe typical CSIP failure progression patterns, in an extensive cycle of laboratory tests [33,34]. It was noted that both EPS and MgO board exhibit distinctly different behavior in compression and tension in all ranges: elastic, plastic, and failure, indicating at the substantial influence of the material bimodularity. An initial FE solution incorporating material bimodularity effect was created and proved to capture CSIP flexural behavior quite well; however, it failed to predict the correct failure pattern in some cases [35]. It was determined that the correct prediction of failure initiation location strongly depends on (1) the accuracy of elastic core response in areas where stress-state transits from compression to shear (e.g. zones over the supports), and (2) how different definitions of hardening in adjacent sample areas influence stress redistribution in the post-yield range.
The main goal of the current work is to propose a reliable and robust computational framework, that can predict failure modes of the MgO board CSIPs under flexural loads, to aid their design process and reduce the number of expensive laboratory tests. To achieve this, a development of the FE model is required. Based on previous work results, it was decided that the most efficient areas of improvements are: (1) introduction of an additional stress-state corresponding to the core shear behavior in the elastic range (only compression and tension were distinguished before), (2) improvement of the post-yield range description accuracy for both core and facings in all stress-states. To accomplish this, a parameter identification study has been performed, based on a comparison between six small-scale experimental flexure tests on CSIP beams with varying lengths and load distributions, and their numerical representations. Afterward, the final version of the FE model has been validated by performing a comprehensive comparison of computational results with laboratory test data for a natural scale panel flexure.

Analyzed panel
The CSIP under consideration is composed of two MgO board facings enclosing a core made of expanded polystyrene (EPS). The layers are bonded together with a one-component polyurethane adhesive. A schematic layout is presented in Fig. 1. All laboratory test specimens came from CSIPs intended for use as a part of a wall assembly, and were provided by the producer. The expanded polystyrene, used as the CSIPs core, is a polymeric closed-cell foam used widely in the building industry. It has low weight, exceptional thermal insulation properties, and very low water absorption properties, making it a popular core material in many types of sandwich panels [19,36]. Its most problematic disadvantages are a low melting temperature and susceptibility to fire. Therefore, only EPS with flame-redundant properties is allowed for use in the construction industry, preferably with a fire protection material cover [19], such as the MgO board. Another crucial issue is the influence of the molding process on the foam density distribution, which has a direct impact on its strength and stiffness properties [36]. Three different densities of EPS were recognized in the tested samples: 15 kg/m 3 , 19 kg/m 3 and 21 kg/m 3 .
While EPS is a well-recognized material, the MgO board is something of a novelty on the market, and its employment as facing material is the main reason why the analyzed panel can be considered as an innovative solution. Available studies on the structural behavior of sandwich panels containing this facing material are quite rare [13,18].
The facings are made of a composite material consisting mainly of magnesium oxide and magnesium chloride, with a glass-fiber mesh used as tensile reinforcement, and often some organic and inorganic additions. The tested board had a thickness of 11 mm and two glass-fiber meshes coated with a thin film of MgO on its surfaces, but in general, the exact composition, geometry, number, and placement of glassfiber meshes varies significantly depending on the producer. It is crucial to note that the resulting quality differs as well, and this factor needs to be taken into consideration in panel design.
Most prominent advantages of the MgO board are superior fire resistance, flexibility, and high compressive and tensile strength [17]. Moreover, it is considered environmentally sustainable, i.e., it reveals multiple "green features": magnesia is a natural material, the board production is an energy-efficient process, they are recyclable, biodegradable, and do not produce any toxic emissions [1,17]. When combined with EPS, the application of the MgO board makes the CSIP, a basis of a structurally efficient and cost-effective modular panel system.

Experimental analysis
A series of flexural tests on CSIP samples of varying geometries and load distributions was performed to obtain experimental data for (1) FE model parameter identification study and (2) for validation of the final version of the model. The experimental study was realized in two stages, corresponding to the stated aims.
Small-scale laboratory tests were planned to produce a range of different flexural responses and failure modes to use as reference data for parameter identification analysis, and the samples' geometries were assumed with this diversity in mind. Two types of load distributions and three beam lengths were tested, resulting in six test setups in total. The longest beams were meant to present typical flexural behavior, and the reduced length samples were prepared to demonstrate the increasing intensity of the core crushing effect.
Full-scale flexure tests were performed on intact wall panels supplied by the producer. Experimental data, obtained in the form of test curves and failure mode observations, were used to validate the FE model after parameter identification study was completed.

Small-scale tests
Simple supported flexure tests were performed on CSIP beams of identical cross-section and three different lengths, L, with the actual support span, L 0 , reduced by 50 mm in all cases (Table 1). Two beams of each length were tested by four-point bending (4PB), with two-point loads and two reaction forces equally spaced from each other (Fig. 2a), and by three-point bending (3PB), with one concentrated load located in mid-span and two reaction forces at the supports (Fig. 2b). Test stands and the overall procedure were based on ASTM C393. The dimensions of the beams were adjusted to produce a range of behaviors necessary for a comprehensive parameter identification investigation. The length of L1 samples corresponded to the longest span available in the test machine and was aimed to produce a global flexural failure while minimizing the core crushing effect. The length of L2 samples was reduced to the half of the original value to induce a behavior dominated by core shear, and the length of L3 samples was reduced even further to cause a strong core crushing effect. The length direction of each sample was parallel to the length direction of a respective source panel, due to the fact, that it is the main direction of flexure for CSIPs used in actual structures. Samples were cut from two different source panels and a variation in their EPS core density was noted: 21 kg/m 3 (EPS21) in the first one, and 19 kg/m 3 (EPS19) in the other. Details of the geometry and numbers of tested samples, n, are presented in Table 1. Experiments were performed in the temperature of (23 ± 2) °C with the air humidity of (50 ± 5) %, on Instron 5569 test-machine with a load cell of 50 kN capacity and a measurement accuracy of ± 0.5%. Steel support and loading cylinders with a diameter of 75 mm were used in both 3PB and 4PB tests. In 4PB tests a spreader beam and a pair of cylinders with a diameter of 25 mm were used additionally to split the load. To avoid local core crushing, steel distribution plates with a width of 60 mm were placed at all points of load application (Fig. 2).
All tests were performed using displacement control with a constant speed of 6 mm/min applied to every sample until failure. Force-displacement curves were recorded continuously as a movement of the traverse with a corresponding value of a resulting reaction force. It is important to note, that the displacement recorded should not be equated with a mid-span deflection of the bottom facing. This has been taken into account throughout the whole analysis, by treating small-scale test deflection as traverse displacement. As the failure initiated, subsequent phases of the failure process were observed and recorded for every sample. Experimental results are presented in comparison with numerical results in Sect. 5.1.

Full-scale tests
Flexure tests were performed in two repetitions on naturalscale CSIPs with the same thickness and layer arrangement as in the small-scale samples. The bending load was realized as four line-loads of equal intensity, distanced from each other by L 0 /4 ( Fig. 3) to create a state corresponding to a uniformly distributed load. Both tested panels have been provided by the producer and had identical geometries. The design of the test procedure and the dimensions of tested samples (full-sized panels) were prepared in accordance with ETAG 016 guidelines. In both cases, the core material was EPS with a density of 15 kg/m 3 (EPS15). Their detailed description is presented in Table 2.
Both test executions were performed in laboratory conditions with an Instron 500 kN load cell and a massive steel test stand (Fig. 4). The stand comprised of two steel bases ensuring simple support conditions with a span length L 0 = 2.2 m, and a loading assembly. All loading and supporting parts were mounted on joints so that all forces could remain normal to the panel's surface throughout the whole test. To prevent local damage, 60 mm wide distribution pads were used in all contact zones between the test stand and the sample (Fig. 3, 4).
The tests were performed using displacement control with a constant speed of 3 mm/min in case of the first panel, and 0.1 mm/min in the second case. The difference in loading rates did not cause any discernible changes in the test results [33]. Visual recording was employed to observe changes in the panel's geometry and to capture the failure mode. A continuous recording of local displacement and strain at several points on top and bottom surfaces of each CSIP was performed at selected locations. Linear variable differential transformer (LVDT) sensors were used to measure the vertical displacement and strain gauges (SG) bonded to the facing surface were applied to obtain longitudinal strain values. Sensor placement is presented in Fig. 3b and the number of readings performed at the most relevant locations is summarized in Table 3. From a total number of 28 SG data series across both tests, 4 were excluded due to a sensor malfunction. The results are presented in comparison with computational outcomes in Sect. 5.2.

Numerical model
A nonlinear, high-fidelity FE model able to provide an accurate prediction of failure modes in CSIPs of different geometries and under various flexural loads was created as the outcome of the present research. The major advantage of the proposed model was introducing the author's procedure into the framework of ABAQUS software [37], allowing to account for the bimodular nature of both EPS and the MgO   board. This enables the model to capture the reduction of core stiffness in areas under compression while maintaining higher stiffness in areas under shear in a straightforward, automated manner. Accounting for the complex interaction between normal and shear stress feature is crucial in designing panels with high core thickness to panel length ratio, prone to failure modes initiated by the core crushing effect.
To effectively support the design process, a macro-scale approach was assumed to make the model a feasible tool applicable even to problems with limited access to material property data. Because the most unfavorable failure conditions at the design stage can be expected in single panels rather than in structural assemblies, the numerical model domain in this research was considered as a stand-alone three-layered sample. The scope of the numerical analysis was determined based on the experimental tests performed. It consists of: (1) a series of small-scale test simulations executed to perform parametric identification of the hardening parameters, and (2) natural-scale test computations performed to validate the FEM solution.

General assumptions
Since the numerical model aims at predicting failure modes, it needs to be able to reproduce effects observed during the course of the experimental study and having a direct influence on the composite's behavior.
Pronounced core deformation caused by concentrated forces and observed during the small-scale tests appeared as a key factor influencing the way in which the composite fails when the core thickness is increased at the design stage. A considerable thickness of the facing material and a brittle nature of its failure indicates that a detailed approach is necessary. Therefore, to present a comprehensive description of how the composite layers deform, a continuum approach was adopted.
Both facings and the core were considered as separate, homogenous constituents. This macro-scale approach does not take into account the presence of a glass-fiber mesh reinforcement and the influence of adhesive permeation into adjoining layers that creates an interfacial transition zone between the core and facings. The lack of data for a detailed description of such meso-scale effects made it necessary to introduce a third constituent: a reinforced band located on the edge surfaces of facing layers where some MgO board material parameter values are being increased in an automated manner, exclusively in regions under tension. The role of this component is to prevent premature failure caused by local strain concentrations near loading points and supports. The constituents are prescribed to their respective sections in a single FE mesh (Fig. 5), meaning that adjacent layers share the same nodes. The numerical model does not allow for slip or delamination behavior, as no such effects were observed in any of the experimental tests prior to failure initiation. Post-failure delamination can be observed in some of the small-scale samples (Figs. 9, 13); however, since it occurred after the samples lost their load-bearing capacity, this effect was not addressed in the numerical model.
All tested samples were modeled adopting a plane stress state conditions. This assumption allows for a considerable reduction of computational power needed for the simulations and is entirely justified by the loading conditions causing only cylindrical bending in all cases. CSIP samples were discretized with 4-node plane stress elements with reduced integration and hourglass control. Loads and supports were modeled with rigid bodies. Loading plates were represented as 60 mm long straight lines, and cylindrical supports as 10 mm long straight segments with reference points placed on perpendicular bisectors, at a distance of 37.5 mm (cylinder radius) to avoid surface-vertex interaction (Fig. 6a,  b). Frictionless contact was postulated between deformable sections and rigid objects (Fig. 5). By taking advantage of the symmetry of the problem, only half of each sample was modeled. Based on mesh density studies, a constant mesh with regular geometry was used in all simulations (Fig. 5): 0.5 × 2 mm in facings (22 elements through layer thickness) and 2 × 2 mm in the core (76 elements through layer thickness). Relative differences in the values of extreme vertical displacements, u z , and normal stresses, σ x , while compared to the results obtained for the mesh twice as dense, did not exceed 3% in the convergence study [34]. The boundary conditions applied in 4PB, 3PB, and full-scale test simulations can be seen in Fig. 6. In case of all 4PB simulations, a preload corresponding to the weight of the spreader beam ( Fig. 2a) was considered additionally. A displacement control was used in all small-scale test simulations, whereas a force control was applied in full-scale test simulations to ensure an equal load intensity in all loading zones.
All small-scale test curves, as presented in Fig. 14, indicate at a nonlinear type of a response in a substantial range before failure load is reached. In L2 and L3 samples, a very clear core crushing effect above the supports was observed both in 3PB and 4PB tests (Figs. 9, 10, 12 and 13). This indicates that the inclusion of material nonlinearity is a vital factor in a correct failure mode prediction. Moreover, significant deformations observed in the soft core imply that geometrically nonlinear effects play an important role as well [25,36]. For these reasons both types of nonlinearity were taken into consideration in the numerical analysis.
Finally, both EPS and MgO board showed bimodular behavior, which means they have different moduli of elasticity in different stress states, most notably in compression and tension. Moreover, other material parameter values, defining inelastic and failure initiation ranges, depend on the stress state as well. In the case of EPS the origin of this phenomenon lies in its microstructure: cellular walls buckle in the elastic range under compression and stretch until failure in tension [32]. For MgO board this behavior resulted from its composite nature: in compression the load is carried by the matrix alone, while in tension the glass-fiber reinforcement activates as well. This behavior has a pronounced influence on the failure mode of the analyzed CSIPs and made creating a numerical model producing a reliable response in different loading conditions a fairly challenging task. The proposed solution to this problem is an author's procedure that identifies stress-state in each integration point in an automated manner (see Sect. 4.2).
The procedure updates stress-state classification changes at the beginning of each increment, so allowing the solver to increase the increment size whenever the solution convergence improves, may result in a potentially significant stressstate shift oversight. For this reason, it is recommended to use small, constant increments with this approach. Moreover, the first increment in an analysis starts with the whole sample being in a default state (SSV = 0, see Sect. 4.2). That speaks for the use of an even smaller load increment at the beginning of the simulation to let the material parameter assignment initiate without causing any significant inaccuracies. All FEA results presented in the current work were computed for an initial increment size of 0.001, and the maximum increment size of 0.01. The use of the procedure itself does not seem to cause noticeable convergence issues. However, whenever the failure initiation criterion nears fulfillment, severe convergence drops occur, causing the analysis to terminate. In order to assure that the solution reaches the failure initiation point before termination, an automatic stabilization option with a damping factor of 1E−9 was enabled. The damping value was assumed small enough not to cause any noticeable changes in the solution prior to failure initiation.

User-defined procedure
To provide a correct response of a material with bimodular properties, several different descriptions of a given material model can be assigned to portions of the numerical sample subjected to different stress states. In the case of simple geometries and loading conditions, e.g. uniaxial compression or tension, appropriate description can be prescribed manually. However, when flexural behavior of CSIPs is considered this approach is highly problematic and unfeasible. For this reason, a procedure that automates the process in the course of the analysis was created.
The user-defined procedure operates on stress state results obtained from the numerical model at the beginning of each increment. It introduces an additional variable using the following set of equations:  (1) to calculate a stress state variable (SSV) at every integration point in a selected portion of the numerical model (the whole CSIP sample in the present research). This results in a continuous distribution of SSV values from -1 to 1. The distribution updates at the start of each subsequent increment, therefore, sufficiently small increments need to be taken. The procedure can be used with any material model that allows to define a dependency from an additional variable (e.g. temperature), and any material parameter can be defined for any number of selected SSV values. If a SSV value at a given integration point is not equal to any of the values specified in the material model definition, then the material property value is being assigned using linear interpolation.
In the present analysis, MgO board material models were defined with two datasets corresponding to compression (SSV = -1) and tension (SSV = 1), and EPS material models were defined with an additional, third dataset, corresponding to shear (SSV = 0). A summary of material parameter dependence from selected SSV values is presented in Table 4.

Material model description
Material models for all layers of the analyzed CSIP were defined in three ranges: linear isotropic elasticity, plasticity with isotropic hardening, and ductile failure initiation. The Drucker-Prager hyperbolic yield criterion was considered sufficiently general and was used to describe both the EPS foam's and the MgO board's plastic response. Since all flexure test samples lost their load-bearing capacity when the initial crack formed, the fulfillment of a basic ductile damage failure initiation criterion [37] worked as an analysis termination condition. As the finding of the failure initiation point was enough to identify the failure mode and the ultimate load, no damage evolution analysis was performed.
Due to substantial scatter of material parameter values of the MgO board [33,34], two variants corresponding to maximum and minimum limits of experimental strength were created; they are referred in the text as "MgO min" and "MgO max". Reinforced facing edge band descriptions were used in these two versions as well. Three different versions of EPS model were defined, corresponding to different densities of the core material used in small-scale samples and full-scale panels as described in Sect. 3.
The summary of material parameter values used in the FEA is presented in Table 4. Fundamental values were obtained from laboratory tests performed in accordance with corresponding standards. Modulus of elasticity, E, and yield stress, σ pl , values of the MgO board were established in procedures based on ASTM C364 (edgewise compression) and EN 12467 (flexure). The same data for EPS21 Table 4 Material parameter values used in FEA E modulus of elasticity, υ Poisson's ratio, σ pl yield stress, E pl modulus of hardening, β angle of friction, ψ dilation angle, p t0 initial hydrostatic tension strength, ε pl,eq equivalent plastic fracture strain, η stress triaxiality factor, * parameter identification study result Material model and EPS19 were obtained from tests described in EN 826 (uniaxial compression) and EN 1607 (uniaxial tension). In the case of EPS15, these properties were calculated from density-based equations proposed in [38,39]. A literature study was also used to assume Poisson's ratios, υ [40,41], as well as angles of friction, β, and dilation angles, ψ [42,43]. More detailed reasoning behind the specified data can be found in research presented in [33][34][35].
The material description was developed further, as a part of the current study, with the addition of SSV value corresponding to the core shear stress-state in the elastic range.
The E values for EPS in the SSV = 0 case were calculated following ASTM C393 using flexure test data from L1 and L2 beams.
The remaining hardening properties, σ pl and E pl , were established in the current work by performing a parameter identification study. Different sets of values were tested until a satisfactory resemblance was obtained between all 6 smallscale test simulation and experimental results. This approach was used to determine values describing the hardening behavior of EPS and the MgO board (marked * in Table 4).
The ductile failure criterion was expressed in terms of the equivalent plastic fracture strain, ε pl,eq , and the stress triaxiality factor, η, defined as a negative ratio of hydrostatic pressure stress to equivalent von Mises stress [37]. Failure initiation criterion parameters were retrieved from numerical models of small-scale samples subjected to loads that caused the failure in corresponding test samples [34].

Results comparison
Three relevant types of results obtained for all samples in the course of experimental tests and numerical simulations are presented in the following section. Stress state variable (SSV) distribution maps, acquired from FEA through the application of the author's procedure, show how the model identifies areas recognized as being under compression (SSV = -1), shear (SSV = 0) or tension (SSV = 1) stress state. Failure modes established in laboratory tests are compared with failure initiation regions recognized with the FE model and presented here in a form of distribution maps of a damage initiation criterion variable (DICV). DICV denotes the degree to which the damage initiation criterion is matched-with numbers 0 for the integration point which does not yield at all and 1 for the case when the criterion is fulfilled. Since damage evolution is not a part of this analysis, the criterion fulfillment is interpreted as the end of each simulation. Finally, force-displacement curves recorded in experimental tests on CSIP samples are compared with the corresponding numerical outcomes. In the case of naturalscale CSIP flexure tests, force-strain curves are additionally examined.

Small-scale sample flexure tests
The final version of the FE model, obtained after completing the parameter identification process, was used to perform numerical simulations of flexure tests for all small-scale CSIP samples. A comparison of computational results with experimental data is presented.
SSV distributions obtained at the moment of failure initiation in FEA of all CSIP beams are presented in Fig. 7. Color-coded maps illustrate the stress state classified in the final increment of the numerical analysis at each given integration point. For the longest beams, the stress state in the top facing was identified as dominated by compression, whereas the bottom facing was recognized as being mainly under tension (cf. Fig. 7a, d). In both facings, local zones of the opposite sign formed near concentrated loads. In the case of samples of a reduced length, the SSV distribution notably changed. The top facing in L2 samples (Fig. 7b, e) can be still recognized as being mainly compressed, however, the SSV distribution in the top facing of L3 samples (Fig. 7c, f) shifted towards almost balanced proportions between compression and tension. The bottom facing in L2 and L3 samples is under compression in the upper section and under tension in the lower section. This indicates that the model identified the longest beams as being under global bending, while in the shorter samples, local bending of both facings became more pronounced. Stress state distribution in the core distinguished compressed areas over the supports and under point loads. In L1 beams (Fig. 7a, d), these areas are limited, however, in 4PB test (Fig. 7a) the compressed area under loading points is visibly more extensive and stretches throughout the whole central region of the beam. The largest part of the core was recognized as being under shear stress. In L2 samples the proportions between compression and shear become more even and in L3 samples compression dominates in the core. The presented results indicate that the author's procedure yielded an anticipated physical response.
FEA results for all CSIP samples are compared with failure modes observed in laboratory tests in Figs. 8, 9, 10, 11, 12 and 13. A pair of numerical results is presented for each case, corresponding to the minimum and maximum strength variants for the MgO board material model. DICV distribution in all samples subjected to 4PB exhibits similar behavior. In all cases, FEA resulted in a tensile failure of the bottom facing, near the support. This type of failure coincides with experimental observations in the case of L2 (Fig. 9) and L3 (Fig. 10) samples. It is worth to mention, that delamination visible in Fig. 9c occurred postfailure and had no influence on the sample's load-bearing capacity. The flexural failure of the longest beam (Fig. 8) occurred in the bottom facing under one of the loading zones, approximately in one-third of the support span. In this case, the computational model predicted a longitudinal strain peak in the correct location until the load intensity reached half of the ultimate value. Then, a plastic yield zone in the core extended over the support area, causing a strain redistribution in the bottom facing. As a consequence, the strain maximum was shifted and a tensile failure occurred at the correct layer but considerably too close to the support.
The difference between results obtained by applying MgO min and MgO max facing material definitions is not significant and is mainly manifested by a formation of increased DICV value regions in the compressed parts of the bottom facing over the support and the top facing under the loading point for MgO min case. This behavior was observed in all 4PB simulations and did not lead to a failure mode change between MgO min and MgO max variants in any of the tests. Additionally, there is a visible difference in the DICV distributions in the core of the L1 beam (Fig. 8). It resulted from a higher bending stiffness of the facings in the MgO max variant, which led to an intensification of plastic strain in the core, translating into an expansion of increased DICV value region.
Simulations of 3PB tests of CSIP samples revealed a more diverse behavior. In the longest samples (Fig. 11) the MgO min variant led to a compressive failure of the top  Fig. 13c was an effect of a post-failure deterioration. A compilation of force-displacement curves for all CSIP samples is presented in Fig. 14. Experimental data were recorded as pairs of the traverse displacement and the corresponding reaction force, and computational results were read at the same location in the numerical model. Individual experimental curves were averaged and resampled together with the numerical ones to acquire It should be noted that all r 2 values are very close to 1, which means that the numerical curves show a very strong resemblance with the experimental ones. Failure initiation conditions were based on the results of 4PB tests of the longest sample and for this reason, failure loads in numerical results coincide with the averaged experimental result (Fig. 14a). The computational value of failure load for a 4PB L2 sample (Fig. 14c) shows also a very good agreement with the experimental result. In the case of 3PB tests of the longest sample (Fig. 14b), on the curve representing results of the computational simulation for MgO max variant, there are additionally marked points indicating shifted positions of DICV extreme value. Ultimate load values for MgO min and MgO max variants create a relatively close range around the averaged experimental value in the case of 3PB test of L1 and L2 samples (Fig. 14b,  d). The numerical model slightly underestimated the failure load values for the L3 samples in both 3PB and 4PB tests (Fig. 14e, f).
Experimental data were converted into facing extreme normal stress with classical sandwich theory (CST) Eq. (2), assuming that bending stiffness of the core can be neglected in the flexural rigidity [7]: where: M y is the bending moment (3PB test: M y = FL 0 /4, 4PB test: M y = FL 0 /6), remaining symbols in accordance with Fig. 1. The lowest relative difference value obtained for the longest 3PB test sample is over 30%. This indicates, that even the longest CSIP beams tested in small-scale experiments do not satisfy CST assumptions and that a more advanced computational model has to be used.
A summary of both numerical and experimental test results is presented in Table 5. It can be seen that failure modes predicted by the FE model agree, completely or partially, with experimental observations in most simulations. Relative error values of ultimate load, δF z u , are within a reasonable range, with a few notable exceptions. All coefficient of determination values, r 2 , are very close to 1, which means, that the shape of force-displacement curves obtained from calculations is very similar to averaged experimental results. The summary shows that the proposed model is able to reproduce the behavior of tested samples in most aspects, however, it can be also seen that there are a few issues that require further attention. Numerical simulation results in the form of DICV maps for both MgO min and MgO max variants (Fig. 16) indicate a tensile failure initiation in the bottom facing near the center of the panel. This is in agreement with the failure mode observed in both laboratory tests. It is important to note that the numerical model showed a second potential tensile failure zone formed in the bottom facing over the support. This observation suggests that if the core thickness of the panel was to be increased, failure mode would presumably shift to the support area, which is a physically anticipated behavior as well.
It was also noted that the values of ε pl,eq of the MgO board under tension, based on small-scale flexure tests, led to an overestimation of the failure load value. Based on the observations made during the laboratory tests, this value was reduced so that damage occurred immediately after the facing material yielded in the tensile zone. The values were In the case of the full-scale flexure test two types of data curves were obtained, both experimentally and numerically. Force-displacement curves were recorded in experimental tests as readings from a load cell and linear displacement sensors located under the bottom facing, while the forcestrain response was additionally recorded from strain gauges on both top and bottom facings (Fig. 3b). For conciseness, computational results were evaluated from the numerical model only at the nodes corresponding to the most relevant locations on the tested panels (Table 3).
A comparison of experimental and numerical force-displacement curves obtained at measurement points located in L 0 /2 and L 0 /4 are presented in Fig. 17. It can be seen that both numerical variants show a very close curve shape similarity with the experimental data, with r 2 values nearly equal 1. The average numerical deflection is slightly underestimated in L 0 /2 and slightly overestimated in L 0 /4, showing, that the predicted numerical response is very similar but not strictly consistent with the experimental observations. The same conclusion can be drawn from the force-strain curve comparison presented in Fig. 18. Results obtained from MgO max computational variant are very close to experimental readings in the top facing, and MgO min variant response is similar to the experimental results in the bottom facing, however, the averaged numerical response underestimates the experimental results in the first case and overestimates them in the second one. A summary of r 2 values for all curves is presented in Table 6 and a comparison  Table 7.

Discussion
A comparison between experimental data and computational results obtained with the use of the proposed numerical model shows a very good qualitative resemblance in all small-and full-scale tests. The outcomes indicate that both the parameter identification process and the FE model validation were successful.
All distributions of SSV values obtained due to incorporation of the author's procedure provided a reasonable interpretation of stress states, sensitive to fundamental changes in CSIP sample's geometry and load distribution. The method produced physically sound results, which were employed for automated control of versatile, stress state based, nonlinear definitions of core and facing material models.
Computational failure mode predictions in the form of DICV maps are generally in agreement with experimental observations. However, few inconsistencies need to be addressed. The simulation of the 4PB test of the CSIP L1 sample identified the failure initiation in the correct layer but positioned it next to the support, in a considerable distance from where it occurred in two laboratory executions. By judging from the strain redistribution in the bottom facing that followed core yielding, it is clear that this inconsistency is caused by the way the core deforms in the plastic region of the model. The problem might lie in a discontinuous transition from plastic compression to plastic shear, pushing the failure zone too close to the support.
Three failure stages observed in the 3PB test simulation of the CSIP L1 beam show that the equivalent plastic fracture strain value used in small-scale tests was too high in the case of this test. Failure mode obtained at stage (I) of failure agrees with the experimental results. It is interesting to notice that the correct failure mode of the full-scale panel occurred for an even smaller value of ε pl,eq , almost immediately after the facing material yielded in the tensile zone. It shows that the transition of this variable from small-scale to full-scale simulations requires special attention and gives ground for further research.
The comparison of force-deflection and force-strain curves for both types of tests shows that the numerical results are very accurate at the point of load application in all small-scale experiments and the L 0 /2 position in natural scale panel simulation, with just a slight discrepancy in the L 0 /4 measurement point. The panel's flexural deformation obtained with FEM is very similar to the experimental Results obtained from the CST equations for the natural-scale panels are in reasonable agreement with both experimental data and numerical outcomes. However, the ultimate stress values obtained for small-scale samples were significantly overestimated in comparison with FEA results. It is an anticipated observation, and it shows that CST should be used in CSIP flexural strength tests with caution.

Conclusions
The comparison studies showed that the proposed numerical model is a robust computational tool, able to capture the critical aspects of analyzed CSIP's behavior based on data obtained from basic, small-scale tests. It is due to the application of the author's procedure, which allows to include bimodular effects in the nonlinear material model description, by distinguishing areas being under different stress states with high accuracy. It enables an automated prescription of varying parameter values in each integration point and updates the distribution with each increment. The study also showed that adding this functionality is essential in a reliable prediction of failure modes of the analyzed CSIP.
Some inconsistencies were observed between the computational results and the test data. The most probable cause is an insufficiently accurate description of the core deformation in the plastic range, caused by the limited small-scale experimental data. A wider variety of measurements, such as facing strain, and visual recording of the core's deformation, could be used for further improvement of the model.
The study also showed that the ultimate stress values calculated for the small-scale samples from CST were significantly overestimated. It implies that this popular engineering method is not a reliable tool in the estimation of the analyzed CSIP's flexural strength. The use of a soft foam core, and the designers' inclination to maximize its thickness for improved thermal insulation and stiffness,    often violates the assumptions of CST and more advanced computational tool needs to be used. The proposed model, obtained by the application of the author's procedure, proved to be a high-fidelity numerical tool, able to aid the design process of analyzed CSIPs and to reduce the extent of necessary laboratory testing. It can be easily adapted for the analysis of other sandwich panels that require a similar level of modeling precision. It also has no limitations preventing it from being used in combination with contact interactions, to account for delamination, or from application in three-dimensional problems with solid only and solid-shell approaches.   Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.