Integration of Peridynamic Theory and OpenSees for Solving Problems in Civil Engineering

Peridynamics (PD) is a powerful method to simulate the discontinuous problems in civil engineering. However, it may take a lot of effort to implement the material constitutive models into PD program for solving a broad range of problems. OpenSees is an open source software which includes a versatile material library and has been widely used by researchers and engineers in civil engineering. In this context, the paper presents a simple but effective approach to integrate PD with OpenSees by using a Client-Server (CS) software integration technique, such that the existing material constitutive models in OpenSees can be directly used by PD. Two applications are presented to verify the new PD-OpenSees platform. The first one is a plate with a pre-crack subject to horizontal load simulated using a three-dimension (3D) multi-yield-surfaces plasticity model, and the second one is a concrete block with a rectangular hole subject to a uniaxial loading condition simulated using a 3D Cap plasticity model. It shows that the generation/ propagation of cracks of the elastoplastic materials (e.g., concrete and soil) can be analyzed by combing PD and three-dimensional plasticity constitutive models without compiling/ linking these complex material models into PD program. Thus, the integrated PDOpenSees platform presented herein is potentially capable to solve a wide range of complex problems in civil engineering.


Introduction
Peridynamics (PD) is a mechanic theory presented by Silling and co-workers [Silling (2000)] that extends classical continuum mechanics to include discrete particles and growing cracks. The continuum theory is based on the partial differential equations, which requires sufficient smooth deformation field. However, the deformation at the region with cracks or displacement discontinuities is singular and the partial derivatives are not available, which leads to inaccuracy or incapability of the continuum theory. The cracks or displacement discontinuities are commonly encountered in the problems of civil engineering, e.g., reinforced concrete components. The PD theory is a practically useful method for these problems. The PD theory follows a set of integral equations [Silling (2000)] which does not require a continuous response fields or continuous material space [Silling, Epton, Weckner et al. (2007); Silling and Lehoucq (2010)]. The internal force acting on a particle is assumed to be only from the particles that are close enough, i.e., inside a so-called material horizon. There are two types of PD, i.e., the bond-based PD and state-based PD. In the bond-based PD, the pairwise force between a particle of interest and a 'pair' particle within its material horizon is calculated based on their relative locations in the reference and current configurations. The interaction force is assumed to be along the direction between the two particles and the magnitude can be calculated using a one dimensional (1D) constitutive model. In the bond-based PD, the Poisson ratio is fixed (i.e., one third or one quarter) for any material, and the plastic incompressibility cannot be considered appropriately [Silling (2000); Gerstle (2015)]. In the state-based PD, a deformation gradient at a particle is calculated based on the information of all particles inside its material horizon, i.e., their relative locations in both undeformed and deformed configurations. A force state may be calculated using 3D material constitutive model based on the deformation. The interaction force vector between two particles does not need to be equal in magnitude and opposite to each other, while they must satisfy the balance of angular momentum. The Poisson ratio can be prescribed for various materials and is not necessary to be fixed [Silling, Epton, Weckner et al. (2007)]. Apparently, the application of the state-based PD is more general than that of the bond-based PD. Some elastoplastic material constitutive models (e.g., Drucker Prager model [Lai, Ren, Fan et al. (2015); Song, Yu and Kang (2018)], Johnson-Holmquist 2 model [Lai, Liu, Li et al. (2018)]) are implemented in the state-based PD for various engineering problems. However, it may still take a lot of effort for researchers to implement new material constitutive models into the PD program. Considering there are abundant material models in the current FEM software, it is interesting to take advantages of these existing models directly and integrate the PD theory and these FEM software platforms in a simple way. OpenSees (the Open System for Earthquake Engineering Simulation) is an open source software developed since 1997 by the pacific earthquake engineering research (PEER) center, which is widely used in simulating the seismic structural or geotechnical responses [Mazzoni, McKenna, Scott et al. (2006);McKenna (2011)]. In the past two decades, OpenSees keeps integrating various material, element, section models and algorithms developed by researchers from all over the world. There are various 1D concrete material models (e.g., the uniaxial Kent-Scott-Park material model [Kent and Park (1971)]), 1D steel models (e.g., the uniaxial bilinear material model), 3D concrete models (e.g., Cap plasticity model [Sandler, Dimaggio and Baladi (1976); Sandler and Rubin (1979);Hofstetter, Simo and Taylor (1993)]) and 3D soil models (e.g., the bounding surface hypoplasticity model [Wang, Dafalias and Shen (1990)] and the multi-yield-surface J2 plasticity model [Gu, Conte, Elgamal et al. (2009)]). OpenSees has included a versatile library and been increasingly widely used by researchers and engineers in civil engineering. A simple but powerful tool command language, Tcl, developed by John Ousterhout and coworkers [Ousterhout and Jones (2009)], is used by OpenSees to create the FE models, analyze the models and get the model responses. However, it is not easy to directly compile and link PD program with OpenSees, since they are both complicated software and are written in Fortran and C++ respectively. In this paper, a simple but effective approach is proposed to integrate PD with OpenSees by using a Client-Server (CS) software integration technique presented by the authors and co-workers [Gu and Ozcelik (2011)]. Based on the CS technique, OpenSees is modified using Tcl to be a "material server" that remains in the memory, handles a material object, waits to receive a strain from client, calculates the corresponding stress based on the material constitutive model and sends it back to client. A socket connection between PD client and OpenSees server is achieved by using a simple Tcl network communication channel (or socket) based on TCP/IP or other network protocols (e.g., UDP). Each material point at a particle in the PD program is a client, which keeps the connection with the material server, sends strain to server and receives stress from the server. Using this approach, the PD program and OpenSees are seamlessly integrated without extra efforts of compiling and linking the two programs. All existing material constitutive models in OpenSees can be directly used by PD. The integration is robust by using the network socket and efficient due to the small data transfer between PD and OpenSees.
Two application examples are presented to verify the PD-OpenSees platform implemented by using the CS technique for solving civil engineering problems. The first one is a plate with a pre-crack subject to horizontal load. The material is simulated using a 3D multiyield-surfaces J2 plasticity model, which is commonly used in geotechnical engineering for simulating behaviors of the clay soil. The second example is a concrete block with a rectangular hole subject to uniaxial loading condition. The 3D Cap plasticity model is used to simulate the concrete behaviors. The examples demonstrate that the generation/propagation of cracks of elastoplastic materials can be simulated efficiently by integrating PD and OpenSees using complex 3D plasticity constitutive models, which may be extremely valuable in solving a wide range of complex problems in civil engineering.

Integration of PD and OpenSees based on CS technique
In this section, the basic theory of PD and two commonly used 3D constitutive models in OpenSees used for the applications are illustrated briefly. An integration method based on the CS technique is reviewed to couple PD and OpenSees.

The PD theory
There are two types of PD, i.e., the bond-based PD and the state-based PD. The equation of motion for both types of PD has the format as, where ρ is the mass density in the reference configuration, u  is the acceleration of particle XA, t is time, b is the body force density (refer to Fig. 1 is a pairwise force that the particle xB exerts on particle xA. As shown in Fig. 1 . The parameter ξ represents the relative position of the two particles in the reference configuration, and η represents their relative displacement in the current configuration, i.e., AB represents the relative position vector in the current configuration. B V is a volume of material particle XB.
The force vector f can be calculated as follows: where s is the stretch (scalar) between the two particles, which is defined as can be calculated by a 1D material constitutive model. A "horizon" δ is used to define the neighborhood of PD particle XA (see Fig. 1).

Figure 2: Deformation of PD particles in state-based PD
In the state-based PD theory, the second term in Eq. (2) is calculated from an integration: where A,B t is called the force density vector as shown in Fig. 2, which can be viewed as the force that particle xB exerts on particle xA, which are calculated using the formula: where KA is the shape tensor of PD particle XA, which is defined in the reference configuration as In which ( ) ω ξ is the influence function. The first Piola-Krichhoff stress PA in Eq. (5) is calculated as where SA is the second Piola-Krichhoff stress and FA is the deformation gradient. The stress SA can be calculated based on a 3D constitutive material model, i.e., in which EA is the deformation of the particle calculated as where F T is the transformation matrix of FA and I is the identity matrix. The key step in the state-based PD is the calculation of the deformation gradient FA used in Eqs. (7) and (9). According to the literatures [Silling, Epton, Weckner et al. (2007); Gerstle (2015)], FA can be calculated by where N is a second order tensor representing the shape of the deformed configuration: , which represents the deformation state.

The material constitutive models in OpenSees
The material constitutive in OpenSees can be described as: where σ , ε and D are the stress tensor, the strain tensor and the tangent stiffness matrix at integration point of FE meshes, respectively. There are various material constitutive models in OpenSees. In this paper, two models are used in the demonstration examples, i.e., a 3D multi-yield surface (MYS) model and a 3D cap plasticity model which are usually used for modeling the soil and concrete materials, respectively. The two material constitutive models are reviewed briefly in the subsections.

Multi-yield-surface (MYS) plasticity model
Multi-yield-surface (MYS) plasticity model is proposed by Iwan [Iwan (1967)], Mroz [Mroz (1967)] and Prevost [Prévost (1977)] for simulating metal and soil materials. It is further developed, enhanced and implemented into OpenSees by Elgamal and co-workers [Elgamal, Yang and Parra (2002); Yang and Elgamal (2002); Yang, Elgamal and Parra (2003)] for study of complicated soil behaviors. Compared with the classic J2 plasticity model which has only one yield surface, the MYS model has a number of similar yield surfaces to represent a better plastic behavior of the soil model under cyclic loading conditions.

The backbone curve and multiple yield surfaces
In geotechnical engineering, the octahedral shear stress strain behavior of soil material can be defined by a backbone curve which is usually determined by experiments, i.e., where G, τ and γ are the low strain shear modulus, the shear stress and the shear strain, respectively. The parameter r γ is a reference shear strain. In the context of MYS, the smooth hyperbolic backbone curve is replaced by a piecewise linear approximation (see Fig. 3). Each line segment corresponds to one yield surface.

Figure 3: The yield surfaces of the MYS model
Each yield surface is defined in the deviatoric stress space as, in which τ , α and K are the deviatoric stress tensor, the back-stress tensor and the size of the yield surface, respectively, as shown in Fig. 3(a). The subscript i denotes the i th yield surface, and NYS is the number of outermost yield surface. Fig. 3(c) shows the yield surfaces at the initial loading time, which are concentric cylindrical surfaces in the deviatoric stress space. The outermost yield surface, i.e., NYS f =0, is the failure surface which defines a geometrical boundary.

Flow rule and hardening law
An associative flow rule is defined in the deviatoric stress space [Gu, Conte, Elgamal et al. (2009)] as follows: where p dε is the plastic strain increment. Q represents the plastic flow direction normal to the yield surface.
: L d = Q τ is the plastic loading function and the symbol 〈〉 represents the MacCauley's brackets, i.e., max(L,0) H is a constant plastic shear modulus. A pure deviatoric kinematic hardening rule is employed to generate the Masing-type hysteretic cyclic response under arbitrary cyclic shear loading conditions ; Gu, Conte, Yang et al. (2011)]. The yield surfaces translate in stress space within the failure envelope [Hill (1998)] and there is no overlapping between any two yield surfaces. The details can be found in the literature [Gu, Conte, Elgamal et al. (2009)].

Cap plasticity model for concrete
Cap plasticity model is defined by a non-smooth, convex surface which consists of three yield surfaces, i.e., a failure envelope, a tensile cut-off region and a hardening cap. It is proposed by FL DiMaggio and IS Sandler in 1971 [DiMaggio andSandler (1971) ;Sandler, Dimaggio and Baladi (1976); Sandler and Rubin (1979)] and further developed by Hofstetter [Hofstetter, Simo and Taylor (1993)].
where κ is a hardening parameter and where T is the maximum hydrostatic tension.

Flow rule and the hardening law:
The flow rule of the model and the consistency condition (i.e., the Kuhn-Tucker condition) are defined as where γ i  is the plastic consistency parameter. The hardening behavior of the Cap model is illustrated by a relationship between the hardening parameter κ and the effective plastic  [Sandler, Dimaggio and Baladi (1976); Sandler and Rubin (1979)].

The integration method
The integration of PD program and OpenSees is achieved based on the Client-Server software technique. The flowchart of PD, OpenSees and their interaction are illustrated in Fig. 6.  As shown in Fig. 6, the procedure of the PD-OpenSees software platform is as follows: Step 1: Build the PD model using PD program and create an OpenSees server.
Step 2: For each particle in PD model, build a corresponding material object in OpenSees server side. A communication channel between PD and OpenSees is established using Tcl socket command. As shown in Fig. 6, the material objects keep in the memory and wait for commands from the PD.
Step 3: In the PD program, for each time step, the shape tensor K, the second order tensor N, the deformation gradient F and the strain tensor E are calculated for each particle using the Eqs. (6), (11), (10) and (9), respectively.
Step 4: PD sends the strain E to OpenSees through the socket.
Step 5: OpenSees receives the strain, calls its material objects to calculate the stress using Eq. (8). The strain, stress and other variables are committed for the calculation of the next time step.
Step 6: OpenSees sends the stress back to PD program, then waits in memory for the next commands (e.g., strain) from PD program.
Step 7: PD program receives the stress, calculates the PK-I stress and the force vector state T using Eqs. (7) and (5), respectively.
Step 8: PD program calculates the acceleration of the PD particle using Eq. (1), updates the velocity and displacement, and then goes to Step 3 for the next time step. The above steps from 1 to 8 are repeatedly performed until the last analysis step is done.

Application examples
Two examples are presented herein to verify the integrated PD-OpenSees platform. The first example is a plate with a pre-crack subject to horizontal tension force, and the material is simulated using the MYS model. The second one is a block with a rectangular hole subject to uniaxial loading condition simulated using the cap plasticity model.

Plate with a pre-crack
The first example presented herein is a plate with a pre-crack subject to horizonal tension force. As shown in Fig. 8, the PD model is discretized into 8910 particles with 2 layers. The length and width of the plate are both 1.89 m, the grid size dx is 0.03 m, the time step dt is 0.000004 s, the critical stretch 0 s is 0.05 and the horizon radius δ is 3.015 times the grid size (i.e., 0.09045 m). The material is simulated using the MYS model with parameters shown in Tab. 1. The computation results are analyzed and compared with the one using linear elastic material model whose parameters are also shown in Tab. 1. The plate is subject to tension forces at two boundaries, as shown in Fig. 8 where xx σ , yy σ , zz σ , xy σ , xz σ and yz σ are the components of the stress σ calculated using Eq. (12).     9 shows the load time history, the deformation and the contour of octahedral shear stress at three representative time steps using the MYS model. In Fig. 9(a), the force increases linearly from zero until it reaches 750 kN at 0.02 second and remains unchanged after then. The stress and deformation of the plate at three representative time steps (i.e., A, B and C corresponding to 0.036 s, 0.048 s and 0.0792 s respectively) are shown in Figs. 9(b)-9(d). For comparison purpose, the load time history, the deformation and the contour of octahedral shear stress at three representative time steps (i.e., A, B and C corresponding to 0.036 s, 0.074 s and 0.0912 s respectively) using the elastic model are shown in Figs. 10(a)-10(d), where the load increases linearly until reaches 5000 kN at 0.02 second and keeps unchanged afterwards. From Fig. 9 and Fig. 10, it is observed that when using MYS model, the stress at the crack tip region grows fast along approximate 45 degrees to the pre-crack, then grows slower after yielding. However, the crack does not obviously extend or propagate along any direction but concentrates in a small region around the crack tip as shown in Fig. 9(d). While using elastic model, the stress increases fast at the crack tip at beginning, similarly with the use of MYS models in Fig. 9, until failure occurs. However, the force at failure is much larger than that for MYS case (i.e., 5000 kN for elastic case instead of 750 kN for MYS case) due to the much larger tension resisting capacity of the elastic material. After then the crack propagates along the pre-crack direction and finally penetrates the plate as shown in Fig. 10(c) and Fig. 10(d).

Figure 11:
Comparison of the crack opening displacement with increasing tension force between MYS model and elastic model Fig. 11 compares the crack opening displacement with increasing tension force between MYS vs. elastic models for the first 5000 loading steps. The original crack opening displacements are 0.06 m for both cases. It is observed that the crack opening begins to propagate when the tension force is 126 kN for MYS model, which is much smaller than the tension force of 1400 kN for elastic model. When the tension force is 750 kN, the crack opening displacement for MYS model reaches 0.0628 m, while that for elastic model remains unchanged, i.e., 0.06 m. This is also because that the elastic material has much larger tension resisting capacity than plastic material.

Figure 12:
The octahedral shear stress-strain behavior of the PD particles when using MYS model The octahedral shear stress-strain behaviors of two representative PD particles (i.e., PD particles m and n) are presented in Fig. 12(a) and Fig. 12(b) respectively. From these figures, it is observed that the calculated octahedral shear stress at the two particles increase along the prescribed hyperbolic backbone curve, which verifies the integrated PD-OpenSees platform presented herein.

A concrete block with a rectangular hole subject to compressive force
The second application is a concrete block with a rectangular hole in the center subject to compressive force. As shown in Fig. 13, the length, width and height of the block are all 0.45m, the side length of the hole inside the block is 0.15 m. The PD model is discretized into 3270 particles with the grid size dx (i.e., the initial particle distance) 0.03 m as shown in Fig. 13. The time step dt is 0.000004 s, the critical stretch 0 s is 0.02, and the horizon radius δ is 0.09045 m. The compressive force is applied at both the top and bottom of the block, and the load time history is shown in Fig. 14(a). The Cap plasticity model is applied in this example and the material parameters are shown in Tab Fig. 14(a)). The stress increases fast at the left and right sides of the hole as shown in Fig. 14(b). After then the brick tends to be compressed in z direction while expanded in x direction due to the Poisson effect (see Figs. 14(c)-14(d)). Vertical cracking occurs at the left and right side of the rectangular hole and grows wider with increasing loads.

Figure 15:
The stress-strain behavior of two representative PD particles Fig. 15 shows the stress-strain behaviors of two representative PD particles p and q. At the beginning of the loading steps, the stress point moves inside the yield functions (i.e., in the elastic region). With increasing loading, the stress point reaches the yield surface of the tensile cut-off region (i.e., 3 0 f = in Fig. 4) at first, and then the yield surface of the failure envelop (i.e., 1 0 f = ). Finally, the stress point reaches the yield surface of the cap ((i.e., 2 0 f = ), which expands the cap outward (i.e., move towards right in the Fig. 15(a) and Fig.   15(b)) with the hardening parameter κ enlarged. It is observed that the results follow the definition of the yield surfaces and hardening law of the Cap plasticity model in Section 2.4, which verifies the integrated PD-OpenSees platform presented herein.

Conclusion
Peridynamics (PD) is a powerful method to solve the discontinuous problems, which is commonly encountered in civil engineering. However, it may take a lot of effort to implement the various material constitutive models into PD program. OpenSees, an open source finite element software, has been widely used in the civil engineering which has versatile material libraries. This paper integrates PD with OpenSees by using a Client-Server (CS) software integration technique. By combining the state-based PD program with 3D plastic constitutive models, the integrated method is capable to solve the complex generation and propagation processes of cracks of elastoplastic material efficiently. Two applications are presented to verify the integrated PD-OpenSees software platform. The first one is a plate with a pre-crack simulated using the MYS model subject to horizontal boundary tension forces. It is observed that the stress at the crack tip region grows fast along approximate 45 degrees to the pre-crack, then grows slower after yielding occurs. However, the crack does not obviously extend or propagate along any direction but concentrates in a small region around the crack tip. This is different from the results using elastic model, in which the crack propagates along the pre-crack direction and finally penetrates the entire plate. The second one is a concrete block with a rectangular hole simulated using the Cap plasticity model subject to vertical compression force. Results show that the brick tends to be compressed in z direction while expanded in x direction with increasing loads. In both application examples, the stress-strain behaviors of representative PD particles follow the definition of the plasticity models in the sense of yield surfaces, flow rule and hardening law, which indicates that the integrated PD-OpenSees software platform presented herein is potentially capable to solve a wide range of problems in civil engineering.