Nonlocal wrinkling instabilities in bilayered systems using peridynamics

Wrinkling instabilities occur when a stiff thin film bonded to an elastic substrate undergoes compression. Regardless of the nature of compression, this phenomenon has been extensively studied through local models based on classical continuum mechanics. However, the experimental behavior is not yet fully understood and the influence of nonlocal effects remains largely unexplored. The objective of this paper is to fill this gap from a computational perspective by investigating nonlocal wrinkling instabilities in a bilayered system. Peridynamics (PD), a nonlocal continuum formulation, serves as a tool to model nonlocal material behavior. This manuscript presents a methodology to precisely predict the critical conditions by employing an eigenvalue analysis. Our results approach the local solution when the nonlocality parameter, the horizon size, approaches zero. An experimentally observed influence of the boundaries on the wave pattern is reproduced with PD simulations which suggests nonlocal material behavior as a physical origin. The results suggest that the level of nonlocality of a material model has quantitative influence on the main wrinkling characteristics, while most trends qualitatively coincide with predictions from the local analytical solution. However, a relation between the film thickness and the critical compression is revealed that is not existent in the local theory. Moreover, an approach to determine the peridynamic material parameters across a material interface is established by introducing an interface weighting factor. This paper, for the first time, shows that adding a nonlocal perspective to the analysis of bilayer wrinkling by using PD can significantly advance our understanding of the phenomenon.

and in human skin during aging [5,6]. The folding pattern in mammalian brains is evidentially linked to intelligence but also to neurological dysfunction [7,8]. The occurrence of wrinkling instabilities is harnessed by engineers, for instance to determine mechanical properties of films [9] or to pattern surfaces at the micro-and nanoscale [10,11]. More recently, the suddenly induced large strains due to buckling have been exploited to enhance the effectiveness of energy harvesters [12,13]. Geometries in which instabilities appear are an elastic half-space [14][15][16], cylindrical or spherical bodies [17,18] and bilayer or multilayer systems [19][20][21], among others. The underlying mechanism of wrinkling instabilities in bilayers is based on a stiff thin film bonded to a compliant elastic foundation. Compression in the film can initiate buckling into a sinusoidal wave pattern. The origins of the compression in the film are manifold, including external compression of the bilayer, relaxation of substrate prestretch, constrained film growth or swelling and substrate shrinkage [22][23][24][25]. Regardless whether the formation of wrinkles is desired or not, there is a need for understanding and controlling the behavior of the system. During the past decades the phenomenon has been explored in extensive analytical, experimental and numerical studies.
The analytical works of Biot [26], considering a beam on an elastic foundation, and Allen [27], concerned with wrinkles in sandwich structures, have laid the ground for closed-form expressions of the critical wrinkling conditions in a bilayer. In [28], Javili and Bakiler have derived a displacement-based solution for growth-induced instabilities which coincides with the Allen solution under certain conditions but also holds for plane strain as well as plane stress conditions. The influence of substrate prestretch has been taken into account in [19,29,30], treating the substrate as a Neo-Hookean material. Holland et al. [19] have additionally distinguished between further origins of compression. Lejeune et al. [21,31] have developed multi-layer models that are capable of modeling a weak attachment of the film to the substrate.
The analytically predicted trends for the influence of the main controlling parameters, e.g. the stiffness ratio of the film and the substrate as well as the film thickness, have been confirmed in experimental studies [6,20,[32][33][34][35][36] as well as numerical simulations [20,25,28,29,[37][38][39][40][41][42][43][44]. To date, the numerical method of choice is the finite element method (FEM). To trigger the instabilities, it is customary to introduce small perturbations to the system through the mesh, boundary conditions or material properties [39]. In order to eliminate the need for subjective perturbations that might affect the resulting wrinkling pattern, in [16,37] the FEM simulation is enhanced by an eigenvalue analysis. This approach facilitates a reliable identification of the critical conditions.
Despite the increasing body of literature, the experimentally observed wrinkling behavior is still not fully explained. Figure 1 depicts the formation of wrinkles in an elastomeric bilayer obtained in experiments by Budday et al. [20]. It can be seen that the wrinkles do not form uniformly along the whole length of the bilayer. Instead, the boundaries seem to have an impact on the buckling. This observation has been left unconsidered in the analysis with classical continuum mechanics (CCM) that is limited to a local conception of the phenomenon. Nonlocal theories, considering long-range forces instead of contact forces only, are typically able to capture such boundary effects. This serves as a motivation to investigate the nonlocal wrinkling behavior in order to shed light on nonlocal effects, e.g. the boundary effect.
Peridynamics (PD) is capable of realizing a nonlocal analysis. Introduced by Silling [45], the method was originally designed to cope with modeling damage where the use of partial differential equations would lead to singularities. To this end, spatial derivatives in the governing equations are replaced by integrals over the neighborhood of a continuum point. The integral terms comprise the interaction Fig. 1 Experimentally observed formation of wrinkles in a bilayer. The buckling is induced by applying an unstretched film on top of a prestretched substrate followed by gradually releasing the prestretch in the substrate. It is clearly noticeable that the wrinkling pattern changes towards the boundaries. The pictures are taken from experiments presented in [20] with kind permission of the authors forces acting across a finite distance and are accountable for the nonlocal character of the method. This aspect carries huge potential to advance the modeling of materials and the application of PD has therefore expanded to fields outside of damage modeling. Exemplary areas of application are multiphysics [46][47][48], multiscale modeling [49][50][51], topology optimization [52,53] and biological systems [54,55]. An extensive overview of recent publications on PD can be found in [56]. The basic version of PD, called bond-based PD, lacks the consideration of volume dilation and a material model is thus restricted to a fixed Poisson ratio [45]. Silling et al. [57] have established a state-based version of PD to overcome this shortcoming. Alternatively, Javili et al. [58] have recently presented the theory of continuum-kinematics-inspired peridynamics. For a detailed description of the computational implementation see [59]. A peridynamic material model is characterized by parameters defining the property of a bond between two points. If dissimilar materials are involved in an interaction, a combination of the bond properties is necessary to compute the interaction force. This issue has been addressed mainly in the context of modeling functionally graded materials using PD. The most common procedure is to average the parameters of the constituent materials [60][61][62]. More elaborate concepts contain weighting parameters following different approaches of weighting. Cheng et al. [63], for instance, used weight functions that are determined by the proportion of the material coefficients, while weighting variables that allow for a more flexible choice were introduced in [64,65]. In the context of modeling diffusion, Diyaroglu et al. [66] approximated the property of an interface bond by taking into account the fraction of the bond length that is associated with the respective material.
The application of nonlocal theories to the buckling analysis of a film-substrate system remains an open field of research. Recently, Ebrahimi [67] has modeled wrinkling of thin films with and without delamination with PD. To the authors' best knowledge, no research has studied the nonlocal wrinkling behavior in bilayers and its influencing parameters. Therefore, the objective of this paper is to employ PD to investigate nonlocal wrinkling instabilities in an elastic bilayer. We aim for new insights into the characteristics of the first instability pattern through considering nonlocal effects and thus to enhance the understanding of the phenomenon. The remainder of this paper is structured as follows. Section 2 provides a brief overview on PD, including an approach on how to define the material parameters at the interface between different constituents. In Sect. 3, first, we introduce the setup leading to the occurrence of wrinkling instabilities. Second, we describe the employed numerical procedure using PD that is enhanced by an eigenvalue analysis. Section 4 comprises the results obtained by the nonlocal analysis of the bifurcation problem in the form of various parameter studies. Similarities as well as discrepancies to the local theory are analyzed. Section 5 concludes the paper by summarizing the main findings.

Peridynamic framework
This section is concerned with the peridynamic formulation that is employed in this paper. In Sect. 2.1, we introduce the basic equations of bond-based PD. Section 2.2 extends the model to heterogeneous materials. In Sect. 2.3, we describe how the equations are discretized for implementation.

Fundamentals
Conceptually speaking, peridynamics incorporates some elements of molecular dynamics into a continuum mechanics framework. Based on continuum mechanics, consider a continuous body B 0 ⊂ R 3 at time t = 0 in the material configuration and its spatial counterpart B t ⊂ R 3 , see Fig. 2. A point X in the material configuration is mapped to the spatial configuration by the nonlinear deformation map y as x = y(X, t) : B 0 × R t → B t . Motivated by molecular dynamics, every point is influenced by the neighboring points within a finite distance. This region of interaction is called the peridynamic horizon H 0 ⊂ B 0 (Lagrangian perspective) with the horizon size δ denoting the radius of the spherical neighborhood of X in the material configuration. In bondbased PD, the interaction between a point X and its neighbor X is considered via the undeformed bond vector = X −X in the material configuration and the deformed bond vector ξ = x − x = y(X ) − y(X) in the spatial configuration.
PD governing equations are continuum field equations but contain a nonlocal integral term. The key governing equation is the balance of linear momentum, in which the internal force acting on a point is determined by an integral of the Fig. 2 Illustration of a continuum body B 0 in the material configuration (left) and its spatial counterpart B t (right). In PD, a point X interacts with its neighbors in a finite neighborhood H 0 defined by the horizon size δ interaction forces over the horizon. In this contribution, we consider the quasistatic case given by with the force density per volume squared p and the external force density per volume in the material configuration b ext 0 . The conservation of the angular momentum is ensured by the nonlocal form of the quasistatic balance of angular momentum which reads A harmonic potential is chosen to model the interaction energy density ψ as With the scalar-valued line measures L = | | and l = |ξ| and the stretch λ, where the term [λ − 1] describes the bond strain. Thus, the elastic material parameter C specifies the resistance against the change of length of a bond. The constitutive law of bond-based PD follows from the differentiation of the interaction energy density with respect to ξ. That is

Material parameters for bilayers
In bond-based PD elastic materials are characterized by the material parameter C. From the constitutive law (4) it follows that the constant C effectively characterizes the stiffness of a bond between two points. This is in contrast to local material parameters in classical continuum mechanics. As a consequence, peridynamic bilayers require special consideration. Across the material interface, two different materials participate in one interaction, as it is illustrated in Fig. 3. In this contribution, we follow an intuitive approach to determine Fig. 3 The neighborhood of a point in a peridynamic bilayer. A bond across the interface connects dissimilar materials the bond stiffness at the interface. Consider a bilayer consisting of material A with C A and material B with C B . The stiffness C I of a bond connecting points from material A and material B is defined as The factor α is referred to as interface weighting factor. It can be chosen as α = 0.5 in order to average C A and C B but also it allows for unequal weighting of the material constants.

Numerical implementation
To set the stage for the computational implementation of the PD framework, the governing equations are discretized. The starting point is the nonlocal balance of linear momentum (1), neglecting body forces for the sake of presentation. The lefthand side is referred to as the residual vector R. Thus, the equation to solve reads The considered body is discretized into a grid of peridynamic points P a described by X a and x a in the material and spatial configurations, respectively. The balance Eq. (6) is evaluated at each point in the sense of collocation. The point-wise contributions are assembled into the global discretized resid- The point-wise residual contains an integral term that is solved by means of numerical integration. To this end, the peridynamic grid points are used simultaneously as quadrature points. The force integral over the horizon translates into a sum of the interaction forces of a point with the neighboring points that are identified by X i and x i in the material and spatial configurations, respectively. The contributions are weighted according to the volume fraction V i of the horizon that is attributed to X i . Applying the numerical integration and inserting the constitutive law (4), the point-wise residual can be expressed as with #N specifying the number of neighbors of P a within its horizon. The discretized bond vector in the deformed configuration is denoted by

Compression-induced wrinkling instabilities in bilayered systems
Analysis of surface wrinkling is a promising application of PD. In Sect. 3.1, we introduce the film-substrate systems in which these wrinkling patterns occur when subject to specific load profiles. In Sect. 3.2, a closed-form solution from literature based on an analytical local analysis is highlighted, which we use to compare to our results. Section 3.3 describes the computational procedure we follow to explore the wrinkling instabilities using PD.

Remark 1
It is important to note that we speak of the effective strainε with the meaning of a geometric quantity measuring the overall difference between the original length and the compressed length of the film in a macroscopic sense. We do not refer to a gradient-like strain that is based on a local view and therefore not present in the PD framework. The overbar is used to emphasize this difference.

Setup
We consider the following setup. A thin stiff film is attached to a deep elastic substrate. The film is subject to a gradually increasing compression in horizontal direction. Once the effective strain in the film reaches a critical valueε crit , the flat surface loses stability and buckles to relax the compressive stresses. Since the film is not freestanding but bonded to the substrate, it is hindered from buckling into a single wave. Instead, short-wavelength buckling is the energetically favored equilibrium state of the system, which appears in form of sinusoidal surface wrinkles. In this paper, two different mechanisms inducing the compression state in the film are considered, which are schematically depicted in Fig. 4. The first mechanism is based on a substrate prestretch prior to film attachment. When the prestretch is released, the substrate relaxes to its initial length. Hence, the substrate is under tension while the film experiences compression triggering the instability. For the second mechanism, the critical conditions originate from a whole-domain compression. As the film is applied to an unstretched substrate, both components are in a state of compression at the onset of instability.

Remark 2
As established above, the computation of the interaction forces is based on the bond elongation with respect to the undeformed bond. The configuration comprising the undeformed bond is referred to as reference configuration. Note that the reference configuration is not necessarily the same for all bonds in the bilayer system. In case of substrate prestretch, the body is enlarged by an additional subdomain, i.e. the film, at t > 0. The reference configurations of bonds in the attached subdomain and bonds connecting the initial domain and the attached subdomain differ from the reference configurations of bonds within the initial domain. However, it is possible to adhere to the commonly used concept of a Lagrangian horizon. That is, the neighbor search is performed within the configuration at t = 0. This can be achieved by approximating the fictitious positions at t = 0 of points in the attached subdomain.

Analytical solution within local theory
To analytically estimate the wrinkling initiation, we consult a linear buckling analysis of a thin film adhered to an infinite half-space derived in [28] within classical continuum mechanics. We use the relations for the critical effective strain ε crit and the critical wavelength λ crit that are expressed in terms of the elastic modulus and Poisson's ratio, which coincide with the classical Allen solution [27] and are given by π t E f E s 1 3 and (8) where E f and E s are the Young's moduli of the film and the substrate, respectively, and ν s is the Poisson ratio of the substrate. Since the derivation relies on the classical beam solution to model the film, it does not account for the Poisson ratio of the film, see [68]. It is pointed out that we consider 2D plane strain which renders an incompressibility limit corresponding to ν = 1. The reduced versions result from inserting ν s = 1 3 which is the fixed Poisson ratio for 2D PD. The closed-form solution is derived for growth-induced buckling and thus assumes a stress-free substrate at the onset of wrinkling. In the literature analytical solutions exist that take into account different origins of compression [19,29,30]. However, the employed nonlinear material models are not well suited for a comparison to the bond-based PD model. Hence, for now, we restrict ourselves to the aforementioned linear analysis.

Numerical analysis using peridynamics
To numerically calculate the formation of wrinkling instabilities, a two-dimensional computational model is developed. The analysis is conducted within the framework of peridynamics using a custom-designed program implemented in MATLAB. The geometry of the model is built according to the schematics in Fig. 4, where W and H are the width and height of the substrate and t is the film thickness. The domain is discretized into a uniform grid of peridynamic points with grid spacing L. In literature, the horizon-over-grid-size is widely chosen as δ/L = 3.01 [69,70]. We adopt this ratio aiming at a good balance between computational effort and accuracy. The prescribed effective strain in horizontal direction is imposed incrementally through Dirichlet-type boundary conditions. Moreover, the domain is supported by a roller constraint at the bottom boundary preventing displacement in the vertical direction. As it is customary in PD, boundary conditions are applied to additional material layers of depth δ along the boundaries of the actual geometry. The two main steps of the numerical procedure, namely the Newton-Raphson scheme and the eigenvalue analysis, are described in the following.

Newton-Raphson scheme
We simulate the behavior of the bilayer by evaluating the deformation problem derived in Sect. 2.3 and given by To solve the system of nonlinear equations, we adopt an iterative Newton-Raphson scheme. The method is based on a linearization of the nonlinear equations that allows for an iterative determination of the vanishing residual vector. The linearization of Eq. (10) at iteration k yields with the tangent stiffness matrix K k at iteration k. For more details on the derivation of K, see [59]. Thus, the unknown deformation increment at iteration k can be computed via The deformation vector x is updated at the end of each iteration by The procedure is repeated until the norm of the residual is within a small given tolerance.

Eigenvalue analysis
The formation of wrinkles mathematically translates into a bifurcation problem. When the effective strain exceeds a critical value, the equilibrium path of the system exhibits bifurcation, meaning there exists more than one equilibrium configuration that satisfies Eq. (10). The point at which the solution path divides into different branches is called the critical point and is indicated by at least one zero eigenvalue of the system matrix. Therefore, an eigenvalue analysis serves as a useful tool to identify the critical point marking the onset of wrinkling, as it was proposed by Javili et al. [37]. We proceed as follows. At every increment of effective strain we compute the five smallest eigenvalues of the matrix K of the last converged solution to check if the critical point is reached. If so, the current effective strain is registered as the critical effective strainε crit . In addition, the eigenvector corresponding to the smallest eigenvalue is determined, giving information about the mode of the system after bifurcation. We update the deformation vector using the obtained eigenvector and the resulting deformation pattern depicts the critical wavelength λ crit .
In close proximity of the bifurcation point, the equilibrium state describing the flat film surface becomes unstable. The system strives to minimize its potential energy. Thus, small perturbations cause the structure to switch to the wrinkled state representing a stable equilibrium with lower potential energy. Due to the nonlocal nature of PD, such perturbations occur "naturally". Peridynamic points that are located at a free surface of a body are not equipped with a full horizon. However, they are described by the same material parameters as a peridynamic point in the bulk possessing the maximum number of neighbors. Therefore, the response at the surface slightly differs compared to that in the bulk. The resulting small deviations in deformation act as perturbations to the system, similar to mesh imperfections that are often manually imposed in FEM simulations. Note that the eigenvalue analysis is carried out accompanying the Newton-Raphson algorithm. If the perturbations have initiated branch switching, this is reflected in the eigenvalue analysis. The smallest eigenvalue does not necessarily pass through zero since the equilibrium configuration follows the stable post-bifurcation path. However, the eigenvalues sharply drop towards zero and reach a plateau, which from a numerical perspective identifies the critical point.

Remark 3
In PD literature, the deviating deformations at a free surface are referred to as "surface effect". This effect is often eliminated by the introduction of correction terms aiming at a good agreement with classical methods concerning the material behavior [71,72]. Note that according to our understanding of PD, a fully nonlocal model implies nonlocal boundaries. Their consequences on the predicted behavior are inherent to the mathematical formulation and we do not treat them as an error that requires correction.

Numerical study on nonlocal wrinkling characteristics
In this section, we conduct a series of numerical studies to investigate the nonlocal wrinkling behavior of the presented bilayer with regard to several influencing parameters. The aim is to demonstrate the ability of the peridynamic computational framework to model the occurring instabilities and to capture nonlocal effects that are entirely elusive within classical models. To this end, we conduct five numerical parameter studies as follows. We consider two different origins of compression. First, in Sect. 4.1 we study the case of compression with substrate prestretch and systematically vary the stiffness ratio, the horizon size, the film thickness and the interface weighting factor. Second, the case of wholedomain compression is considered in Sect. 4.2. We carry out a nonlocality study by decreasing the horizon size, which allows for a comparison to a computation using FEM.

Compression due to substrate prestretch
In the first part, we model instabilities driven by the relaxation of the prestretched substrate. The numerical model consists of a substrate of an initial width of W = 2 and initial height of H = 0.5 as well as a film of thickness t that is varied throughout the study but never exceeds 8 % of the substrate height. The substrate is prestretched by 30 %. Unless otherwise stated, an interface weighting factor α = 0.5 is used, specifying material A as the film and material B as the substrate. As it is inherent to two-dimensional bond-based PD, the Poisson ratio of all modeled materials is ν = 1/3.

Influence of the stiffness ratio
We study the influence of the stiffness ratio C r = C f /C s , with the elastic coefficients C f and C s of the film and the substrate, respectively, by assuming different values between 50 and 1000 at a constant film thickness of t = 0.01 and a constant horizon size of δ = 0.0301. The resulting wrinkling instabilities are analyzed in terms of their main features, i.e. the critical effective strain in the filmε crit and the critical wavelength λ crit . Figure 5 (left) compares the numerically determined critical effective strainε crit as a function of C r with the predictions of the local analytical solution. The data suggest two findings. First, the critical effective strain decreases with increasing stiffness ratio. This is expected as the stiffer the film, the less resistance it experiences from the compliant substrate and therefore less effective strain is required to trigger buckling. Second, our results show that the peridynamic simulations predict the same qualitative trend as the local theory. It is important to bear in mind that quantitative differences are expected since we compare a local (linear) approach to a nonlocal (nonlinear) continuum model. However, the deviations should decrease for C r → ∞ and δ → 0, which is further elaborated on in Sect. 4.1.2. Our study shows that the nonlocal model predicts a higher critical strain than the local model. It can be concluded that nonlocal models prove to be more stable than local models. This finding is in agreement with reports from literature for other cases of instabilities, e.g. shear band evolution and localization instability, as it is described in early contributions to this field, see [73] and [74], among others. Figure 5 (right) depicts the plots of the smallest eigenvalues of the system matrix over the effective film strain for a selection of stiffness ratios. The eigenvalues are computed at each increment throughout relaxation of the substrate. As the compressive effective strain increases, the eigenvalues eventually drop towards zero since the instability point is approached. With decreasing stiffness ratio this can be observed at a later stage of relaxation. Figure 6 (top left) illustrates the influence of the variation of C r on the critical wavelength λ crit . The analytical results are included for comparison. The wavelength is calculated via the number of folds N with λ crit = W /N to be able to compare it against the analytical wavelength derived for growth-induced wrinkling. In line with the theoretical trend, the wavelength grows with augmenting stiffness ratio. The discrete jumps in wavelength occur due to the prescribed displacement at the left and right boundaries that allow for either half or full folds. Thus, with a fixed computational domain of finite width the continuous course of the wavelength over the stiffness ratio cannot be reproduced numerically. Figure 6 (top right) gathers a set of representative examples of the wrinkling pattern for different stiffness ratios. The respective deformed configurations of the domain are depicted in Fig. 6 (bottom). All shown plots result from the eigenvectors corresponding to the smallest eigenvalue at the critical point that are imposed as deformation to obtain the resulting wave pattern. It can be observed in all patterns that the waves close to the boundaries are less pronounced than those in the center. A similar wrinkling behavior has been captured by experiments, as shown in Fig. 1. Our results indicate that nonlocal material behavior might be responsible for this boundary effect. Note that in nonlocal formulations the application of boundary conditions is more complex than in local formulations. We have verified that the boundary effects also occur for different boundary conditions application strategies including a variation of the material type and size of the additional boundary region.

Influence of the horizon size
To study the influence of the horizon size δ on the wrinkling characteristics, we carry out a nonlocality study. This is achieved through varying δ while fixing the number of neighbors within the horizon by adjusting L according to δ/L = 3.01. In this manner, the observed effects can be attributed to the changing horizon sizes. With this study we pursue two objectives. The first goal is to analyze how the level of nonlocality of the material model affects the critical effective strain. Second, we aim at verifying our model through a comparison with the analytical solution for the classical model. It is known from the literature that the PD theory converges to the local solution of CCM for δ → 0. Therefore, we expect the numerical results to approach the analytical solution with decreasing δ. The numerical study is conducted with four different decreasing values for δ (0.0602, 0.0301, 0.01505, 0.007525) at a constant film thickness t = 0.02. We test three different stiffness ratios C r (50,100,200). We note that we are restricted in the choice of L and therefore also δ. Since we use a uniform grid and keep the film thickness t constant, it must be possible to resolve t by L. Figure 7 illustrates the critical effective strainε crit as a function of the horizon size δ. For every modeled stiffness ratio the critical effective strain increases with increasing horizon size. It can be concluded that larger horizons, i.e. more nonlocal models, require larger effective strains to trigger the instabilities. Moreover, diminishing the horizon size Fig. 7 Critical effective strainε crit versus horizon size δ for different stiffness ratios C r leads to smaller differences between the nonlocal numerical results and the local analytical predictions. As expected, in the limit of δ → 0 the PD results approach the local solution. This confirms that our simulation model is capable of quantitatively predicting the critical conditions for surface wrinkles.
The corresponding wrinkling patterns are depicted in Fig. 8 for the example associated with C r = 50. In contrast to the perfectly sinusoidal local analytical solution, the PD solutions are noticeably influenced by boundary effects. These effects do not completely vanish for small horizon sizes. This might be explained by the high sensitivity of the instability problem to small numerical disturbances especially for small horizon sizes. However, the wavelength slightly decreases with shrinking horizon sizes, showing the trend, approaching the classical solution.

Remark 4
In contrast to a standard mesh convergence study for FEM, there are two types of convergence studies in PD due to the two independent parameters of grid spacing L and horizon size δ. One option is to fix δ while decreasing L resulting in more neighbors within the horizon of each point. The discretized nonlocal solution is expected to converge to the continuum nonlocal solution. The second option is to keep a sufficient number of neighbors constant by fixing δ/L while decreasing δ. In the limit of δ → 0 the nonlocal solution is expected to converge to the local solution. In this paper, we have conducted the second option in order to investigate the effect of nonlocality. In Fig. 7, the convergence of our PD solution to the local solution is shown. A comprehensive convergence study also following the first option has been carried out in our previous study [59], demonstrating that both types of convergence are achieved.

Remark 5
The degree of nonlocality depends on both the considered length scale as well as the material at hand. At a certain length scale, e.g. at atomistic scale, all materials behave nonlocally. Hence nonlocality is material-dependent in the sense that it becomes more significant for some materials at a given length scale than for others. Furthermore, the extent of nonlocality too can vary significantly between different materials. The horizon size is a material parameter that captures this extent. At atomistic scale, it can be compared to the cutoff radius commonly used in molecular dynamics Fig. 9 Critical effective strainε crit versus film thickness t for different stiffness ratios C r that can also vary significantly depending on the material and interatomic potential type.

Influence of the film thickness
The next parameter of interest is the film thickness t. To investigate its influence on the critical effective strain ε crit , we consider four different values of t (0.01, 0.02, 0.03, 0.04). The study is conducted for six different stiffness ratios (50, 75, 100, 125, 150, 200) with a fixed horizon size δ = 0.0301. As is shown in Fig. 9, the numerical analysis predicts a decrease of critical effective strain with increasing film thickness. It is interesting to note that this relation is not covered by the analytical solution within the local theory which renders a constant critical effective strain independent of the film thickness. This observation suggests that the influence of t onε crit is introduced by nonlocality.
In addition, it can be found from Fig. 9 that the nonlocal model renders a higher critical strain than the local solution which is in accordance with our conclusion from Sect. 4.1.1 and 4.1.2. With increasing film thickness the difference to the analytical solution shrinks. Since the PD model implies two length scale parameters, i.e. δ and t, the relative length scale δ/t is an additional measure of nonlocality. Therefore, the approach to the local solution with increasing t is explained by a decrease of δ/t resulting in a more local model.
To further clarify the relation between the critical strain and the film thickness, we carry out a nonlocality study for three different film thicknesses (0.01, 0.02, 0.04). That is, we test four different horizon sizes (0.0602, 0.0301, 0.01505, 0.007525) at a constant stiffness ratio C r = 100. We remark that since we adjust L according to δ/L = 3.01 and use uniform grids, we are limited to a maximum horizon size δ = 0.0301 for t = 0.01. Figure 10 compares the critical effective strain obtained numerically to the local analytical It can be observed that the differences between the numerical results for different film thicknesses decrease with decreasing δ and approach the local analytical result. Consequently, the data demonstrate that the influence of the film thickness on the critical effective strain can be considered a nonlocal effect. We again point out that the film thickness is the only length scale parameter inherent to the local theory. Since PD additionally includes the horizon size, an interplay of the two length scale parameters might cause the observed effect.

Influence of the interface weighting factor
In Sect. 2.2, we introduced the interface weighting factor α to determine the elastic coefficient C I across a material interface. So far, results have been obtained with α = 0.5 yielding an intuitive averaging of the material parameters. In this section, to gain a better understanding of α and its physical interpretation, we perform a set of numerical studies by considering values in the range of 0 ≤ α ≤ 1. Figure 11 summarizes the results for the film thickness t = 0.02, horizon size δ = 0.0301 and three different stiffness ratios C r (100, 200, 500). The left graph shows an increase in critical wavelength with increasing α. Since a higher stiffness of the interfacial bonds leads to a higher overall stiffness of the interface layers, this trend is in line with the conclusion in Sect. 4.1.1. The results of the local analytical solution are included in the graph for the sake of completeness. It can be observed that the deviations of the numerical to the analytical results are larger for smaller values of α. However, we point out that the deviations cannot be explained solely by the variation of the interfacial stiffness as they might also result from nonlocality as well as the dependence of the wavelength on the computational domain as described in Sect. 4.1.1. Analyzing the right graph in Fig. 11, we make three observations. First, in the proxim-ity of α = 0 increasing α causes a sharp increase of the critical effective strain. Second, once a maximum is reached, the critical effective strain continuously decreases. The latter effect can again be attributed to the higher overall stiffness of the interface region. However, the stiffness of the interface bonds also affects the strength of the attachment between the film and the substrate, resulting in the first observation. For α approaching zero, the interface bonds are too weak to fully couple the film and the substrate. Consequently, less effective strain is required for the film to overcome the resistance of the substrate and buckle. Figure 12 further illustrates this effect by means of deformation plots for different values of α. Third, the local analytical solution consistently predicts smaller values than the numerical results. This agrees with the findings from Sect. 4.1.2, where we established that a nonlocal material model yields a higher critical effective strain than a local model.

Whole-domain compression
In this section, wrinkling instabilities induced by a wholedomain compression are considered. We conduct a nonlocality study analog to Sect. 4.1.2. The aim is to investigate the effect of nonlocality for the origin of compression at hand. Furthermore, we intend to compare the obtained numerical results not only against analytical but also numerical results based on the local theory. We expect to see differences between the nonlocal PD solution and the local results that diminish with decreasing level of nonlocality. The dimensions of the numerical model are W = 3, H = 0.5 and t = 0.02. The numerical analysis in the framework of CCM is carried out with FEM enhanced by an eigenvalue analysis using MATLAB. The domain is discretized into 30000 bi-quadratic quadrilateral elements. The peridynamic study comprises four different horizon sizes δ (0.0602, 0.0301, 0.01505, 0.007525) while δ/L = 3.01 is kept constant, resulting in different degrees of nonlocality. The material parameters are chosen as C r = 100, α = 0.5 and ν = 1/3 pursuant to bond-based PD. The FEM model employs a Neo-Hookean material model. Figure 13 shows the results of the nonlocality study in terms of the critical effective strain. In accordance with the findings from Sect. 4.1.2, the graph suggests that nonlocal effects cause an increase in critical effective strain. Following our expectations, the PD results approach both the numerical and the analytical local solution. We attribute the small deviations of the FEM solution from the analytical solution to the fact that the analytical solution is based on linear elasticity. The results of this study demonstrate that our methodology using PD is capable of modeling the effect of nonlocality by varying δ while approximating the local theory in the limit of δ → 0.  Figure 14 illustrates the results of the eigenvalue analysis for the FEM simulation as well as the PD simulation with δ = 0.007525. The purpose is to elaborate on two main aspects that distinguish an eigenvalue analysis with PD from an eigenvalue analysis with FEM. First, although in both cases the smallest eigenvalues approach zero when reaching the critical effective strain, in the PD analysis none of the eigenvalues eventually hits zero. Recall that the eigenvalues correspond to the system matrix of the equilibrium state at the current increment. Since there are virtually no perturbations present in the FEM framework, the system is in an unstable equilibrium, which manifests itself in an eigenvalue passing through zero. Due to small perturbations that are naturally inherent to the PD framework, the system switches to the undulated state. Since this equilibrium state is stable, no negative eigenvalues can be observed. Second, the wrin- Fig. 14 The results of the eigenvalue analysis using FEM (top) and PD (bottom). The course of the five smallest eigenvalues over the effective strain are depicted in the left diagrams. The effective strain increments are refined in the proximity of the critical point, where the eigenvalues experience a sharp drop. The zoom boxes highlight the differences: With FEM, the smallest eigenvalue passes through zero. With PD, the smallest eigenvalue reaches a plateau slightly above zero. The plots on the right depict the wrinkling patterns resulting from the eigenvectors corresponding to the smallest eigenvalues kling pattern resulting from the FEM analysis is a perfectly sinusoidal wave. The PD analysis, however, renders a sinusoidal wave that does not have a constant amplitude. It can be observed that the amplitude decreases towards the boundaries of the domain. Therefore, in contrast to FEM, PD provides an understanding for the boundary effect witnessed in experiments. The nonlocal approach uncovers an influence of the boundaries on the wrinkling that is not accessible for local models.

Conclusion
In the past, the buckling behavior of a bilayer has been the subject of a multitude of studies mainly based on classical continuum mechanics. Here, as a first attempt to add a nonlocal perspective to the topic, we have presented a numerical study of wrinkling instabilities in a film-substrate system using peridynamics.
We have developed a computational model that is capable of predicting the main characteristics of the wrinkling pattern. The numerical procedure includes an accompanying eigenvalue analysis that allows to precisely capture the critical conditions of the instability. We have presented the results of several parameter studies considering two different origins of compression, namely substrate prestretch and whole-domain compression. Throughout the studies, we have repeatedly validated our computational model by showing that our results converge to the local results for a decreasing horizon size. The comparison of our nonlocal results regarding the influence of different controlling parameters to CCM suggest both similarities and differences. First, the influence of the stiffness ratio follows the same trend as in the local theory. The nonlocality manifests itself in quantitative differences. Second, an effect of the film thickness on the critical effective strain that is not covered by local predictions is found. Third, the eigenvalues of the stiffness matrix do not pass through zero at the critical point, as it is observed with FEM. This can be attributed to nonlocal surface effects that act as perturbations and cause the system to switch to the stable wrinkled state. Another key finding of our study is that the nonlocal model is able to explain the experimental observation of deviating wrinkling behavior at the boundaries of the bilayer. The PD simulations show weaker waves towards the boundaries, while this effect is not reproduced in FEM simulations. This suggests nonlocality as a physical origin of the experimentally found behavior. Moreover, we have established an approach to determine the elastic coefficient across a material interface. Our study of the consequences on the material response reveals a link between the stiffness of interfacial bonds and the strength of the film attachment.
In summary, this paper shows that a nonlocal analysis of instability patterns in bilayers can provide new insights into the nature of the phenomenon. PD has been found to be a suitable tool to realize the nonlocal approach. In future research, the work will be extended to an investigation of nonlocal effects on secondary wrinkling instabilities in the post-buckling regime. Furthermore, modeling materials with different levels of compressibility will be achieved by employing continuum-kinematics-inspired peridynamics [58].