Micropolar Peridynamic Stress-Based Model for Concrete Structures

1. Introduction

Enhanced computational and modeling tools are required to study new construction materials, such as concrete. Classical Continuum Mechanics (CCM) and Finite Element Methods (FEM) are advanced methods used to model the behavior of material structures. However, when loads induce crack formation and propagation, damage evolves, violating the assumptions of CCM, and singularities start to appear. To address the shortcomings of continuum models, a nonlocal model called the peridynamic model was proposed in 2000 by Silling [1].

The governing equations in the peridynamic model do not assume spatial differentiability and allow discontinuities to naturally arise as a part of the solution [1]. The bond-based peridynamic model was applied to model concrete structures in ([2,3,4,5,6,7,8]); however, it cannot model materials with Poisson’s ratio different from 1/4 in three dimensions [9]. The micropolar peridynamic model was proposed to improve the bond-based peridynamic model, which allowed the moment and force densities to interact among particles within a finite limit, known as the material horizon ([9,10,11]). In addition, a more general model, called the state-based peridynamic model [12] was developed. State-based models have been applied to concrete structures in ([13,14,15]). Peridynamic state-based models can be categorized into ordinary state-based and nonordinary state-based models. The force state is assumed to be parallel to the deformed state in the former model, whereas the shape and stress tensors are used in the force state in the latter model [12]. These state-based models are more complicated and computationally more expensive than the bond-based model.

Micropolar elasticity, or micropolar continuum, is an extension of CCM; it was first developed by the Cosserat brothers at the beginning of the twentieth century [16]. Micropolar elasticity introduces a characteristic length in its constitutive equations[17]. Micropolar Peridynamics (MPD) was first developed by Gerstle et al. and was based on Euler-Bernoulli beams comprising shear force and moment densities, as well as axial force densities with displacements and rotations ([9,10,11]). Based on the micropolar continuum and the state-based peridynamic model, Debasish Roy established a state-based micropolar peridynamic theory for linear elastic solids in 2015 [18]. Subsequently, Diana et al. formulated a generalized micropolar peridynamic model with shear deformability and proposed failure criteria based on deformation [19]. Using this model and based on a stretch and shear failure criteria, Diana et al. simulated fractures in rocks [20]. Based on the micropotential energy that defines orthotropy, Diana and Ballarini developed an orthotropic micropolar peridynamic model using the stretch and fracture energy failure criteria [21]. Guo et al. also developed a Cosserat peridynamic model, in which the rotation was independent of translational displacement [22]. Other types of micropolar peridynamic models that incorporate viscoelasticity, nonuniform material horizon, have also been proposed [23] and [24]. By combining the peridynamic differential operator and micropolar peridynamics, Lei et al. developed a micropolar damage model for concrete structures[25].

The classical interpretation of stress, formalized by Augustine-Louis Cauchy in the early 19th century [26], assumes a local contact force acting on an infinitesimal area element within a continuum. The Cauchy stress and its Lagrangian counterparts, the Piola-Kirchhoff tensors [27], remain the standard for macroscopic structural analysis; its reliance on spatial derivatives makes it mathematically ill-defined at material discontinuities. To maintain thermodynamic consistency and facilitate a bridge between nonlocal interactions and classical continuum mechanics, the concept of a Peridynamic Stress Tensor was subsequently derived. The peridynamic stress tensor was first introduced by Silling in 2000 [1] and later modified by Lehoucq et al. in 2008 [28]. Fallah et al. proposed a computational approach to calculate the peridynamic stress tensor and compared their calculations with the stress tensor using FEM [29]. Ballarini et al. analyzed singular stress fields using the peridynamic stress tensor and found that as the material horizon decreased, the peridynamic solution converged into the analytical solution [30]. Parallel to these continuum advancements, the rise of molecular dynamics demanded a measure of stress that could bridge the gap between discrete atomic trajectories and macroscopic observables. This is addressed by virial stress, which originates from Clausius’s work in the 19th century, but was adapted for atomic systems by Irving and Kirkwood [31]. Li and his collaborators showed that the peridynamic stress is the first Piola-Kirchhoff virial stress, and by deriving a simpler formula of the peridynamic stress, they developed an expression that can be easily implemented in numerical computations [32]. This framework was also used in [33] to obtain a peridynamic model based on deformations.

Accurate simulation of failure often hinges on the selection of an appropriate damage criterion. Traditionally, peridynamic bond failure is governed by critical stretch or critical bond energy density criteria. Recent advances have shifted towards stress-based failure criteria within the peridynamic framework. The stress-based failure criterion using micropolar peridynamics was used in [34] and [35] where State-Based peridynamics was used to obtain the peridynamic stress tensor and in [36] and in [37], where the CCM version of the stress tensor was used. By establishing connections between the macroscopic and mesoscopic definition of the peridynamic stress tensor, Song et al. developed a framework that integrates the stress-state-related Mohr-Coulomb failure criterion for concrete structures [38]; however, their approach is only valid for a Poisson’s ratio of 1/3 for plane stress problems.

The maximum distance of the interactions between material points is defined by the peridynamic horizon, which is an important parameter that imparts non-locality. Choosing a proper horizon involves a trade-off between accuracy, computational efficiency, and the ability to reproduce physical behaviors such as crack propagation and branching. Although many peridynamic studies treated the horizon as a non-physical parameter, in [9,39,40] expressions for the peridynamic horizon were developed as functions of material properties such as the Fracture Energy, the Modulus of Elasticity, and the allowable maximum stretch.

In this study, the micropolar peridynamic stress tensor was postulated for plane stress problems, in which the peridynamic horizon was expressed as a function of Poisson’s ratio, the material strength and the fracture toughness.

In order to obtain direct relationships between the analytical peridynamic stress tensor and material parameters, Lehoucq’s definition of the peridynamic stress tensor was employed. Nevertheless, due to progressive link failure, the calculation of peridynamic stress is time-consuming; therefore, calibration parameters that depend on the state of stress and stretch at each point were proposed. In addition, an efficient numerical scheme was proposed to calculate the micropolar peridynamic stress tensor, which was incorporated into an explicit dynamic relaxation solver to model plane stress concrete problems. In these problems, failure criteria based on the maximum stretch and stress were implemented.

The remainder of this paper is organized as follows. Section 2 describes the basis of micropolar peridynamics and the micropolar peridynamic stress tensor, where closed-form solutions for linear elastic problems were described and relationships between the conventional stress-strain and the micropolar peridynamic stress-strain tensor were found. This Section also explains a numerical technique and introduces a damage model for concrete structures. Section 3 provides numerical examples of the model showing comparisons with physical tests. Conclusions and a discussion are given in Section 4.

2. Materials and Methods

Micropolar peridynamics is summarized herein on the basis of previously proposed definitions in [1,9], and the Micropolar Peridynamic Stress Tensor is postulated based upon the definition in [28]. Two-dimensional linear elastic examples and applications to linear elastic fracture mechanics are provided. A numerical scheme is explained and a proposed damage model for concrete structures is described.

2.1. The Micropolar Peridynamic Stress Tensor

The bond-based peridynamic model postulates a vector function in units of force per unit volume squared $f_{ij}$ that represents the interaction between particles i and j inside a material horizon $\delta$ [1]. The sum of all internal forces per unit volume acting on a particle i is expressed as follows:

$${\int }_{V_j}{f}_{ij}d{V}_j+b_i={\rho }_i{\displaystyle \frac{{\partial }^2{u}_i}{\partial t^2}},$$

(1)

where $b_i$ is the body force acting on particle i; ${\displaystyle \frac{{\partial }^2{u}_i}{\partial t^2}}$ and ${\rho }_i$ are the acceleration and density of particle i , respectively. This integral was performed on all particles j within $\delta$. In the micropolar peridynamic model, a moment equation was added [9].

$${\int }_{V_j}{m}_{ij}d{V}_j+{\widehat{m}}_i=0,$$

(2)

where the moment per unit volume squared is $m_{ij}$ and ${\widehat{m}}_i$ is a moment density function [9]. These moments and forces are functions of the relative positions r, the relative displacements ($u_j-u_i$, $v_j-v_i$) and the rotations of particles i and j (${\theta }_i$,${\theta }_j$). Forces $f_{ij}$ and moments $m_{ij}$ between particles can be expressed using the following matrix at the linear microelastic level as follows:

$$\left[\begin{array}{c}\widehat{f_x^i}\\ \widehat{f_y^i}\\ \widehat{m_z^i}\\ \widehat{f_x^j}\\ \widehat{f_y^j}\\ \widehat{m_z^j}\end{array}\right]=\left[\begin{array}{cccccc}c/r& 0& 0& -c/r& 0& 0\\ 0& 12d/r^3& 6d/r^2& 0& -12d/r^3& 6d/r^2\\ 0& 6d/r^2& 4d/r& 0& -6d/r^2& 2d/r\\ -c/r& 0& 0& c/r& 0& 0\\ 0& -12d/r^3& -6d/r^2& 0& 12d/r^3& -6d/r^2\\ 0& 6d/r^2& 4d/r& 0& -6d/r^2& 4d/r\end{array}\right]\left[\begin{array}{c}\widehat{u_i}\\ \widehat{v_i}\\ \widehat{{\theta }_z^i}\\ \widehat{u_j}\\ \widehat{v_j}\\ \widehat{{\theta }_z^j}\end{array}\right].$$

(3)

The relations between peridynamic constants $(c,d)$ and conventional linear elastic constants $(E,\nu )$ can be found in [10].

Cauchy’s first equation of motion can be described as ${\sigma }_{k,l}+b_k=\rho {\ddot{u}}_k$, where ${\sigma }_{k,l}$ is the divergence of the stress tensor, $b_k$ is the body force, and ${\ddot{u}}_k$ is the acceleration of each point in the domain. Compared with Equation (1), the following equality must hold [28]:

$${\sigma }_{k,l}={\int }_{V^{\prime }}{f}_{k}d{V}^{\prime },$$

(4)

where ${\sigma }_{k,l}$ is the divergence of the stress tensor, $f_k$ is the force between particles i and particle j , and $d{V}^{\prime }$ is the volume of particle j . Thus, the peridynamic stress tensor at point x with normal n was defined as [28]:

$${\sigma }_{kl}\left(x_l\right)={\displaystyle \frac{1}{2}}{\int }_S{\int }_z{\int }_y{(y+z)}^2{f}_k(x_l+y{m}_l,x_l-z{m}_l)m_{l}dydzd\mathsf{\Omega },$$

(5)

where $m_l$ is a unit vector in the link direction; y and z are scalars; and $d\mathsf{\Omega }$ is a differential of an angle over a sphere S with radius $\delta$. Based on this definition, the peridynamic stress tensor at $x_l$ can be obtained by calculating the interaction between two volumes with centers located at $x_l$ with volumes $d{A}_{y}d{y}_l$ and $d{A}_{z}d{z}_l$ respectively [28]. For two-dimensional models, the stress tensor is given by

$${\sigma }_{kl}\left(x_l\right)={\displaystyle \frac{1}{2}}{\int }_L{\int }_z{\int }_y(y+z)f_k(x_l+y{m}_l,x_l-z{m}_l)m_{l}dydzd\mathsf{\Omega },$$

(6)

where $d\mathsf{\Omega }$ is the differential of an angle on the arc of length L. Based on Equation (5), for a prescribed and positive material horizon, the stress tensor in peridynamic theory can also be expressed using spherical and cylindrical coordinates as follows:

$${\sigma }_{kl}\left(x_l\right)={\displaystyle \frac{1}{2}}{\int }_0^{\pi }{\int }_0^{2\pi }{\int }_0^{\delta }{\int }_0^{\delta -y}{(y+z)}^2{f}_k(p_l,q_l)m_{l}sin\theta dzdyd\varphi d\theta ,$$

(7)

where $p_l=x_l+y{m}_l$ and $q_l=x_l-z{m}_l$. Using cylindrical coordinates, the two-dimensional peridynamic stress tensor becomes

$${\sigma }_{kl}\left(x_l\right)={\displaystyle \frac{1}{2}}{\int }_0^{2\pi }{\int }_0^{\delta }{\int }_0^{\delta -y}(y+z)f_k(p_l,q_l)m_{l}dzdyd\theta ,$$

(8)

where $f_k$ is the force between the particles i and j and $m_l$, is the unit vector with its origin located at i, and has the orientation of the link between these particles.

Using these definitions, the stress tensor using micropolar peridynamics and the conventional stress tensor were compared at the linear microelastic level. Assuming constant strains, the displacement field $u_k$ can be expressed as follows.

$$u_k={ϵ}_{kl}{\xi }_k,$$

(9)

where ${ϵ}_{lk}$ is the conventional strain tensor, ${\xi }_k$ is the relative position, expressed as ${\xi }_1=rcos\theta$, ${\xi }_2=rcos\theta sin\varphi$, ${\xi }_3=rsin\theta sin\varphi$ for three-dimensional models and ${\xi }_1=rcos\theta$, ${\xi }_2=rsin\theta$ for two-dimensional models. By differentiating Equation (9), the following relationships are obtained $u_{k,k}={ϵ}_{kk}$, and ${\gamma }_{kl}=2{ϵ}_{kk}=u_{k,l}+u_{l,k}$. Furthermore, rotations at the microelastic level are expressed as ${\omega }_{kl}={\displaystyle \frac{1}{2}}(u_{k,l}-u_{l,k})=0$. Thus, rotations are eliminated, and the following transformation matrices are used in equations (7) and (8), respectively.

$$\left[T\right]=\left[\begin{array}{ccc}sin\theta sin\varphi & cos\theta & sin\theta cos\varphi \\ -cos\theta sin\varphi & sin\theta & -cos\theta cos\varphi \\ -cos\varphi & 0& sin\varphi \end{array}\right],\left[T\right]=\left[\begin{array}{cc}cos\theta & sin\theta \\ -sin\theta & cos\theta \end{array}\right].$$

(10)

The following expressions for the stress tensor using micropolar peridynamics for the three-dimensional case (11) and the two-dimensional case (12) were obtained after completing the integrations in equations (7) and (8) with the displacement field in Equation (9), assuming a uniform displacement field for a linear isotropic solid.

$${\sigma }_{kl}=\left\{{\displaystyle \frac{\pi c{\delta }^4}{30}}-{\displaystyle \frac{\pi d{\delta }^2}{15}}\right\}{ϵ}_{mm}{\delta }_{kl}+2\left\{{\displaystyle \frac{\pi c{\delta }^4}{30}}+{\displaystyle \frac{\pi d{\delta }^2}{10}}\right\}{ϵ}_{kl},$$

(11)

$${\sigma }_{kl}=\left\{{\displaystyle \frac{\pi c{\delta }^3}{24}}-{\displaystyle \frac{3\pi d\delta }{2}}\right\}{ϵ}_{mm}{\delta }_{kl}+2\left\{{\displaystyle \frac{\pi c{\delta }^3}{24}}+{\displaystyle \frac{3\pi d\delta }{2}}\right\}{ϵ}_{kl}.$$

(12)

The second term within brackets in equations (11) and (12) denotes the shear modulus G, which can also be obtained by rotating a pure shear state by ${45}^{\circ }$, ${\sigma }_{max}=\tau$, ${\sigma }_{min}=-\tau$. equations (11) and (12) can be obtained by differentiating the micropolar peridynamic strain energy density as follows:

$${\sigma }_{kl}={\displaystyle \frac{\partial U}{\partial {ϵ}_{kl}}}.$$

(13)

The stress tensor in equations (11) and (12) is valid only at the linear microelastic level with zero microrotations. For the case where microrotations have different values (${\theta }_i={\varphi }_1$, ${\theta }_j={\varphi }_2$), the shear stress component (${\sigma }_{kl},k\ne l$) for two-dimensional isotropic plane stress solids becomes

$${\sigma }_{kl}=\left\{{\displaystyle \frac{\pi c{\delta }^3}{12}}+3\pi d\delta \right\}{ϵ}_{kl}-3\pi d\delta ({\varphi }_1+{\varphi }_2),$$

(14)

and

$${\sigma }_{lk}=\left\{{\displaystyle \frac{\pi c{\delta }^3}{12}}+3\pi d\delta \right\}{ϵ}_{lk}+3\pi d\delta ({\varphi }_1+{\varphi }_2).$$

(15)

Moreover, the couple stress components $m_{kl}$ for plane stress problems are

$$m_{kl}=\pi d{\delta }^2({\varphi }_1+2{\varphi }_2)$$

and

$$m_{lk}=\pi d{\delta }^2(2{\varphi }_1+{\varphi }_2).$$

For values of $d=0$ or $\nu ={\displaystyle \frac{1}{3}}$, ${\sigma }_{lk}={\sigma }_{kl}$.

By assuming linear elasticity, isotropy, and known displacements, closed-form solutions for the stress distribution using micropolar peridynamics were obtained in the following examples.

2.1.1. A Bar Subjected to Uniaxial Stress and Uniaxial Strain

A bar was subjected to a uniaxial load, where the stress tensor was obtained using Equations (7) and (8). The displacement field for particle i is given by $u_1={ϵ}_{11}{x}_1$, $u_2={ϵ}_{22}{x}_2$ for the two-dimensional model and $u_1={ϵ}_{11}{x}_1$, $u_2={ϵ}_{22}{x}_2$, $u_3={ϵ}_{33}{x}_3$ for the three-dimensional model. Similarly, for particle j, $u_1={ϵ}_{11}(x_1+(y+z)m_1)$, $u_2={ϵ}_{22}(x_2+(y+z)m_2)$, with $m_1=cos\theta$, $m_2=sin\theta$ for the two-dimensional model and $u_1={ϵ}_{11}(x_1+(y+z)m_1)$, $u_2={ϵ}_{22}(x_2+(y+z)m_2)$, $u_3={ϵ}_{33}(x_3+(y+z)m_3)$, with $m_1=cos\theta ,m_2=cos\theta sin\varphi$, and $m_3=sin\theta sin\varphi$ for three-dimensional models. For the uniaxial stress case ${ϵ}_{33}=-\nu {ϵ}_{11}$ and ${ϵ}_{22}=-\nu {ϵ}_{11}$, ${\sigma }_{22}={\sigma }_{33}=0$, ${\sigma }_{12}={\sigma }_{13}={\sigma }_{23}=0$. For uniaxial strain ${ϵ}_{33}={ϵ}_{22}=0$, and after calculating the integrals, the uniaxial stress can be expressed as

$$\left[\sigma \right]=E\left[\begin{array}{cc}{ϵ}_{11}& 0\\ 0& 0\end{array}\right],\left[\sigma \right]=E\left[\begin{array}{ccc}{ϵ}_{11}& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right],$$

using both models. The result for the uniaxial strain obtained using the three-dimensional model is as follows.

$$\left[\sigma \right]={\displaystyle \frac{E}{(1-2\nu )(1+\nu )}}\left[\begin{array}{ccc}(1-\nu ){ϵ}_{11}& 0& 0\\ 0& \nu {ϵ}_{11}& 0\\ 0& 0& \nu {ϵ}_{11}\end{array}\right],$$

and the two-dimensional model is

$$\left[\sigma \right]={\displaystyle \frac{E}{(1-{\nu }^2)}}\left[\begin{array}{cc}{ϵ}_{11}& 0\\ 0& \nu {ϵ}_{11}\end{array}\right].$$

Because of the assumption of uniaxial strain and plane stress conditions, the stresses obtained are exactly the same as in classical elasticity.

2.1.2. A Plate Subjected to Plane Stress and Plane Strain

Stresses on a two-dimensional plate were calculated for plane stress and plane strain conditions assuming a uniform displacement field given by

$$u_l={ϵ}_{kk}{x}_l,$$

(16)

for particle i, where $k=l$ and $k=l=1,2$ for the plane stress model and $k=l=1,2,3$ for the three-dimensional model. Similarly, for particle j,

$$u_l={ϵ}_{kk}(x_l+(y+z)m_l),$$

(17)

where $m_1=cos\theta$, $m_2=sin\theta$ for the two-dimensional model and $m_1=cos\theta ,m_2=cos\theta sin\varphi$, and $m_3=sin\theta sin\varphi$ for the three-dimensional model. Thus, displacements for the plane stress case expressed in equations (16) and (17) using ${ϵ}_{33}={\displaystyle \frac{\nu }{(\nu -1)}}\left\{{ϵ}_{11}+{ϵ}_{22}\right\}$, in which ${\sigma }_{33}=0$. The stress tensor for this case is given by

$$\left[\sigma \right]={\displaystyle \frac{E}{(1-{\nu }^2)}}\left[\begin{array}{ccc}{ϵ}_{11}+\nu {ϵ}_{22}& 0& 0\\ 0& {ϵ}_{22}+\nu {ϵ}_{11}& 0\\ 0& 0& 0\end{array}\right].$$

The stress tensor for the plane strain case using ${ϵ}_{33}=0$ is given by

$$\left[\sigma \right]={\displaystyle \frac{E}{(1-2\nu )(1+\nu )}}\left[\begin{array}{ccc}(1-\nu ){ϵ}_{11}+\nu {ϵ}_{22}& 0& 0\\ 0& (1-\nu ){ϵ}_{22}+\nu {ϵ}_{11}& 0\\ 0& 0& \nu ({ϵ}_{11}+{ϵ}_{22})\\ \end{array}\right],$$

The three-dimensional model was used for this calculation. In both cases, ${\sigma }_{12}={\sigma }_{21}={\sigma }_{32}={\sigma }_{23}={\sigma }_{13}={\sigma }_{31}=0$ and $m_{12}=m_{21}=m_{32}=m_{23}=m_{13}=m_{31}=0$, since ${\displaystyle \frac{\partial u_i}{\partial x_j}}=0$ for $i\ne j$.

These results show that the micropolar peridynamic stress tensor matches the stress tensor in CCM. Using this methodology, expressions of the micropolar peridynamic stress are obtained using Mode I displacement fields.

2.1.3. Application to Linear Elastic Fracture Mechanics

In this subsection, stresses near the crack-tip are analyzed using micropolar peridynamics and compared with Linear Elastic Fracture Mechanics (LEFM). To obtain the micropolar peridynamic stress, known displacements were assumed and boundary conditions were verified. Numerical integration with adaptive quadrature was used, in which microrotations were assumed to be zero to avoid singularities. In addition, integration was evaluated from $0<\theta <\pi$ to avoid discontinuity in the notch. Assuming sharp notches with infinitesimally small widths, the displacements (u, v) for Mode I are given by [41]:

$$u={\displaystyle \frac{K_I}{2G}}\sqrt{{\displaystyle \frac{R}{2\pi }}}cos\left({\displaystyle \frac{\theta }{2}}\right)\left[\kappa -1+2{sin}^2\left({\displaystyle \frac{\theta }{2}}\right)\right],$$

(18)

$$v={\displaystyle \frac{K_I}{2G}}\sqrt{{\displaystyle \frac{R}{2\pi }}}sin\left({\displaystyle \frac{\theta }{2}}\right)\left[\kappa +1-2{cos}^2\left({\displaystyle \frac{\theta }{2}}\right)\right],$$

(19)

where $K_I$ is the stress intensity factor for Mode I, G is the shear modulus, $\nu$ is the Poisson’s ratio and $\kappa =(3-\nu )/(1+\nu )$. Polar coordinates were used for integration, in which particles i are located at a distance R from the crack-tip, and particles j at a distance $R^{\prime }$ as shown in Figure 1. $\mathsf{\Delta }=\varphi -\theta$, $\alpha ={180}^o-\varphi +\theta$, $R^{\prime }=R^2+r^2-2Rrcos\alpha$, $\beta =arcsin\left(r/R^{\prime }sin\alpha \right)$, $x_o^j=x_o^i+(y+z)cos\varphi$, $y_o^j=y_o^i+(y+z)sin\varphi$. Thus, employing equations (18) and (19), the displacements for particles i and j are given by the following:

$$u^i={\displaystyle \frac{K_I}{2G}}\sqrt{{\displaystyle \frac{R}{2\pi }}}cos\left({\displaystyle \frac{\theta }{2}}\right)\left[\kappa -1+2{sin}^2\left({\displaystyle \frac{\theta }{2}}\right)\right],$$

$$u^j={\displaystyle \frac{K_I}{2G}}\sqrt{{\displaystyle \frac{R^{\prime }}{2\pi }}}cos\left({\displaystyle \frac{\theta +\beta }{2}}\right)\left[\kappa -1+2{sin}^2\left({\displaystyle \frac{\theta +\beta }{2}}\right)\right],$$

$$v^i={\displaystyle \frac{K_I}{2G}}\sqrt{{\displaystyle \frac{R}{2\pi }}}sin\left({\displaystyle \frac{\theta }{2}}\right)\left[\kappa +1-2{cos}^2\left({\displaystyle \frac{\theta }{2}}\right)\right],$$

$$v^j={\displaystyle \frac{K_I}{2G}}\sqrt{{\displaystyle \frac{R^{\prime }}{2\pi }}}sin\left({\displaystyle \frac{\theta +\beta }{2}}\right)\left[\kappa +1-2{cos}^2\left({\displaystyle \frac{\theta +\beta }{2}}\right)\right].$$

Observe that at the tip of the crack, as $R\to 0$, $R^{\prime }\to r$, $\beta \to 0$, $\theta \to \varphi$ (Figure 1), stresses at the tip of the crack in the x and y directions are given by:

$${\sigma }_{xx}={\displaystyle \frac{16K_I(27{\nu }^2-263\nu -10)}{525({\nu }^2-1)}}\sqrt{{\displaystyle \frac{2}{{\pi }^3\delta }}},$$

(20)

$${\tau }_{xy}={\displaystyle \frac{16K_I(153{\nu }^2-732\nu +235)}{525(1-{\nu }^2)}}\sqrt{{\displaystyle \frac{2}{{\pi }^3\delta }}},$$

(21)

$${\tau }_{yx}={\displaystyle \frac{8K_I(219{\nu }^2-1336\nu +405)}{525({\nu }^2-1)}}\sqrt{{\displaystyle \frac{2}{{\pi }^3\delta }}},$$

(22)

$${\sigma }_{yy}={\displaystyle \frac{32K_I(39{\nu }^2+79\nu -100)}{525({\nu }^2-1)}}\sqrt{{\displaystyle \frac{2}{{\pi }^3\delta }}}.$$

(23)

As an example, assume that we have a single edge notch plate subjected to a tensile stress $\sigma$ of 2 MPa. The notch value a is 30 mm and the width W is 150 mm. The material parameters are given by $E=$18,000 MPa and $\nu$ = 0.20. The stress intensity factor for this case is obtained with the following formula:

$$K_I=Y\sigma \sqrt{\pi a},$$

(24)

$$Y={\displaystyle \frac{2W}{\pi a}}tan\left({\displaystyle \frac{\pi a}{2W}}\right){\displaystyle \frac{\left[0.752+2.02\frac{a}{W}+0.37{\left(1-sin\left(\frac{\pi }{a}\right)\right)}^3\right]}{cos\left(\frac{\pi a}{2W}\right).}}$$

(25)

the value of $K_I$ obtained was $26.98$ MPa-mm^1/2, and the stress ${\sigma }_{yy}$ using Equation (23) is $20.4$ MPa. Figure 2 shows a plot of micropolar peridynamic stress ${\sigma }_{yy}$ using numerical integration compared to a plot of the same stress using LEFM. In addition, the stress value of $20.4$ MPa is shown at the tip of the crack.

On the basis of MPD and the micropolar peridynamic stress, a numerical model was developed, which is explained in the following section.

2.2. Numerical Model

A single-processor computer code was developed and, in order to maximize its efficiency, matrix functions were maximized, and loops and logical blocks were minimized; several arrays of stored displacements and material properties were distributed on grids that have one unit of measure. In this model, a material horizon of three times the grid spacing was chosen. Consequently, the stress tensor was calculated using a matrix of material points j with a size of 7 x 7 at each point i Figure 3. This scheme uses a dynamic relaxation method, in which displacement solutions are expressed with the following equation.

$$\left[M\right]\left[\ddot{u}\right]+\left[C\right]\left[\dot{u}\right]+\left[K\right]\left[u\right]=\left[f^{ext}\right],$$

(26)

Where $\left[M\right]$ is the mass matrix, $\left[C\right]$ is the damping matrix, $\left[K\right]$ is the stiffness matrix, $\left[u\right]$ are the displacements and $\left[f^{ext}\right]$ the external forces. By choosing an adequate damping matrix, the transient solution of Equation (26) vanishes.

The geometry and model parameters were distributed in a $n\times m$ array, where n is the number of points in x and m is the number of points in the y direction. The rotational stiffness and stiffness matrices were obtained as follows:

$$k_x={\left[k_x\right]}_{n\times m},k_y={\left[k_y\right]}_{n\times m},k_z={\left[k_z\right]}_{n\times m}.$$

The frequencies were calculated using matrix operations as follows:

$${\omega }_x=\sqrt{{\displaystyle \frac{k_x}{m}}},{\omega }_y=\sqrt{{\displaystyle \frac{k_y}{m}}}.$$

The mass moment of inertia was defined at each point such that the rotational frequency coincided with the translational frequency,

$$I={\displaystyle \frac{k_z}{{\omega }_x^2}},{\omega }_z=\sqrt{{\displaystyle \frac{k_z}{I}}},$$

and the time differential $\mathsf{\Delta }t$ was chosen such that

$${\omega }_{max}=max({\omega }_x,{\omega }_y,{\omega }_z),\mathsf{\Delta }t<{\displaystyle \frac{1}{{\omega }_{max}}}.$$

The properties and displacements of the material were ordered in 7 x 7 arrays to determine the stresses at the node i (Figure 3). To accelerate the calculation of the stress tensor at all points in the domain, some vectors were previously calculated by computing integrals by Lagrange interpolation. The following 49 x 1 distance vectors $\left[{\xi }_1\right]$ and $\left[{\xi }_2\right]$: $\left[{\xi }_1\right]=[-3,..,x_j-x_i,..,3]{\mathsf{\Delta }}_s$, $\left[{\xi }_2\right]=[-3,..,y_j-y_i,..,3]{\mathsf{\Delta }}_s$ were defined. $x_j,y_j,x_i,y_i$ are the coordinates of the particles j and i, respectively, and ${\mathsf{\Delta }}_s$ is the grid spacing. Let the matrix $\left[a\right]$ of 49 x 49 be defined for each element as follows:

$$a_{k,i+7(j-1)}={\xi }_1^{i-1}\left(k\right){\xi }_2^{j-1}\left(k\right).$$

(27)

and, the matrix $\left[c^p\right]$ as

$$c_{ij}^p=a_{ij}^{-1},$$

(28)

for $i=1..7$ , $j=1..7$ and $k=1..49$ where $a_{ij}$ is obtained with Equation (27).

Let $\tilde{x}=(y+z)cos\theta$ and $\tilde{y}=(y+z)sin\theta$, so that $\tilde{x}$ and $\tilde{y}$ are within a radius circle of $(y+z)$. Let a vector 1 x 49 $\left[p\right]$ be defined as

$$\left[p\right]=[{\tilde{x}}^0{\tilde{y}}^0{\tilde{x}}^0{\tilde{y}}^1...{\tilde{x}}^6{\tilde{y}}^6],$$

A coefficient vector $\left[c^{uv}\right]$ was obtained using Equation (28):

$$\left[c^{uv}\right]=\left[c^p\right]\left[p\right].$$

Using Equation (3) and the transformation matrix in Equation (10), the peridynamic link stiffness is given by:

$$\left[k\right]={\left[T\right]}^T\left[\widehat{k}\right]\left[T\right],$$

so that the peridynamic force for stresses in the $xx$ direction and displacements in the x direction can be expressed as:

$$\left[f_{xx}^u\right]=k_{11}\left[c^{uv}\right]cos\left(\theta \right).$$

Similarly, for $xy$ stresses and displacements in the y direction:

$$\left[f_{xy}^v\right]=k_{22}\left[c^{uv}\right]cos\left(\theta \right).$$

Thus, material parameters c for $xx$ stresses and displacements in the x direction are calculated using the following equation:

$$\left[c_{xx}^u\right]={\int }_0^{2\pi }{\int }_0^{\delta }{\int }_0^{\delta -y}\left[f_{xx}^u\right](y+z)dzdyd\theta ,$$

(29)

similarly, $xy$ stresses and displacements of parameters d and in the y direction:

$$\left[d_{xy}^v\right]={\int }_0^{2\pi }{\int }_0^{\delta }{\int }_0^{\delta -y}\left[f_{xy}^v\right](y+z)dzdyd\theta ,$$

(30)

For each material point i with coordinates $(p,q)$, displacements, rotations, and material properties for all particles j with coordinates $(r,s)$ near each particle i, can be represented using the following tensors:

$$u_{pqrs}^s,v_{pqrs}^s,{\theta }_{pqrs}^s,c_{pqrs}^s,d_{pqrs}^s,$$

for $u,v,\theta$ displacements and $c,d$ material properties. Where $-3\le r\le 3$, and $-3\le s\le 3$, in which the boundaries were verified to obtain real physically reasonable values. These tensors were reduced to third-order tensors using 49 nonzero indices instead of r and s.

$$u_{pqk}^s,v_{pqk}^s,{\theta }_{pqk}^s,c_{pqk}^s,d_{pqk}^s.$$

The material parameters c and d can be represented as the following first-order tensors:

$$c_k^{uxx},c_k^{uyy},c_k^{uxy},c_k^{uyx},c_k^{vxx},c_k^{vyy},c_k^{vxy},c_k^{vyx},$$

$$d_k^{uxx},d_k^{uyy},d_k^{uxy},d_k^{uyx},d_k^{vxx},d_k^{vyy},d_k^{vxy},d_k^{vyx},$$

and the following third-order tensors for $xx$ stresses are computed for displacements in the x direction (u):

$$s_{pqk}^{uxx}=c_{pqk}^s{c}_k^{uxx}+d_{pqk}^s{d}_k^{uxx},(\mathrm{No}\mathrm{contraction};k\mathrm{is}\mathrm{now}\mathrm{a}\mathrm{free}\mathrm{index}),$$

and the same stresses in the y direction (v):

$$s_{pqk}^{vxx}=c_{pqk}^s{c}_k^{vxx}+d_{pqk}^s{d}_k^{vxx},(\mathrm{No}\mathrm{contraction};k\mathrm{is}\mathrm{now}\mathrm{a}\mathrm{free}\mathrm{index}).$$

Displacements at each point $(p,q)$ are represented as first-order tensors $u_k$ and $v_k$, so the stresses are computed using the following equations.

$${\sigma }_{pq}^{xx}=s_{pqk}^{uxx}{u}_k+s_{pqk}^{vxx}{v}_k$$

(31)

$${\sigma }_{pq}^{yy}=s_{pqk}^{uyy}{u}_k+s_{pqk}^{vyy}{v}_k$$

(32)

$${\sigma }_{pq}^{xy}=s_{pqk}^{uxy}{u}_k+s_{pqk}^{vxy}{v}_k$$

(33)

$${\sigma }_{pq}^{yx}=s_{pqk}^{uyx}{u}_k+s_{pqk}^{vyx}{v}_k.$$

(34)

At each point $(p,q)$, ${\sigma }_{xx}$, ${\sigma }_{yy}$, ${\sigma }_{xy}$, ${\sigma }_{yx}$ are represented using matrices for two-dimensional plane stress problems such as:

$$\left[\sigma \right]=\left[\begin{array}{cc}{\sigma }_{xx}& {\sigma }_{xy}\\ {\sigma }_{yx}& {\sigma }_{yy}\end{array}\right],$$

(35)

where $\left[\sigma \right]$ is not symmetric for most micropolar peridynamic problems. One way to ensure symmetry is by defining ${\overline{\sigma }}_{xy}=({\sigma }_{xy}+{\sigma }_{yx})/2$; The principal stress can be calculated with Equation (35), and different failure criteria can be applied to plane stress problems.

For comparison purposes, a two-dimensional FEM single-processor computer code was also developed. In this model, quadrilateral elements of one unit of measure were used with three by three Gauss points for each element. Examples using these schemes are explained in the following subsections.

2.2.1. Circular Hole Plate Subjected to Tensile Stress

A concrete plate with a circular hole in the center subjected to tensile stress was simulated with the Micropolar Peridynamic Model and the Finite Element Model.

The geometry and boundary conditions are described in Figure 4. The properties of the material are the elastic modulus $E=18,000$ MPa and the Poisson ratio $\nu =0.20$, the thickness of 50 mm, the mesh size of $150\times 300$ = $45,000$ particles and a material horizon of $\delta =3$ mm. For the FEM code $150\times 300$ = $45,000$, quadrilateral elements were also used. The applied tensile stress $\sigma$ was $8.0$ MPa.

Figure 5 show plots of the stress field in MPa (${\sigma }_{yy}$) using the Finite Element Model and the Micropolar Peridynamic Model, respectively. Figure 6 shows the stress plots at the right edge of the circular hole ($y=150$mm) using the last two models and the analytical solution (CCM). These plots show pretty good agreement between MPD and CCM stress.

2.2.2. Single Edge Notched Plate Subjected to Tension

A single edge notched concrete plate subjected to plane stress loading was simulated using the models described above. The geometry of the model and the boundary conditions are described in Figure 7. The Modulus of Elasticity is $E=18,000$ MPa and the Poisson’s ratio $\nu =0.20$. The thickness of the plate is 50 mm and the material horizon is 3 mm. The number of particles and quadrilateral elements was 150x300 for the finite element and the micropolar peridynamic model, respectively.

Figure 8 show plots of the stress field in MPa (${\sigma }_{yy}$) with the use of FEM and MPD models, and Figure 9 shows a comparison of the same stress in the middle of the plate, using MPD and FEM.

2.3. Damage Model for Concrete

A concrete damage model was implemented using the single processor code described above. Based on Figure 10 and an explicit dynamic relaxation method, an implicit linear elastic model was used to find the load required to break the first bond. Subsequently, an explicit nonlinear explicit dynamic relaxation model was used to find solutions of Equation (26). A nonlinear bond-based damage function for concrete was obtained by calculating the forces between the particles i and j using a microelastic damage function. Nonlinear microelastic damage models were previously developed in [9] where the connection force changed if the maximum stretch exceeded some prescribed values. Figure 11 shows an example of a microelastic damage model for concrete structures, where softening behavior was considered. The force function is a function of the stretch s between the particles i and j as well as the stress $\sigma$ at node i. In this model, the link remains linearly elastic as long as its stretch s exceeds the compressive stretch limit $s_c$ or does not exceed the tensile stretch limit $s_t$. The tensile peridynamic force remains constant after this value and before $s_c$, and the compressive force is constant. If the value of ${\alpha }_t{\sigma }_t$ exceeds another specified value, the force at node i decreases to zero. Similarly, if ${\alpha }_c{\sigma }_c$ decreases to some other value, the force on the particle i also falls to zero. In this case, ${\sigma }_t$ and ${\sigma }_c$ are the maximum tensile stress and the compressive strength at the point i, respectively, and ${\alpha }_c$, ${\alpha }_t$ are the calibration parameters. These factors $\alpha$ need to be greater than one because once one link is broken, the peridynamic stress should not consider the force of that broken link. However, considering the decrease in stress due to broken links, more computing time will be required for subsequent calculations. Thus, the numerical value of the peridynamic stress with the calibration factors was approximately equal to the real value.

For each simulation, the mesh size, applied forces, material parameters, and dynamic relaxation parameters were chosen. In an implicit linear module, displacements and rotations $\left[d\right]$ were obtained by inverting the global stiffness matrix $\left[K\right]$:

$$\left[u\right]={\left[K\right]}^{-1}\left[f^{ext}\right],$$

where $\left[f^{ext}\right]$ is the applied force. The initial displacements in x and y were obtained using the expression mentioned above ($\left[u\right]$,$\left[v\right]$) with rotations and stretches. The critical link was localized and the load factor was calculated to break this link. The load was increased using a prescribed load increment factor, and the peridynamic forces and moments were calculated. In addition, the stretch intervals were recalculated using this increase in load. Once the initial loads and displacements are obtained, the stress tensor is calculated at each point. The process of computing the stress tensor is described in Section 2.2

Using the Mohr failure criterion, if the maximum or minimum stresses (${\sigma }_1$,${\sigma }_3$) fall outside the failure surface (Figure 12), all peridynamic links that attach the particle i drop to zero (Figure 11).

Once a peridynamic link fails, the peridynamic stress changes. For efficiency purposes, the stiffness of all the links attached to the particle i is not changed; therefore, the calibration factors ${\alpha }_c$ and ${\alpha }_t$ should take into account the decrease in stress due to progressive link failure.

As mentioned in Section 1, many peridynamic studies treat the material horizon as a free parameter. This parameter is generally expressed as a function of the modulus of Elasticity, the stretch, the maximum stress, and the fracture energy ([10,39]), in which the material horizon for a typical concrete was on the order of 800 mm.

Using the equations described in Section 2.1, the material horizon can be expressed as a function of the fracture toughness, the maximum tensile stress and the Poisson’s ratio. Therefore, if failure occurs ($K_I$ = $K_{Ic}$), the material horizon $\delta$ can be expressed, using the maximum tensile stress of equations (20) to (23).

$$\delta =f_p{\left({\displaystyle \frac{K_{Ic}}{{\sigma }_t}}\right)}^2,$$

(36)

where $f_p$ is a function of $\pi$ and the Poisson’s ratio $\nu$, for example, for a Poisson’s ratio of 0.22, the value of $f_p$ is 0.05618899.

Application to specific problems are discussed in the following Section.

3. Results

In this section, plane stress simulations are depicted using plain concrete problems and compared with physical results from the laboratory.

3.1. Double Edge Notched Specimen Subjected to Uniaxial Tension

A double-edge-notched-plain concrete specimen subjected to tensile stress was analyzed using the method described in Section 2.2. These plates were modeled using extended FEM in [42], and some experiments were carried out by [43] in 1991. Using the material properties reported in [42], the concrete modulus of elasticity was 18.0 GPa, ${\sigma }_c$ = 30 MPa, and ${\sigma }_t$ = 3.2 MPa, geometry and boundary conditions are shown in Figure 13, and the following equation was used to obtain $K_{Ic}$([44]):

$$K_I=f\left({\displaystyle \frac{a}{W}}\right)\sigma \sqrt{\pi a},$$

(37)

$$f\left({\displaystyle \frac{a}{W}}\right)={\displaystyle \frac{1.122-0.561(a/W)-0.205{(a/W)}^2-0.190{(a/W)}^4}{\sqrt{1-(a/W)}}},$$

where W is the semi-width of the plate, a is the depth of the notch, and $\sigma$ is the applied maximum stress of 3.2 MPa. Thus, for a thickness of 50 mm, the value of $K_{Ic}$ was 0.45 MPa·m^1/2. To model the notches, the material properties were significantly weaker than the properties of the concrete material. Using Equation (36), where an assumed value of the Poisson ratio of $0.22$ was incorporated, the value of the material horizon $\delta$ is 0.0013 m, which is approximately 1.5 mm. According to Section 2.2, $\delta$ is three times the grid spacing; therefore, the model grid spacing is 0.5 mm, so a mesh size of $120\times 250$ was used. The tension ($s_t$) and compression ($s_c$) stretches were 0.00018 and 0.0018, respectively. After calibrating the values of $\alpha$, the values of ${\alpha }_c$ and ${\alpha }_t$ were 5.0 and 2.5, respectively.

Figure 14 shows the crack trajectories of this plate using MPD and Figure 15 shows the stress-elongation curves obtained using this model, compared to the experiment reported in [43].

3.2. Four Point Single Edge Notched Beam Subjected to Flexure

A single edge-notched four-point loaded plain concrete beam was analyzed with the method described in Section 2.2. Notched beams with four-point loads have been modeled using both the finite element method and peridynamics in [42,45]. The properties of concrete were $E=$35.0 GPa, ${\sigma }_t$ = 3 MPa, and ${\sigma }_c$ = 45 MPa. Using a Poisson’s ratio of 0.2 and assuming a typical fracture toughness value of 1.4 MPa·m^1/2, the value of $\delta$ according to Equation (36) is 3.3 mm, so the grid spacing for this model was 1.0 mm. Figure 16 shows the configuration of the problem. After calibrating the values of $\alpha$, the values of ${\alpha }_c$ and ${\alpha }_t$ were 2.0 and 1.5, respectively.

Figure 17 and Figure 18 show the crack trajectory and the load deflection curve of this model using the scheme shown above. This path coincides with those found in [42] and[45] and is in agreement with the experimental results reported by [42].

3.3. Three-Point Single Edge-Notched Beam Subjected to Flexure

Based on the methodology described above, a three-point plain concrete beam was simulated using micropolar peridynamics and compared with laboratory results. Geometry and boundary conditions are shown in Figure 19. According to laboratory tests, the concrete modulus of elasticity of 21.1 GPa, ${\sigma }_c$ = 20 MPa, ${\sigma }_t$ = 2.0 MPa and with the help of the following equation ([41]):

$$K_I={\displaystyle \frac{P}{B\sqrt{W}}}f\left({\displaystyle \frac{a}{W}}\right),$$

(38)

$$f\left({\displaystyle \frac{a}{W}}\right)={\displaystyle \frac{3S}{W}}\sqrt{{\displaystyle \frac{a}{W}}}{\displaystyle \frac{\left[1.99-\frac{a}{W}\left(1-\frac{a}{W}\right)\left(2.15-3.93\frac{a}{W}+2.7{\left(\frac{a}{W}\right)}^2\right)\right],}{2\left(1+2\frac{a}{W}\right)\sqrt{{\left(1-\frac{a}{W}\right)}^3}}}$$

(39)

where W is the width of the beam, S is the span, a is the depth of the notch, B is the thickness of the beam (150 mm) and P is the maximum applied load (11.7 kN), the value of $K_{Ic}$ was 0.51 MPa·m^1/2. Employing Equation (36) with a Poisson’s ratio of 0.22, the value of the material horizon $\delta$ is 0.0033 m, which is approximately 3 mm, and in this case the model grid spacing is 1 mm. The tension ($s_t$) and compression ($s_c$) stretches were 0.0001 and 0.001, respectively, and the values of ${\alpha }_c$ and ${\alpha }_t$ were 5.0 and 2.0, respectively.

Figure 20 shows the crack path of this beam using the aforementioned model and Figure 21 shows the crack of the three point loaded beam after a physical test. This path coincided with those previously reported [42] and in the laboratory. Figure 22 shows the load-displacement curve obtained in the laboratory compared to the one obtained with the model.

4. Discussion

In this study, a novel micropolar peridynamic model was proposed based on stress and stretch failure criteria to model concrete structures. The main conclusions are provided as follows:

1.

On the basis of the definition of the peridynamic stress tensor and the original micropolar peridynamic model, a one-to-one relationship between the conventional constitutive stress-strain tensor and the micropolar peridynamic stress-strain tensors for linear elastic materials was obtained.

2.

This framework was able to evaluate the stress for discontinuous displacement fields. For linear elastic problems, the results agree with the CCM solutions.

3.

Using displacement fields from LEFM, stresses at the crack-tip were expressed as functions of the stress intensity factor, the material horizon, and the Poisson’s ratio.

4.

Based on the proposed expressions of the micropolar peridynamic stress at the crack-tip, the material horizon was expressed as a function of three key parameters: fracture toughness, maximum tensile stress, and the Poisson’s ratio.

5.

A numerical matrix-based scheme was implemented to model concrete structures based on an explicit dynamic relaxation model. A nonlinear bond-based damage function was implemented for concrete structures, where the force that connects the particles changes if the maximum stretch exceeds some prescribed values. In contrast to other bond-based models, the link forces in this model drop to zero if the maximum tensile stress or the maximum compressive stress exceeded some threshold values for tension and compression.

6.

To verify the performance of the model, several concrete structures subjected to plane stress were analyzed and compared with the physical results. The tension of a double edge notched specimen and the flexure of a single edge notched for a three- and four point loaded plain concrete beam were simulated using the proposed numerical scheme. Compared with the crack trajectories and experimental curves obtained from the physical tests, the results agree well with the model.

In the future, large deformation, three-dimensional models and connection with finite element models will be examined using the proposed framework.