Equivalent Double-Layer Model for a Laminated Glass: A Proof of Concept Open Access

Ships and offshore structures are exposed to large environmental loads from waves, wind, ice, and currents [1]. Typically, these are steel structures, but alternative materials have also been used increasingly over the last decades. Glass is used in marine structures to ensure observations from the maritime environment during shipping and offshore operations, but as an enabling element allowing immersion of passengers to the maritime environment when enjoying life. The fact that glass is a brittle material, however, challenges its use when exposed to random loads from the maritime environment. One way to improve the brittleness of glass structures is to build the windows from laminated glass. Laminated glass has at least two glass panes joined by a polymeric interlayer. This lamination increases safety because the interlayer retains the glass fragments in place in case of an accidental failure, preventing the complete collapse of the window [2] (see Fig. 1 [3]). Typical applications for laminated glass are shown in Fig. 2 [4].

In the design of ships and offshore structures, the challenge is, however, that the structures are (a) large and complex (e.g., see Fig. 3) and (b) the time to design these, often prototype structures, is limited. The lamination layer is often considerably thinner than the glass panes it connects. It is also often much weaker in terms of stiffness. Thus, the laminated glass does not behave as a monolithic entity. Consequently, modeling large and complex structures cannot be easily done with commercial finite element (FE) codes as the kinematics of typical elements are often limited to the first-order shear deformation theory (FSDT), i.e., the so-called Reissner–Mindlin elements. Instead, continuum shells with additional degrees-of-freedom are needed for each layer of the laminated glass. Hence, the solution time increases. Combined with thickness optimization, in which different design alternatives are tested for optimal design, the model must be remeshed every iteration. This methodology adds to the computational burden and thus reduces the efficiency of structural design.

Therefore, different theories for laminated plates seem tempting. For example, Reissner [5] formulated the differential equations for sandwich plates, including the effect of transverse shear and normal stress deformation (that are constant through thickness) where the through-thickness behavior is continuous. Reddy [6] presented a generalized higher-order shear deformation theory for laminated plates, including the transverse shear and normal stress deformation (that are not constant through thickness) where the through-thickness behavior is also continuous. Carrera [7] reviewed zig-zag theories for laminated plates where the through-thickness behavior of the element is described by functions that may be continuous or discontinuous over the thickness of the laminate. While such theories are promising, they are not available in commercial FE packages. Another way to model the laminate glass is to use the equivalent thickness method [8]. One replaces the multilayer plate with one layer whose thickness has been derived based on stiffness or stress equivalency. This approach allows using a single-layer shell element in any commercial FE package. The downside is that it requires at least two separate analyses to obtain correct deflection and stress. Furthermore, it requires a shape parameter that depends on the window shape, size, and boundary conditions. Hence, such an approach may not be efficient in thickness optimization or modeling large complex structures. On the other hand, Liang et al. [9] proposed an equivalent shell element for laminated glass, showing great potential, but unfortunately, it is not readily available in any commercial FE packages.

The purpose of this article is to propose a methodology for the strength analysis of laminated glass structures by using any commercial FE codes that contain FSDT shell elements and constraint equation option in modeling. We test this methodology in four cases, which include the asymmetry of the laminate, shear transfer between the layers, and geometrical nonlinearity by von Karman strains. We validate the approach by detailed continuum shell models in which different layers are explicitly modeled.

The idea is to relax the through-thickness linearity assumption of strains and further the stresses in the following manner: (1) The FSDT-induced strain field is piecewise constant through the thickness. Each layer has constant normal strains, and the variation of these piecewise constant strains is described by linear distribution over the mid-planes of each layer (see Fig. 5). (2) The strain field that is linearly varying over each layer is introduced to the modified FSDT field. This layerwise local strain field is of pure bending type and thus equals zero at the mid-plane of each layer. The strain fields (1) and (2) are added for the full strain field. At the monolithic limit, this results in the classical FSDT strain field and at the limit of zero out-of-plane shear stiffness of the interlayer at the dual layer, i.e., $Df=D1+D2+..Dn$⁠.

Romanoff et al. [11] presented a structural analysis approach for a web-core sandwich plate by dividing it into two separate FE meshes to alleviate computational burden. Since the structure of the sandwich plate is analogous to that of laminated glass, the same approach is adopted. Consider a laminated glass consisting of two glass layers and one interlayer in between. These three layers are discretized as two different element sets that are overlayed with each other. These are called the laminate and stabilizing layers, which share the same $x$- and $y$-coordinates for nodes but are separated infinitesimally in the $z$-direction as shown in Fig. 6. The laminate layers consist of the two glass panes and the one interlayer. Its purpose is to include the shear transfer through the interlayer. The stabilizing layers consist only of the individual bending layers of the two glass panes. Its purpose is to carry more load when the interlayer shear modulus weakens. Thus, the two layers act as planer spring fields in series, sharing the same out-of-plane displacement, i.e., deflection, $wstab=wlam$⁠. Without the stabilizing layer, the displacement of the structure would grow to infinite as the interlayer shear modulus gets smaller.

The transverse shear stiffness of the stabilizing layer consists only of the glass layers in Eq. (6) (i.e., $k=1,3$⁠).

Finally, the boundary conditions and the external loads are only applied to the laminate layer. The coinciding nodes in the $x$- and $y$-coordinates of the two layers have shared $z$-translation. That is, the nodes on top of each other have identical deflection. See Fig. 7. Note that $x$- or $y$-translation is not shared.

The amount each layer carries the load depends on the shear modulus of the interlayer. When the shear modulus is near zero, i.e., no shear transfer, the stabilizing layer carries all the load. On the other hand, the laminate layer starts to carry more load as the shear modulus increases. We may observe this behavior by plotting the bending moment for each layer. For example, see Fig. 8 where a square laminate glass with a side length of 1000 mm and a thicknesses of $(t1,tint,t2)=(5,1,5)$ mm is subjected to uniformly distributed load under simply supported boundary conditions in a linear analysis. When the shear modulus of the interlayer is 0 MPa, i.e., no shear transfer, the stress resultants in the laminate layer are zero. When the shear modulus is 24,476 MPa (calculated with 70 GPa elastic modulus and 0.43 Poisson’s ratio), i.e., full shear transfer, the laminate layer carries much greater portion of the load than the stabilizing layer.

In linear analysis, the bending loads do not produce normal forces; $N11=N22=NN12=0$⁠. For example, Fig. 9 shows the normal stress distribution through the thickness due to bending moments. In no shear transfer case, the laminate layer has no stresses. On the other hand, in full shear transfer, the stresses in the stabilizing layer become very small.

We first present (1) a 3D solid-shell FE model of laminated glass and then (2) a shell 2D FE model where the theory presented in Sec. 2.2 is implemented. We use the 3D model to validate the proposed 2D modeling technique. That is, we have two insulating glass unit (IGU) models: one with 3D solid-shell elements and one with 2D shell elements. The latter may have a coarser mesh; hence, the total nodal degrees-of-freedom are smaller, which is computationally more attractive. We use ansys mechanical apdl in all the analyses.

The glass panes and the interlayer are modeled with solid-shell elements (SOLSH190 in ansys). They are used to model thin to moderately thick shell structures. The element has a 3D topology and eight nodes with three translation degrees-of-freedom each. The advantage of solid-shell to shell elements is that through-thickness strain and stress profiles are available. The advantage over a full solid element is that a solid shell is less prone to shear locking in bending dominant problems. Hence, we can use only one element through the thickness for each layer with solid-shell elements. However, one must avoid a mesh that is too coarse because, with a thin interlayer, the aspect ratio (width to thickness) of the element increases rapidly when the mesh size increases. Too high aspect ratios can cause locking problems, which requires a mesh convergence analysis. The mesh converge analysis has been done for the same model by Heiskari et al. [12]. They found a mesh size of 50 mm sufficient there and used it in this study.

One cannot bond or clamp a glass pane so tightly that using simply supported (all edge translations fixed) is realistic when large deformations are present. The validation in Ref. [13] showed that such boundary conditions heavily underestimate the deflections. Therefore, we use boundary conditions that allow the edges to slide inward. To achieve the same condition with the solid-shell elements, we can place an additional set of shell elements on the top and bottom edges of the solid-shell elements representing the support, e.g., SikaFlex bonding. Then, one can apply the geometric boundary conditions to these support elements’ “free” edges. See Fig. 10 for visualization and Fig. 11 for meshed rectangular laminated glass model. Note that even when all the translations are fixed in the supports, the laminated glass can still slide in the in-plane direction.

For the shell model, we use four-node structural SHELL181 elements. Each node has six degrees-of-freedom: three rotations and translations. The FSDT governs the element behavior. We use gens section type for the shell, which is a pre-integrated shell section type. Therefore, we must only specify and assign the $A$⁠, $B$⁠, $D$⁠, and $T$ matrices from Eq. (3) through Eq. (6) to the elements. These can be given in ansys apdl with SSPA, SSPB, SSPD, and SSPE commands, respectively. A constraint equation (ce command in ansys) constrains the nodes coinciding in the xy-coordinates. ansys does not accept a zero $A$ matrix, so small values are used instead for the stabilizing layer. Similarly, the edges of the stabilizing layer must be bounded in an in-plane direction to prevent rigid body motion or large pivots in nonlinear analysis. The edges of the laminate layers are only bound in the $z$-direction. The nodes at lines going through the laminate layer at $a/2$ and $b/2$ are bounded in the $y$-direction and $x$-direction, respectively, to prevent the rigid body motion of the laminate layer. See Fig. 12 for the meshed model and explanation of the boundary conditions. The external loads only apply to the laminate layer.

In this section, we first compare the laminated glass FE models. The deflections are plotted for interlayer shear modulus values ranging from 0.001 to 10,000. The analyses are conducted on a linear and nonlinear basis. Then, we compare the stresses between the models. For brevity, we only consider rectangular configurations: (1) $a=2500$ mm, $b=1500$ mm, $t1=t2=11$⁠, $tint=1.52$ mm, and (2) $a=b=2000$ mm, $t1=11$ mm, $t2=7$ mm, $tint=1.52$ mm. The former is a rectangle with an asymmetric thickness configuration, while the latter is a square with a symmetric thickness configuration. The applied load is 2.5 kPa and 10.0 kPa for linear and nonlinear analyses, respectively. The 2.5 kPa load is a typical design load for windows on the upper decks according to the classification rules. 10 kPa load is chosen for nonlinear analysis for demonstrating the accuracy of the proposed model also in the geometrically nonlinear region. Finally, all the numerical data used to plot the graphs are presented in Appendix.

See Fig. 13 for maximum deflection results and relative and absolute accuracies obtained from linear analysis for varying interlayer shear modulus values and 2.5 kPa applied load. For $2500×1500$ size, the relative and absolute accuracies are 0.7–2.6% and 0.03–0.11 mm, respectively. For $2000×2000$ size, the relative and absolute accuracies are 0.6–1.7% and 0.04–0.10 mm, respectively. Hence, the models agree with each other.

Next, we change the analysis type to nonlinear and the load to 10 kPa to make the large deflection effect (von Kármán strains) more pronounced. The sizes and thickness configurations are the same as shown in Fig. 13. See Fig. 14 for the deflection results and relative and absolute accuracies. For the $2500×1500$ size, the relative and absolute accuracies are $−$1.0–2.4% and 0.02–0.19 mm, respectively. For the $2000×2000$ size, the relative and absolute accuracies are $−$1.3–1.2% and 0.01–0.35 mm, respectively. Hence, the models also agree with each other on a nonlinear basis.

The maximum normal stresses for each interlayer stiffness are recorded from the 3D solid-shell model. In the 2D double-shell model, one must calculate the maximum stress from the maximum normal and bending stress resultants using Eqs. (24)–(26). Observing only the stress in the bottom glass pane in a symmetric thickness configuration is sufficient. We must observe the stress on both glass panes in an asymmetric thickness configuration. See Fig. 15 for the normal stress results and the corresponding absolute and relative accuracies for symmetric rectangular and asymmetric square laminates on a linear basis. The accuracy is generally around 1% at maximum, but there is a large spike at 45 MPa of 8.4%. This spike may be because of some numerical problem in the solution. However, the spike corresponds to 0.11 MPa in absolute accuracy, which is insignificant.

See Fig. 16 for the normal stress results and the corresponding absolute and relative accuracies for symmetric rectangular and asymmetric square laminates on a nonlinear basis. The accuracy is generally worse with very small interlayer stiffness values. The accuracy increases with higher interlayer stiffness values. The accuracy varies between $−4$% and 6%. These correspond to $−$1 MPa and 2 MPa, respectively. Considering that glass has a design strength of 40 MPa in marine structures with a safety factor of 3 or 4 (according to classification rules, e.g., [14,15]), such accuracy is acceptable.

There is a need for an efficient and readily available finite element modeling technique for laminated glass plates in marine structures. One can model the individual windows in these structures locally using full 3D solid elements (e.g., Ref. [12]), but that becomes computationally too expensive when considering large global models (e.g., Fig. 3). The global ship model can be used to identify critical locations, i.e., where the stresses or deflections of structural members are high. Because global analysis is typically linear and the mesh is coarse, these critical locations can be analyzed more locally (e.g., in partial global models or individually) in detail. For windows, simplified 2D methods exist (e.g., Ref. [8]), but their shortcomings prevent efficient analysis. Therefore, this article proposed an FE modeling concept for laminated glass with two glass panes joined by an interlayer, which is a step toward including the laminated windows efficiently and accurately in a global model. It is suitable for any FE software with a shell element and constraint equation option.

The four case studies showed that such a concept could accurately estimate deflection and stresses efficiently, even at large deflections. The accuracy in nonlinear analysis varied from $−$1.2% to 2.3% and $−$4% to 6% in terms of maximum deflection and maximum stress, respectively. These results are important as the geometric nonlinearity can have a significant positive influence on the stresses and deflections of the glass panes, which can potentially offer structural weight savings by reducing the glass pane thickness [13]. The double-shell element model has fewer degrees -of-freedom and can have coarser mesh than a full 3D solid model of laminated glass. Hence, the computational burden is reduced, which is advantageous for global FE models. However, it is still unreasonable to expect to run geometrically nonlinear analysis of the full ship model. A long-term goal is to be able to analyze, for example, larger glass structures such as the AquaDome in a nonlinear analysis (Fig. 3). This article only proved that the theory works. How one implements this in their model determines how effective the modeling is.

This article is aligned with the classification rules and marine standards in that the hermetically sealed gas-filled cavity in IGU type of windows is neglected. Hence, one only analyzes the glass panes directly in contact with the environmental loads. This approach is a simplification and results in a conservative response of the deflection and stresses [12]. In the future, we also want to implement the present modeling concept in IGUs. Furthermore, now the proposed model works well as a local model. However, the connection of the local model to the global model requires further studies.

The corresponding author was supported by the School of Engineering at Aalto University and Meyer Turku Oy. Further, the research presented in this article has received funding from Business Finland (Grant No. 3409/31/2022). The work is part of the carbon-neutral lightweight ship structures using advanced design, production, and life-cycle Services (CaNeLis) project within the scope of Climate-Neutral Cruise Ship (NEcOLEAP) and Sustainable Manufacturing Finland roadmaps. All financial support was greatly appreciated.

There are no conflicts of interest.

The authors attest that all data for this study are included in the paper.

The numerical data used to plot results are presented here in a shortened form for better readability. The tables have been reduced to 20 rows (from 140 rows) that capture the curves of the presented figures. Tables 1 and 2 relate to Fig. 13, Tables 3 and 4 relate to Fig. 14, Tables 5–8 relate to Fig. 15, and Tables 9–12 relate to Fig. 16.