Peridynamic model for visco-hyperelastic material deformation in different strain rates

This study presents a peridynamic (PD) constitutive model for visco-hyperelastic materials under homogenous deformation. The constitutive visco-hyperelastic model is developed in terms of Yeoh strain energy density function and Prony series. The material parameters in the model are identified by optimizing the classical stress–strain relation and tension test data for different strain rates. The peridynamic visco-hyperelastic force density function is proposed in terms of the peridynamic integral and the Yeoh strain energy density. The time-dependent behaviour for different strain rates is captured by numerical time integration representing the material parameters. The explicit form of peridynamic equation of motion is then constructed to analyse the deformation of visco-hyperelastic membranes. The numerical results concern the deformation and damage prediction for a polyurea membrane and membrane-type acoustic metamaterial with inclusions under homogenous loading. Different surface defects are considered in the simulation. The peridynamic predictions are verified by comparing with finite element analysis results.

Constitutive models for visco-hyperelastic materials have been recently developed and can be divided into quasi-static and time-dependent models [7].
For quasi-static cases, the visco-hyperelastic materials exhibit hyperelastic behaviour. Hyperelastic behaviour is represented by empirical continuum approaches [8,9]. These approaches require a vast amount of material parameters which usually leads to several deformation tests implemented for parameter identification. The physical interpretation of these parameters is often unclear [3]. On the other hand, some models are based on statistics of molecular chains network which are also been proposed with idealized assumptions [10,11]. Moreover, Eremeyev and Naumenko [12] developed a relationship between effective work adhesion and peel force for thin hyperelastic films undergoing large deformation.
Time-dependent models associated with viscoelasticity phenomenon have been mostly formulated by hereditary integral approaches or experimental methods. In these models, the effect of strain history on the current stress state is reflected by fading-memory function by considering different strain rates [13][14][15]. Some other models are based on the constitutive relation between strain energy density and strain invariants to include the effect of time-dependent behaviour. The choice for appropriate strain energy function is determined by material parameters identification [3]. By using these models, the strain rate dependency can be considered to predict the large deformation behaviour of visco-hyperelastic materials [16].
However, the classical continuum mechanics (CCM) for visco-hyperelastic materials faces mathematical challenges in the presence of defects. This is due to the fact that the equation of motion in CCM, including wellknown finite element method [17,18], involves spatial derivatives for displacement components. In addition to CCM, there exist variational approaches with and without considering damage [19][20][21]. As an alternative approach, the peridynamic theory, a new continuum mechanics formulation introduced by Silling [22,23], removes the aforementioned challenges of CCM using the integral representation of internal forces. Therefore, it is suitable for damage prediction in structures. According to dell'Isola et al. [24,25], the origins of peridynamics go back to Gabrio Piola, who has made many contributions in mathematics and physics, especially nonlocal continuum theory, including nonlocal internal interactions, strong form of variational principle and higher gradient mechanics. Besides, PD theory is not limited to linear elastic material [26]. Madenci [27] provided the PD integrals for strain invariants to investigate the PD form of the strain energy density for both linear elastic and hyperelastic materials. It was concluded that the deformation response based on new strain energy density function is similar to that of bond-based peridynamics. Bang and Madenci [28] developed the PD strain energy density form for neo-Hookean-type hyperelastic membrane under different loading conditions. Madenci and Oterkus [29] proposed the ordinary state-based PD constitutive model for viscoelastic deformation in terms of Prony series. They captured the relaxation behaviour of viscoelastic material under mechanical and thermal loads. Madenci and Oterkus [30] and Mitchell [31] presented the ordinary state-based PD model for plastic deformation. Taylor [32] and Forster et al. [33] also considered the PD model for viscoplastic materials.
In this study, PD model for visco-hyperelastic materials is developed. This study first presents a Yeoh-type visco-hyperelastic constitutive model which works in a broad range of strain rates. The hyperelastic part of Cauchy stress captures the quasi-static behaviour using the Yeoh strain energy density, while the viscoelastic part of Cauchy stress represents the rate-dependent behaviour by using hereditary integral. This constitutive model is verified by experimental data for high-damping rubber [34] and polyurea [35]. Secondly, based on PD integrals for the first strain invariant of right Cauchy-Green strain tensor [27], the derivation of PD force density is presented by using Yeoh visco-hyperelastic material parameters and Cauchy stress definition. In order to capture the time-dependent properties of viscoelastic materials, the PD viscoelastic force density is described in terms of Prony series. The developed PD model is verified by comparing with finite element predictions. The developed PD model is used to simulate polyurea membrane with a defect in the form of a hole and a pre-existing crack. Finally, the PD model is used to predict membrane-type acoustic metamaterial morphology.

Visco-hyperelastic constitutive model
A visco-hyperelastic material can be described by an elastic spring A in parallel with a Maxwell element B as shown in Fig. 1. Therefore, stress can be separated into two parts as (1) Fig. 1 Rheological model of the visco-hyperelastic material [5] Meanwhile, there is a common assumption that captures the effects of both deformation and time in viscohyperelastic constitutive model [36] σ(ε, t) = σ 0 (ε)g(ε, t) where σ and ε are stress and strain tensor, respectively. The parameter σ 0 (ε) represents the instantaneous stress-strain relationship of the material. The function g(ε, t) is generally expressed by Prony series which varies during the deformation [5] where θ i 's are relaxation times corresponding to time dependent behaviour and g ∞ and g i are dimensionless constants which reflect the contribution of equilibrium and viscous response, respectively. Those constants should satisfy the following conditions [5]: For the time-dependent material behaviour, a visco-hyperelastic constitutive relation can be expressed as When the deformation loading is quite slow (t = ∞) which is close to the equilibrium response of the material, the spring A is the only source of the elastic stress component which becomes When the loading is quite fast with high strain rate (t = 0), the damper does not have enough time to adapt with the applied deformation. Thus, the deformation is endured by adjoining the springs. The instantaneous response can be represented by the spring A and the adjoining spring in part B as For a visco-hyperelastic material, the stress components can be represented as

Hyperelasticity
For an isotropic incompressible hyperelastic material, the Cauchy stress tensor can be expressed as [37] where p h is the hydrostatic pressure for incompressible material, I is the second-order identity tensor, B is the left Cauchy-Green deformation tensor with invariants of in which H is the displacement gradient tensor and the symbol "tr" represents trace of the matrix. W represents the strain energy density function which can be defined by the invariants of B.
In the view of the compressible properties of the rubber, the strain energy density function of Yeoh model has been proposed to investigate the rubbery mechanism which provides a better fit to the experimental data [38] as where N represents the number of terms in the model, nd represents the number of dimensions and c n0 represents material parameters. The 3-D Yeoh model can be written as and the 2-D Yeoh model can be written as In this study, Yeoh strain energy density is represented by using two terms, N = 2. For 3-D uniaxial tension condition, the deformation gradient tensor F and the left Cauchy-Green deformation tensor B can be expressed as where λ represents the principle stretch. The Cauchy stress for the equilibrium response of hyperelastic material (part A) can be obtained by using Eqs. (9) and (11b) as

Visco-hyperelasticity
Considering the time-dependent behaviour in viscoelastic relaxation tests, the viscoelastic constitutive relation is generally expressed in the form of a convolution integral [5] in which t represents time and τ represents time integral variable. The viscous part of Cauchy stress with hydrostatic pressure p v becomes [5] σ by using the 3-D form of Yeoh model for viscous strain energy density function with N = 2 Considering one-term Prony Series for viscous part by using i = N = 1, a single relaxation time can be implemented in Eq. (16). Single relaxation time can be related to strain rate as [16] [5]. Therefore, Eq. (18) becomes The viscous part of Cauchy stress becomes the following formula by plugging Eqs. (16) and (17) into Eq. (15) The viscous part of Cauchy stress for uniaxial loading can be expressed by using the deformation tensor for uniaxial loading provided in Eq. (12) The viscous strain energy density function in Eq. (17) for i = 1 can be written as and By plugging Eq. (22) in Eq. (21) the viscous part of the Cauchy stress components can be written as Cauchy stress definition for visco-hyperelastic material can be obtained by summing Eqs. (13) and (24) as By setting uniaxial loading condition, σ 22 = σ 33 = 0, Eq. (25b) becomes The hydrostatic pressure can be obtained from Eq. (26) as By plugging Eq. (27) into Eq. (25a), the stress component in the loading direction for visco-hyperelastic material can be obtained as Equation (28) can be integrated in terms of deformation instead of time for a constant strain rate,ε 11 . The principle stretch at time t can be expressed as λ =ε 11 t + 1. An infinitesimal deformation variable ζ which satisfies ζ =ε 11 τ +1 where 1 ≤ ζ ≤ λ is introduced to substitute the integral variable τ in Eq. (28). Therefore, the relation between λ and τ becomes: Equation (28) can be rewritten as Eq. (29b) In this case, the integral variable can be replaced by using dζ = d (ε 11 τ + 1) =ε 11 dτ . The upper limit of the integral can be changed to λ while the lower limit can be changed to 1. Finally, Eq. (29b) can be rewritten as Here, the magnitude of time rate of the right Cauchy-Green deformation tensor ( Ċ ) can be chosen asε [5] For uniaxial deformation with a constant strain rate, Cauchy-Green deformation tensor becomes and the absolute value for the first-order derivative of Cauchy-Green deformation tensor can be expressed as The relaxation time in Eq. (19) can be rewritten by using Eqs. (31)(32) as

Identification of constitutive model parameters
The data of uniaxial compression test on high-damping rubber (HDR) material presented by Amin et al. [34] are used for parameter identification. The data include stress-strain relations for equilibrium and low-to-moderate strain rate cases. The experimental equilibrium case is used to determine the material parameters associated with spring element A as c 10 and c 20 . The experimental nonequilibrium case with strain rate 0.24 s −1 is used to determine the material parameters associated with the Maxwell element B asc 1 10 ,c 1 20 , ϕ ref and η. The Table 1 Predicted visco-hyperelastic parameters from experimental data by Amin et al. [34] Hyperelastic parameters Viscoelastic parameters  parameter identification is described in Appendix A. Table 1 shows the material parameters from monotonic compression test for HDR material. Figure 2 represents the comparison of predicted and experimental stretch and stress relations in different strain rates for HDR material. As it can be seen from the figure, there is a good agreement between the proposed visco-hyperelastic model results and experimental data. Amin et al. [34] reported that the behaviour of HDR did not show significant change for strain rates higher than 0.88 s −1 .
The data of uniaxial compression test from low-to-moderate strain rates performed by Roland et al. [35] are used for parameter identification for polyurea material. The lowest strain rate, 0.15 s −1 , is used to determine the material parameters associated with spring A, c 10 and c 20 . The experimental nonequilibrium case with strain rate 14 s −1 is used to determine the material parameters associated with Maxwell element B,c 10 ,c 20 , ϕ ref and η. Table 2 shows the material parameters for polyurea material. Figure 3 represents the comparisons of predicted and experimental stretch and stress relations in different strain rates for polyurea material. As it can be seen from the figure, there is a good agreement between the proposed model results and experimental data even in higher strain rates, illustrating the applicability of the proposed model for the wide broadband of strain rates.

Peridynamic theory
Peridynamic theory concerns the nonlocal representation of a physical field at a material point. Each material point has nonlocal interaction with other material points in a certain horizon, and the strength of interaction is specified by a weight function. As derived by Siling et al. [23], peridynamic form of equation of motion at any instant time t can be defined as [39] ρü ( where ρ is the density of the material, V is the volume of material point x and b is the body force density vector.
According to the PD theory, the displacement field at material point x is influenced by the collective interaction of all the family members x located in its interaction domain, H x , which is defined by its horizon, δ, as shown in Fig. 4 [22].
The spatial integral in Eq. (34) can be evaluated by using a meshless scheme and can be constructed in a discretized form by replacing the integration with a finite summation as [40] (35) where N denotes the number of family members of material point x (k) and V ( j) represents the volume of each material point x ( j) .The parameters t (k)( j) and t ( j)(k) are the force density vectors being parallel to the relative position vector in the deformed state as shown in Fig. 4. As derived by Madenci and Oterkus [40], determination of force density vectors requires the explicit PD form of strain energy density function. If the PD form of the strain energy density function is known, the PD force density vector can be obtained as Force density and kinematics of PD material points [39] According to Madenci [27], the displacement gradient tensor for a homogeneous deformation can be expressed as where the symbol ⊗ represents the dyadic product, I represents the identity matrix and w x − x represents the weight function. The weight function is chosen as The parameter m in Eq. (37) is defined as The parameter m can be evaluated as m = V δ 2 with V representing the volume as V = 4πδ 3 /3 for 3D, V = π hδ 2 for 2D and V = 2 Aδ 2 for 1D [27]. The parameter h denotes the thickness of the membrane, and A denotes the cross-sectional area of the bar. Similarly, the deformation gradient tensor can be obtained as [27] in which the parameter m is specified as m = 4πδ 5 /3 for 3-D, m = π hδ 4 for 2-D and m = 2 Aδ 4 for 1-D models.
The PD form of the first invariant I 1 can be represented as [27] 4 Visco-hyperelastic Peridynamic formula

Hyperelastic response based on Yeoh model
The force density-stretch relation for a material point is defined [40] Equation (42) can also be rewritten by considering the PD form of the first invariant as By using Eq. (41), it can be obtained as By using Eq. (11a) and according to the hyperelastic membrane with two-term Yeoh model presented in Sect. 2.1 with N = 2, it can be obtained as Therefore, Eq. (43) becomes

Visco-hyperelastic response based on Yeoh model
The correspondence principle is used to establish the force density for a linear visco-hyperelastic isotropic material. If the body experiences a viscoelastic deformation, the viscous force density then can be expressed in terms of time integral as Eqs. (14) and (43) t in which the subscript i represents the terms of Prony series and the parameter θ i represents the relaxation time of i th Prony series. W v i is the Yeoh model for viscous strain energy density given in Eq. (17). In accordance with Eq. (45), the viscous property can be expressed by two-term Yeoh model for strain energy density as For one-term Prony series, Eq. (49) can be simplified as As for homogenous deformation of a 2-D membrane by using nd = 2, tr(I) = 2 and parameter m = π hδ 4 for 2-D model, the viscous force density in Eq. (50) can be rewritten as in which the first invariant I 1 can be calculated by Eq. (10) with the homogenous deformation loading λ where the corresponding classical deformation gradient tensor is The instantaneous peridynamic stretch λ ins can be expressed as Therefore, the time integral in Eq. (51) then can be rewritten as Defining dτ as the viscous stretch which represents the interaction between material points x and x, the viscous stretch is discretized in order to perform the numerical time integral. Therefore, the time integral in Eq. (54) can be rewritten as For the time step t n+1 , the viscous stretch can be written as Under the assumption that λ ins varies linearly with time the first term and the second term of Eq. (56) can be evaluated as With the Euler approximation and t θ → 0 the viscous stretch γ (t n+1 ) subsequently can be represented as Considering the single relaxation time (one-term Prony series) for viscous part, θ can be represented by Eq.
The total bond force density for visco-hyperelastic membrane can be expressed by using Eqs. (47) and (61) as In PD theory, the failure prediction is relatively simpler with respect to traditional methods which is mainly based on bond breakage. The interaction between two material points x (k) and x ( j) can be permanently terminated, leading to a modification of PD equation of motion by introducing the parameter μ (k)( j) [41] as where In the case of visco-hyperelastic deformation, the force density between two material points, x (k) and x ( j) , has time-dependent behaviour. Thus, the visco-hyperelastic micropotential between these points can be expressed by the summation of the instantaneous hyperelastic and viscous components of the micropotential. Here, 0 + represents the initial time of deformation, t 0 + (k)( j) represents the instantaneous hyperelastic force density and t (k)( j) represents the visco-hyperelastic force density. The micropotential w (k)( j) between two material points x (k) and x ( j) can be calculated as [29,30] in which 0 + represents the initial moment of the deformation, t 0 + (k)( j) represents the instantaneous hyperelastic force density and the stretch is defined as The energy release rate G (k)( j) for each interaction can be related to the micropotential which is required to eliminate the interaction between two material points x (k) and x ( j) that are across the unit crack surface, A = h x, with x denoting the material points spacing and h denoting the thickness as [29,30] The critical energy release rate for each material point is assumed to be distributed equally to each interaction of this material point. The critical energy release rate for each interaction G c (k)( j) can be expressed as [29,30] in which N c represents the number of interactions creating a unit crack surface. For 2-D model with a horizon size of δ = 3.015 x, the total number of interactions passing through a unit crack surface can be calculated as N c = 36 [42,43]. Therefore, the failure criteria for each interaction can be expressed as

Numerical procedure
In this study, the solution of PD equation of motion requires spatial integration which is performed by using a meshfree method. Adaptive dynamic relaxation (ADR) method is implemented for quasi-static analysis [44]. Euler explicit time integration scheme is used for dynamic analysis. The solution concerning calculation of viscoelastic term can be achieved by the following procedure: (1) At the load step [n + 1], the viscous stretch between material points x (k) and x ( j) can be calculated by using Eq. (60) as in whichc 1 10 andc 1 20 represent material parameters for viscous property and I 1 can be calculated for homogenous deformation loading λ by using Eq. (52). (3) Calculate the weight function between material points x (k) and x ( j) by using Eq. (38) as: (4) Calculate the viscous force density between material points x (k) and x ( j) for the load step [n + 1] by using Eqs. (51), (71) and (72) as in which m = π hδ 4 for 2-D [27]. (5) The corresponding hyperelastic force density between material points x (k) and x ( j) can be calculated from Eq. (47) as in which I 1 can be calculated for homogenous deformation loading λ by using Eq. (52) and c 10 , c 20 represent material parameters for hyperelastic property.  (6) The total force density between material points x (k) and x ( j) for the load step [n + 1] can be calculated as The micropotential w (k)( j) for the load step [n + 1] between two material points x (k) and x ( j) can be calculated by using the force density from Eq. (64) as leading to the failure parameter for the interaction between two material points x (k) and x ( j) by using Eq. (68), and local damage for the material point x (k) by using Eq. (69).  Table 2. The hyperelastic membrane is tested without a defect, with a defect in the form of a hole and a pre-existing crack. The deformed configuration of a membrane-type acoustic metamaterial (MAM) is subsequently verified. Secondly, the visco-hyperelastic membrane under situations of without defect, with a defect in the form of a hole and a pre-existing crack and MAM-type structure are validated. The visco-hyperelastic parameters based on Yeoh model and Prony series, c 10 , c 20 ,c 1 10 ,c 1 20 , ϕ ref and η, are given in Table 2. The membranes are discretized into uniform grids. The discretization is chosen as 100 × 100 for the membrane without defect and with a defect of a hole and a pre-existing crack. In terms of the membrane-type acoustic material (MAM), the discretization is chosen as 120 × 120 in order to capture the interaction between multiple materials. In view of the convergence for PD predictions, the horizon in this study is specified as δ = 3.015 x.
The displacement boundary condition is imposed through the boundary layers. The crack propagation problems are investigated by using dynamic analysis with a time step size of t = 10 −4 s. All other cases are investigated for quasi-static analysis by using ADR method.
In order to verify the accuracy of PD model, the membrane is also simulated by FEM. The acceptable finite element discretization is achieved by 100 × 100. The element type PLANE 182 with Yeoh model and Prony series are utilized in the ANSYS model.

Hyperelastic membrane without a defect
The hyperelastic membrane is subjected to the homogenous principle stretch λ as shown in Fig. 5. The principle stretch λ is varied from 1.5 to 4 with an increment of 0.5. Figure 6 represents the PD and ANSYS horizontal displacement predictions at the specified point (x = 0.25 m, y = 0.25 m). As can be seen from this figure, the comparison of PD and ANSYS horizontal displacement predictions indicates a very good agreement.
The PD stress prediction at specified point (x = −0.165 m, y = 0.125 m) is presented in Fig. 7 for both elastic and hyperelastic membranes. PD stress can be calculated by using the definition provided by Zimmermann [45] as in which t represents the force density vector and represents the volume of the body. The force density for elastic membrane can be calculated by using bond-based theory as [40] t e (k)( j) = c · s · y ( j) −y (k) and the force density for hyperelastic membrane can be calculated by using force density given in Eq. (74). It is obvious that the PD stress result for the hyperelastic membrane is greater than the result for the elastic membrane. As can be seen from Fig. 7, the difference between PD stress is 0.31 × 10 6 Pa for λ = 1.5 and 9.12 × 10 6 Pa for λ = 4. The increasing deviation between the PD stress predictions indicates the hyperelastic material property of the membrane.

Hyperelastic membrane with a hole
The hyperelastic membrane with a hole of radius R = 0.08 m is subjected to homogenous principle stretch as shown in Fig. 8. Figures 9, 10 represent the displacement components for λ = 4 and λ = 7 in the deformed configuration. As it can be seen from the figures, the shape of the hole remains circular. The contour plots of PD deformation predictions agree well with the ANSYS predictions, validating the PD model even for large deformations.

Hyperelastic membrane with a pre-existing crack
In the presence of a pre-existing crack, the hyperelastic membrane in Fig. 11 includes a vertical crack with a length of l = 0.1m in the centreline of the membrane. In the situation not allowing failure, the membrane is subjected to homogeneous principle stretch as λ = 1.5 and λ = 4. Figures 12, 13 represent the displacement components in the deformed configuration. As it can be seen from the results, PD predictions compare well with ANSYS results. There is a slight difference near the crack tips between the PD and ANSYS results. This is due to more refined mesh in ANSYS model near the crack tips. In Fig. 13, as expected, the crack opening increases with increased loading.
When failure is allowed, the critical energy release rate is chosen as G c = 6.69 N/m for the hyperelastic material. The homogenous loading stretch is specified asλ = 205 (1/s). The crack propagation at the end of t = 0.022 s, t = 0.023 s, t = 0.024 s and t = 0.025 s is presented in Fig. 14. The failure expands in both x-and y-directions simultaneously. As the initial crack is aligned with y-direction, the top and bottom part of failure zone continuously expands as crack branches.

Hyperelastic membrane with rigid inclusions
The membrane-type acoustic metamaterial (MAM) is generally constructed by elastic membrane and rigid inclusions. The solid inclusions are bonded into the membrane. As shown in Fig. 15a, four inclusions and a hyperelastic membrane are considered to construct the MAM. The radius of solid inclusions is r 0 = 0.0625, and their centres are located at (−2.5r 0 , 0), (2.5r 0 , 0), (0, −1.5r 0 ) and (0, 1.5r 0 ). Since the hyperelastic membrane has a Young's modulus of E = 1.1976 MPa, the material parameters for solid inclusions are chosen as E 2 = 144 GPa [46] and ν = 0.3. Note that if the material point is inside the rigid inclusion, the force density is calculated based on Eq. (79); on the other hand, if the material point is inside the membrane, the force density is calculated based on Eq. (74).
In the absence of a pre-existing crack, the MAM is subjected to homogenous principle stretch as λ = 1.5. Figure 16 represents the displacement components in the deformed configuration. As can be seen from this figure, the PD predictions are in good agreement with those of ANSYS. The predictions capture the effect interaction of the solid inclusions with the hyperelastic membrane.
In the presence of a crack, MAM includes a horizontal crack with length of l = 0.1m as shown in Fig. 15b. The PD predictions of displacement contour plots for λ = 1.5 are presented in Fig. 17. Comparing with displacement contours of PD and ANSYS results, the effect of surrounding inclusions is accurately captured.
When failure is allowed, the critical energy release rate is chosen as G c = 6.69 N /m for hyperelastictype membrane. The homogenous loading stretch is specified asλ = 1 (1/µs) with dt = 0.1 µs. Figure 18 represents the crack propagation at the end of 150 µs, 170 µs, 180 µs and 190 μs. The maximum damage values are observed, especially close to the tips of initial crack. The failure zone expands in both x-and ydirections. The crack propagates around the edge of the solid inclusions since the critical energy release rate for solid inclusions is much greater than for the hyperelastic membrane.

Visco-hyperelastic membrane without a defect
In the absence of defect, the visco-hyperelastic membrane is subjected to the homogenous principle stretch λ which is varied from 1.25 to 2.5 with an increment of 0.25. The visco-hyperelastic behaviour is evaluated at t = 2.5 µs, and strain rate is specified asε 11 = 0.15. Figure19 presents the PD stress at (x = −0.165 m, y = 0.125 m) for both visco-hyperelastic and hyperelastic membrane. As it can be seen from Fig. 19, PD stresses for hyperelastic membrane are bigger than those of visco-hyperelastic membrane. The difference between PD stress is 0.24 × 10 5 Pa for λ = 1.25 and 1.20 × 10 6 Pa for λ = 2.5. The stiffness property of visco-hyperelastic membrane degrades with time, leading to the increasing deviation between PD stress of hyperelastic and visco-hyperelastic membranes.

Visco-hyperelastic membrane with a hole
In the presence of a hole, the visco-hyperelastic membrane is subjected to the applied loading stretch of λ = 2.5 and λ = 3.0. The constant strain rate is specified asε 11 = 0.15. Figures 20, 21 present the displacement components in the deformed configuration. As it can be seen from the figures, the shape of the hole remains circular. The contour plots of PD deformation predictions agree well with the ANSYS predictions. 6.7 Visco-hyperelastic membrane with a pre-existing crack In the presence of pre-existing crack, the visco-hyperelastic membrane is subjected to homogeneous principle stretch as λ = 1.5 and λ = 2.5. When failure is allowed, the critical energy release rate is chosen as G c = 6.69 N/m for the viscohyperelastic material. The homogenous loading stretch is specified as λ = 159.4. Figure 24 presents the crack propagation at the end of 0.025 s, 0.029 s, 0.030 s and 0.031 s. Since the crack is aligned with the y-direction, the crack initially propagates in y-direction. Later, the crack propagates in both x-and y-directions as it can be seen in Fig. 24c, d. The damage pattern is also presented in deformed configuration in Fig. 25. It is observed that a blunted zone occurs at the tips of the crack as shown in Fig. 25.

Visco-hyperelastic membrane with rigid inclusions
For the visco-hyperelastic MAM material, the following cases are considered. In the absence of a pre-existing crack, the membrane is subjected to homogenous principle stretch as λ = 1.5. As it can be seen from Fig. 26, PD predictions are in good agreement with ANSYS results.
In the presence of a pre-existing crack, the visco-hyperelastic membrane includes a horizontal crack with a length of l = 0.1 in the centreline of the membrane. In this case, the membrane is subjected to homogeneous principle stretch as λ = 1.2. Figure 27 represents the contour plots of displacement fields. As can be seen from this figure, both PD and ANSYS predictions agree very well. When failure is allowed, the critical energy release rate is chosen as G c = 6.69 N /m for visco-hyperelastictype membrane. The homogenous loading stretch is specified asλ = 1 (1/μs) with dt = 0.1 µs. Figure 28 presents the crack propagation at the end of 150 µs, 170 µs, 180 µs and 190 µs. By comparing with Fig. 18, the visco-hyperelastic membrane shows less damage since the stresses for visco-hyperelastic membrane are smaller than those of hyperelastic membrane.

Conclusion
The novelty of this study is to develop a PD model for visco-hyperelastic material. It specially concerns Yeohtype membrane under homogenous stretch loading conditions. The Yeoh strain energy density is expressed by the invariants of Cauchy-Green tensor. The material parameters in the strain energy function are identified by curve fitting of tension test data. The PD form of the visco-hyperelastic force density function is represented in terms of the PD integral. The viscous stretch term for viscoelastic properties is derived based on Yeoh strain energy density. The developed PD model is utilized to predict the deformation and damage for a polyurea membrane with a defect in the form of a hole and a pre-existing crack. The deformation and damage for a membrane-type acoustic metamaterial with inclusions are also simulated by using the developed PD model. The numerical results capture the quasi-static deformation and dynamic fracture propagation of the polyurea membrane and membrane-type acoustic metamaterial under homogenous loading. The developed PD model is verified by obtaining good agreement with ANSYS predictions.

Appendix A
The model parameters identification for high-damping rubber and polyurea in Section 3.2 is based on the basic optimal curve-fitting method [47]. This nonlinear curve-fitting problem is solved in least-squares sense named as Isqcurvefit in MATLAB curve-fitting toolbox, using the function below where F represents the objective function, a represents the material parameters that need to be identified, x data represents the given input data and y data represents the observed output. To find the material parameters, the objective function F (a, x data ) is specified as the constitutive stressstretch relation from Eq. (30) as where λ data represents the stretch data from the compressive tests for different strain rates. Therefore, the least-squares function in Eq. (A.1) can be specified as where σ 11,data represents the stress data from the compressive tests for different strain rates. The integration in Eq. (A.2) is solved by numerical integration function trapz in MATLAB. This numerical function performs numerical integration via the trapezoidal method. This method approximates the integration over an interval by breaking the area down into trapezoids with more easily computable areas. For integration with N + 1 evenly spaced points, the approximation can be defined as Therefore, the integration in Eq. (A.2) then can be evaluated by using the following equation: where the integration parameter N is defined as 100.
In least-squares method, the better solution of parameter identification can be obtained by defining initial value, upper and lower bounds for each parameter. The initial values, lower bound and upper bounds for c 10 , c 20 ,c 1 10 ,c 1 20 , ϕ ref and η for both high-damping rubber and polyurea are given in Table 3.