Parametric Sensitivity of Shear Correction Factors for Multiwall Corrugated Structures

1. Introduction

Lightweight layered and cellular structures are widely used in engineering whenever high stiffness-to-weight and controllable mechanical response are required. Typical examples include sandwich panels, rib-stiffened shells, corrugated cores, honeycomb configurations, and architected metamaterials. In many applications, the global structural response is governed not only by material properties but also by geometry-driven mechanisms such as local bending of webs, shear deformation within compliant cores, and stress redistribution around voids. As a result, accurate analysis and design require modelling frameworks that can connect microstructural geometry to macroscopic stiffness and strength while remaining computationally efficient for parametric studies and optimization tasks. [1,2,3,4,5]

A common strategy in this context is to replace a detailed heterogeneous configuration with a homogenized or equivalent continuum model. For layered and periodic structures this often takes the form of an equivalent orthotropic plate or beam, typically derived using energy-based homogenization or computational representative volume elements (RVE). [6,7,8,9,10,11,12,13] While such reduced-order models offer major computational advantages, their accuracy depends critically on how they represent the internal deformation mechanisms of the original structure. Among these mechanisms, transverse shear is particularly challenging: it is negligible for thin homogeneous plates but becomes decisive for thick laminates, soft-core sandwich structures, and periodic cores with voids or inclined webs. [14,15,16,17,18,19,20]

In classical plate theory, shear deformation is neglected (Kirchhoff–Love model), leading to good predictions for slender plates dominated by bending. In contrast, in the Timoshenko–Mindlin family of models, transverse shear is included and a key parameter appears: the shear correction factor $k_s$. This factor compensates for the fact that the assumed through-thickness shear strain (often constant) does not match the true shear stress distribution, which depends on cross-sectional geometry and material heterogeneity. In homogeneous rectangular beams and plates, reference values of $k_s$ are available (e.g., the well-known factors for rectangular cross-sections). However, for strongly heterogeneous, perforated, or architected cross-sections, $k_s$ cannot be assumed universal and must be identified consistently. [21,22,23,24,25,26]

The need for reliable shear correction is especially evident in periodic cores and multi-layered configurations. For corrugated or honeycomb-like microstructures, shear stresses may localize in thin load paths, bypass void regions, and be redistributed by stiff faces. In such systems, the effective shear response is governed by an interplay of geometry (period, height, thickness, topology) and phase stiffness contrast. Consequently, inaccurate shear correction can lead to systematic errors in predicted deflections, stress resultants, and even in derived design conclusions when simplified plate models are used for parametric design. [27,28,29]

A second challenge is that shear correction factors are frequently introduced in reduced-order models in a way that is not energy-consistent, especially when one attempts to reuse simple textbook factors for cross-sections with complex heterogeneity. Energy-consistent identification is crucial because the effective shear stiffness must reproduce not only the global deflection but also the correct shear energy stored in the structure under transverse loading. For heterogeneous layered systems, energy-based equivalence provides a robust path to define the effective shear stiffness $\left(GA{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right.$ and to derive $k_s$ relative to an explicit reference model. [10,26]

A practical difficulty remains: even if an energy-consistent definition is adopted, computing $\left(GA{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right.$ repeatedly for many design variants can be expensive, particularly if full-scale finite element simulations are required for each geometric configuration. This becomes a bottleneck when the goal is not a single calibration but a systematic assessment of how sensitive $k_s$ is to changes in geometry and material properties across a realistic design space. Such sensitivity knowledge is important in at least three contexts: (i) robust design under manufacturing tolerances and paper variability, (ii) development of fast surrogate models for design optimization, and (iii) deciding which parameters must be measured and controlled to ensure predictive reliability of reduced-order models. [30,31,32]

Corrugated board is a representative and industrially relevant example of such a geometry-dominated, heterogeneous structure. It consists of flat paper liners bonded to one or more sinusoidally shaped fluting layers, forming 3-ply (single-wall) and 5-ply (double-wall) assemblies commonly used in packaging. Corrugated board is simultaneously lightweight, recyclable, and mechanically efficient, yet its response is strongly influenced by flute geometry, paper anisotropy, and moisture-dependent stiffness. Because its internal architecture contains voids and inclined webs, transverse shear deformation may contribute significantly in certain test configurations and structural applications, especially at short spans and under out-of-plane loading. [33,34,35,36,37]

A substantial body of research addresses bending stiffness, edge crush strength, and constitutive modelling of corrugated board under various loading conditions. Reduced-order modelling approaches, including homogenized orthotropic plates, are widely used to predict global stiffness and to support package design. [38,39] Nevertheless, compared to bending-dominated metrics, the identification and systematic characterization of shear correction remains less developed, particularly for multi-wall boards and for broad variations in geometry and material stiffness that arise across flute families and paper grades.

In addition, sensitivity studies reported in the literature often focus on local effects (small perturbations around a baseline configuration) or on a limited set of parameters. While such analyses can be informative, they do not necessarily capture the full range of coupled effects present in corrugated structures. When a response is governed by nonlinear geometric interactions, void topology, and stiffness contrast, the ranking of influential parameters may change across the domain. In this case, global sensitivity analysis, based on space-filling sampling strategies such as Latin Hypercube Sampling (LHS), provides a more appropriate framework to evaluate parameter influence, detect interactions, and identify robust trends [31,40].

Motivated by these gaps, the present work investigates the sensitivity of the shear correction factor $k_s$ for corrugated board across a broad multi-parameter space, focusing on both 3-ply and 5-ply configurations. The study adopts a pixel-based, energy-consistent procedure for computing $\left(GA{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right.$ directly from the cross-sectional geometry and spatially varying shear modulus fields, avoiding ad-hoc assumptions and enabling high-throughput evaluation [41,42,43,44,45,46,47]. The parameter space is explored using LHS, and sensitivity is quantified in a normalized, unit-free manner using measures that directly express percentage-to-percentage influence, enabling an interpretable ranking of drivers across variables with different physical dimensions [10,35].

The main contributions of this paper are: (a) a global (nonlocal) sensitivity assessment of $k_s$ for corrugated paperboard performed over a wide admissible parameter space, rather than relying on local perturbations around a single reference design. This reveals dominant drivers and interactions relevant for real product families and manufacturing variability; (b) a unified framework for 3-ply and 5-ply boards, explicitly allowing distinct geometries and stiffnesses of two flutings in the double-wall case, which is essential for capturing multi-layer coupling effects that cannot be inferred from single-wall trends; (c) energy-consistent, pixel-based evaluation of $\left(GA{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right.$ enabling high-throughput computation of $k_s$ across thousands of sampled designs, offering a practical route toward surrogate modelling, robust design, and reliable shear-corrected reduced-order models; (d) Normalized, unit-free sensitivity indices (%/%), based on log-elasticities and complementary global measures, providing directly interpretable statements such as “a 1% change in parameter $x$produces an average $S$ % change in $k_s$” across the explored domain.

These outcomes are important because they help transform $k_s$from a loosely chosen tuning factor into a geometry- and material-dependent quantity with quantified uncertainty and sensitivity, improving the reliability of reduced-order models used in packaging mechanics and other periodic-core structural applications.

2. Materials and Methods

2.1. Scope of the Study and Investigated Configurations

This study quantifies the sensitivity of the transverse shear correction factor $k_s$ to the key geometric and material parameters of corrugated paperboard. Two structural families are considered: (i) 3-ply (single-wall): top liner + sinusoidal fluting + bottom liner; (ii) 5-ply (double-wall): top liner + fluting 1 + middle liner + fluting 2 + bottom liner.

For the double-wall configuration, the two flutings may exhibit different periods, heights and thicknesses, enabling a systematic evaluation of geometric interactions. The analysis is conducted over a wide admissible parameter space using global sampling (Latin Hypercube), which allows identification of dominant parameters and parameter couplings beyond local perturbations around a single reference design.

2.2. Definition of the Shear Correction Factor

The shear correction factor is defined by comparing the effective transverse shear stiffness of the heterogeneous cross-section with the shear stiffness of a flat reference model:

$$k_s={\displaystyle \frac{\left(GA{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right.}{G^{\ast }{A}^{\ast }}}.$$

(1)

The reference area is taken as

$$A^{\ast }=b H,$$

(2)

where $b$ denotes the out-of-plane width (unit width is adopted in computations) and $H$ is the adopted reference height corresponding to the total thickness of the considered corrugated board configuration.

To ensure a stable normalization across broad geometric and material variations, a thickness-weighted reference shear modulus is used:

$$G^{\ast }=\sum _{k=1}^{n_{\mathcal{l}}}{w}_k G_k,w_k={\displaystyle \frac{t_k}{\sum _{j=1}^{n_{\mathcal{l}}}{t}_j}},$$

(3)

where $t_k$ and $G_k$ are the thickness and shear modulus of the $k$-th layer (liners and flutings). This definition avoids imposing a priori that only one phase carries shear and ensures consistent comparison between 3-ply and 5-ply configurations.

2.3. Parametric Geometry and Pixel-Based Representation

2.3.1. Sinusoidal Fluting Geometry

The fluting centerline is defined by:

$$z\left(x\right)=z_0+{\displaystyle \frac{h}{2}}\mathrm{c}\mathrm{o}\mathrm{s}\left({\displaystyle \frac{2\pi }{p}}x+\phi \right),$$

(4)

with period $p$, height $h$(peak-to-peak), vertical offset $z_0$, and phase shift $\phi$. Finite thickness $t_f$ is introduced by constructing a constant-thickness band around the centerline.

For 5-ply boards, two flutings are described independently by $\left(p_1,h_1,t_{f1}\right)$ and $\left(p_2,h_2,t_{f2}\right)$ and are placed into two non-overlapping core bands between (top, middle) and (middle, bottom) liners, respectively.

2.3.2. Pixel Image Generation and Phase Encoding

The cross-section is represented by a grayscale pixel image generated in MATLAB. The following encoding is used:

• background (void): intensity $1.0$,
• liners: intensity $0.5$,
• fluting(s): intensity $0.0$.

The physical pixel size $\mathrm{p}\mathrm{i}\mathrm{x}$ (mm/pixel) controls geometric resolution. To prevent loss of thin layers and boundary artefacts, the admissible domain is constrained such that:

$$t_{\mathrm{m}\mathrm{i}\mathrm{n}}\ge 2 \mathrm{p}\mathrm{i}\mathrm{x},$$

(5)

where $t_{\mathrm{m}\mathrm{i}\mathrm{n}}$ is the minimum thickness among all considered liners and flutings. The same generator is employed both for methodological figures and for the numerical evaluation of $k_s$, ensuring full consistency between geometry visualization and computation.

2.4. Pixel-Based Identification of $\left(GA{)}_{eff}\right.$

2.4.1. Energy-Consistent Formulation

The effective transverse shear stiffness is identified by enforcing equivalence between the shear strain energy of a flat Timoshenko model and the heterogeneous cross-section under the same shear resultant $V$. For a representative segment of length $L$, the flat-model energy reads:

$$U_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{t}}={\displaystyle \frac{1}{2}} {\displaystyle \frac{V^2}{\left(GA{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right.}} L.$$

(6)

For a heterogeneous cross-section, the shear energy may be expressed using a normalized 2D shear-stress shape function $\varphi (x,z)$ defined such that:

$${\tau }_{xz}(x,z)=V \varphi (x,z),{\int }_A\varphi (x,z) \mathrm{d}A=1,$$

(7)

which guarantees that integration of ${\tau }_{xz}$ over the cross-section yields the shear resultant $V$. Under linear elasticity, the heterogeneous shear energy becomes:

$$U_{\mathrm{p}\mathrm{i}\mathrm{x}}={\displaystyle \frac{1}{2}}{\int }_A{\displaystyle \frac{{\tau }_{xz}^2(x,z)}{G(x,z)}} \mathrm{d}A={\displaystyle \frac{1}{2}}{V}^2{\int }_A{\displaystyle \frac{{\varphi }^2(x,z)}{G(x,z)}} \mathrm{d}A.$$

(8)

By imposing $U_{\mathrm{f}\mathrm{l}\mathrm{a}\mathrm{t}}=U_{\mathrm{p}\mathrm{i}\mathrm{x}}$, the effective shear stiffness is obtained as:

$${\left(GA\right)}_{\mathrm{e}\mathrm{f}\mathrm{f}}={\left({\int }_A{\displaystyle \frac{{\varphi }^2(x,z)}{G(x,z)}} \mathrm{d}A\right)}^{-1}.$$

(9)

Finally, $k_s$ follows from Eq. (1).

2.4.2. Pixel Discretization and Material Mapping

The cross-section is discretized into pixels (or subpixels) of area $\mathsf{\Delta }A$. Void pixels are excluded from the mechanical domain. Each material pixel $i$ is assigned a shear modulus $G_i$ based on its phase label (liner vs fluting; in 5-ply cases, fluting 1 and fluting 2 may have different $G$). This yields a discrete material field $G(x,z)\to G_i$ and a corresponding set of integration weights $\mathsf{\Delta }{A}_i$.

2.4.3. Construction of the 2D Shear-Stress Shape Function $\varphi (x,z)$

To avoid introducing additional modelling assumptions beyond geometry and material phase distribution, the shear-stress shape function $\varphi (x,z)$ is constructed directly from the pixelized cross-section using a geometry-consistent 2D formulation. The approach relies on the classical relationship between shear stresses and the first moment of area (static moment) while retaining full 2D resolution required by corrugated geometries.

First, the area and centroid of the material domain are computed from pixel quadrature:

$$A={\int }_{A}1 \mathrm{d}A\approx \sum _{i=1}^{N_p}\mathsf{\Delta }{A}_i, \bar{z}={\displaystyle \frac{1}{A}}{\int }_{A}z \mathrm{d}A\approx {\displaystyle \frac{1}{A}}\sum _{i=1}^{N_p}{z}_i \mathsf{\Delta }{A}_i.$$

(10)

Next, the second moment of area about the centroidal axis is evaluated:

$$I={\int }_A{(z-\bar{z})}^2 \mathrm{d}A\approx \sum _{i=1}^{N_p}{(z_i-\bar{z})}^2 \mathsf{\Delta }{A}_i.$$

(11)

To define the shear-stress distribution across the cross-section, a 2D static-moment field $Q(x,z)$ is introduced. For a point $\left(x,z\right)$ in the material domain, $Q(x,z)$ represents the first moment of area of the portion of the cross-section located above that point (with respect to the centroidal axis). In discrete form, for every material pixel $i$ at coordinates $\left(x_i,z_i\right)$, we define:

$$Q_i\equiv Q(x_i,z_i)={\int }_{A(z\ge z_i)}\left(z-\bar{z}\right) \mathrm{d}A\approx \sum _{\begin{array}{c}j=1\\ z_j\ge z_i\end{array}}^{N_p}(z_j-\bar{z}) \mathsf{\Delta }{A}_j.$$

(12)

The local thickness measure entering the shear-stress relation is defined from the pixel domain as an effective local width $b(x,z)$, computed as the material measure in the horizontal direction associated with the given pixel row. In the discretized setting, this can be expressed as:

$$b(z_i)\approx {\displaystyle \frac{1}{\mathsf{\Delta }z}}\sum _{\begin{array}{c}j=1\\ z_j\in [z_i-\Delta z/2,z_i+\Delta z/2]\end{array}}^{N_p}\mathsf{\Delta }{A}_j,$$

(13)

which corresponds to the length of material present at the level $z=z_i$ (for unit out-of-plane width). This definition naturally captures the fact that corrugated sections may be discontinuous along $x$ at a given height due to void regions.

With these ingredients, the unnormalized shear-stress shape field is computed as:

$${\tilde{\varphi }}_i={\displaystyle \frac{Q_i}{I b\left(z_i\right)}}.$$

(14)

Finally, the field is normalized to satisfy the resultant condition in Eq. (7):

$${\varphi }_i={\displaystyle \frac{{\tilde{\varphi }}_i}{\sum _{k=1}^{N_p}{\tilde{\varphi }}_k \mathsf{\Delta }{A}_k}}.$$

(1)

Equations (12)–(15) define a fully 2D, geometry-consistent shape function computed directly from the pixel representation. This construction avoids prescribing shear localization a priori and allows the subsequent energy evaluation to reflect both geometric complexity (corrugation) and spatially varying shear modulus $G(x,z)$.

2.4.4. Discrete Evaluation of $\left(GA{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right.$ and $k_s$

With ${\varphi }_i$and $G_i$defined on material pixels, Eq. (9) becomes:

$${\left(GA\right)}_{\mathrm{e}\mathrm{f}\mathrm{f}}={\left(\sum _{i=1}^{N_p}{\displaystyle \frac{{\varphi }_i^2}{G_i}} \mathsf{\Delta }{A}_i\right)}^{-1}.$$

(16)

The shear correction factor is then computed from Eq (1):

All computations are performed for a unit width $b=1$, and the reported values of $k_s$ are therefore independent of scaling in the out-of-plane direction.

2.5. Global Sensitivity Analysis Using Latin Hypercube Sampling

2.5.1. Parameter Vectors

For the 3-ply configuration, the parameter vector is:

$${\mathbf{x}}_3=\left[p,\hspace{0.25em} h,\hspace{0.25em} t_f,\hspace{0.25em} t_{L,t},\hspace{0.25em} t_{L,b},\hspace{0.25em} G_f,\hspace{0.25em} G_L\right].$$

(17)

For the 5-ply configuration:

$${\mathbf{x}}_5=\left[p_1,h_1,t_{f1},\hspace{0.25em} p_2,h_2,t_{f2},\hspace{0.25em} t_{L,t},t_{L,m},t_{L,b},\hspace{0.25em} G_{f1},G_{f2},G_L\right].$$

(18)

This parametrization enables separation of geometric and material influences and explicitly allows different flute families in the double-wall board.

2.5.2. Parameter Ranges and Feasibility Constraints

The admissible ranges are selected to cover representative corrugated geometries and material variability while ensuring stable discretization:

3-ply:

• $p\in [4.0, 8.5]$mm,
• $h\in [1.0, 4.5]$mm,
• $t_f\in [0.10, 0.35]$mm,
• $t_{L,t},t_{L,b}\in [0.10, 0.40]$mm,
• $G_f\in [0.15, 1.50]$GPa,
• $G_L\in [0.25, 2.50]$GPa.

5-ply:

• $p_1,p_2\in [4.0, 9.5]$mm,
• $h_1,h_2\in [1.0, 5.0]$mm,
• $t_{f1},t_{f2}\in [0.10, 0.35]$mm,
• $t_{L,t},t_{L,m},t_{L,b}\in [0.10, 0.40]$mm,
• $G_{f1},G_{f2}\in [0.15, 1.50]$GPa,
• $G_L\in [0.25, 2.50]$GPa.

All samples must satisfy: (i) no overlap of phases (enforced by construction in the image generator), and (ii) the resolution condition in Eq. (5).

2.5.3. LHS Sampling Strategy

A Latin Hypercube Sampling (LHS) design of size $N$is generated for each configuration family. Geometric variables are sampled from uniform distributions within their bounds. To represent multiplicative uncertainty in stiffness, shear moduli are sampled from a log-uniform distribution:

$$\log G\sim \mathcal{U}\left(\log G_{\mathrm{m}\mathrm{i}\mathrm{n}},\log G_{\mathrm{m}\mathrm{a}\mathrm{x}}\right).$$

(19)

For each sampled vector ${\mathbf{x}}^{\left(n\right)}$, the computational pipeline is:

1.

generate pixel geometry,

2.

assign shear modulus field $G_i$,

3.

compute the 2D shape function ${\varphi }_i$using Eqs. (12)–(15),

4.

evaluate $\left(GA{)}_{\mathrm{e}\mathrm{f}\mathrm{f}}\right.$by Eq. (16),

5.

compute $k_s$by Eq. (1).

This produces the dataset $\left\{({\mathbf{x}}^{\left(n\right)},k_s^{\left(n\right)}){\}}_{n=1}^N\right.$ used in subsequent analysis.

2.5.4. Normalized global sensitivity metrics

To quantify sensitivity in a unit-free and comparable manner across variables, two complementary normalized indices are adopted. First, we use log-elasticities, defined as

$$S_j={\displaystyle \frac{\partial \ln k_s}{\partial \ln x_j}},$$

(20)

which directly represent the percentage change of $k_s$ per percentage change of the input variable $x_j$. In practice, $S_j$ is estimated by fitting a linear model in logarithmic space,

$$\ln k_s={\beta }_0+\sum _{j=1}^d{\beta }_j\ln x_j+\epsilon ,$$

(21)

so that $S_j\approx {\beta }_j$. Second, we compute partial rank correlation coefficients (PRCC), which quantify the influence of each variable while controlling for the remaining inputs and are robust to monotonic nonlinearities. Both indices are evaluated on the dataset generated by Latin Hypercube Sampling and are reported as normalized, dimensionless sensitivity measures.

2.6. Computational Implementation and Reproducibility

All calculations are performed in MATLAB. Geometry generation and pixel-based integration are fully automated to enable high-throughput evaluation of $k_s$ across LHS samples. For stochastic designs, a fixed pseudo-random seed is used to ensure reproducibility. A mesh-resolution robustness check is conducted on a representative subset of samples to confirm that the computed $k_s$ values and sensitivity trends are not affected by pixel discretization artefacts.

3. Results

3.1. Dataset Overview and Numerical Stability

Global datasets for 3-ply and 5-ply corrugated boards were generated using Latin Hypercube Sampling (LHS) within the prescribed geometric and material parameter ranges summarized in Table 1. For each LHS realization, the corresponding geometry was constructed directly in the pixel domain, and the shear correction factor $k_s$ was computed using the pixel-based energy-consistent method described in Section 2. After enforcing feasibility and minimum-resolution constraints, the final datasets consisted of $N_3$ valid samples for the 3-ply configuration and $N_5$ samples for the 5-ply configuration.

$p$$h$$t_f$$p_2$$h_2$$t_{f2}$$t_{L,t}$$t_{L,m}$$t_{L,b}$$G_f$$G_L$Representative examples of the pixel-based geometries generated within these ranges are shown in Figure 1, illustrating the fundamental structural differences between the single-flute 3-ply board and the double-flute 5-ply configuration.

To verify the numerical robustness of the pixel discretization, a resolution sensitivity study was performed on a representative subset of samples by repeating the computations using a refined pixel size. The relative changes in the shear correction factor, $\mathsf{\Delta }{k}_s/k_s$, are summarized in Table 2.

The limited magnitude of these variations confirms that the global trends discussed in the following sections are not governed by discretization artefacts but reflect genuine geometric and material effects.

3.2. Distribution of the Shear Correction Factor $k_s$

Across the investigated parameter space, the computed shear correction factors exhibit a pronounced spread, reflecting the strong dependence of transverse shear behavior on both geometry and material contrast. The empirical distributions obtained from the LHS datasets for the 3-ply and 5-ply configurations are shown in Figure 2.

While the 3-ply boards exhibit a relatively compact distribution, the 5-ply boards display a markedly wider dispersion and heavier tails, as evident in Figure 2(b). A quantitative summary of these distributions is provided in Table 3.

3.3. Global Sensitivity: Normalized Effects (%/%)

To quantify parameter influence in a unit-free and directly interpretable manner, two complementary global sensitivity measures were employed: (i) log-elasticities (Eq. (20)) and (ii) Partial Rank Correlation Coefficients (PRCC).

The resulting normalized sensitivity indices for the 3-ply and 5-ply configurations are reported in Table 4 and Table 5, respectively.

$p$$h$$t_f$$t_{L,t}$$t_{L,b}$$G_f$$G_L$

$p_1$$h_1$$t_{f1}$$p_2$$h_2$$t_{f2}$$t_{L,t}$$t_{L,m}$$t_{L,b}$$G_f$$G_L$For a concise visual ranking of the dominant drivers, the log-elasticities are displayed in the tornado plots shown in Figure 3.

3.4. Response Surfaces for the Dominant Geometric Drivers

For the 3-ply configuration, the response surface $k_s(p,h)$evaluated at median values of the remaining parameters is shown in Figure 4.

Analogous response surfaces for the 5-ply configuration are presented in Figure 5, where each fluting family is examined separately.

3.5. Material Contrast Effects and Key Interactions

Material contrast between liners and fluting layers, typically expressed via shear-modulus ratios, significantly affects the shear correction factor. The dependence of $k_s$ on modulus contrast is illustrated in Figure 6.

Interaction effects between material contrast and geometry are further examined using binned plots shown in Figure 7.

3.6. Comparison Between 3-ply and 5-ply Sensitivity Patterns

A direct comparison of normalized sensitivities for the 3-ply and 5-ply configurations is presented in Figure 8, allowing universal and configuration-specific drivers of $k_s$ to be identified.

4. Discussion

The present study provides a comprehensive, global assessment of the shear correction factor $k_s$ for multiwall corrugated structures, combining large-scale parametric sampling with an energy-consistent, pixel-based formulation. By systematically analyzing both 3-ply and 5-ply configurations, the results allow clear identification of universal drivers, configuration-specific effects, and key geometric–material interactions governing transverse shear behavior.

4.1. Reliability of the Dataset and Numerical Robustness

A prerequisite for meaningful sensitivity analysis is numerical robustness. The adopted parameter ranges and sampling strategy, summarized in Table 1, ensure broad coverage of realistic corrugated-board configurations while maintaining physical feasibility. The pixel-based geometries illustrated in Figure 1 confirm that the discretization captures essential geometric features for both single- and double-flute systems.

The resolution sensitivity study reported in Table 2 demonstrates that the computed shear correction factors are only weakly affected by pixel refinement. The small median and upper-percentile deviations indicate that the global trends identified in subsequent analyses are intrinsic to the structural parameters rather than artefacts of discretization. This robustness is particularly important for the 5-ply case, where increased geometric complexity could otherwise amplify numerical sensitivity.

4.2. Global Variability of the Shear Correction Factor

The distributions of $k_s$obtained from Latin Hypercube Sampling, shown in Figure 2, reveal substantial variability across the investigated parameter space. The statistical descriptors provided in Table 3 confirm that this variability is significantly larger for 5-ply boards than for 3-ply boards. This finding reflects the additional deformation mechanisms introduced by the second fluting and the middle liner, which increase both geometric freedom and interaction effects.

The broader dispersion and heavier tails observed in Figure 2(b) indicate that extreme combinations of geometry and material contrast can lead to shear correction factors far from classical reference values. This result directly challenges the common practice of adopting a single, constant shear correction factor for corrugated structures, especially for multiwall configurations.

4.3. Dominant Parameters and Normalized Sensitivity Patterns

The normalized sensitivity analysis provides a quantitative ranking of parameter influence in a unit-free form. The results reported in Table 4 (3-ply) and Table 5 (5-ply) show that geometric parameters associated with fluting shape consistently dominate the response. This conclusion is reinforced by the tornado plots in Figure 3, where flute height and period emerge as the primary drivers of $k_s$ across the global parameter space.

While log-elasticities and PRCC largely agree in their ranking, discrepancies for selected parameters suggest the presence of nonlinear or non-monotonic effects. Such behavior is expected in corrugated structures, where local changes in geometry can alter load paths and energy distribution in a non-proportional manner. Importantly, the sensitivity patterns differ between 3-ply and 5-ply boards, indicating that conclusions drawn from single-wall systems cannot be directly transferred to multiwall configurations.

4.4. Physical Interpretation via Response Surfaces

For the 3-ply configuration, the response surface shown in Figure 4 reveals strong nonlinear coupling between flute period and height. Regions of high $k_s$ correspond to geometries that promote more uniform shear strain distribution across the cross-section, whereas low values are associated with slender or widely spaced flutes that concentrate shear deformation.

The 5-ply response surfaces in Figure 5 further demonstrate that each fluting family contributes distinctly to the global shear response. Fixing one flute while varying the other reveals asymmetric influence patterns, highlighting the role of inter-layer interactions. These results confirm that the effective shear correction factor in multiwall boards cannot be interpreted as a simple superposition of individual flute contributions.

4.5. Role of Material Contrast and Interaction Effects

Beyond geometry, material contrast between liners and fluting layers plays a critical role in shaping transverse shear behavior. The scatter plots in Figure 6 show that $k_s$ varies systematically with shear-modulus ratios such as $G_f/G_L$, but the relationship is strongly modulated by geometry. This observation underscores that material contrast alone is insufficient to predict shear correction behavior without considering geometric context.

The interaction plots in Figure 7 provide further insight by demonstrating how stiffness contrast amplifies or suppresses geometric effects. For certain modulus ratios, changes in flute height lead to pronounced variations in $k_s$, while for others the response becomes comparatively insensitive. Such interactions explain the partial discrepancies observed between log-elasticities and PRCC in Table 4 and Table 5, and they emphasize the need for multidimensional analysis when assessing corrugated-board shear behavior.

4.6. Comparison Between 3-ply and 5-ply Systems

A direct comparison of normalized sensitivities, presented in Figure 8, highlights both universal and configuration-specific drivers of the shear correction factor. Flute geometry emerges as a universal determinant of $k_s$ across both systems, confirming its fundamental role in transverse shear deformation.

However, the 5-ply configuration exhibits additional influential parameters associated with the middle liner and the second fluting layer. These configuration-specific effects lead to higher variability and more complex interaction patterns, as already suggested by the broader distributions in Figure 2 and the response surfaces in Figure 5. Consequently, design rules or simplified models calibrated for 3-ply boards may be inadequate or even misleading when applied to multiwall corrugated structures.

4.7. Implications for Modeling and Design

Taken together, the results presented in Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8 and Table 1, Table 2, Table 3, Table 4 and Table 5 demonstrate that the shear correction factor for corrugated boards is not a fixed material constant but a highly configuration-dependent quantity. The strong sensitivity to geometry, material contrast, and their interactions implies that simplified assumptions, such as adopting a universal shear correction factor, can lead to significant modeling errors, particularly for multiwall systems.

The proposed energy-consistent, pixel-based framework provides a systematic means of quantifying these effects and can serve as a foundation for reduced-order models or design-oriented correlations. Importantly, the distinction between universal and configuration-specific drivers identified in Figure 8 offers guidance for developing simplified yet physically grounded design rules for corrugated structures of varying complexity.

5. Conclusions

This study presented a comprehensive global sensitivity analysis of the shear correction factor $k_s$ for 3-ply and 5-ply corrugated board structures using a pixel-based, energy-consistent numerical framework combined with Latin Hypercube Sampling. Based on the results reported in Table 1, Table 2, Table 3, Table 4 and Table 5 and Figure 1, Figure 2, Figure 3, Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8, the following conclusions can be drawn:

1.

The resolution sensitivity study confirmed that the computed shear correction factors are only weakly affected by pixel refinement. The small relative variations observed for both 3-ply and 5-ply configurations demonstrate that the identified trends are intrinsic to the structural parameters rather than discretization artefacts, supporting the reliability of the proposed methodology.

2.

The distributions of $k_s$reveal substantial variability across the investigated parameter space, with significantly broader ranges for 5-ply boards compared to 3-ply boards. This increased dispersion reflects the additional deformation mechanisms and interaction effects introduced by multiple flutings and intermediate liners, indicating that $k_s$ cannot be treated as a universal constant.

3.

Normalized sensitivity indices consistently identify fluting geometry—particularly flute height and period—as the primary drivers of the shear correction factor in both configurations. These parameters exert a stronger influence on $k_s$ than material properties alone, underscoring the central role of geometric design in governing transverse shear behavior.

4.

Response surfaces and interaction plots demonstrate pronounced nonlinear coupling between geometric parameters as well as between geometry and material contrast. In multiwall configurations, the contribution of each fluting family is not additive, and interactions between layers play a decisive role in shaping the effective shear response.

5.

While shear-modulus ratios between fluting and liners significantly affect $k_s$, their influence is strongly dependent on the underlying geometry. Material contrast can either amplify or suppress geometric effects, explaining observed differences between sensitivity measures and reinforcing the need for multidimensional analysis.

6.

Although certain drivers of $k_s$, such as fluting geometry, are universal, multiwall (5-ply) boards exhibit additional configuration-specific sensitivities associated with the middle liner and the second fluting layer. Consequently, conclusions drawn from single-wall systems cannot be directly extrapolated to more complex corrugated structures.

7.

The results demonstrate that adopting a single, fixed shear correction factor—commonly used in simplified beam and plate models—may lead to significant inaccuracies, particularly for multiwall corrugated boards. The presented framework provides a physically grounded basis for developing configuration-dependent shear correction factors suitable for advanced numerical modeling and design-oriented reduced-order formulations.

Overall, this work establishes a systematic methodology for quantifying the sensitivity of shear correction factors in corrugated structures and provides new insight into the interplay between geometry and material contrast. The findings offer a foundation for improved mechanical modeling of multiwall corrugated boards and motivate the development of refined design rules that account explicitly for structural configuration.