Peridynamic Model for a Mindlin Plate Resting on a Winkler Elastic Foundation

In this study, a peridynamic model is presented for a Mindlin plate resting on a Winkler elastic foundation. In order to achieve static and quasi-static loading conditions, direct solution of the peridynamic equations is utilised by directly assigning inertia terms to zero rather than using widely adapted adaptive dynamic relaxation approach. The formulation is verified by comparing against a finite element solution for transverse loading condition without considering damage and comparing against a previous study for pure bending of a Mindlin plate with a central crack made of polymethyl methacrylate material having negligibly small elastic foundation stiffness. Finally, the fracture behaviour of a pre-cracked Mindlin plate rested on a Winkler foundation subjected to transverse loading representing a floating ice floe interacting with sloping structures. Similar fracture patterns observed in field observations were successfully captured by peridynamics.


Introduction
In many engineering applications including marine, civil and transport engineering, analysis of structures resting on an elastic foundation is an important problem of interest [1]. To represent the elastic foundation, Winkler and Pasternak formulations are widely utilised. In Winkler formulation, the elastic foundation is represented by the distribution of springs to resist the lateral deflection of the structure resting on the elastic foundation. On the other hand, Pasternak formulation can capture the shear interaction between springs [1].
Although there are numerous studies in the literature considering elastic foundation problem, only few of them investigated the behaviour of an existing crack inside a structure resting on an elastic foundation. Amongst these, Matysiak and Pauk [2] performed stress analysis near a crack tip in an elastic layer resting on a Winkler foundation by using the method of Fourier transforms and dual integral equations. Farjoo et al. [3] investigated rolling contact fatigue cracks in railway tracks and used a simplified finite element model (FEM) and extended finite element method (XFEM). They observed that the elastic foundation leads to an additional bending stress which increases the crack growth rate significantly. In another study, Attar et al. [1] investigated the free vibration of a shear deformable beam with multiple open edge cracks using lattice spring model. Finally, Nobili et al. [4] presented a full-field solution for the linear elastostatic problem of a homogeneous infinite Kirchhoff plate with a semi-infinite rectilinear crack resting on a two-parameter elastic foundation. They calculated stress intensity factors for both symmetric and skew-symmetric loading conditions.
In this study, an alternative approach, peridynamics [5], is used for the analysis of a Mindlin plate resting on a Winkler type elastic foundation. Peridynamics was originally introduced to overcome the limitations of classical continuum mechanics. The equations of motion in peridynamics are in the form of integro-differential equations, and they do not contain any spatial derivatives. Therefore, these equations are valid regardless of discontinuities. Peridynamics has been successfully used to analyse different material systems and geometrical configurations [6][7][8][9][10][11][12][13][14]. An extensive literature survey on peridynamics is given in Madenci and Oterkus [15] and Javili et al. [16]. Aforementioned benefits of peridynamics have attracted interest in solving solid mechanics problems particularly those involving damage and fracture. Majority of such attempts deal with full 3D models or 2D plane stress/plane strain models. There are relatively few peridynamic models considering structures resisting transverse deformation with one dimension (e.g. the thickness) significantly smaller than the other two (e.g. aircraft fuselage, ship hull, pressure vessel) including Silling and Bobaru [17] for 2D membranes, Taylor and Steigmann [18], O'Grady and Foster [19], Diyaroglu et al. [20] and Reddy et al. [21] for plates and flat shells, and Chowdhury et al. [22] for shells. This study is an extension of the Mindlin plate formulation developed by Diyaroglu et al. [20]. A similar approach was presented by Di Paola et al. [23] for non-local modelling of a beam on an elastic foundation. The current formulation is capable of analysing Mindlin plates resting on an elastic Winkler foundation with damage prediction capability. Moreover, the direct solution approach (Bobaru et al. [24], Breitenfeld et al. [25]) is presented to obtain the solution in static conditions rather than using widely adapted adaptive dynamic relaxation (ADR) scheme [26]. Finally, several verification and demonstration cases including a Mindlin plate with or without an initial crack subjected to transverse loading or pure bending loading conditions are presented to validate the current formulation and demonstrate its capabilities.

Peridynamic Theory
Peridynamic (PD) theory was introduced by Silling [5] as an alternative continuum mechanics formulation to classical continuum mechanics (CCM). As opposed to the localised concept of CCM, PD theory is based on non-local interactions between material points. Therefore, material points which are far from each other but within their interaction range, called horizon, can interact with each other (see Fig. 1). Material points, x ′ , inside the horizon, H x , of the material point, x, can be considered as the family members of the material point, x. Moreover, PD theory uses displacements rather than derivatives of displacements. Therefore, the equation of motion of a material point is expressed as an integro-differential equation and this equation is valid anywhere in the structure regardless of discontinuities such as cracks. In Eq. (1), u is the displacement vector, b is the body load, ρ is the mass density and dV x 0 denotes the volume of the material point x ′ . 'Dot' symbol represents the time derivative, and f is the peridynamic force density vector representing the force that the material point x ′ exerts on the material point x. The relative position and relative displacements of the material points x ′ and x are defined, respectively, as and In the original peridynamic formulation, i.e. bond based peridynamics, peridynamic force for an elastic and isotropic material can be expressed as where f(| ξ + η| , ξ) is a scalar-valued function which depends on the bond stretch, s, and the bond constant, c, as The bond stretch can be defined as The bond constant, c, can be specified in terms of elastic modulus, E, and horizon size δas where h is the thickness and A is the cross-sectional area. In PD theory, the material damage is included as part of the constitutive relationship by introducing a failure parameter, so that if the stretch exceeds a critical stretch value, failure parameter reduces the peridynamic force value to zero. In other words, the peridynamic bond between two initially interacting material points is broken.
In order to solve the PD equation of motion given in Eq. (1), the meshless method is widely used. Therefore, Eq. (1) can be rewritten in a discrete form as ρ :: where k is the main material point, j is the family member and N is the number of material points inside the horizon of the material point k.

Peridynamic Mindlin Plate Formulation
The peridynamic formulation presented in the previous section is for material points having translational degrees of freedom only. If rotational degrees of freedom are desired to be included to represent Mindlin plate formulation in peridynamics, appropriate changes to the original PD formulation should be made as explained in Diyaroglu et al. [20]. In Mindlin plate formulation, each material point has three degrees of freedom including transverse deflection, w and rotation of planes around x-axis, ϕ y and y-axis, ϕ x (see Fig. 2). As presented in Diyaroglu et al. [20], the transverse shear angle and curvature can be respectively expressed in peridynamic form as Fig. 2 Initial and deformed configuration of a Mindlin plate [20] and where θ is the peridynamic bond orientation with respect to x-axis. Moreover, the peridynamic equations of motion for the material point k can be derived using the principle of virtual work as ρh :: :: ρ h 3 12 :: Using transverse shear angle and curvature equations given in Eqs. (9) and (10), Eqs. (11)-(13) can be rewritten as ρh :: :: :: where and k s represents the shear correction factor. To describe mode-I and mode-III type fracture modes, Diyaroglu et al. [20] defined critical curvature and critical shear angle parameters, respectively, as where G Ic and G IIIc represent mode-I and mode-III critical energy release rates, respectively.

Direct Solution of the Peridynamic Mindlin Plate Formulation
In peridynamics, the static solution can be obtained by using adaptive dynamic relaxation (ADR) [26] or direct approach [24]. In ADR, artificial damping is introduced to the system and the solution converges to the static solution after a certain number of iterations. In the direct approach, the inertia term is specified to zero and a matrix solution is required. Therefore, the PD force function can be expressed in terms of the second-order micromodulus tensor C as [5] f In the case of PD Mindlin plate formulation, micromodulus tensor, C can be defined as a Jacobian matrix which is a matrix of all first-order partial derivatives of a vector-valued function. Therefore, for the force vector function f which is a function of shear angle φ and curvature κ, the micromodulus tensor can be defined as: where f z , m ϕ x and m ϕ y represent force or moment functions between material points arising from transverse shear deformation and bending. Utilising peridynamic equations given in Eqs. (11)-(13), force and moment functions can be obtained as Therefore, using Eq. (23) micromodulus tensor C takes the form of The force and moment functions between material points j and k can be rewritten by substituting Eqs. (9) and (10) into Eq. (28) as After reorganising Eq. (29), the following matrix expression of force and moment functions can be obtained as For static and quasi-static problems, the acceleration terms, ϕ y can be omitted from the equation of motion as

Peridynamic Mindlin Plate Resting on an Elastic Foundation
In this study, the Winkler foundation is considered as the elastic foundation and coupled with PD Mindlin formulation presented in Section 4. Winkler foundation was originally introduced by Winkler [27] for modelling the soil-structure interactions. Winkler method assumes that vertical translation of the soil, w, at a point depends only upon the contact pressure, p, acting at that point in the idealised elastic foundation and a proportionality constant, k, as The proportionality constant, k, is commonly referred to as the modulus of subgrade reaction or the coefficient of subgrade reaction. This model was first used to analyse the deflections and resultant stresses in railroad tracks. In the following years, it has been applied to many different soil/fluid-structure interaction problems, and it is known as die Winkler model (Fig. 3).
In order to combine the Winkler foundation with PD Mindlin plate matrix formulation, Winkler foundation formulation can be written in matrix form as where k is the spring stiffness and h is the thickness of the plate. It is assumed that the Winkler foundation only affects transverse deflection.

Numerical Results
As part of the numerical results, simple static loading conditions are considered first to compare the PD predictions with the finite element analysis (FEA) results using ANSYS, a commercial FEA software. Next, a plate with a central crack under pure bending resting on a Winkler foundation with very small spring stiffness is considered as a validation case to compare against results obtained in Diyaroglu et al. [20]. Then, fracture behaviour of a precracked ice sheet floating on water under transverse loading condition is investigated. . This problem was first introduced by Lu et al. [28] to simulate displacement distribution for a finite size ice floe interacting with sloping structures. As it was stated by Lu et al. [28], there is no analytical closed-form solution to calculate the deflection and stress distribution of a finite plate with free edges under evenly distributed edge pressure within a half circular area. Therefore, a numerical solution is adopted in order to verify PD results. The length of the square plate is L = 0.43 m with a thickness of h = 0.01 m. The radius of the loading area is R = 0.2 L. The Young's modulus of the plate is specified as  The peridynamic solutions for transverse displacement and rotations are compared with finite element solutions obtained by using ANSYS shell element, which is suitable for thick/thin shell structures. As depicted in Figs. 6 and 7, PD and FE solutions are in good agreement with each other and this verifies that the PD direct solution correctly captures the deformation behaviour of the Mindlin plate rested on an elastic foundation.

Pre-cracked Mindlin Plate Rested on a Winkler Foundation Subjected to Pure Bending Conditions
In the second example, a verification study is considered as in Diyaroglu et al. [20] to investigate the behaviour of a pre-existing crack in a Mindlin plate. A square plate The Young's modulus of the plate is specified as E = 3.227 GPa, and the shear modulus is G = 1.21 GPa. The distance between material points is Δx = 0.01 m. The horizon size is chosen as δ = 3.015Δx. The stiffness of the Winkler foundation is set to be a very small value, k = 10 −9 N/m, in order to represent the original example of Diyaroglu et al. [20] which is free from the elastic foundation. The material is chosen as polymethyl methacrylate (PMMA), which shows brittle fracture behaviour. Mode-I fracture toughness of this material is given as 1.33 MPa ffiffiffiffi m p [29]. In order to show simple mode-I crack growth, a bending moment loading is applied to a single row of points at horizontal boundary regions of the plate. Small increments of resultant body load of Δ e b x = 0.05 N/m are induced in order to obtain a stable crack growth. Crack starts to grow approximately at e b x = 284 N/m as shown in Fig. 9, and a similar crack pattern is obtained as in Diyaroglu et al. [20].

Pre-cracked Mindlin Plate Rested on a Winkler Foundation Subjected to Transverse Loading
This case represents a further development of the example presented in Section 6.1. Ice floe length is specified as L = 3.01 m. Load area radius representing the sloping structure load is set to R = 0.086 m. The thickness of the plate is h = 0.01 m (Fig. 4). Ice is modelled as an isotropic material with Young's modulus of E = 5.5 GPa and shear modulus of G = 2.0625 GPa. The distance between material points is Δx = 0.01935 m. The horizon size is chosen as δ = 3.015Δx. Winkler foundation stiffness k is set to k = 10055.25 Pa/m which roughly approximates water behaviour. Ice is a complex material and can exhibit either ductile or brittle fracture properties depending on the conditions [30]. For this example case, sea ice is considered a brittle material as considered in Lu et al. [28]. Mode-I fracture toughness of sea ice is given as 0.12 MPa ffiffiffiffi m p [30]. To the authors' knowledge, there is no available value for mode-III fracture toughness of sea ice in the current literature. We assumed mode-III toughness to be 7 times greater than mode-I by comparing the ratios to other brittle materials such as PMMA.
In order to generate initial damage, a small initial crack was introduced into the model. The size of the initial crack was set to 6Δx and crack orientation was perpendicular to the free edge.
Initial load is set to 0, and then small increments of resultant body load, Δ b b z = 0.1 N/m 2 , are induced in order to obtain a stable crack growth.  According to Lu et al. [28] and Nevel [31], the specified plate size in this problem case can be considered as a semi-infinite plate. According to observations in the field made by Kerr [32], the failure mechanism of a semi-infinite ice plate subjected to a force P at the free edge proceeds as follows. First, a radial crack forms, which starts under the load and propagates normal to the free boundary. This is followed by the formation of a circumferential crack that causes final failure. This behaviour is clearly captured by peridynamics and shown in Fig. 10 where Fig. 10 a shows initial crack before the plate is loaded, Fig. 10

Final Remarks
In this study, a new peridynamic model is presented for a Mindlin plate resting on a Winkler type elastic foundation. The formulation is validated by comparing against FEA results for a transverse loading condition for a plate without a crack. For a pure bending loading condition applied to a plate with a central crack free from the elastic foundation provided a similar crack pattern that was obtained in an earlier study. Finally, a pre-cracked ice sheet floating on water under transverse loading conditions was investigated. As observed in field observations, peridynamic results showed that first a radial crack forms and propagates normal to the free boundary followed by a circumferential crack.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.