A meshless method for the nonlinear von Kármán plate with multiple folds of complex shape

We present a meshless discretisation method for the solution of the non-linear equations of the von Kármán plate containing folds. The plate has Mindlin–Reissner kinematics where the rotations are independent of the derivatives of the normal deflection, hence discretised with different shape functions. While in cracks displacements are discontinuous, in folds, rotations are discontinuous. To introduce a discontinuity in the rotations, we use an enriched weight function previously derived by the authors for cracks (Barbieri et al. in Int J Numer Methods Eng 90(2):177–195, 2012). With this approach, there is no need to introduce additional degrees of freedom for the folds, nor the mesh needs to follow the folding lines. Instead, the folds can be arbitrarily oriented and have endpoints either on the boundary or internal to the plate. Also, the geometry of the folds can be straight or have kinks. The results show that the method can reproduce the sharp edges of the folding lines, for various folding configurations and compare satisfactorily with analytical formulas for buckling or load-displacement curves from reference solutions.

However, it is worth noticing that most of the simulations of origami or folded structures consider plates as infinitely stiff (or mechanisms), and the deformations are essentially rigid motions [12,13] or bar-hinge models [14].
Instead, for deformable plates, one has to solve the governing equations for plates, either for small or large normal deflections w. To introduce the folding lines, one needs to enforce a discontinuity on the rotations of the plate. For this purpose, various researchers introduced internal hinge lines. Numerous works attempted analytical or semianalytical solutions for the buckling and vibrations of linear elastic plates with internal hinges. In [15], Xiang and coworkers used the Levy's method (separation of variables) with an assumed trigonometric solution and a domain decomposition method to study the buckling of rectangular plates with through-width internal hinges. The idea is to divide the domain into two parts, provide different solutions for each sub-domain, and then patch the solutions together by imposing the continuity of the displacements, moments and shear forces, but not for the rotations. The subdivision is possible because, in these examples, the folding lines are straight and intersect the boundary of the plate.
This method was successfully used for linear elastic plates. Examples include buckling loads [16], natural vibrations [17] for a first-order shear deformation plate theory. In [18], Grossi used a combination of the Ritz method with Lagrangian multipliers for anisotropic plates and also with curved internal hinge lines [19]. Nonetheless, these methods become less general and cumbersome for von Kármán plates with multiple folding lines that do not allow a domain decomposition. In these cases, numerical methods are easier to use and more versatile.
Several examples of studies of folded plates can be found in the literature, especially in Civil Engineering. It suffices to think that a folded plate is a common architectural building block, for example used for roofs or windows [20]. However, in these cases, the fold is often created as a stress-free ridge in the initial mesh. We will not consider these cases in this paper, as we focus on stress-free flat plates that develop ridge-like deformation a posteriori as the result of applied loads.
Notable examples in the meshfree literature include [21,22], where the fold geometry is introduced a priori, by joining different plates through boundary conditions. A meshfree method discretises each plate, and the folded plate is an assembly of individual plates.
In this work, instead, we present a method that can handle folding lines that are internal to the plate and does not need to consider a folded plate as an assembly of several plates. In fact, in some cases, it is not possible to separate the domain, for example for complicated folding lines, or folds having a finite length.
The sketch in Fig. 1 presents the idea of this paper: a fold is a line of discontinuity (indicated with symbol · ) for the rotation fields, like cracks are discontinuities for the displacement fields. We hereby define as infinite folds those lines extending fully across the domain and finite as those folds with extremes within the interior of the domain. In the Mindlin-Reissner kinematics, rotations are discretised with independent shape functions from the displacements. Therefore, introducing a discontinuity in the rotations will generate a folding line. For this reason, the same machinery developed for cracks applies to folding. In this paper, we use a method for cracks in 2D for meshfree methods developed by one of the authors [23].
The paper proceeds as follows: firstly, in Sect. 2 we present the von Kármán theory of plates with large deflections; then, in Sect. 3 we reformulate the problem in weak form; in Sect. 4 we introduce the meshless discretisation and the nonlinear governing equations. In Sect. 5 we show several examples of plates containing single and multiple folds. In particular, we compare the deformations of plate with and without folds. We present two verification examples for a plate with one fold, with the derivation of the analytical solutions pre-(a) A crack with a discontinuity u in displacements.

(b)
A fold with a discontinuity ϕ in rotations.

Fig. 1
Difference in kinematics between cracks and folds sented in "Appendix C". We also compare the deflection with a reference solution for a simply supported plate with uniform distributed load and for the Euler's critical buckling load. Finally, in Sect. 6 we draw the conclusions, and point out to the reader the limitations of the method and possible directions for improvement.

The von Kármán plate
We define as a plate the following cylinder Ω 0 in the reference configuration where S 0 is the mid-plane of the plate and h is the thickness (Fig. 2). Defining l = max (L, b) as a characteristic length of S 0 For thick plates, ≥ 10%. For these values of , the Mindlin-Reissner kinematics applies where u 0 , v 0 , w 0 are the mid-plane displacements and ϕ 1 , ϕ 2 are the mid-plane rotations. Differently from to the classical von Kármán theory, we consider mid-plane rotations independent from the gradients of w 0 . We then make the hypothesis of small in-plane strains and moderate rotations, which allows to retain the geometrical non-linearities only for the normal deflection w 0 [24]. The Green-Lagrange strain tensor E is, in Voigt notation, where E 0 is the Green-Lagrange strain tensor of the midplane and κ is the curvature of the mid-plane. The constitutive model for the von Kármán plate is linear isotropic elastic material, with E being the Young's modulus and ν is the Poisson ratio. The Second Piola-Kirchhoff stress is given by with D is the stiffness tensor.

Weak form of the equations
The problem is then formulated as: find the displacement where H 1 (Ω 0 ) is the space of vectorial functions in L 2 (Ω 0 ) that are square-integrable along with their first derivatives; Π [u] is the potential energy functional, defined as where U [u] is the the internal energy functional where W is the strain energy function; P [u] is the external work functional (in absence of body forces) (9) and λ is the load level; the boundary Γ 0 t ⊂ ∂Ω 0 is where the traction λt 0 is prescribed, witht 0 being a unit vector and λ the load level; Γ 0 u ⊂ ∂Ω 0 is the boundary where the displacementū is prescribed, with α being the load level for the applied displacements and γ is a penalty factor. Equation (6) can be reformulated as Integrating over the thickness, the variation (11) becomes where with S given by Eq. (5). A shear correction factor of 5/6 multiplies the transverse shear forces N 23 and N 13 to take into account the through-thickness parabolic variation of the transverse shear stress, as predicted by the three-dimensional theory, because the Mindlin-Reissner kinematics predicts a constant transverse shear stress [24]. Also, , the subscript (·) N and (·) T respectively indicate the normal and tangential direction of L 0 t in the reference configuration, while the subscript (·) Z means the through-thickness direction; the bar( ·) indicates applied displacements and rotations ; the symbols N and M indicate respectively the stress resultants and the moment resultants over the thickness oft 0 .

The discretization of the weak form
We will use a meshfree setting, namely the Reproducing Kernel Particle Method (RKPM) [25], combined with the intrinsic enrichment presented in [23,26,27] to introduce discontinuities in the rotations.
The mid-plane S 0 is discretized with a cover of N overlapping spheres Q I ⊂ S 0 of variable radii r I , I = 1, . . . , N , such that S 0 ⊂ N I =1 Q I . We call nodes the centres of these spheres X I , and we consider H as an average measure of the distance between two neighbouring nodes. Because of the overlap, H < r I , I = 1, . . . , N .
Following a Bubnov-Galerkin method, we approximate test and trial functions as a linear combination of compact support RKPM shape functions. The where U I , V I , W I , Φ 1 I , Φ 2 I are nodal coefficients (not nodal values) and φ I (X) shape functions centred in X I for the mid-plane displacements and ψ I (X) are the shape functions for the rotations.
The shape functions are computed as where the weighting function ω is defined as The vector P(X) denotes a complete basis of the subspace of polynomials of degree k, and r is the average of the radii r I . In the following, we use the empirical rule for the support radii r = k + 1 to avoid singularity of the moment matrix. The moment matrix can be inverted quickly using an iterative algorithm based on the Sherman-Morrison formula [28] which provides explicit equations for M −1 and proved to reduce sensibly the computational costs associated with Eqs. (24) and (21). The same formula (21) stands for ψ I (X).

Weight function enrichment
To introduce a fold, we need to introduce an enrichment in ψ I (X). We use an weight function enrichment, where andē is the enrichment function. The enrichment function is briefly recalled in "Appendix A" for the convenience of the reader. In Eq. (26), the enrichmentē depends on the discretization node X I in the following manner: imagine to extend the fold line to reach the boundaries of S 0 ; this imaginary line splits S 0 into two parts, a positive one S + 0 where S > 0 and a negative part S − 0 , with S being the axis perpendicular to segment s (Fig. 3). Applying function e in Eq. (44) to the nodes X I ∈ S − 0 would simply zero their weight functions. The simply supported plate as in [30] To avoid such occurrence, the weight functions are modified with the complement to 1 of ē e(X, The enrichment in Eq. (27) is applied only to those nodes whose support is cut by the folding line: as showed in Fig. 4, where d is defined in Eq. (43).
Let us now assume that there exist n folding lines, denoted by s 1 , s 2 , . . . s n . Letē 1 ,ē 2 , . . .ē n be the respective enrichments. Then, these enrichments apply iteratively to the ω, as shown in Fig. 5: where ω 0 e = ω. The greatest advantage of the proposed enrichment is that Eq. (29) applies also for intersecting folding lines. However, care must be taken in choosing an appropriate r I for each node, or an appropriate nodal spacing H , to avoid singularity of the moment matrix (Eq. (24)).

Discretised equations of motion
Substituting Eqs. (16), (17), (18), (19) and (20) into Eqs. (13) and (15), and eliminating the variation of the nodal values, we arrive at where F (i) is the internal force vector and β and where N (·,·) are the applied in-plane normal forces per unitlength, M (·,·) are the applied moments per unit-length and Q (·) are the applied shear forces per unit-length. The subscripts refer to the X and Y axis in the reference configuration and n X , n Y are the components of the normal unit vector of curve L 0 t . Finally, Whenever the load level λ or α is not provided, such as in tracing bifurcation diagrams, one more equation is necessary.
The closure of Eq. (30) is the arc-length constraint The simply supported plate as in [30] with fourfolds where l is the arc-length parameter, and Δd = d − d 0 , where d 0 is a previously converged state. We solve the nonlinear equations with a quasi-Newton trust-region algorithm for sparse problems [29], with a provided sparsity pattern for the Jacobian to speed up the iterations. Both Eqs. (30) and (37) are normalized, so that the error in each iteration is expressed as a percentage. The convergence is usually achieved within machine precision in 4 or 5 iterations.

Simply supported plate without and with folds
We first solved a classic benchmark example (Fig. 6) for the von Kármán plate, with the aim of identifying the optimal number of nodes and the order of the basis in Eq. (23). All the quantities are appropriately normalised to obtain a dimensionless problem: L = 1 m, E = 1 Pa, ν = 0.316, h/L = 0.1. All four edges are simply supported, u 0 = v 0 = w 0 = 0. A normal uniform pressure p is applied on the mid-plane.
The direction of the load is fixed, i.e. it does not follow the deformed normal direction of the plate. "Appendix B" reports a detailed analysis on the number of nodes, order of quadrature, number of quadrature cells and order of polynomial basis for this problem. For simplicity, we used structured discretisations (Fig. 8), although the method is applicable also to unstructured discretisations. In this section we report the main findings. It turns out that a good convergence of 15 nodes per side of the plate is obtained for various thickness-to-length ratios if a full thirdorder polynomial basis is used. Contrarily to finite elements, the use of such a high order basis can be easily obtained with an RKPM approximation. In addition, we use 60 first-order rule Gaussian quadrature cells per length. Figure 7 shows a comparison with a semi-analytical (yet approximate) solution derived from a Fourier series expansion of the vertical deflection w 0 [30]. A third-order basis also seems to alleviate shear-locking issues as the thickness of the plate becomes small. Figure 7 shows no degradation of the solution even for h/L = 1 %. Therefore, in the following examples, unless otherwise stated, we will use a third-order basis and the same nodal spacing. To better capture the kinematics of a plate containing multiple folds and to avoid singularities of the moments matrix, as a precaution, we use 25 nodes per side of the plate and 100 quadrature cells per side of the plate (Fig. 8).
Next, we keep the same plate and boundary conditions, but we add fourfolds as in Fig. 9. The effect of the folds is clearly visible in Fig. 10: for the same load p, the vertical deflection is smaller than the one for the plate with no folds. This stiffening agrees with the intuition: a flat thin sheet of paper will bend under its own weight. But, if the sheet of paper contains folds, it will bend less (or not bend at all) (Fig. 11).

Verification of the method
We verified the method against analytical solutions obtained from the Euler-Bernoulli beam Eq. (45). We assumed that a fold exists in the middle of the beam. We considered two cases: an unloaded simply supported beam and a beam loaded only by a bending moment. The first one is with the beam at its rest state, and consists in a roof-like deformation (Eq. (50) and Fig. 12a) with no curvature; the second solution (Eq. (53) and Fig. 12b) has a uniform non-zero curvature with a cusp at the fold line. For both cases, we computed the L 2 norm of the error: with w H 0 given by Eq. (18) and w being the exact solution given in "Appendix C". We also evaluated the error for the rotations with ϕ H 1 given by Eq. (19) and θ is the exact solution given in "Appendix C".
The order of reproducibility of the shape functions is 1 and the order of quadrature is 1. The number of quadrature cells per shorter side of the plate is always 4 times the number of nodes (per shorter side of the plate). We kept the same H in both axes: for example, if n is the number of nodes on the shorter side (in this case b, y-axis), the number of nodes in the x-axis is L/b n.
All the quantities are appropriately normalised to obtain a dimensionless problem: The dihedral angle of the fold (the amplitude of the discontinuity in the rotations) was set at ϕ 0 = π/10. Figure 13 shows the percentage error for different nodal spacing H : even with a small number of nodes (3) along the y-axis (H = b/2) the percentage error is 0.1965% for the displacements and 0.1699% for the rotations. Figure 14 shows the error in the displacements and rotations for the second case (Eq. (53)) considered for the verification. Also in this case, the errors are very small: for

Rectangular plate with folds in the middle
Next, we show more explicitly the effect of the enrichment in a plate by inserting one or more folds in the middle. We construct an ad-hoc example as in Fig. 15.  Figure 16a shows the deformation of the plate without a fold: the rotation field ϕ 1 (X , Y ) is smooth and continuous, and the deformation shows no kinks. Instead, with a fold in the middle, the rotation field is discontinuous (Fig. 16b) and the deformation shows a ridge. The enrichment allows the introduction of multiple and kinked folds, and Fig. 17 show the same example with two inclined parallel infinite folds separated by L/8 and a central fold composed by two lines with a kink in the middle.

Buckling of plates with folds
This example is the plate shown in Fig. 18. The plate contains fourfolds, one is a central horizontal fold running through the length, and the remaining three are vertical, long half the width of the plate, staggered and separated by L/4. On these last threefolds, the transverse deflection is fixed. The plate contains an initial imperfection on the right and left edges, as an initial applied rotation α = 5 • . A compressive load λ is applied on the right edge. The presence of the folds increases the critical Euler's load, since the effective length is the distance between the vertical fold lines L/4 P c = π 2 E I with I being the minimum area moment of inertia of the cross section. Figure 19 shows the load-displacement curve, with the load approaching asymptotically the critical load. The deformation stages are displayed in Fig. 20, where the fold lines create a shiny reflection in the deformed mid-plane surface.

Conclusions, limitations and future directions
This paper demonstrated a meshless approach to the simulations of folds in nonlinear von Kármán plates. The folding lines in this paper can be infinite (extending throughout the mid-plane) or finite (with end-points internal to the mid-plane) folds. We showed examples of plates containing multiple folds of straight or kinked shape. The main conclusion of the article is that with the same methods developed for cracks, it is possible to simulate folding of plates. In fact, folds are discontinuities in the rotations, just like cracks are discontinuities in the displacements. Therefore, if one has a method for reproducing discontinuous shape functions, it is possible to use such method to discretise the rotation fields.
For this reason, we believe that this paper can bridge the fields of computational fracture mechanics and computational structural instabilities. Imagine all the vast literature of numerical methods for fracture mechanics that can now be used to simulate folding of structures, which includes the realm of origami mechanics.
This approach, of course, is possible if the shape functions approximating the rotations are different from the ones used for displacements. This is the case of thick plates, which possess Mindlin-Reissner kinematics. For thin plates, therefore, a possible future direction could be the introduction of discontinuities directly into the derivatives of the shape functions. Examples include intrinsic [37] or extrinsic enrichments [38,39].
Using Mindlin-Reissner kinematics for thin plates might generate shear locking issues. However, owing to the merits of the RKPM approximation, the use of a higher-order reproducibility allowed good accuracy even to thin plates, with aspect ratio h/L = 0.01. Of course, we do not claim that this is a solution to the shear locking issue: a more thorough analysis of this problem is outside the scopes of the paper. We refer to the works of [40][41][42] for a detailed analysis of the shear locking issue in meshfree methods.

A Enrichment function
Consider a segment s ⊂ S 0 with endpoints (X 1 , Y 1 ) and (X 2 , Y 2 ), midpoint (X m , Y m ) and length L s . Define axis T collinear with s and S perpendicular to it, with origin in (X m , Y m ) (Fig. 21). The one-dimensional signed distance function of the endpoints along T is The positive part of d (1) s (T ) is The distance function of s is and showed in Fig. 22a. The enrichment function is the normal derivative of the distance function e = 1 2 and it is showed in Fig. 22b.

B Analysis on the number of nodes, quadrature rules and reproducibility order
As anticipated in Sect. 5.1, a good comparison with a semianalytical solution is obtained when using 15 nodes per side of the plate and a first-order quadrature rule for various thickness-to-length ratios, if a full third-order polynomial basis is used. In fact, Figs. 23, 24, 25 and 26 show the degradation of the accuracy with decreasing thickness ratio h/L for RKPM with linear reproducibility order. Figures 27 and 28 show instead that, for small thickness ratios, for example h/L = 0.01, quadratic and cubic order converge faster with the increase in number of nodes: in particular, the cubic order reaches convergence even with 15 nodes per side of the plate. A similar trend is seen for higher thickness ratios, although they are not reported here for the sake of conciseness.     Finally, Fig. 29 shows that when using RKPM with a third order reproducibility, a Gaussian quadrature rule of order 1 is sufficient, provided that the number of cells per side of the plate is 4 times the number of nodes per side of the plate.

C Analytical solutions of the Euler-Bernoulli beam with a fold
Let us consider the Euler-Bernoulli beam where (·) = d/dX . To ease the notation, we removed the bar (·) and assumed that all the quantities are dimensionless. For example, where θ is the rotation, M is the moment, S is the shear force. Let us assume that a fold exists at X = X 0 and consider two cases, shown in Fig. 30.
For an unloaded simply supported beam (Fig. 30a) with a fold of dihedral angle ϕ 0 , the boundary conditions are w(0) = 0 (49a) w (0) = ϕ 0 /2 (49b) The solution of Eq. (47) is and rotation which is a roof-like deformation (Fig. 31a), with H(X ) being the Heaviside function and X = X + |X | 2 (52) the ramp function. For a simply supported beam with an applied moment M 0 at both ends and X 0 = 1/2, the solution of Eq. (47) is shown in Fig. 31b. The discontinuity in the rotations at X = 1/2 has magnitude