*Article* **Estimation of the Compressive Strength of Corrugated Board Boxes with Shifted Creases on the Flaps**

**Damian Mrówczy ´nski 1 , Tomasz Garbowski <sup>2</sup> and Anna Knitter-Pi ˛atkowska 3, \***


**Abstract:** In the modern world, all manufacturers strive for the optimal design of their products. This general trend is recently also observed in the corrugated board packaging industry. Colorful prints on displays, perforations in shelf-ready-packaging and various types of ventilation holes in trays, although extremely important for ergonomic or functional reasons, weaken the strength of the box. To meet the requirements of customers and recipients, packaging manufacturers outdo each other with new ideas for the construction of their products. Often the aesthetic qualities of the product become more important than the attention to maintaining the standards of the load capacity of the packaging (which, apart from their attention-grabbing functions, are also intended to protect transported products). A particular flaps design (both top and bottom) and its influence on the strength of the box are investigated in this study. An updated analytical–numerical approach is used here to predict the strength of packaging with various flap offsets. Experimental results indicated a significant decrease in the static load-bearing capacity of packaging in the case of shifted flap creases. The simulation model proposed in our previous work has been modified and updated to take into account this effect. The results obtained by the model presented in this paper are in satisfactory agreement with the experimental data.

**Keywords:** corrugated board; box strength estimation; packaging flaps; crease line shifting

## **1. Introduction**

The relentless increase in consumption all around the contemporary world is reflected in the significant growth in the production of various goods. This, in turn, entails the necessity of their packing, safe storing and transportation to any destination. Due to growing ecological awareness and concern for the environment, the perfect choice is undoubtedly corrugated cardboard boxes. The undeniable facts are that they are recyclable, easy for disposal, ecological, durable under appropriate conditions and easy to store in a flat form after manufacturing. Among their numerous advantages, one cannot fail to mention the easy imprint of brand names on them. This is highly useful in cases of shelfready packaging (SRP) or retail-ready packaging (RRP) when, after being transported to the site, the packaged products are placed directly on the shelves. Upon opening the cardboard boxes along the specially designed and made perforations, products are ready to purchase. Such a solution is a huge time-saver for large companies.

In the case of individual recipients of merchandise, especially when shopping online (which nowadays is a significant part of the sales market), a very important factor is the possibility of smoothly returning purchased products if the consumer is not satisfied, for a range of reasons. Retailers that offer reusable packaging to send back purchased goods are very competitive on the market. Again, corrugated cardboard boxes are perfect in such situations. They are easy to open thanks to well thought-out perforations and, after

**Citation:** Mrówczy ´nski, D.; Garbowski, T.; Knitter-Pi ˛atkowska, A. Estimation of the Compressive Strength of Corrugated Board Boxes with Shifted Creases on the Flaps. *Materials* **2021**, *14*, 5181. https:// doi.org/10.3390/ma14185181

Academic Editor: Dimitrios Tzetzis

Received: 29 July 2021 Accepted: 6 September 2021 Published: 9 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

re-sealing with the built-in adhesive strip, are ready to send back. However, it must not be forgotten that the packaging must have sufficient durability to survive the return transport.

Therefore, in view of the above, scientific research while applying analytical as well as numerical methods and/or laboratory tests has been an inherent part of a separate branch of industry, i.e., the production of corrugated cardboard packaging, for many years.

The proper mechanical strength of the paperboard or corrugated cardboard boxes is directly connected with two characteristic in-plane directions of orthotropy. Machine direction (MD) is perpendicular to the main axis of the fluting and parallel to the paperboard fiber alignment, whilst cross direction (CD) is parallel to the fluting. In order to examine the strength of corrugated cardboard boxes, one can perform some fundamental physical tests, i.e., compressive, tensile or bursting strength tests, which, in practical terms, are the most significant. The most prevalent are the box compression test (BCT) and the edge crush test (ECT) for corrugated cardboard.

A significant impact on the load-bearing capacity of packages is undoubtedly the various perforations, openings and flap locations on corrugated cardboard boxes. The first two issues have been meticulously discussed by Garbowski et al. in [1] and [2], respectively. In the present study, the influence of the flap locations on the strength of corrugated cardboard boxes, as another article in a series, is discussed. The conducting of physical experiments usually involves a great deal of time and cost. Therefore, recently, other methods of testing corrugated boxes have emerged to determine their strength by physical testing only.

Alternatively, the compressive strength of boxes can be assessed based on formulae that have been presented in numerous literatures. Their adoption, thanks to their simplicity, results in quick and easy solutions for practical applications. Moreover, no additional experiments are necessary. The parameters that are introduced in these formulae can be systemized into three groups: paper, board and box parameters [3]. In the first group one can specify: the ring crush test (RCT), Concora liner test (CLT), liner type, weights of liner and fluting, corrugation ratio and a constant related to fluting. In the second one: thickness, flexural stiffnesses in MD and CD, ECT and moisture content. Finally, in the third: dimensions and perimeter of the box, applied load ratio, stacking time, buckling ratio and printed ratio. Nearly 70 years ago, the paper (RCT, flute constant) and box (perimeter, box constant) parameters were applied for the prediction of boxes' compressive strength in the formula presented by Kellicutt and Landt [4]. The dependence of critical force on paper parameters (CLT, type of liner) and cardboard box dimensions in the BCT was presented in [5].

Generally applicable in the packaging industry is the procedure proposed by McKee et al. [6], in which the parameters of the paperboard (ECT, flexural stiffnesses) and the box perimeter were introduced. Nevertheless, the provided formula is applicable only for comparatively simple boxes. Throughout the years, many scientists endeavored to broaden the applicability of the McKee's analytical formulae. Allerby et al. [7] modified the constants and exponents in the above-mentioned approach. Schrampfer et al. [8], in turn, amended McKee's approach by extending the possibility of implementing a broader range of cutting methods and equipment. Batelka and Smith [9] enhanced the relationship with the dimensions of the box and Urbanik and Frank [10] introduced the Poisson's ratio as well. The arbitrary chosen constant value as a parameter in the McKee's formula limited its applicability to simple standard boxes. Moreover, Garbowski et al. [1,2,11] examined this approach for more sophisticated cases and modified the McKee's formula. One cannot forget that the compression strength of corrugated paperboard boxes [12] depends on many factors, such as moisture content of the box [13,14], the presence of openings, ventilation holes and perforations [1,2,15], storage time, stacking conditions [16] and numerous others.

An alternative option to compute the strength of the boxes is to implement, the well-known in engineering, finite element method (FEM). It has been involved in a lot of research, including the problems of numerical analysis with regard to the transverse shear stiffness of corrugated cardboards [17–21] as well as buckling and post-buckling phenomena [22]. The method that efficiently allows one to simplify the examined models is homogenization [23–27]. The result of this procedure is one single layer described by the effective properties of the composite, rather than building the layers made out of different materials. The advantage of this approach to the problem is a significant saving of calculation time while maintaining the appropriate accuracy of the results. The approach based on strain energy, applicable to sandwich panels in the issue of homogenization, was presented by Hohe [28]. For this purpose, a representative element of the heterogeneous and homogenized elements was proposed. Another method, using a periodic homogenization technique considered by Buannic et al. [29], allows not only for an equivalent membrane and the pure bending characteristics of period plates but also, in a modified version, includes the transfer of shear effect in the analysis. The FEM was applied by Biancolini [30] for the examination of a micromechanical part of the considered plate. In the aftermath of application, the energy equivalence between the model and the equivalent plate as well as the stiffness properties of the sandwich plate were obtained. In turn, Abbès and Guo [31] analyzed the plate, which was decomposed into two beams in the directions of the plate, which allowed them to find the torsion rigidity of the orthotropic sandwich plates. The method of treating the quasi-static equilibrium of a material subjected to deformation with hardening was proposed in [32]. Therefore, the experimental data obtained in the dynamic case of deformation could be compared with the data calculated for the quasi-static case. The laboratory tests, properly chosen and scheduled, were performed right on the composite. Layered elements, on which effective parameters can be measured directly, are an alternative method for homogenization. This very approach is proposed in the present research.

An operation during which fold and perforation lines are introduced is defined as creasing. One cannot neglect its impact on the load-bearing capacity of corrugated paperboard. Undeniably, those lines reduce the mechanical strength of the manufactured corrugated paperboard boxes, hence the results of extensive research can be found in the literature. The comparison between the experimental and FEM numerical results, performed in order to examine the creasing influence on the local strength of corrugated paperboard, was discussed by Thakkar et al. [33]. The impact of creasing and subsequent folding on the mechanical properties of laminated paperboard has been picked up by Beex and Peerlings [34], who performed physical as well as numerical experiments, whilst Giampieri et al. [35], to acquire the mechanical response of creased paperboard after folding, used a constitutive model. Domaneschi et al. [36] and Awais et al. [37] proposed an essential (from a practical point of view) solution for the packaging industry, basing it on the FEM simulations of paperboard creasing. Experimental, as well as numerical, studies on the influence of the creasing process during press forming on the paperboard mechanical properties were conducted by Leminen et al. [38].

The particular top and bottom flaps design, which is directly related to the flap creases, and their influence on the strength of corrugated cardboard boxes is investigated in this study. An updated analytical–numerical approach is used to predict the strength of the packaging with various flap offsets. Experimental results pointed out a significant decrease in the static load-bearing capacity of packaging in the case of shifted flap creases. The simulation model, proposed in the previous works of the authors [1,11], has been modified and updated to take this effect into account as well. The results obtained during the analysis of the numerical model proposed in the paper are in adequate agreement with the experimental data. This approach by which the prediction of the strength of boxes with offset flaps is analyzed is, to our knowledge, very pioneering and constitutes an innovative contribution to the development of the field related to the prediction of the load capacity of corrugated cardboard packaging.

#### **2. Materials and Methods**

#### *2.1. Corrugated Board Packaging with Shifted Flaps*

In previous works, the authors analyzed packages with perforations [1] and openings [2]. Here, the focus is on packages with offset flaps. Such packaging is becoming standard in retail-ready packaging that is also used for shipping. The shifting of the crease line (see Figure 1) makes the flaps more adjustable after closing. Unfortunately, the load-bearing capacity of packaging with shifted flaps significantly diminishes.

**Figure 1.** Box with offset flaps.

The drop in box strength results from a certain sequence of loading, in which the edges of the two shifted (elongated) walls of the package are loaded first, while the other two are only loaded after buckling and/or crushing of the first two (see Figure 2).

**Figure 2.** Sequence of loading of the package vertical walls: (**a**) edges loaded in the first step; (**b**) edges loaded in the second step.

In order to derive a simplified calculation algorithm for estimating the strength of boxes with offset flaps, a series of tests was first performed in the laboratory for various boxes made of different corrugated cardboard. All studies were carried out on the BCT press [39] (see Figure 3). In order to be able to perform computer predictions of the packaging load capacity, it is required in the first step to identify the material parameters of the corrugated board, then to select the appropriate material model and finally to build a numerical or analytical model that takes into account the geometry of the analyzed box.

**Figure 3.** BCT tests: (**a**) BCT press; (**b**) packaging with the shifted offsets on the flaps.

The following sections describe the laboratory testing of corrugated board, the constitutive modeling of corrugated cardboard, a numerical simulation model and a simple analytical algorithm for estimating the load capacity of corrugated cardboard packaging.

#### *2.2. Laboratory Testing of Corrugated Board*

Laboratory tests of the corrugated cardboard were performed to determine its stiffness and strength. The four most commonly used tests are: edge crush test, shear stiffness testing, torsional stiffness test and 4-point bending test. The edge crush test (ECT) measures the compressive strength of a corrugated board sample. This test is performed for relatively stocky specimens, so that the failure mechanism is the crushing of the sample, not the loss of stability. The ECT value is often used to determine the load capacity of the corrugated cardboard package in analytical [6], analytical–numerical [1,2,11] or purely numerical [40,41] approaches.

The shear stiffness test (SST) is used to measure the shear stiffness of a sample by applying two equal forces at opposite corners. The measurement of displacements and reaction forces on the supports enables the required stiffness to be calculated. The SST is characterized by a high sensitivity to crushing the sample, resulting in processes such as die-cutting and laminating. The torsional stiffness test (TST) consists of twisting the sample by 10 degrees in both directions and is performed to determine the torsional stiffness. Only the linear part of the bending moment/angle of rotation diagram is being considered for this purpose. The obtained TST values are valid even for highly crushed, broken and flaccid samples.

The bending stiffness test (BNT) is used to determine the bending stiffness in the 4-point bending test. The static scheme of the tested sample allows a constant bending moment and a shear force equal to zero between the internal supports to be obtained, which provides more accurate measurement of the bending stiffness value. On the other hand, the presence of a shear force between internal and external supports makes it possible to take into account the effect of the shear stiffness as well.

#### *2.3. Corrugated Board: Material Model and Constitutive Parameters*

Since paperboard is an orthotropic material, many material parameters are needed for its correct mathematical description. Therefore, more laboratory tests should be carried out. In papermaking laboratories one can determine visual, functional and mechanical properties of paperboard or corrugated board. The most popular mechanical tests include, for example: (a) short span compression test (SCT) of paperboard; (b) tensile test of paperboard; (c) resistance to bursting of paperboard or corrugated board; (d) edge crush test (ECT) of corrugated board; (e) flat crush test (FCT) of single walled corrugated board; (f) corrugated board bending stiffness (4-point bending test).

Some of these tests can be directly used for linear elastic material model calibration, namely the plane strain Young's modulus in two perpendicular directions, Kirchhoff's modulus and Poisson's ratio. The modulus of elasticity (i.e., Young's modulus) is a quantity

well known to designers and engineers, but less common in paper specifications in the cardboard packaging industry. Traditionally, the stiffness modulus can be determined while performing a uniaxial tensile test of a sample. As paperboard is an orthotropic material, more tests are required to determine all elastic parameters (see Figure 4).

**Figure 4.** The load–displacement curves in MD, CD and 45 deg.

Determining the elastic parameters is an important step in the box load-bearing capacity estimation procedure, thus the brief introduction to some basic definitions, the constitutive description of the paperboard and the method of calibrating material constants will be presented in the subsequent sections. For orthotropic materials in a plane stress state, the relationship between elastic strains and stresses can be written as:

$$
\begin{bmatrix}
\varepsilon\_{11} \\
\varepsilon\_{22} \\
2\varepsilon\_{12}
\end{bmatrix} = \begin{bmatrix}
1/E\_1 & -\nu\_{21}/E\_2 & 0 \\
0 & 0 & 1/G\_{12}
\end{bmatrix} \begin{bmatrix}
\sigma\_{11} \\
\sigma\_{22} \\
\sigma\_{12}
\end{bmatrix} \tag{1}
$$

<sup>ଵ</sup> <sup>ଶ</sup> ଵଶ ଵଶ ଶଵ where *E*<sup>1</sup> is Young's modulus in the Machine Direction (MD); *E*<sup>2</sup> is Young's modulus in the Cross Direction (CD); *G*<sup>12</sup> is Kirchhoff's modulus and *ν*12, *ν*<sup>21</sup> are Poisson's coefficients. Due to the symmetry of the material compliance/stiffness matrix, the relationship between the Poisson's coefficients is as follows:

$$\frac{\nu\_{12}}{E\_1} = \frac{\nu\_{21}}{E\_2}.\tag{2}$$

The Hill model [42] can be successfully employed to describe the behavior of the paper in an inelastic phase. Implementation of the Hill model requires the definition of the elastic domain described by the plastic yield function and the description of the material hardening:

$$\dot{f}(\sigma,\kappa) = \sqrt{a\_1\sigma\_{11}^2 + a\_2\sigma\_{22}^2 - a\_{12}\sigma\_{11}\sigma\_{22} + 3a\_3\sigma\_{12}^2} - \sigma\_0(\kappa) \le 0,\tag{3}$$

<sup>ଵ</sup> = <sup>ଶ</sup> = ଵଶ = <sup>ଷ</sup> = 1 ሺሻ ଵଶ ଶ ଶ ଶ where √ ∗ is an effective stress *σe f f* , which can be reduced to classical Huber-Mises criterion for isotropic materials if *a*<sup>1</sup> = *a*<sup>2</sup> = *a*<sup>12</sup> = *a*<sup>3</sup> = 1; *σ*0(*κ*) is a yield stress function; *κ* is a hardening parameter, usually related to effective plastic strains; *σij* are the stresses in main orthotropic directions; *a<sup>i</sup>* and *a*<sup>12</sup> are called anisotropic parameters, which can be determined from simple tensile tests in the main orthotropic directions:

$$a\_1 = \frac{\sigma\_0^2}{\sigma\_{10}^2}, \quad a\_2 = \frac{\sigma\_0^2}{\sigma\_{20}^2}, \quad a\_3 = \frac{\sigma\_0^2}{3\sigma\_{120}^2}.\tag{4}$$

ଶ ଵଶ ଵଶ where: *σ*<sup>0</sup> is the initial yield stress in the reference direction; *σ*<sup>10</sup> is the yield stress in first direction (e.g., MD); *σ*<sup>20</sup> is the yield stress in second direction; *σ*<sup>120</sup> is the yield stress in shearing.

ଵଶ = <sup>ଵ</sup> <sup>ଶ</sup> 3<sup>ଷ</sup> − 4ସହ,

The remaining parameter *a*<sup>12</sup> can be determined from the equation:

$$a\_{12} = a\_1 + a\_2 + \mathfrak{Z}a\_3 - 4a\_{45'} \tag{5}$$

where *a*<sup>45</sup> is the anisotropic parameter determined from a tensile test in an angled direction of 45 deg. As for most materials, only the values *σ*<sup>10</sup> , *σ*<sup>20</sup> and *σ*<sup>120</sup> are known, in practical applications for the coefficient *a*<sup>12</sup> usually a simplified relationship is assumed, e.g.,: ସହ ଵ ଶ ଵଶ ଵଶ

$$a\_{12} = \frac{\sigma\_0^2}{\sigma\_{10}\sigma\_{20}}.\tag{6}$$

It is a known fact that paperboard behaves differently under tension and compression. Therefore, the chosen plasticity criterion (which is symmetric in case of tension and compression) is not appropriate for this type of material. However, for simple strength calculations with a stress state dominated by compression, this model is a sufficient approximation. For the correct analysis of the structure in the complex stress state, one of the more sophisticated constitutive models should be used, e.g., [43–48].

#### *2.4. Numerical Predictive Model*

The numerical model of the box was built in the Abaqus Unified FEA software (2020, Dassault Systèmes Simulia Corp., Providence, RI, USA.) [49]. Two types of models had to be created: (i) the non-offset packaging and (ii) the package with flaps offset. In order to simplify the computations and save the computing time, only 1/8 part of the box was modeled instead of the whole packaging (see Figure 5). The material used in the model was linear elastic orthotropic model with Hill plasticity.

**Figure 5.** Scheme of the 1/8 part of the package.

To obtain the appropriate behavior of the numerical model, symmetry boundary conditions were defined on each edge (see Figure 6). For the packaging model without offset, only one computation step was defined, in which the displacement was applied on both edges. In the case of the package with offset flaps, in the first step only the offset edge was loaded and in the second step the load was then applied to the non-offset edge. The 4-node quadrilaterals shell elements with full integration, named S4, were used for all computations. For different dimensions of packaging, different values of mesh size were assumed. For example, for the package dimensions 500 × 500 × 500 mm the approximate global size of the element was 12 mm, which ultimately gave 882 elements, 946 nodes and 5676 degrees of freedom. To add the initial deformations (resulting from imperfections) of the box vertical walls, a buckling analysis was performed before the main calculations. The first buckling mode of the model found in this way was later introduced in the next step in the form of scaled imperfections.

**Figure 6.** Boundary conditions for the case of: (**a**) the non-offset package; (**b**) the package with offset flaps (first step); (**c**) the package with offset flaps (second step).

#### *2.5. Analytical Predictive Model*

 = <sup>ଶ</sup> ெ = <sup>ଵ</sup> The simplified procedure for estimating the compressive strength of a corrugated cardboard box with offset flaps proposed here is based on an analytical model. The algorithm exploits the basic constitutive parameters of a single box wall, namely: *ECTCD* compressive strength in CD, *ECD*= *E*2—compressive stiffness of corrugated boards in CD and *EMD* = *E*1—compressive stiffness of corrugated boards in MD. Since in some cases the instability of a single wall may occur before plasticization, it is also necessary to determine the critical load for an orthotropic rectangular plate, e.g., from the formula [1,2,11]:

$$P\_{cr}^{i} = \frac{\pi^2}{B\_l^2} \frac{t\_l^3}{12} \sqrt{E\_{\text{CD}}^i E\_{MD}^i} \left(\frac{mB\_l}{H} + \frac{H}{mB\_l}\right)^2 \tag{7}$$

 where *B<sup>i</sup>* is the width of the *i*-th panel; *t<sup>i</sup>* is the *i*-th panel thickness; *H* is the box height; *m* is the number of half-waves for which *P i cr* reaches the minimum.

The analysis of strength estimation of a box with shifted flaps, as already discussed in the previous section, consists of two stages, in which the higher walls (i.e., the shifted ones) are loaded first (see Figure 2a), while the lower walls are loaded only if preliminary crushing and/or buckling of the first two walls occurs (see Figure 2b). Therefore, the overall load capacity of the packaging is the sum of the load capacity of two pairs of opposite walls of the box, namely:

= <sup>ଵ</sup> <sup>ଶ</sup>

$$BCT = \alpha BCT\_1 + BCT\_2 \tag{8}$$

where

$$BCT\_1 = 2kECT^r \left(P\_{cr}^1\right)^{1-r} \gamma\_1 \gamma\_2 B\_1 \tag{9}$$
  $\text{ifted walls, while}$ 

,

<sup>ଵ</sup> = 2 is the load capacity of the shifted walls, while

$$BCT\_2 = 2kECT^r \left(P\_{cr}^2\right)^{1-r} \gamma\_3 \gamma\_4 B\_2 \tag{10}$$

is the load capacity of lower walls.

 ∈ ሺ0,1ሻ <sup>ଵ</sup> <sup>ଶ</sup> In Equations (9) and (10) *k* is a certain constant and *r* is an exponent, *r* ∈ (0, 1), and *γ<sup>i</sup>* are the reduction coefficients. *B*<sup>1</sup> and *B*<sup>2</sup> are base dimensions, which are shown in Figure 7. The *α* coefficient reduces the value of the first term in Equation (8) due to the initial failure

and/or buckling of the walls loaded in the first step (see Figure 8). This factor can be calculated using the formula below:

$$\mathfrak{a} = 1 - \frac{\mathfrak{u}\_{off} - \mathfrak{u}\_0}{\mathfrak{u}\_{max} - \mathfrak{u}\_0} \, \mathrm{} \tag{11}$$

where *uo f f* is an offset of higher walls; *umax* = *H* is assumed to be equal to the height of the box; *u*<sup>0</sup> is the vertical deformation corresponding to the maximum load. The latter can be calculated from Hooke's law considering the stiffness in the CD direction, *ECD*; single box wall height, *H* (see Figure 7); shifted wall width, *B*1; board thickness, *t*; the compressive strength, *BCT*<sup>1</sup> (see Figure 8). Thus, finally we obtain: ௫ = <sup>ଵ</sup> <sup>ଵ</sup>

$$
\mu\_0 = \frac{BCT\_1}{2tB\_1E\_{CD}}H.\tag{12}
$$

**Figure 7.** Box dimension symbols.

<sup>ଷ</sup> = min ቈ൬<sup>ଶ</sup> య **Figure 8.** Force-displacement visualization of the proposed method.

 ସ The reduction factors *γ<sup>i</sup>* are always less than one and depend on the ratio of the box dimensions and the exponents *r<sup>i</sup>* . The *γ*<sup>1</sup> factor in Equation (9) reads:

൰

, 1,

$$\gamma\_1 = \min \left[ \left( \frac{B\_1}{H} \right)^{r\_1}, 1 \right], \tag{13}$$

while *γ*2:

$$\gamma\_2 = \min\left[ \left( \frac{B\_1}{B\_2} \right)^{r\_2}, 1 \right]. \tag{14}$$

<sup>ଵ</sup> <sup>ଶ</sup> <sup>ଷ</sup>

Similarly, the coefficient *γ*<sup>3</sup> in Equation (10) is:

$$\gamma\_3 = \min \left[ \left( \frac{B\_2}{H} \right)^{r\_3}, 1 \right], \tag{15}$$

while *γ*4:

$$\gamma\_4 = \min \left[ \left( \frac{B\_2}{B\_1} \right)^{r\_4}, 1 \right]. \tag{16}$$

All unknown factors in Equations (8)–(10), namely constant *k* and exponent *r*, and the four exponents *r<sup>i</sup>* in Equations (13)–(16), can be found by calibration with experimental data. The calibration procedure will be presented in the following section.

#### *2.6. Calibration Procedure*

The main goal of this study is to propose a reliable analytical model for the quick estimation of the load capacity of offset packaging. Therefore, the calibration of the coefficients in the analytical equations is particularly important. Unfortunately, the limited number of laboratory results creates a risk that the analytical model will be valid only for a small set. In order to extend the applicability of the proposed model, a calibration procedure consisting of two stages was engaged: (i) in the first step, special attention was paid to the correct mapping of experimental results into a numerical model; (ii) in the second one, the already tuned numerical model was used to generate much larger sets of cases, which were then utilized to identify the sought parameters in the analytical model.

In the first step, the only unknowns are the initial imperfections. Therefore, a very simple strategy is used, in which the numerical model is calibrated with experimental data by appropriate scaling of the initial deformations of the vertical box walls. In the second step, the coefficients in Equations (9) and (10) are identified in the assumed order: first the constant *k* as well as exponents *r* and *r*1, then *r*<sup>2</sup> and *r*3. In both cases, simple techniques were used to minimize the discrepancy between analytical model prediction and numerical results with the use of the least squares method.

#### **3. Results**

#### *3.1. Corrugated Board: Material Testing*

In order to correctly determine the properties of the material, it was necessary to examine samples of corrugated board in several typical laboratory tests. For this purpose, a FEMat BSE device (FEMat Sp z o.o., Poznan, Poland) [50] was used. In total, seven different types of corrugated cardboard with a grammage of 350 to 965 g/m<sup>2</sup> were tested. Since cardboard is a very heterogeneous material, at least 10 samples in each test were examined for each grade in order to obtain statistically reliable results. In Table 1, the sample results for the BC-780 grade are summarized. The first column represents a test number, the second column shows the sample thickness and in the third to ninth columns the results obtained from different tests in both orthotropy directions are demonstrated (all test symbols are explained in the previous section).

Figure 9 demonstrates the force-displacement curves from all tests of the corrugated cardboard. Since both shape of the curve and the calculated shear stiffness (SST) in the machine and cross directions are almost identical, only the values in the MD are shown. In Table 2, the mean values of the tests for all seven grades are presented. The first column represents grades that were used in the packaging for which the box compression test was carried out (details will be discussed in the next section). In the second to ninth columns, the measured stiffnesses obtained from the BSE device are shown.


**Table 1.** Test values for BC-780 corrugated cardboard grade.

**Figure 9.** Force-displacement curves for BC-780 corrugated cardboard in various tests: (**a**) ECT; (**b**) BNT–MD; (**c**) BNT–CD; (**d**) SST; (**e**) TST–MD; (**f**) TST–CD.



*3.2. Box Compression Test (BCT)*

In the next step, the load capacity of the packaging was checked. For this purpose, the FEMat BCT-20T20 compact press (FEMat Sp. Z o.o., Poznan, Poland) [38] was exploited (see Figure 3a). A total number of 18 samples of various dimensions and materials were prepared. The analysis was carried out for two types of packaging: without and with an offset. In Table 3, the results obtained with the box compression test are presented. In the first column, corrugated cardboard grades are shown. The second, third and fourth columns show the dimensions of the package (see Figure 7). For offset packaging, the edge of the *B*<sup>1</sup> dimension is the offset edge. The fifth column represents the value of the load capacity of the package without offset. Columns six and seven are the BCT values for the offset package: the sixth column is the value of the first extreme and the seventh column is the value of the second extreme. ଵ


**Table 3.** Main dimensions and BCT values of various corrugated cardboard packaging. 

In Figure 10, the force-displacement diagrams for boxes with dimensions 300 × 200 × 200 mm, with and without offset, made of BC-780 and EB-965 corrugated cardboard are shown.

**Figure 10.** Selected measurements from a BCT press for grades: (**a**) BC-780; (**b**) EB-965.

#### *3.3. Prediction Results of the Numerical Model*

Having the geometry of all the tested boxes and the material properties of the corrugated cardboards, it was possible to build numerical models and calibrate the only one remaining component: the initial imperfections. These are especially important in the geometrically nonlinear FE analysis. To introduce preliminary deformations into the model, first a buckling analysis was carried out to find the first preferred buckling mode, which was then introduced as a deformed shape of the load-bearing panels of the box.

The influence of the imperfection size on the load capacity of the box 300 × 200 × 200 mm made of BC-780 is shown in Figure 11.

**Figure 11.** Influence of the imperfections on the load capacity of the BC-780 box.

After a successful calibration procedure, the results obtained with the numerical model are summarized in Table 4, which also shows the differences between the calculated and the measured values of the BCT.

**Table 4.** Comparison of measured and numerically determined BCT values for various corrugated cardboard packaging.


#### *3.4. Prediction Results of the Analytical Model*

As already discussed, the main step was to calibrate the coefficients in the analytical formulas for the load capacity estimation of corrugated board packaging. For this purpose, synthetically generated results were utilized. Thanks to the use of numerical results, the range of packaging dimensions was much wider, which resulted in a greater number of analyzed cases and therefore made the calibration more reliable.

ሾ<sup>ଷ</sup> , <sup>ସ</sup> ሿ <sup>ଷ</sup> <sup>ସ</sup> Table 5 shows all coefficients found in the minimization process used inEquations (9) and (10), while in Figure 12 the discrepancy function in two-dimensional space [*r*3,*r*4] is shown. It can be seen that in the selected range of parameters *r*<sup>3</sup> and *r*<sup>4</sup> there is only one local minimum, which is also the global minimum (see Figure 12).

**Table 5.** Coefficients values.


− − −

− − −

<sup>ଷ</sup> <sup>ସ</sup>

− −

− −

Figure 13 shows the estimation errors obtained from the analytical model in the calibration procedure for all offset and non-offset boxes. Table 6 presents a comparison of the results obtained from the tuned analytical model with the experimental results.

Figures 14 and 15 show the distribution of the prediction error in the design space, which are the main dimensions of the box (L, B, H). It can be seen that the greatest error occurs with boxes that are short and long.

**Figure 14.** The prediction error distribution obtained using the analytical model for the first extreme: (**a**) H = 200 mm; (**b**) H = 500 mm; (**c**) H = 800 mm.

**Figure 15.** The prediction error distribution obtained using the analytical model for the second extreme: (**a**) H = 200 mm; (**b**) H = 500 mm; (**c**) H = 800 mm.

#### **4. Discussion**

Since corrugated board is an orthotropic and non-homogeneous material, a large number of tests were required for the correct characterization of its mechanical parameters. This means that when testing both corrugated cardboard and boxes made of such material, one can expect a large dispersion of test results. This is related to the heterogeneity of the paper itself, as well as the corrugated cardboard, and the inaccuracy of the assemblies of the tested packaging. Thus, the number of tests in the case of boxes should not be less than five, as is the case with testing the corrugated cardboard samples.

Among the many mechanical tests of corrugated board available in the papermaking laboratory, the most important define not only the static edge crush resistance but also the flexural and torsional stiffness of the specimen. In this study, the BSE system [50], which allows the examination of five physical parameters of cardboard (for three of them in both directions of orthotropy), was exploited. Based on the results of all laboratory tests (see Table 1 and Figure 9), homogenized parameters describing the elastic and plastic behavior of the particular corrugated cardboard were obtained.

In order to diversify the set of BCT laboratory results, various corrugated cardboards (a total of seven types) that are used for box production and nine different dimensions of the packaging structure, in two variants (without and with offsets), were tested. The results are presented in Table 3, where among the dimensions of the boxes and the symbols of the corrugated board one can also find the BCT results for two cases: (a) without offsets in column five and (b) with offsets in columns six and seven. The sixth column of Table 3 presents the force value for which two offset walls have been crushed. This is clearly seen in Figure 10: the first peak in the blue result plots. Column seven of Table 3 shows the maximum force value obtained in the BCT test.

Both predictive models take into account the behavior shown in Figure 10, which is characteristic for the offset boxes. The numerical model is loaded sequentially, first on the walls with an offset and then when the displacement of the upper surface exceeds the given offset, the two remaining walls are also loaded. As already mentioned, after calibrating the material model and for the given geometry, the only unknown was the size of the imperfections of the vertical walls. These parameters were in each case adjusted so that the estimates agreed with the laboratory results. The effect of the applied imperfection in a specific case (box EB-780) is shown in Figure 11.

This phenomenon is treated slightly differently in the analytical model, which is based solely on the geometry of the box, its strength in CD and both stiffnesses in MD and CD. In this case, the imperfections are embedded in the predictive model through the critical load term in Equations (9) and (10), while the sequential crushing of the shifted and non-shifted faces is captured by independently determining two values and scaling the maximum force in the first peak by the factor α (see Equation (8)). This allows the degraded resistance of walls with an offset and the resistance of walls without an offset to be taken into account in the second peak. The tuning exponents found by the minimization procedure (shown in Table 5) reached the optimal values of 0.5 or 1.0, while the constant *k* and exponent *r* reached values of 0.75 and 0.55, respectively.

The use of data synthetically generated by the calibrated numerical model allowed a much greater accuracy of the tuned parameters in the analytical model to be obtained. This was mainly due to a larger range of results numerically generated for various geometric dimensions of boxes that could not be physically produced and tested in the BCT press. Figure 13 shows the distribution of the prediction error of the load capacity of various corrugated board packages without and with an offset. The largest discrepancies occur for packages with a relatively large proportion of dimensions (see Figures 14 and 15). However, the average error in both cases does not exceed 7%. Overall, the proposed predictive analytical model can capture the first peak in any experimentally tested sample fairly correctly, and the error in most cases is less than 7%. The greatest differences can be observed for samples B-400-1 and B-400-2, where the error was 8% and 9%, respectively. Similar conclusions can be drawn when predicting a second peak. In most cases, the error did not exceed 5%; only for the B-400-2 sample did it reach 9%.

In general, the application of analytical models existing in the literature, e.g., those proposed in [4–10] or even more the recent models presented by Garbowski et al. [1,2,11], does not allow one to predict the strength of boxes with shifted flaps. The reason is that these models do not take into account the sequential crush of the package walls. Particular attention should also be paid to modeling with purely numerical models, because special techniques for sequential loading of the walls with appropriate imperfections should be considered as well. The results presented in Tables 4 and 6 show the precision with which both the numerical model and the proposed analytical model reflect the laboratory results for selected constructions of corrugated cardboard boxes. The results obtained from both models do not differ by more than 10% from the experimental results.

#### **5. Conclusions**

This article presents numerical and analytical models for predicting the strength of boxes with displaced flaps. The obtained results are in accordance with the conducted laboratory tests. In both models, the mechanical parameters of the corrugated board obtained from the selected laboratory tests were implemented. Both models are based on a sequential approach for the loading of the vertical walls of a box; the walls with an offset are loaded first, then the walls without an offset. At the moment of loading the walls without offset, the two walls loaded in the first step are already partially damaged. Therefore, this type of packing is characterized by much lower load-bearing capacity than packages with flaps without an offset. Thanks to the methodology presented in this paper and utilization of such predictive tools, it is possible not only to design packaging more consciously, but also to deliver and optimally use the material for their manufacturing, and thus improve the sustainable economy of the production plant.

**Author Contributions:** Conceptualization, T.G.; methodology, T.G.; software, T.G. and D.M.; validation, D.M., A.K.-P. and T.G.; formal analysis, D.M.; investigation, D.M. and A.K.-P.; resources, D.M.; data curation, D.M.; writing—original draft preparation, A.K.-P., D.M. and T.G.; writing—review and editing, A.K.-P. and T.G.; visualization, D.M.; supervision, T.G.; project administration, T.G.; funding acquisition, A.K.-P. and T.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The APC was funded by the Ministry of Science and Higher Education, Poland, the statutory funding at Poznan University of Life Sciences, grant number 506.569.05.00 and the statutory funding at Poznan University of Technology, grant number 0411/SBAD/0004.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Acknowledgments:** Special thanks to the FEMat Sp. z o. o. company (Poznan, Poland) (www. fematsystems.pl—accessed on 21 July 2021) for providing the commercial software.

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

## **References**


**Zbigniew Pozorski 1, \* , Jolanta Pozorska 2 , Ireneusz Kreja <sup>3</sup> and Łukasz Smakosz 3**


**Abstract:** This paper deals with the local loss of stability (wrinkling) problem of a thin facing of a sandwich panel. Classical solutions to the problem of a facing instability resting on a homogeneous and isotropic substructure (a core) are compared. The relations between strain energy components associated with different forms of core deformations are discussed. Next, a new solution for the orthotropic core is presented in detail, which is consistent with the classic solution for the isotropic core. Selected numerical examples confirm the correctness of the analytical formulas. In the last part, parametric analyses are carried out to illustrate the sensitivity of wrinkling stress to a change in the material parameters of the core. These analyses illustrate the possibility of using the equations derived in the article for the variability of Poisson's ratio from −1 to 1 and for material parameters strongly deviating from isotropy.

**Keywords:** sandwich panels; local instability; strain energy; wrinkling; orthotropic core

#### **1. Introduction**

In a typical sandwich element, the two facings are joined to each other by a relatively thick but deformable core. The deformations and stresses in the sandwich panel are caused by the acting loads (wind, snow, self-weight, live load), but they are also largely due to thermal loads. As a result of these interactions, the facing can be compressed, and because it is connected to a susceptible substructure (a core), it very often experiences local loss of stability (wrinkling).

Wrinkling is undoubtedly one of the most common damage mechanisms of a sandwich element. For this reason, the correct estimation of the stress leading to the loss of facing stability is a key issue that has been undertaken by many researchers using various approaches: analytical, numerical, experimental, or mixed (or some combination of these approaches). Numerical methods allow for solving many complex problems, and the performed experiments allow for the verification of the obtained results. Nevertheless, analytical solutions should also be treated as very valuable, even if they are obtained with significant simplifications. Simple formulas are easy for engineering application and allow for a very quick (and continuous) assessment of the sensitivity of the solution to a change in design parameters.

With full awareness of the new challenges related to the sandwich structures (anisotropy [1], influence of extreme excitations [2], new production technologies [3], and many others), this work is an attempt to take a deeper look at the known classical solutions to the local instability problem [4–6]. The presented solution for an orthotropic core is based on the work of [7], in which sandwich columns under compression were considered, and the solution was presented in the form of hyperbolic functions. It also

**Citation:** Pozorski, Z.; Pozorska, J.; Kreja, I.; Smakosz, Ł. On Wrinkling in Sandwich Panels with an Orthotropic Core. *Materials* **2021**, *14*, 5043. https://doi.org/10.3390/ma14175043

Academic Editor: Antonio Gloria

Received: 8 July 2021 Accepted: 1 September 2021 Published: 3 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

clearly refers to the classic solution for an isotropic core [6], where a facing and a core were assumed as infinite and the differential equation written for the facing was used.

The above-mentioned classic approaches to the problem of facing instability are constantly being used and extended to more and more complex issues. The analytical model that leads to wrinkling of the orthotropic face layer supported by a transversely isotropic core was presented in [8]. Wrinkling of a composite-facing sandwich panel under biaxial loading was discussed in [9]. Article [10] presents the solution to the symmetrical face sheet wrinkling problem using the energy method. The approach focused on a 3D case of wrinkling of orthotropic face sheets was presented in [11]. The analytical approach to the problem of anisotropic facing instability was presented in the works [12,13]. Wrinkling in sandwich structures with a functionally graded core was discussed in [14]. The papers [15–17] are examples of work on wrinkling, in which the core was modeled using higher-order theories, which allowed, among others, to take into account the influence of the core transverse compressibility.

This paper is divided into three parts. In the first one, we present some relations between the classical solutions to the analyzed problem of wrinkling. We believe that they will shed a slightly different light on known solutions. This applies to the conditions of reaching the critical stress, the influence of the Poisson ratio on the wrinkling stress, and the relationship between strain energy components. In the second part, the solution for the orthotropic core is derived and discussed, and we focus on the interpretation of the solution and the question of the conditions for obtaining it. In the third part, a parametric analysis of the solution for the orthotropic core is presented, illustrating the sensitivity of the solution (especially the wrinkling stress) to a change in some material parameters. In our opinion, this is essential for the optimal design of layered structures. By assuming certain constraints on material parameters, we can specify a solution with the maximum value or the minimum sensitivity.

#### **2. Formulation of the Problem**

We are considering a sandwich panel consisting of two thin facings and a thick but deformable core. Due to the bending of the composite panel, considerable compressive stresses may be generated in its facing, resulting in a local loss of stability. The instability has the form of wrinkling. In general, due to the variety of support and load conditions, the problem can be very complex; however, in practical civil engineering problems, a facing is usually compressed unidirectionally [18].

The wrinkling phenomenon may be considered as a compression effect of a thin facing (treated as a beam or plate) supported by a continuous elastic core (Figure 1). The facing in tension is ignored because the deformation of the core quickly disappears as the distance from the compressed facing increases. It is convenient to assume that the compressed facing is infinitely long and the core extends to infinity on one side of the facing. The wrinkling is associated with short waves of buckling of the facing. Figure 1 shows a fragment of the deformed facing supported by the core.

It is assumed that the face layer is in a uniform stress and strain state. The deformations of the facing, which are infinite and periodic, induce strain and stress in the core. Core deformations quickly decay as the variable *z* increases, and the rate of this decay depends on the assumed displacement field.

The core and facing materials are homogeneous. Suppose the core is isotropic or orthotropic with one of the orthotropic axes coinciding with the direction of the compression. The facing material could be orthotropic if its axes were aligned with the material axes of the core. These are quite strong assumptions, but they give analytical results that are relatively easy to interpret.

**Figure 1.** Assumed shape of a wrinkling.

#### **3. Classical Solutions of the Wrinkling Problem**

#### *3.1. Energy Method—Linear Decay Function*

Following the proposition of Hoff and Mautner [4], the core is affected only in a small zone with depth *h* (smaller than the thickness of the core). The shape of the face deformation is assumed in the sinusoidal form (Figure 1), and the core deformation field vanishes linearly with coordinate *z*:

$$w\_C = w\_F \frac{(h-z)}{h} = \mathcal{W} \frac{(h-z)}{h} \sin \frac{\pi \chi}{l},\tag{1}$$

 = ி ሺି௭ሻ = ሺି௭ሻ sin గ௫ where *wC*, *wF*, and *W* denote the vertical displacement of the core, face, and the displacement amplitude, respectively. The term *l* is a half wavelength of the wrinkles. Comparing the sum of strain energy of the core and the facing (per half wavelength) with the external work done by an applied work, the expression on the compressive stress in the facing is obtained:

$$
\sigma\_{\rm x} = \frac{E\_{\rm C}l^2}{\pi^2 t\_F h} + \frac{hG\_{\rm C}}{3t\_F} + \frac{\pi^2 E\_F}{12} \left(\frac{t\_F}{l}\right)^2. \tag{2}
$$

Symbols *E<sup>C</sup>* and *G<sup>C</sup>* denote the modulus of elasticity and shear modulus of the isotropic core material, respectively. The thickness of the facing is *tF*, whereas the modulus of elasticity of the isotropic facing material is *EF*.

௫ = ா మ గమ௧ಷ + ீ + <sup>గ</sup> <sup>మ</sup>ா<sup>ಷ</sup> ଵଶ ቀ ௧ಷ ቁ ଶ The minimum value of the compressive stress (2) corresponds to the critical (wrinkling) stress, and it can be found by using derivatives of *σ<sup>x</sup>* with respect to *h* and *l*:

ଷ௧ಷ

$$
\sigma\_w = \sqrt[3]{\frac{3}{4}} \cdot \sqrt[3]{E\_\mathbb{C} \mathbf{G}\_\mathbb{C} E\_F} \cong 0.909 \cdot \sqrt[3]{E\_\mathbb{C} \mathbf{G}\_\mathbb{C} E\_F}.\tag{3}
$$

*σ* It is worth noting that reaching the wrinkling stress corresponds to a situation in which each term on the right-hand side of Equation (2) is equal to each other.

#### ଷ య య *3.2. Energy Method—Exponential Decay Function*

<sup>௪</sup> = ට ସ ∙ ඥி ≅ 0.909 ∙ ඥி య Plantema [5] assumed the displacement field of the exponential form

$$w\_{\mathbb{C}} = w\_{\mathbb{F}}e^{-kz} = \mathcal{W}e^{-kz} \sin \frac{\pi \mathfrak{x}}{l} \,\tag{4}$$

where *k* ≥ 0 is an auxiliary constant (with the unit inverse to the unit of variable *z*). The use of the strain energy of the core makes it possible to represent the compressive stress of the facing as

$$
\sigma\_{\rm x} = \frac{E\_{\rm C}kl^2}{2\pi^2 t\_F} + \frac{G\_{\rm C}}{2kt\_F} + \frac{\pi^2 E\_F}{12} \left(\frac{t\_F}{l}\right)^2. \tag{5}
$$

The wrinkling stress is obtained from the conditions of zeroing the derivatives of *σ<sup>x</sup>* with respect to *k* and *l*. As we can see, a slightly different assumption of the displacement field leads to a different result. First of all, the condition for reaching the extreme (minimum) stress *σ<sup>x</sup>* is different. Again, when the critical stress is reached, each term of Equation (5) has the same value. The wrinkling stresses (3) and (6) are independent of the Poisson ratio of the core material.

$$
\sigma\_w = \frac{3}{2 \cdot \sqrt[3]{6}} \cdot \sqrt[3]{E\_\mathrm{C} G\_\mathrm{C} E\_F} \cong 0.825 \cdot \sqrt[3]{E\_\mathrm{C} G\_\mathrm{C} E\_F} \tag{6}
$$

#### *3.3. Differential Equation Method*

The solution based on the differential equation method was presented by Allen [6]. Stresses in the elastic isotropic medium can be defined using the Airy stress function *F* (*x*,*z*). The strain compatibility in the *x*–*z* plane leads to the bi-harmonic differential equation.

$$\frac{\partial^4 F}{\partial z^4} + 2\frac{\partial^4 F}{\partial x^2 \partial z^2} + \frac{\partial^4 F}{\partial x^4} = 0. \tag{7}$$

Equation (7) is satisfied by the function

$$F(x,z) = A \sin \frac{\pi x}{l} (1 - Bz)e^{-\frac{\pi z}{l}},\tag{8}$$

where *A* and *B* are constants. Constant *B* can be found by using the condition that the *x*-displacements and strains at the surface of the core (*z* = 0) are equal to zero. Constant *A* can be expressed by the amplitude *W* of the *z*-displacement at *z* = 0. By using Allen's method, nearly the entire mechanical field is obtained, which depends on *x* and *z* variables. If the state of plane stress is assumed, then displacements *u*, *v*, and *w*, strains *εx*, *εy*, *εz*, and *γxz*, and stresses *σx*, *σy*, and *τxz* are non-zero.

The equilibrium differential equation for the facing has the form

$$B\_F \frac{d^4 w}{dx^4} + P \frac{d^2 w}{dx^2} = \sigma\_z. \tag{9}$$

where the stress *σ<sup>z</sup>* is the effect of the interaction between the facing and the core (see Figure 1). The symbol *B<sup>F</sup>* denotes the face bending stiffness per unit width. For the beam theory (as used here), *B<sup>F</sup>* = *EFt* 3 *F* /12; in the case of the plate theory, *B<sup>F</sup>* = *EFt* 3 *F* /12 1 − *ν* 2 *F* . Using the function of the facing displacement

$$w = \mathcal{W}\sin\frac{\pi\mathbf{x}}{l},\tag{10}$$

and the parameter *m* = *l*/*tF*, the compressive stress in the facing can be expressed as

$$
\sigma\_{\chi} = \frac{\pi^2 E\_F}{12m^2} + \frac{a}{\pi}m = \sigma\_1 + \sigma\_2. \tag{11}
$$

where

$$a = \frac{2E\_{\mathbb{C}}}{(1 + \nu\_{\mathbb{C}}) \cdot (3 - \nu\_{\mathbb{C}})} \tag{12}$$

is the material constant. The two terms of the solution for (11) are denoted as *σ*<sup>1</sup> and *σ*2, respectively.

From the condition for the extreme, d*σx*/d*m* = 0, we can find *m* = *π* · <sup>√</sup><sup>3</sup> *<sup>E</sup>F*/6*<sup>a</sup>* and the minimum critical (wrinkling) stress:

$$
\sigma\_{\mathcal{W}} = \sqrt[3]{\frac{9}{2(1+\nu\_{\mathcal{C}}) \cdot (3-\nu\_{\mathcal{C}})^2}} \cdot \sqrt[3]{E\_{\mathcal{C}} \mathcal{G}\_{\mathcal{C}} E\_F} = r \cdot \sqrt[3]{E\_{\mathcal{C}} \mathcal{G}\_{\mathcal{C}} E\_F}.\tag{13}
$$

If we assume the facing stiffness as for the plate, *B<sup>F</sup>* = *EFt* 3 *F* /12 1 − *ν* 2 *F* , the modulus *E<sup>F</sup>* should be replaced by *EF*/ 1 − *ν* 2 *F* . ௪ = ට ଶሺଵାఔሻ∙ሺଷ ି ఔሻ మ ∙ ඥி = ∙ ඥி ி = ிி ଷ /12ሺ1 − ி ଶ ሻ

య

ଽ

య

It is interesting that for the minimum value of *σ<sup>x</sup>* (11), the second expression (*σ*2) is exactly two times higher than the first (*σ*1) [19]. From some literature sources, e.g., [6] p. 159, Figure 8.3, it can be drawn incorrectly that both of these values are equal. The value of the first root (*r*) depends only on the Poisson ratio of the core material *νC*, but for the typical range of this parameter, the root *r* reaches the value from 0.780 to 0.794. It is also worth noting that as the Poisson ratio tends to −1, the critical stress would increase to infinity, although this is a rather theoretical situation. ி/ሺ1 − ி ଶ ሻ *σ σ σ ν* −

#### *3.4. Comparison of Classical Solutions*

#### 3.4.1. Influence of the Poisson Ratio

Let us return first to Allen's solution. The result (13) was obtained for a plane stress state. Assuming a plane strain state, the procedure is analogous; however, the functions of stresses, strains, and displacements are different. Equation (11) is valid, but:

$$a = \frac{2E\_{\mathbb{C}}(1 - \nu\_{\mathbb{C}})}{(1 + \nu\_{\mathbb{C}}) \cdot (3 - 4\nu\_{\mathbb{C}})},\tag{14}$$

*σ* =∙ ඥி/6 <sup>య</sup>

య

*σ σ*

$$
\sigma\_{\mathcal{W}} = \sqrt[3]{\frac{9\left(1-\nu\_{\rm C}\right)^{2}}{2\left(1+\mathsf{V}\_{\rm C}\right)\cdot\left(3-4\nu\_{\rm C}\right)^{2}}} \cdot \sqrt[3]{\mathcal{E}\_{\rm C}\mathcal{G}\_{\rm C}\mathcal{E}\_{\rm F}} = s \cdot \sqrt[3]{\mathcal{E}\_{\rm C}\mathcal{G}\_{\rm C}\mathcal{E}\_{\rm F}}.\tag{15}
$$

Of course, the value of *s* in (15) is different than *r* in (13). To compare Allen's solutions in the case of the plane stress and plane strain states, see Figure 2.

**Figure 2.** Comparison of Allen's solution in the case of the plane stress (*r*) and plane strain (*s*) states.

− It should come as no surprise that for *v<sup>C</sup>* = 0, the coefficients *r* and *s* are identical and equal to 0.794. For negative values of *vC*, the coefficients *r* and *s* take similar values that are much higher than 0.794. For *v<sup>C</sup>* tending to −1, the values of *r* and *s*, and hence the critical stress values, tend to infinity. In the range of *v<sup>C</sup>* (−1; +0.5), the critical stresses in the plane strain state are higher than in the plane stress state, but the greatest differences between *r* and *s* appear for *v<sup>C</sup>* close to 0.5. This is obvious because in a plane state of stress, the material has the potential to deform in the *y*-direction (perpendicular to the plane), which facilitates the deformation of the facing. The plane strain condition limits the deformation (in the *y*-direction) and makes it difficult to buckle the facing. The greater the Poisson ratio, the greater the significance of this effect. For some order, let us remind you that the Poisson

ratio of the core material does not affect the critical stresses in the case of the solutions given by Hoff and Mautner (3) and Plantema (6).

#### 3.4.2. Assumptions and Strain Energy Considerations

The Hoff–Mautner and Plantema solutions are based on an energy approach. The assumption of a specific displacement field turns out to be very effective and quickly leads to a solution. However, it is worth noting that in contrast to Allen's solution, the assumed displacement fields ((1) or (4)) result in non-fulfillment of most of the differential equilibrium equations of a solid (mass forces were omitted in Equation (16)):

$$
\sigma\_{\text{ji},\text{j}} = 0.\tag{16}
$$

Let us return to the solution presented by Allen [6]. In the case of a plane stress state, constant *B* is

$$B = -\frac{\pi}{2l}(1 + \nu\_{\mathbb{C}}),\tag{17}$$

and the stresses in the core are expressed as the corresponding derivatives of the function *F* (*x*,*z*). By using commonly known physical and geometric relationships, we determine the fields of strain and displacement. Therefore, we can calculate the appropriate components of the strain energy of the core, obtaining, respectively:

$$\frac{1}{2} \int\_{0}^{\infty} \int\_{0}^{l} \sigma\_{\text{X}} \varepsilon\_{\text{X}} d\mathbf{x} dz = \frac{1}{16} \frac{A^2}{E\_{\mathbb{C}}} \frac{\pi^3}{l^2} (1 + \nu\_{\mathbb{C}}) \left(1 - \nu\_{\mathbb{C}}^2\right) \tag{18}$$

$$\frac{1}{2} \int\_{0}^{\infty} \int\_{0}^{l} \sigma\_{z} \varepsilon\_{z} d\mathbf{x} dz = \frac{1}{16} \frac{A^{2}}{E\_{\mathbb{C}}} \frac{\pi^{3}}{l^{2}} (1 + \nu\_{\mathbb{C}}) \left(13 - 4\nu\_{\mathbb{C}} - \nu\_{\mathbb{C}}^{2}\right) \tag{19}$$

$$\frac{1}{2} \int\_{0}^{\infty} \int\_{0}^{l} \pi\_{\text{xz}} \gamma\_{\text{xz}} d\mathbf{x} d\mathbf{z} = \frac{1}{16} \frac{A^2}{E\_{\mathbb{C}}} \frac{\pi^3}{l^2} (1 + \nu\_{\mathbb{C}}) \left(10 - 4\nu\_{\mathbb{C}} + 2\nu\_{\mathbb{C}}^2\right). \tag{20}$$

For *v<sup>C</sup>* = 0, the ratio of energies expressed in (18)–(20) is 1:13:10. Let us recall that under the condition of loss of stability, the elastic energy in the facing is half of the elastic energy in the core. Commenting on the relations between the energies in the core, we can say that the share of energy (18) resulting from the deformation of the core along the *x*-direction (*εx*) is small, which can justify the omission of this term in classical energy methods. For the sake of order, we note that the fulfillment of the condition of loss of local stability for each of the previously discussed classical energy methods means that the energy components on the left side of Equations (19) and (20) are equal to each other, and the integral (18) is equal to zero.

An additional point requires clarification. Each of the presented classical methods differs in the final result, but only because of different assumptions and not because of the method itself. For example, this is easily demonstrated by using Allen's mechanical fields for the energy approach. Then, it turns out that the obtained expression for the critical stress is identical to (13).

For an illustration of the assumptions made in Allen's solution, Figure 3 presents the displacement fields *w*(*x*, *z*) (perpendicular to the facing i.e., along the *z*-axis) and *u*(*x*, *z*) (along the *x*-axis). The values are given assuming the constant *A* = 1, see (8). The isotropic material of the core *E<sup>C</sup>* = 4 MPa, *v<sup>C</sup>* = 0.05 and a facing with a thickness of *t<sup>F</sup>* = 0.5 mm made of an isotropic material *E<sup>F</sup>* = 210 GPa were assumed. The range of the *x*-axis corresponds to 2*l* = 74 mm, while the *z*-coordinates are given in millimeters. In Figure 3a, for *z* = 0, we can see a full sinusoid with an extreme equal to 3.248 <sup>×</sup> <sup>10</sup>−<sup>5</sup> , which disappears with increasing *z*. The amplitude of the sinusoid decreases 10 times for *z* = 36 mm. In Figure 3b, according to the assumption, for *z* = 0, the horizontal displacements are equal to zero. The variability of the function *u*(*x*, *z*) in the *x*-direction is described by the cosine function, the extreme value of 4.253 <sup>×</sup> <sup>10</sup>−<sup>6</sup> is reached for *z* = 12 mm; at a distance of *z* = 60 mm, the function value is 10 times smaller than the extreme. The rapid disappearance of displacements

with the increment of the *z*-coordinate, and the values of *w*(*x*, *z*) being one order greater than *u*(*x*, *z*), are both noteworthy, as they, among other things, justify the omission of longitudinal deformations in classical energy methods.

−

−

**Figure 3.** Solution obtained by using the differential equation method: (**a**) vertical displacement *w*(*x*, *z*), (**b**) horizontal displacement *u*(*x*, *z*).

#### **4. Solution for the Orthotropic Core**

#### *4.1. Differential Equation*

A certain solution to the problem of facing wrinkling resting on an orthotropic elastic substructure and loaded on the edge (in the facing plane) was presented in [7]. A similar approach was used in [8]. The following is a detailed solution, which is an extension of [6], formally based on [7,8], but it differs in some nuances. Efforts were made to present the solution precisely in order to also discuss the conditions for obtaining this solution.

Suppose we have an orthotropic core, in which orthotropic axes coincide with the axes of the element. The facing is compressed uniaxially, and the load direction is according to the material axes of the core. Such a situation is very common in practice [20]. The constitutive relation for the orthotropic core material is:


In the case of a 2D problem, relation (21) can be simplified to:

$$\begin{cases} \varepsilon\_{\mathbf{x}} = a\_{\mathbf{x}\mathbf{x}} \sigma\_{\mathbf{x}} - a\_{\mathbf{x}\mathbf{z}} \sigma\_{\mathbf{z}}\\ \varepsilon\_{\mathbf{z}} = -a\_{\mathbf{x}\mathbf{z}} \sigma\_{\mathbf{x}} + a\_{\mathbf{z}\mathbf{z}} \sigma\_{\mathbf{z}}\\ \varepsilon\_{\mathbf{x}\mathbf{z}} = (1/2 \mathbf{G}\_{\mathbf{x}\mathbf{z}}) \tau\_{\mathbf{x}\mathbf{z}} \end{cases} \right\}. \tag{22}$$

In the case of plane stress state, material constants *axx*, *azz*, and *axz* are:

$$\begin{aligned} a\_{xx} &= 1/E\_x\\ a\_{zz} &= 1/E\_z\\ a\_{xz} &= \nu\_{xz}/E\_x \end{aligned} \qquad\qquad\qquad\qquad\qquad\tag{23}$$

whereas for the plane strain state, we have:

$$\begin{cases} a\_{xx} = \frac{1 - \nu\_{xy}\nu\_{yx}}{E\_x} \\ a\_{zz} = \frac{1 - \nu\_{yz}\nu\_{zy}}{E\_z} \\ a\_{xz} = \frac{\nu\_{xz} + \nu\_{xy}\nu\_{yz}}{E\_x} \end{cases} \tag{24}$$

The compatibility of strains in the *x–z* plane requires:

$$
\frac{
\partial^2 \varepsilon\_x
}{
\partial z^2
} + \frac{
\partial^2 \varepsilon\_z
}{
\partial x^2
} - 2 \frac{
\partial^2 \varepsilon\_{xz}
}{
\partial x \partial z
} = 0.
\tag{25}
$$

After introducing the Airy stress function *F* (*x*, *y*) such that

$$
\sigma\_{\mathbf{x}} = \frac{\partial^2 F}{\partial \mathbf{z}^2}, \sigma\_{\mathbf{z}} = \frac{\partial^2 F}{\partial \mathbf{x}^2}, \tau\_{\mathbf{x}\mathbf{z}} = -\frac{\partial^2 F}{\partial \mathbf{x}\partial \mathbf{z}'} \tag{26}
$$

condition (25) takes the following form:

$$a\_{zz}\frac{\partial^4 F}{\partial x^4} + 2\left(\frac{1}{2G\_{xz}} - a\_{xz}\right)\frac{\partial^4 F}{\partial x^2 \partial z^2} + a\_{xx}\frac{\partial^4 F}{\partial z^4} = 0. \tag{27}$$

By using substitution

$$\eta = \epsilon z = \left(\sqrt[4]{a\_{zz}/a\_{xx}}\right) z \tag{28}$$

we obtain

 

$$\frac{\partial^4 F(\mathbf{x}, \boldsymbol{\eta})}{\partial \mathbf{x}^4} + 2\kappa \frac{\partial^4 F(\mathbf{x}, \boldsymbol{\eta})}{\partial \mathbf{x}^2 \partial \boldsymbol{\eta}^2} + \frac{\partial^4 F(\mathbf{x}, \boldsymbol{\eta})}{\partial \boldsymbol{\eta}^4} = 0,\tag{29}$$

where *κ* is a dimensionless quantity and depends only on the material parameters of the core,

$$\kappa = \frac{1}{\sqrt{a\_{\text{xx}} a\_{\text{zz}}}} \left( \frac{1}{2G\_{\text{xz}}} - a\_{\text{xz}} \right). \tag{30}$$

For an isotropic material, *κ* = 1.

#### *4.2. Solution of the Differential Equation*

To find a solution of (29), we separate variables:

$$F(\mathbf{x}, \eta) = G(\mathbf{x})H(\eta) \tag{31}$$

and assume the sinusoidal form of function *G* (*A*<sup>1</sup> is a constant)

$$G(\mathbf{x}) = A\_1 \sin \frac{\pi \mathbf{x}}{l} \tag{32}$$

which leads to

$$\frac{d^4H}{d\kappa^4} - 2\kappa \left(\frac{\pi}{l}\right)^2 \frac{d^2H}{d\eta^2} + \left(\frac{\pi}{l}\right)^4 H = 0. \tag{33}$$

By assuming that the function *H*(*η*) = *e λη* is a general solution of Equation (33), we obtain a solution in the form of a linear combination of this function:

$$H(\eta) = \mathbb{C}\_1 e^{\lambda\_1 \eta} + \mathbb{C}\_2 e^{\lambda\_2 \eta} + \mathbb{C}\_3 e^{\lambda\_3 \eta} + \mathbb{C}\_4 e^{\lambda\_4 \eta},\tag{34}$$

where

$$\begin{aligned} \lambda\_1 &= +\frac{\pi}{l}\sqrt{\kappa - \sqrt{\kappa^2 - 1}} \\ \lambda\_3 &= +\frac{\pi}{l}\sqrt{\kappa + \sqrt{\kappa^2 - 1}} \end{aligned} \quad \begin{aligned} \lambda\_2 &= -\frac{\pi}{l}\sqrt{\kappa - \sqrt{\kappa^2 - 1}} \\ \lambda\_4 &= -\frac{\pi}{l}\sqrt{\kappa + \sqrt{\kappa^2 - 1}} \end{aligned} \right\}. \tag{35}$$

The positive solutions for *λ* have to disappear to allow an exponential decrease in the stresses in the thickness direction *z*. Therefore, *C*<sup>1</sup> = 0, *C*<sup>3</sup> = 0, and

$$F(\mathbf{x}, \mathbf{y}) = \begin{bmatrix} \mathbf{C}\_2 e^{-\frac{\pi}{\mathsf{T}} \sqrt{\mathbf{x} - \sqrt{\mathbf{x}^2 - 1} \mathsf{c} \mathsf{c}}} + \mathbf{C}\_4 e^{-\frac{\pi}{\mathsf{T}} \sqrt{\mathbf{x} + \sqrt{\mathbf{x}^2 - 1} \mathsf{c} \mathsf{c}}} \\\\ \mathbf{B}\_1 e^{-\frac{\pi}{\mathsf{T}} \sqrt{\mathbf{x} - \sqrt{\mathbf{x}^2 - 1} \mathsf{c} \mathsf{c}}} + \mathbf{B}\_2 e^{-\frac{\pi}{\mathsf{T}} \sqrt{\mathbf{x} + \sqrt{\mathbf{x}^2 - 1} \mathsf{c} \mathsf{c}}} \end{bmatrix} \mathbf{A}\_1 \sin\frac{\pi \mathbf{x}}{\mathsf{T}} \end{bmatrix} \tag{36}$$

where *B*<sup>1</sup> = *C*<sup>2</sup> *A*<sup>1</sup> and *B*<sup>2</sup> = *C*<sup>4</sup> *A*<sup>1</sup> are constants. These constants can be calculated with the following boundary conditions:

$$
\varepsilon\_x(z=0) = 0,\tag{37}
$$

which reflects the observation that the face material is typically much stiffer than the core material, and

$$
\sigma\_z(z=0) = A \sin \frac{\pi \chi}{l},\tag{38}
$$

because the stress at the interface in the *z*-direction is distributed as a sine wave with a certain amplitude *A* corresponding to the assumed wave deformation. From the assumed boundary conditions, we obtain:

$$\begin{cases} \begin{array}{c} B\_1 = -A \left( \frac{l}{\pi} \right)^2 \frac{a\_{xx} + a\_{xx} \epsilon^2 \left( \kappa + \sqrt{\kappa^2 - 1} \right)}{2a\_{xx} \epsilon^2 \sqrt{\kappa^2 - 1}}\\ \end{array} \\\ B\_2 = A \left( \frac{l}{\pi} \right)^2 \frac{a\_{xx} + a\_{xx} \epsilon^2 \left( \kappa - \sqrt{\kappa^2 - 1} \right)}{2a\_{xx} \epsilon^2 \sqrt{\kappa^2 - 1}} \end{array} . \tag{39}$$

It is easy to note that *B*<sup>1</sup> + *B*<sup>2</sup> = −*A l π* 2

The equilibrium differential equation for the facing has the same form as (9). By integrating *ε<sup>z</sup>* (22), we can find the following expression for the facing displacement *w<sup>F</sup>* = *w*(*z* = 0)

.

$$w(z=0) = \frac{\pi}{l} \sin \frac{\pi \mathbf{x}}{l} \left[ a\_{zx} B\_1 \varepsilon \sqrt{\kappa - \sqrt{\kappa^2 - 1}} + a\_{zx} B\_2 \varepsilon \sqrt{\kappa + \sqrt{\kappa^2 - 1}} + a\_{zz} B\_1 \frac{1}{\varepsilon \sqrt{\kappa - \sqrt{\kappa^2 - 1}}} + a\_{zz} B\_2 \frac{1}{\varepsilon \sqrt{\kappa + \sqrt{\kappa^2 - 1}}} \right]. \tag{40}$$

By performing some additional algebraic transformations as suggested by Allen [6], one can arrive at the analogy of (11) with

$$a = \frac{\sqrt{\kappa^2 - 1}}{\sqrt{\kappa + \sqrt{\kappa^2 - 1}} - \sqrt{\kappa - \sqrt{\kappa^2 - 1}}} \frac{2a\_{\text{xx}}\mathfrak{c}}{2a\_{\text{xx}}a\_{\text{xz}}\mathfrak{c}^2 - a\_{\text{xz}}^2 + a\_{\text{xx}}a\_{\text{zz}}(2\kappa + 1)}.\tag{41}$$

From the condition for the extreme d*σx*/d*m* = 0, we have as before

$$m = \pi \cdot \sqrt[3]{\frac{E\_F}{6a}}\tag{42}$$

and the minimum critical (wrinkling) stress is obtained (the solution is consistent with (13)): *κ*

$$
\sigma\_w = \frac{3}{2\sqrt[3]{6}} \cdot \sqrt[3]{a^2 E\_F}.\tag{43}
$$

It is easy to prove that again, *σ*<sup>2</sup> = 2*σ*<sup>1</sup> (see also Figure 4).

*σ*

**Figure 4.** Wrinkling stress as a function of *m* parameter.

Based on the quick analysis of Equation (41), it can be concluded that for *a*, which has the nature (and the measurement unit) of the stiffness modulus (of the core), the condition *κ* > 1 must be satisfied to reach real values. With smaller values of *κ*, the roots in Equation (41) are complex numbers. However, it is somewhat surprising that despite the complex roots in (41), the value of *a*,

$$a = \sqrt{\frac{\kappa + 1}{2}} \frac{2a\_{\text{xx}}\epsilon}{2a\_{\text{xx}}a\_{\text{xz}}\epsilon^2 - a\_{\text{xz}}^2 + a\_{\text{xx}}a\_{\text{zz}}(2\kappa + 1)},\tag{44}$$

is real if the condition *κ* > −1 is satisfied. In order for parameter *m* to be positive, the denominator in expression (41) must be positive (the nominator is positive). When the denominator in (41) approaches 0<sup>+</sup> , *a* and consequently also *σ<sup>w</sup>* tend to infinity.

Let us take a moment to analyze the value of *κ* in a plane stress state. According to (30), we have

$$\kappa = \sqrt{E\_{\rm x} E\_{\rm z}} \left( \frac{1}{2G\_{\rm xz}} - \frac{\nu\_{\rm xz}}{E\_{\rm x}} \right). \tag{45}$$

Modules *Ex*, *Ez*, and *Gxz* must be positive. In this situation, if *νxz* is negative, then *κ* will always be positive. If *νxz* = 0.5, then *κ* is positive when *E<sup>x</sup>* > *Gxz*; if *νxz* = 1, then *κ* is positive when *E<sup>x</sup>* > 2*Gxz*. Let us recall that in the case of orthotropic materials, the condition

for the stability of the material behavior is not only the positive values of the *Ex*, *Ez*, and *Gxz* modules, but also, among others [21],

$$|\nu\_{\rm xz}| < \sqrt{E\_{\rm x}/E\_{\rm z}}.\tag{46}$$

#### **5. Examples**

#### *5.1. Analytical Solutions*

The first example concerns the facing with a thickness of *t<sup>F</sup>* = 0.5 mm made of an isotropic material (*E<sup>F</sup>* = 210 GPa) placed on an isotropic core (*E<sup>C</sup>* = 4 MPa, *ν<sup>C</sup>* = 0.05; therefore *G<sup>C</sup>* = *EC*/2(1 + *ν<sup>C</sup>* ) = 1.905 MPa). According to the approach of Hoff–Mautner (3), Plantema (6), and Allen (13), we will obtain the following wrinkling stresses, respectively: *σ H-M* = 106.32 MPa, *σ <sup>P</sup>* = 96.49 MPa, and *σ <sup>A</sup>* = 92.36 MPa.

Now let us consider the same facing (*t<sup>F</sup>* = 0.5 mm *E<sup>F</sup>* = 210 GPa) supported by an orthotropic substructure (*E<sup>x</sup>* = 10 MPa, *E<sup>z</sup>* = 4 MPa, *νxz* = 0.05, *Gxz* = 3 MPa).

We assume a plane state of stress in the *x–z* plane and look for the critical stress that will cause the wrinkling of the facing. According to (28), (30), and (41), we will get *ǫ* = 1.257, *κ* = 1.022, and *a* = 3.256 MPa. The wrinkling stress is achieved for *m* = 69.336 (see (42)) and according to (43), *σ<sup>w</sup>* = 107.8 MPa. Figure 4 shows the dependence of the critical stress (solid line) on the *m* parameter. The blue and brown lines show both stress components (11).

#### *5.2. Numerical Solutions*

Numerical analysis of the instability problems of all kinds of structures is an intriguing and fascinating task, but it is not easy. First of all, it should be realized that numerical models are often much more complex than analytical models. This is due to the fact that commercial software (using, for example, the finite element method) allows for a relatively quick creation of spatial models. However, the problem is that the appropriate model class requires boundary conditions corresponding to this model. Therefore, these conditions are usually different than in the analytical model, which makes it difficult to compare the solutions. This issue was pointed out by numerous researchers [22–24]. This problem also arises when it comes to determining the critical stresses in a thin facing resting on a susceptible substructure.

The numerical analysis of the discussed issue was prepared using ABAQUS, which is a software suite for finite element analysis and computer-aided engineering. The problem was solved using two different classes of numerical models: 2D and 3D. A detailed description of the 3D model is presented below. The results obtained for the 2D model are presented at the end of the subsection.

The three-dimensional model was created in order to fully analyze the phenomenon of loss of stability in conditions close to the plane stress state. Of course, we also tried to make the numerical model as close as possible to the analytical model. The model space is not infinite, but the dimensions have been defined so that the displacements, strains, and stresses at the edge of the model are relatively small; the core body was 1.2 m long and 0.3 m high. The core thickness was 0.05 m, which should provide a freedom of deformation along the *y*-axis. A facing strip 0.5 mm thick and 0.7 m long rests on such a substructure. The geometry of the system is shown in Figure 5.

**Figure 5.** Numerical model of the problem of compression of a thin facing resting on a susceptible substructure.

In order to compress the facing, an area of 0.10 m × 0.05 m was determined on its two ends, to which a uniform pressure *p* was applied in the *x*-direction (tangent to the facing, opposite at the two ends). The load application area is distant from the edge of the substructure (0.25 m). The decision was made to apply the load distributed over the surface because the attempt to load the system in the form of a linear load applied to the edge of the facing caused too much local disturbance. As the system had to be supported, after several attempts, it was decided to support the bottom surface of the substructure. The displacement conditions *u<sup>y</sup>* = 0, *u<sup>z</sup>* = 0 were assumed on the entire bottom surface, and additionally, *u<sup>x</sup>* = 0 was assumed in the middle of this surface, which is shown in Figure 5. This support had a small influence on the behavior of the system, while ensuring its necessary stabilization. When trying to limit the displacements on the sides of the substructure, it turned out that these limitations affect the behavior of the system and cause stress disturbances. The attempt to define the boundary conditions identical to those in the analytical model was unsuccessful. The assumption that the horizontal displacement of the facing equals zero made it practically impossible to induce the appropriate stress state in this facing. The description of the model shows that despite all efforts, the numerical model has some deviations from the theoretical model in which the core (substructure) is an infinite elastic half-space.

*ν ν ν ν ν* The material parameters of the 3D numerical model corresponded to the analytical model. The facing material was assumed to be isotropic elastic (*E<sup>F</sup>* = 210 GPa, *ν<sup>F</sup>* = 0.3). In the case of an isotropic core, it was assumed *E<sup>C</sup>* = 4 MPa, *ν<sup>C</sup>* = 0.05. When the case with the orthotropic core was analyzed, its parameters were defined as follows: *E<sup>x</sup>* = 10 MPa, *E<sup>y</sup>* = 10 MPa, *E<sup>z</sup>* = 4 MPa, *νxy* = *νxz* = *νyz* = 0.05, *Gxy* = *Gxz* = *Gyz* = 3 MPa. The 3D model uses C3D8 solid elements (core) and S4 shell elements (facing), in which there is no reduced integration. Interaction between the facing and the core was defined using a TIE connection, which causes the displacements of nodes of one surface to be identical to the displacements of nodes on the other surface. The size of the finite element mesh was constant and equal to 0.01 m. It is worth mentioning that the problem of facing wrinkling is mesh-dependent. The mesh should be dense enough to allow deformation of the core and facing. Thus, the mesh size is dependent on the finite element itself (a shape function) as well as the properties of the facing and core materials. In the case of S4 finite elements (doubly curved general purpose shell, finite membrane strains), it is sufficient if there are two finite elements per half-wavelength *l*. In our case, the half-wavelength was in the order of 0.035–0.038 m; therefore, the size of the finite elements turned out to be small enough (0.01 m).

Since the phenomenon of face wrinkling is associated with the local deformation of the compressed face, a geometrically nonlinear static analysis and the Riks method were used. Due to the symmetry of the problem, it turned out to be beneficial to introduce into the model's initial imperfections as a linear combination of buckling modes of the structure. The buckling modes were solved independently. The size of the introduced imperfections was very small. The sum of four modes multiplied by 0.00001 was introduced, which

meant that the positions of the model nodes were disturbed about 0.025 mm. This means that the amplitude of the imperfection was 5% of the facing thickness.

The load applied to the model could increase to the value of 2000 kPa, which corresponds to a compressive force of 10 kN and a compressive stress in the facing of 400 MPa. Obviously, such a load value was never realized because the facing had previously buckled. The applied load level was determined on the basis of the LPF (Load Proportionality Factor) value.

Another interesting challenge of numerical analysis is the question of recognizing when a structure loses stability and when it does not. Unfortunately, as in real conditions, and unlike in analytical solutions, in a numerical solution, there is usually no unambiguous parameter indicating the state of the system (stable–unstable). Wrinkles in the compressed facing appear very quickly, which is illustrated in Figure 6a (only the facing was presented). A certain determinant of instability may be the appearance of a nonlinear relationship between LPF and arc length factor (Figure 6b), indicating the nonlinear nature of the process [25]. One should also pay attention to the difference between the compressive stress in the facing in the *x*-direction calculated on the basis of the currently applied force divided by the facing cross-section area and the stress obtained in the FE model that takes into account nonlinear effects, i.e., local deformations. The comparison of these stresses for the first eight load increments is presented in Table 1. The stress values estimated at the FE nodes are much higher due to the effect of the load acting on the distance resulting from the deformation of the facing. Since the material was originally assumed to be perfectly elastic, the stresses in the model can be very high.

*σ* **Figure 6.** Numerical solution of the compression of the elastic facing resting on the susceptible core: (**a**) *σ<sup>x</sup>* stress at the top of the facing (fifth load increment) (**b**) LPF–arc length relation.



*σ σ*

This situation, which is complex for evaluation, definitely changes after assuming that the facing material is perfectly elastic–plastic. Assuming the yield point *f<sup>y</sup>* = 270 MPa (the value is consistent with the characteristics of typical steel sheets used for the production of sandwich panels), in the ninth load increment, for LPF = 0.239, the LPF–arc length diagram breaks down (Figure 7), which corresponds to the theoretical facing compression stress 0.239 × 400 = 95.6 MPa.

**Figure 7.** Numerical solution of the compression of the elastic–plastic facing resting on the susceptible core: LPF–arc length relation.

A very similar relationship can be observed in the analysis of the problem with the orthotropic core. A breakdown of the LPF–arc length relationship occurred for LPF = 0.284, which corresponds to the theoretical stress 113.6 MPa. Of course, the stresses in the nodes of the model are different and reach the yield point of the material.

The obtained numerical results (95.6 MPa and 113.6 MPa) are close to the theoretical values (92.36 MPa and 107.8 MPa). Introduction of the yield stress for the facing material facilitates the interpretation of the numerical results. It is also worth paying attention to the fact that for the seventh or eighth load increment (Table 1), the stresses in the core reach the values close to the strength of typical core materials.

A number of numerical analyses were also carried out using the 2D model. The geometry and boundary conditions of this model corresponded to the geometry and boundary conditions of the 3D model. The main difference between the models was that they used plane (not spatial) finite elements: CPS4 for the core and B23 beam elements for the facing. It turned out that the 2D model behaves very similarly to the 3D model. Among other things, there are similar difficulties in interpreting the moment of loss of stability. This situation changes after assuming that the facing material is perfectly elastic–plastic. For the isotropic core, the LPF–arc length relationship is very similar to the relation presented in Figure 7, but the breakdown occurs for LPF = 0.227, which corresponds to the facing compression stress 0.227 × 400 = 90.8 MPa. For the orthotropic core, the extreme LPF value is 0.259, which corresponds to the stress of 103.6 MPa. These values are similar to the analytical results, although they are slightly lower than in the case of the 3D model.

#### **6. Parametric Analysis**

*Description of the Models*

Using the derived formulas, the influence of the material parameters of the orthotropic core on the value of the wrinkling stress was calculated and illustrated (Figure 8). Modulus *E<sup>x</sup>* = 10 MPa was assumed as constant. The modules *E<sup>z</sup>* and *Gxz* are variable. Moreover, each of the graphs corresponds to a different value of the Poisson ratio *νxz*, namely −1.0, 0.0, 0.5, and 1.0. For additional illustration of the problem, the graphs of the parameter *κ* are also presented in Figure 8.

The basic conclusions from the analysis of the graphs are quite obvious and consistent with the case of the isotropic core: the greater the stiffness of the core, the higher the wrinkling stress. It gets more interesting when *νxz* = 1, because with large *E<sup>z</sup>* and *Gxz* the parameter *κ* approaches −1 and the parameter *a* increases strongly. In the case when *νxz* = −1, the parameter *κ* takes values in the typical range (positive values), but with large values of *E<sup>x</sup>* and *Gxz*, the parameter *a* (44) reaches much higher values and grows faster than the parameter *m* drops (42). It should be emphasized that for the presented range of variability, *κ* > −1 and *a* is a positive value. *ν* − *κ*

**Figure 8.** *Cont*.

*κ ν* − *ν ν ν* **Figure 8.** Influence of core material parameters on the value of wrinkling stress and the parameter *κ*: (**a**) *νxz* = −1.0, (**b**) *νxz* = 0, (**c**) *νxz* = 0.5, (**d**) *νxz* = 1.0.

#### **7. Conclusions**

*ν κ* − *ν* − *κ κ* − The first part of the article contained a short survey of the classical solutions to the problem of instability of a facing resting on a homogeneous and isotropic substructure (core). It was presented how the assumptions concerning the displacement field affect the solution of the problem. Next, the dependence of the solution [6] on the value of the Poisson ratio was presented, and strain energy analyses were carried out to investigate the relationships between the individual components of the deformation energy of the core. In the second part of the paper, the derivation of the formula for the critical stress in the case of uniaxial compression of the thin facing resting on the orthotropic core was presented. The conditions for the existence of the solution were discussed, which in principle are met for a wide range of variability of material parameters. The numerical example confirming the compliance of the selected analytical solution with the numerical one was also presented. The article discussed the applied models in detail and explained the difficulties associated with determining the load and support boundary conditions. The presented numerical model has not been experimentally verified, although a similar model was verified in [19] for a core with the Poisson ratio ranging from 0 to 0.3. The third part of the article presented the results of the parametric analysis, i.e., the effect of changing the material parameters of the orthotropic core on the wrinkling stress. This type of analysis can be of great importance in the optimal design of sandwich systems where local loss of stability plays a significant role. The developed solution can be easily introduced into the optimization procedure.

The presented work confirms that the further development of analytical methods in solving the discussed problem is advisable and important, both from a scientific and

engineering point of view. Undoubted benefits also come from the possibility of numerical analysis of the issue under discussion. The applied FE models revealed that due to the local loss of stability, the stresses in the facing locally increase to the yield point, and the stresses in the core reach values similar to the strength of the core material. This means that if we want to accurately understand the stress state in the facing and the core, a relatively simple and attractive analytical approach should be supplemented with a numerical solution.

**Author Contributions:** Conceptualization, Z.P. and J.P.; methodology, Z.P.; software, Ł.S.; validation, I.K. and Ł.S.; formal analysis, J.P.; investigation, Z.P.; resources, Z.P. and I.K.; data curation, J.P.; writing—original draft preparation, Z.P.; writing—review and editing, I.K. and Ł.S.; visualization, Ł.S.; supervision, I.K.; project administration, Z.P.; funding acquisition, Z.P. and J.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This APC was funded by Poznan University of Technology, grant number 0411/SBAD/0004 and Czestochowa University of Technology.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

