Finite element model based on refined plate theories for laminated glass units

Laminated glass units exhibit complex response as a result of different mechanical behavior and properties of glass and polymer foil. We aim to develop a finite element model for elastic laminated glass plates based on the refined plate theory by Mau. For a geometrically nonlinear description of the behavior of units, each layer behaves according to the Reissner-Mindlin kinematics, complemented with membrane effects and the von K\'{a}rm\'{a}n assumptions. Nodal Lagrange multipliers enforce the compatibility of independent layers in this approach. We have derived the discretized model by the energy-minimization arguments, assuming that the unknown fields are approximated by bi-linear functions at the element level, and solved the resulting system by the Newton method with consistent linearization. We have demonstrated through verification and validation examples that the proposed formulation is reliable and accurately reproduces the behavior of laminated glass units. This study represents a first step to the development of a comprehensive, mechanics-based model for laminated glass systems that is suitable for implementation in common engineering finite element solvers.


Introduction
Laminated glass has been developed to increase safety and security of fragile thin glass sheets. When the load is transferred in both in-plane directions, the unit cannot be considered as a beam and a more accurate description is needed. Thus, the goal of this paper is an extension of our earlier finite element model of laminated glass beams  towards laminated glass plates. For this reason, the overview of available literature below will focus on laminated glass plates, whereas the discussion related to beams can be found in Zemanová et al. (2014). Let us summarize that the laminated glass units behave in an unusual manner due to several reasons. First, the polymer foil has the elastic modulus which is by orders of magnitude smaller than the elastic modulus of the glass. Such heterogeneity in elastic properties of individual layers complicate the modeling of laminated glass units. Second, the polymer foil is highly viscoelastic with a great temperature dependency. Third, geometric nonlinearity should be included in the laminated glass behavior, which is caused by the small thickness of plates compared to the other dimensions.
A number of experimental measurements, e.g. (Vallabhan et al., 1992;Behr et al., 1993;Bennison et al., 1999), studied influence of temperature and load duration on a structural behavior of a laminated glass unit since the polymer foil is a time/temperature-dependent material, whose performance determines the coupling along glass layers. The basic method to approximate the mechanical response of laminated glass units involves the utilization of one of the limiting cases: layered case as an assembly of independent glass layers and monolithic case with thickness equal to the combined thickness of glass layers and interlayers 1 determined under the assumption of geometric linearity. 2 In Behr et al. (1985), the recommendation for the choice of the limiting case can be found, which takes into account only the temperature dependence, whereas the extension of the previous experimental stress analysis by load duration influence is summarized in Behr et al. (1993). Several approaches to the more precise analytical or numerical analysis of laminated glass structures have been proposed in order to ensure a reliable and efficient design.
Single-layer approaches, useful especially in the preliminary phases of the design procedure, approximate the behavior of unit by an equivalent homogeneous plate and can be found in a draft version of European standard for glass CEN/TC 129 (2013a). Benninson et al. (2008) extended the concept of effective thickness to the case of laminated glass beams and determined the thickness of equivalent monolithic glass unit as a function of load duration and temperature history.  proposed new relations for the deflection-and stresseffective thicknesses of geometrically linear laminated glass plates using variational principles and approximation of shape function for deformation of laminated glass plate. Practical expressions for the design of laminated glass was summarized in . The effects of geometric nonlinearity are not considered in this approach although they can be significant for laminated glass units with a small value of the shear modulus of polymer foil.
Analytical methods provide closed-form solutions only for laminated beams with very specific boundary conditions and when excluding the effects of geometric nonlinearity; for plates only series-based solutions are available. Foraboschi (2012) described a mechanical behavior of laminated glass plates by a system of differential equations derived under the Kirchhoff assumptions for two thin glass plates with the same thickness. He considered only shear stresses and strains in the polymer interlayer and solve the resulting system by the Fourier series expansion for the load and displacements. His approach is valid for a simply supported rectangular laminated glass plate subjected to lateral uniformly distributed static loading. An extension of this model to other sandwich plates with stiff interlayer can be found in Foraboschi (2013). These analytical approaches are rather difficult to be generalized to units with multiple layers and are mainly useful for the verification of numerical models.
We find the numerical analysis of the discretized formulation of the problem as the most advantageous for the modeling of nonlinear behavior. The most common tool for structural analysis is the finite element method, nevertheless a few approaches based on finite difference method can be found in literature, e.g. Vallabhan et al. (1993) and Aşık (2003). The first possibility for the simulation of behavior of laminated glass units utilizes three-dimensional solid elements for individual layers and their interaction, cf. (Duser et al., 1999). This approach requires a large number of solid finite elements and therefore leads to expensive computations, since the thickness of the laminates, especially the polymer foil, is small compared to the other dimensions. To avoid fully resolved three-dimensional simulations, we used a refined plate theory in our approach and considered the layer-wise linear variation of in-plane displacements in the thickness direction. Most of the commercial finite element codes do not have such a feature, but rely on classical laminate formulations. In this framework, the laminated glass is modeled as a homogeneous unit with effective stiffness determined from its layout and standard plate (shell) elements are employed in the analysis. Because of a great mismatch of the material properties of PVB and glass, this approach is not appropriate as it is too coarse to correctly represent the inter-layer interactions.

Model assumptions
The proposed numerical model for laminated glass plates is inspired by a specific class of refined plate theories (Mau, 1973;Šejnoha, 1996;Kruis et al., 2002). In this framework, each layer is treated as a shear deformable plate with independent kinematics. Interaction between individual layers is captured by the Lagrange multipliers (with a physical meaning of nodal forces), which result from the fact that no slippage is supposed to occur between glass plate and polymer foil.
As for constitutive modeling, glass typically responds in a linear elastic manner, whereas the polymer foil behaves in a viscoelastic manner. However, in many cases even the foil is modeled as being linear elastic by means of an effective shear modulus related to temperature and duration of loading. The hypothesis that the viscoelastic material is approximated by the linear elastic one, assuming a known temperature and duration of loading, may be constraining. On the other hand, a viscoelastic model of interlayer is more difficult to manage, and will be discussed independently and in more details in forthcoming publications.
As for the effects of geometric nonlinearity, we have reported that they should not be neglected when modeling laminated glass structures, and also confirmed this statement by numerical examples in Zemanová et al. (2014), based on the fully nonlinear Reissner beam formulation. Before extending the model to the laminated plates, we compared first the response of the fully nonlinear formulation with a simpler model based on the von Kármán assumptions, in which only deflections are large whereas the in-plane displacements and rotations remain small. Since the response of these two models was found to be very similar, see Zemanová (2014)[Section 7.1], our strategy will be to accommodate the nonlinear effects due to vertical deflections only in membrane strains.
The following nomenclature is used in the text. Scalar quantities are denoted by lightface letters, e.g. a, and the bold letters are reserved for matrices, e.g. a or A. A T standardly stands for the matrix transpose and A −1 for the matrix inverse. The subscript in parentheses, e.g. a (i) , is used to emphasize that the variable a is associated with the i-th layer.

Model formulation 3.1 Kinematics
Being inspired by Pica et al. (1980), each layer is modeled as the Reissner-Mindlin plate with nonlinear membrane effects, by adopting the following kinematic assumptions: • a straight line segment, normal to the underformed mid-surface, remains straight but not necessarily normal to the deformed mid-surface, • vertical displacements do not vary along the thickness of the plate.
x y (2) Figure 1: Kinematics of laminated plate assuming the layer-wise linear variation of in-plane displacements in the thickness direction Thus, the nonzero displacement components in each layer u (i) , v (i) , and w (i) are parametrized as, see Figure 1, where i = 1, 2, 3 refers to individual layers, u 0 , v 0 , and w 0 denote the displacements of the mid-surface, ϕ x and ϕ y represent the straight line segment rotations about the corresponding axes, and z (i) is measured in the local coordinate system relatively to the mid-plane of the i-th layer. For the layer surfaces to be connected, the corresponding displacements must satisfy the compatibility conditions with i = 1, 2.
The nonzero components of the Green-Lagrange strain tensor, e.g. (Jirásek and Bažant, 2002, Section 24.1), at the i-th layer can be divided into the membrane effects and the transverse shear strains Under the von Kármán assumptions (Timoshenko and Woinowsky-Krieger, 1959), the derivatives of u (i) and v (i) with respect to all coordinates are small and these parts of strain tensor become and are expressed in terms of the mid-plane displacements, rotations, and pseudo-curvatures as Observe that the shear strain components and pseudo-curvatures are the same as in the geometrically linear case, but the membrane strain components include nonlinear contributions depending on the deflections of the layer mid-surface The matrix S and the differential operators ∂ and ∇ are provided by

Constitutive relations and specific internal forces
Having parametrized the strains with quantities defined on the mid-surface, we proceed to stresses. We assume that the normal stress in the z direction is negligible compared to the other stresses. It follows from the generalized Hooke law for isotropic bodies that the membrane stress components and the shear stress components are given by, cf. (Pica et al., 1980), where is the membrane stiffness matrix, E (i) and G (i) are Young's modulus of elasticity and shear modulus of the i-th layer, and I is the identity matrix. Notice that, since the Reissner-Mindlin kinematics were utilized, the transverse shear stresses are approximated as constant in the z (i) direction.
The specific internal forces energetically conjugate to the generalized strains (5), are obtained by suitable weighted stress averages in the thickness direction. As for the membrane strains ε (i) m , these are associated with the specific normal forces n (i) in the form The quantity corresponding to the pseudo-curvature κ (i) is the matrix of specific bending moments m (i) provided by The specific shear forces q (i) , conjugate to the shear strains γ (i) , are obtained in an analogous way but need to be corrected due to the assumed constant shear stresses by a constant k (1) = k (3) = 5/6 for top and bottom layers and k (2) = 1 for interlayer. 3 In Eqs. (11)-(13), D s denote the matrices of normal, bending, and shear stiffness of the i-th layer, respectively.
Thus in the expressions for the specific internal forces, the geometric nonlinearity due to von Kármán assumptions appears in the normal forces only, and the specific moments and shear forces remain in the same form as in the geometrically linear case.

Energy functional
The governing equations of the discretized model are derived by the energy-minimization arguments. To this purpose, we express the internal energy of the i-th layer in the form where Ω (i) denotes the mid-surface of the i-th layer. The potential energy of external loading acting at the i-th layer reads and m (i) referring to the intensity of the distributed in-plane, out-of-plane and moment loads. The total energy of the i-th layer is expressed as the sum of internal and external components in order to assemble the total energy of the system from the contributions of individual layers The true fields u then follow from the minimization of functional Π, subject to the prescribed boundary conditions and the compatibility conditions (2).

Numerical implementation 4.1 Discretization
A standard discretization is employed, in which we approximate the unknown fields by bi-linear functions at the element level, see Appendix A for explicit expressions. As a result of this step, the internal and external energies are approximated as where r (i) e is the vector of generalized nodal displacements of the e-th element and the i-th layer and f ext,e stores the generalized nodal forces resulting from the distributed loading in (15), see again Appendix A for more details.

Resulting system of governing equations
In order to obtain the nodal degrees of freedom corresponding to external loading, we rely again on variational arguments, and start with the discretized energy functional to be minimized under the kinematic constraints and compatibility conditions. This produces the Lagrangian function in the form where r stores arbitrary kinematically admissible generalized nodal displacements, and λ is the vector of admissible Lagrange multipliers (with physical meaning of forces between layers due to compatibility). The block of matrix C, implementing the tying conditions between two adjacent layers at the j-th node, now assumes the form (i = 1, 2) so that the condition Cr = 0 enforces the compatibility equations (2) discretely at all interfacial nodes. Notice that, due to the assumption that the rotations ϕ remain small, the second term in (22) is linear in r, which significantly simplifies the derivation as well as implementation of the numerical model, when compared to the formulation presented in Zemanová et al. (2014). The Karush-Kuhn-Tucker optimality conditions, e.g. (Bonnans et al., 2003, Chapter 14), associated with (22) read where r * and λ * are vectors of true nodal degrees of freedom and Lagrange multipliers. Due to the presence of quadratic terms in the normal strain (5a), the resulting system of nonlinear equations must be resolved by the Newton scheme, by iteratively searching for the correction k+1 δr to the known nodal displacements k r in the form To this purpose, we linearize the term and plug it into the optimality condition (24a). As a result, we receive a linear system in the form in which Due to the fact that each layer is treated separately, the matrices appearing in Eq. (27) possess the block structure and are obtained by the assembly of the contributions of individual elements, cf. (Bittnar anď Sejnoha, 1996) The explicit expression for the internal forces and the tangent matrices are elaborated in detail in Appendix A.

Implementation issues
The resulting version of the algorithm is shown in Algorithm 1, where the convergence of the iteration is checked by the equilibrium residual.
Algorithm 1: Conceptual implementation of the Newton method for von Kármán plate Data: initial displacement 0 r, tolerance Result: r * , λ * k ← 0, 0 λ ← 0, assemble k f int and C while ( k η > ) do assemble k K solve for ( k+1 δr, k+1 λ) from Eq. (27) k+1 r ← k r + k+1 δr assemble k+1 f int k ← k + 1 r * ← k r, λ * ← k λ Note that the finite element model of geometrically nonlinear laminated glass plate was implemented in MATLAB-based application LaPla. In the post-processing step, the stress components are evaluated from the discretized constitutive equations (9) at the element level and extrapolated to nodes using standard smoothing procedures, e.g. (Hinton and Campbell, 1974). An interested reader is referred to Zemanová (2014) for a brief description of MATLABbased program LaPla.

Examples and discussion
We validated and verified the developed finite element model, Section 5.1, and compared it against the semi-analytical and effective thickness solutions, Sections 5.2 and 5.3.

Experimental validation and numerical verification
We considered a laminated glass unit simply supported along all four edges as an example problem for the deflection and stress analysis. Vallabhan et al. (1993) performed a series of tests on such laminated units. In-plane dimensions of the unit were 1.5 m × 1.5 m and the unit was subjected to a transverse loading, gradually increasing to 6.9 kPa in 60 s. The temperature during a series of tests varied between 21 • C and 27 • C. The maximum lateral displacements and strains at several locations at the top and bottom of the plate were measured and converted to the principal stresses. Table 1 summarizes the thicknesses and material properties of individual layers. The material parameters of glass were taken from Vallabhan et al. (1993), because these values were assumed for the conversion of measured strains to stresses. The Poisson ratio for interlayer is according to draft European standard by CEN/TC 129 (2013b). The value of the shear modulus of PVB was found by matching the measured central deflection of laminated glass units with the values obtained by the numerical models. Vallabhan et al. (1993) obtained an agreement between experiment and finite-difference solver using the value of 689.5 kPa in their study, whereas we matched experimental data with a smaller value of 400 kPa. Our finite element model is stiffer due to the explicit inclusion of interlayer thickness and this fact is in accordance with the results of the fully three-dimensional finite element analysis by Duser et al. (1999), and so is the decreased optimal value of the shear modulus. Ly 2 ], on the increasing intensity of distributed load f z appears in Figures 2-5. We reproduced the experimental data from Vallabhan et al. (1993) and determined the principal stresses only from the membrane components in order to compare the principal stresses from the model directly to experimental values. The response of the finite element model was computed for two values of interlayer shear modulus, 400 kPa and 689.5 kPa. In addition, we included the prediction of the monolithic (thickness of two glass panes and interlayer together) and the layered glass plate (thickness of two independent glass panes), corresponding to the lower and upper bound for deflections, but not for the principal stresses at the center of the unit. Due to the effects of geometric nonlinearity, the location of the extreme of principal stresses is not at the center of the plate, but moves with an increasing load. Such phenomenon is best visible in Figure 3 for the minimum principal stresses, where the measured values become practically constant at the load range of 3-5 kPa and then even slightly increase. In Figure 2, we present the comparison of experimental data for the central deflection with the response of the nonlinear model of laminated glass plates. The proposed nonlinear model represents the overall behavior of the unit well. We reached the best agreement for the values of central deflections corresponding to the intensity of the distributed load between 3-5 kPa, which suggests that the constant value of the shear modulus G (2) is optimally adjusted to this interval. For lower values of f z , the experimental values of central deflections are smaller and closer to the monolithic limit, for f z > 5 kPa, the behavior of the unit is closer to the layered limit.
We attribute this phenomenon to the time-dependent behavior of the shear modulus of the interlayer under the constant temperature. The load duration, increasing with the value of the load intensity, corresponds to decreasing values of the effective shear modulus. Therefore, the effective value of the shear modulus for f z < 3 kPa is higher than the optimized one, which results in lower values of the central deflections. The effective shear modulus of PVB is lower for f z > 5 kPa, resulting in higher values of the central deflections. Thus, more accurate results can be obtained by viscoelastic analysis, but that requires the specification of the shear relaxation modulus of the interlayer, which is not available in the present case. Figure 3 shows the comparisons of the results for the maximum principal stresses at the center of the bottom surface of the bottom layer and for the minimum principal stresses at the center of the top surface of the top layer. We found excellent agreement between experimental and computed values for the maximum principal stresses and good agreement for the minimum principal stresses. Similar conclusions also follow from the comparison of extreme stresses at the top and bottom layers in a quarter of the plate, Figures 4 and 5, in which we obtained very good agreement for tensile stresses and good agreement for compression stresses between the model predictions and experimental values.
In Figure 6, we plotted the distributions of normal stresses through the thickness of the previous simply supported laminated glass plate at the center [ Lx 2 ,  The geometrically nonlinear behavior was assumed and the load intensity was set to 5 kPa. The normal stresses in the interlayer are negligible when compared with the normal stresses in the glass layers. The values of normal stresses at the mid-surfaces are generally nonzero, due to the presence of membrane stresses in the layers, and the load is equally distributed into the two glass layers, as visible from the identical slopes of the normal stresses. After validating the finite element model, we also present its verification against the finite difference solution by Aşık (2003). The input data for the problem are the same as in Vallabhan et al. (1993), except for slightly adjusted geometric parameters: the span of the plate is 1.6 m × 1.6 m and thicknesses of glass units equal to 5 mm. As shown in Figure 8, the distributions of extreme principal stresses on the top and bottom surfaces of the structure match very well. The largest difference in maximum values is around 8.5%, which can be caused by only single-digit precision of the finite difference results for the top surface (the error due to this inaccuracy can be up to 10%). The differences in maximum values for the bottom surface are under 2.5%. Such an agreement is reasonable, given that the results were obtained by different numerical method and that the contours were constructed from different interpolations. Finally, we also checked the convergence of response of the proposed model upon uniform mesh refinement. The dependence of the central deflections and the principal stress at the center and at the whole bottom layer appears in Figure 7. Both the central deflections and central extreme stresses exhibit a rapid convergence, while the accuracy of the largest principle stress found in the plate shows significant mesh-dependency. We attribute this phenomenon to the fact that, as a result of geometric nonlinearity, the extreme stresses become localized to a small area in the vicinity of the corner, Figure 8, where an adapted finite element mesh would be more suitable. Nevertheless, the discretization by 50 × 50 elements appears to be sufficiently accurate for all the three quantities of interest.

Verification against semi-analytical model
We analyzed again a simply supported square plate loaded with a uniformly distributed pressure 750 Pa. The dimensions of the plate were 3 m × 3 m, Table 2 summarize the thicknesses and material parameters of layers. The value of the shear modulus of the polymer interlayer was varied between 10 −5 MPa and 10 5 MPa. 5 -10 −5 -10 5 Poisson's ratio ν (i) [-] 0.22 0.49 The details of the semi-analytical model for a simply-supported rectangular three-layered plate can be found in Foraboschi (2012). In short, the glass plies were assumed to behave as thin Kirchhoff plates under small deflections. The interlayer was supposed to be thin, compliant and to exhibit a state of pure shear. The resulting system was solved by the Fourier series expansion. The comparison of maximum deflection and normal stress for the semi-analytical and the finite element model appears in Figure 9. We achieved an excellent agreement for deflections.
The differences in values for the semi-analytical solution and geometrical linear finite element results stay under 2.5% up to the shear modulus of 10 MPa. For larger values of the shear modulus, the deviations increase up to 10%, because the semi-analytical model in Foraboschi (2012) is not suitable for such relatively stiff interlayers, for which the model does not approach the true monolithic response. The results for geometrically linear and nonlinear finite element model are almost identical for the shear modulus of interlayer greater than 1 MPa, but the differences increase up to 40% for small values of the shear modulus (close to the layered limit). The maximum deflection at the center of plate for the geometrically nonlinear prediction ranges from 1 600 to 1 200 of the span, which implies that the effects of geometric nonlinearity can be significant in practically relevant cases.
Analogous conclusions follow for the normal stress distributions, Figure 9. The agreement of solutions is good, the maximum error for geometrically linear model is about 8.5%, which corresponds with conclusions in Foraboschi (2012). For the geometrically linear and nonlinear solution, the monolithic limit is almost the same, but the layered limit again shows important differences.

Comparison with effective thickness approach
One of the popular simplified approaches to the structural design of laminated glass units are effective thickness methods, e.g. Benninson et al. (2008) or . We compared the solution of our finite element model for laminated glass plates with the effective thickness approach using the enhanced effective thickness (EET) method of , which gives more accurate results for deflection and stress calculations for different boundary conditions than the approach of Benninson et al. (2008). We analyzed two rectangular plates loaded with a uniformly distributed pressure 750 Pa with the size 3 m × 2 m. Table 3 summarize the thicknesses and material parameters of layers. The value of shear modulus of the polymeric interlayer was varied between 0.01 MPa and 10 MPa. Four different boundary conditions were considered: four, three or two edges simply supported, and two opposite edges simply supported and one fixed, cf. Figure 10. We compared the effective thicknesses computed according to  with the values obtained from the finite element modeling, Figure 10. The effective thickness method provides the deflection-and stress-effective thickness, which is the constant thickness of the equivalent monolithic homogeneous plate that has the same maximal deflection or maximal stress as laminated glass unit under the same boundary and load conditions. The effective thickness approach does not take into account the geometric nonlinearity, but in this example, the effect of geometric nonlinearity is minor. Therefore, we considered only the geometrically linear response of the finite element model. We found the effective thicknesses for the finite element model (FEM) by matching the maximum value of deflections and normal stresses of the Figure 10: Comparison of the effective thickness obtained with the EET approach and computed from the response of the proposed model (FEM) for rectangular plates (a) with four edges simply supported, (b) with three edges simply supported, (c) with two opposite edges simply supported, (d) with two opposite edges simply supported and one fixed laminated glass unit to the response of the equivalent monolithic homogeneous plate. Due to the symmetry of the problem, we used only a quarter or a half of the plate for computations and discretized it with 40 × 60 or 40 × 120 elements per a layer, respectively.
It can be observed that the enhanced effective thickness approach gives a very good approximations for the response of rectangular laminated glass subjected to a uniformly distributed transverse loading, which in practice coincide with the finite element solution for the plate with two or three edges simply supported. On the other hand, the enhanced effective method predicts effective thicknesses slightly greater than these from the numerical simulations for the other two boundary conditions, so that it underestimates the maximal deflections and stresses.

Conclusions
The assemblies of multiple elastic elements bonded by a compliant interlayer are frequently encountered in modern constructions. The main contribution of this work is the development of a finite element model describing the laminated glass units based on the Mau refined theory for plates, which avoids fully resolved three-dimensional simulations, and its detailed verification and validation. We showed in the presented work that such an approach circumvents the limited accuracy of similar models available in most finite element systems, which employ a single set of kinematic variables and average the mechanical response through the thickness of the plate.
Although this paper was focused exclusively on the laminated glass units with two glass plates and an interlayer, the numerical models for multi-layered sandwich structures developed in this work are general and may find applications beyond glass structures. As a particular example, consider structural insulating panels consisting of a layer of polymeric foam sandwiched between two layers of structural board or laminated hybrid-glass units.
This study represents a first step to the development of a comprehensive, mechanics-based model for laminated glass systems, which is suitable for implementation in common engineering finite element solvers. As shown by the presented results, the proposed numerical methods are well-suited for the modeling of laminated glass beams and plates, mainly because of the reduction of computational cost and accurate representation of the structural member behavior. The other advantage of the proposed numerical framework is its conceptual simplicity, which will enable us to easily include additional refinements of the model, such as a more accurate description of the time/temperature-dependence of the interlayer or delamination.
where the first three quantities correspond to the degrees of freedom contributing to the midsurface membrane strains, pseudo-curvatures, and transverse shear strains, respectively, and the last matrix r (i) K,e stores the degrees of freedom relevant to the nonlinear membrane strains, Eq. (7).
The matrices of basis function are formed in the standard way The matrices needed to evaluate the generalized strains then follow directly from Eq. (5) in the form We used 2 × 2 Gauss quadrature with four integration points [x g , y g ] 4 g=1 to evaluate the normal and the bending terms, and 1 × 1 quadrature at the center for the shear terms as a remedy to shear locking. As a result, we obtained , and α g and α 0 are weight functions for four-or one-point numerical integration.
Differentiating the terms f The treatment of the "von Kármán term" is more involved, due to the nonlinear interaction of quantities r K,x,e,c (x g , y g )B (i) K,y,e T (x g , y g ) + B (i) K,y,e,c (x g , y g )B (i) K,x,e T (x g , y g )    n (i) e (x g , y g ).
Note that, for example, the matrix K (i) nK,e is localized to rows corresponding to the degrees of freedom collected in r