Effect of seismic wave propagation in massed medium on rate-dependent anisotropic damage growth in concrete gravity dams

Abstract

Seismic modeling of massive structures requires special caution, as wave propagation effects significantly affect the responses. This becomes more crucial when the path-dependent behavior of the material is considered. The coexistence of these conditions renders numerical earthquake analysis of concrete dams challenging. Herein, a finite element model for a comprehensive nonlinear seismic simulation of concrete gravity dams, including realistic soil–structure interactions, is introduced. A semi-infinite medium is formulated based on the domain reduction method in conjunction with standard viscous boundaries. Accurate representation of radiation damping in a half-space medium and wave propagation effects in a massed foundation are verified using an analytical solution of vertically propagating shear waves in a viscoelastic half-space domain. A rigorous nonlinear finite element model requires a precise description of the material response. Hence, a microplane-based anisotropic damage–plastic model of concrete is formulated to reproduce irreversible deformations and tensorial degeneration of concrete in a coupled and rate-dependent manner. Finally, the Koyna concrete gravity dam is analyzed based on different assumptions of foundation, concrete response, and reservoir conditions. Comparison between responses obtained based on conventional assumptions with the results of the presented comprehensive model indicates the significance of considering radiation damping and employing a rigorous constitutive material model, which is pursued for the presented model.

Graphical abstract

Keywords

soil–structure interaction / massed foundation / radiation damping / anisotropic damage

Introduction

Concrete dams are one of the most important infrastructures in every country. Owing to their social and economic effects, their seismic safety assessment is a primary task in dam engineering. The dam–foundation–reservoir system involves a complex interaction as it comprises three domains with different behaviors. Hence, for an accurate simulation, some important factors should be considered in each domain. An appropriate representation of the nonlinear concrete response, geometrical nonlinearity of the body, and damping mechanism should be considered for the appropriate modeling of the dam. To provide a suitable representation of the far-field domain, various factors including foundation mass, radiation damping, seismic wave propagation effect, appropriate earthquake input mechanism, and absorbing boundary conditions are required. The compressibility of reservoir water, absorption of seismic waves by sedimentary materials at the bottom of the reservoir, and unbounded nature of the reservoir domain should be considered in the fluid domain.

Various attempts for considering such details in the seismic analysis of the dam–foundation–reservoir system have been performed [1–7]. An appropriate foundation model can be provided using a viscoelastic half-space domain, and the finite element modeling of such a domain presents some major challenges. The unbounded foundation should be truncated at a reasonable distance and an artificial boundary condition should be imposed at the truncated regions to absorb the outgoing waves; otherwise, false reflected waves are produced within the domain. The artificial boundary in combination with an appropriate earthquake input mechanism should reproduce free-field effects correctly. Hence, the implementation of the procedure is not trivial. In addition, truncated boundaries should absorb the outgoing waves with arbitrary angles and arbitrary frequencies. Owing to these challenges, most researchers used the massless foundation assumption, which was first proposed by Clough [8]. Owing to its simplicity and favorability in commercial software, this approach is widely employed in engineering practice. In this approach, the mass of the foundation is disregarded and only its flexibility is considered. Therefore, seismic waves propagate with infinite velocity in the medium, i.e., the free-field earthquake record imposed at the truncated region of the foundation reaches the base of the dam instantly. Accordingly, in the massless foundation assumption, only a finite domain of foundation with zero mass density is modeled and an earthquake free-field record is applied at the truncated boundary. Using this technique, no absorbing boundary condition is required.

Regardless of simplicity, studies have shown that the massless foundation assumption is erroneous. Tan and Chopra [9,10] and Chopra [11] concluded that the resulting stresses and displacements of a dam using a massless foundation model are overestimated. Consequently, finite element studies have been performed to include the foundation mass in analyzing the dam–foundation–reservoir system. Numerous studies have been conducted to develop effective yet simple artificial boundaries to represent the truncated far-field domain. Absorbing boundary conditions [12,13], perfectly matched layers [14–16], infinite elements [17,18] and free-field boundary conditions [19,20] are most typical artificial boundary conditions used in the finite element analysis of the dam–foundation–reservoir system. In addition to the absorbing boundary condition, the earthquake input mechanism significantly affects the reliable massed foundation model. Léger and Boughoufalah [21] investigated four typical earthquake input mechanisms and reported that some of them were unrealistic. Wang et al. [22] employed viscous-spring boundary conditions to analyze a concrete arch dam with a massed foundation and verified their results using a computer program for three-dimensional earthquake analysis of concrete dams, named EACD [23]. Bielak and Christiano [24] introduced an effective method for the earthquake input mechanism in soil–structure interaction problems based on the domain reduction method [25]. The method was developed for the seismic analysis of large-scale problems involving seismic sources, propagation paths, and local site effects. It is a useful tool for soil–structure interaction problems, in which earthquake excitation is directly applied inside the computational domain. It has been validated by Yoshimura et al. [26] and introduced as an advanced technique for considering both earthquake input mechanism and absorbing boundary condition appropriately.

Meanwhile, the appropriate seismic simulation of a concrete dam requires an authentic representation of the nonlinear concrete response. As the seismic loading of such massive structures, even in moderate or small magnitudes, is accompanied by concrete cracking, disregarding nonlinear behavior will yield unreliable results. Many researchers have attempted to consider nonlinearities due to concrete cracking in various manners, including the use of linear elastic fracture mechanics [27–29], nonlinear fracture mechanics [30–32], the fictitious crack model [33], the crack band model [34], and the smeared crack method [35].

Two major approaches for resembling fracture processes are discrete and continuous approaches. Damaged regions are of measure zero in the former, whose topology is defined either explicitly or implicitly. Remeshing techniques [36–38] and the extended finite element method [39–41] represent cracks and localized bands explicitly, whereas the cracking node method [42] and peridynamics [43] employ implicit representations. Meanwhile, the finite thickness of bands can be defined by regularizing the partial differential equations of the problem [44–46]. This enables a continuous modeling via continuum damage mechanics [47,48], either with [49] or without the plasticity of a material [50]. Meanwhile, the phase field model [51,52] introduces a length-scale into the model using a scalar field known as the phase field [53]. Employing an additional differential equation, crack propagation is tracked naturally and without defining supplementary criteria [54]. This renders phase-field models suitable for simulating multiple interacting cracks, merging of cracks, and crack branching [55]. Therefore, phase-field models have garnered extensive attention in the field of fracture mechanics. Recent publications regarding these models include the phase-field modeling of dynamic cracking [56], fracture in porous media [57,58], hydraulic fracture [59], and fracture based on neural networks [60].

Among continuum damage models, tracking nonlinear concrete response using the damage–plastic formulation is advantageous compared with the other models, as both stiffness degradation and permanent deformation are reproduced satisfactorily. However, a rigorous plastic damage model must account for different tensile and compressive behaviors of concrete, closing and reopening of cracks, volumetric expansion due to compression, and damaged-induced anisotropy. In this study, the nonlinear behavior of concrete was simulated using the damage–plastic model of Daneshyar and Ghaemian [61]. The model can reproduce damaged-induced anisotropy, which results in distinct degenerations in different directions, in conjunction with irreversible deformations due to excessive straining of concrete. Induced anisotropy is defined using the history of directional damage growth by defining two distinct tensors in the space of tensile and compressive damage. Owing to the different descriptions of damage in the two states, the dissimilar behaviors of concrete in tension and compression as well as microcrack closing and reopening are obtained.

As the direct finite element modeling of massed half-space medium must be performed well to impose the appropriate earthquake input and absorbing boundaries simultaneously, radiation damping is primarily disregarded in the seismic safety evaluation of concrete dams. Radiation damping is either considered using the boundary element technique in the frequency domain, or completely disregarded in the favor of nonlinear modeling. In this study, a finite element model based on the domain reduction method in conjunction with a standard viscous boundary was used for the direct nonlinear analysis of a dam-massed foundation–reservoir system, and the nonlinear response of concrete was monitored using a novel microplane-based rate-dependent anisotropic damage–plastic model. To facilitate understanding, the domain reduction method (DRM) formulation is described in Section 2. To verify the model, a numerical analysis was performed in the elastic half-space domain, and the result was compared with analytical data. Another verification was performed by comparing the analytical deconvolution results against the response of the viscoelastic half-space domain. In Section 3, the fundamental assumptions and formulation of the anisotropic damage-plastic model of concrete are reviewed, and the verification of the model using available tests is presented. Finally, a dam–foundation–reservoir system is analyzed based on different assumptions of massless/massed foundation, linear/nonlinear concrete response, and empty/full reservoir conditions. Remarks and conclusions are presented as well.

DRM

The DRM is a two-step procedure that was developed for the seismic analysis of large-scale problems. The key point of the method is the separation of the far-field (seismic source) from the near-field (local site) by changing the variables of the system of equations. In the first step, a free-field analysis is performed for a half-space domain that includes the seismic source. In this step, the effective equivalent forces of a stripe of finite elements are calculated based on the free-field response. In the second step, the effective forces are applied within a reduced domain that includes the local site and structure.

Formulation

The formula of the original DRM [25] was employed in this study; however, the material damping, which was disregarded in the original work, was considered. The general framework of the dam-massed half-space foundation, including the seismic source, is illustrated in Fig. 1.

The half-space domain was truncated at {\mathrm{\Gamma }}^{+} to reduce the computational cost (Fig. 2). The computational domain was segmented into \mathrm{\Omega } and {\mathrm{\Omega }}^{+} by the virtual boundary \mathrm{\Gamma } (Fig. 3). In these figures, u_{\mathrm{i}}, u_{\mathrm{e}}, and u_{\mathrm{b}} are the nodal displacement vectors of the internal domain \mathrm{\Omega }, external domain {\mathrm{\Omega }}^{+}, and boundary \mathrm{\Gamma }, respectively

Finite element discretization of the Navier equation in \mathrm{\Omega } yields

(1)\left[\begin{array}{cc}{M}_{\mathrm{i}\mathrm{i}}& M_{\mathrm{i}\mathrm{b}}\\ M_{\mathrm{b}\mathrm{i}}& M_{\mathrm{b}\mathrm{b}}\end{array}\right]\left\{\begin{array}{c}{\ddot{u}}_{\mathrm{i}}\\ {\ddot{u}}_{\mathrm{b}}\end{array}\right\}+\left[\begin{array}{cc}{C}_{\mathrm{i}\mathrm{i}}& C_{\mathrm{i}\mathrm{b}}\\ C_{\mathrm{b}\mathrm{i}}& C_{\mathrm{b}\mathrm{b}}\end{array}\right]\left\{\begin{array}{c}{\dot{u}}_{\mathrm{i}}\\ {\dot{u}}_{\mathrm{b}}\end{array}\right\}+\left[\begin{array}{cc}{K}_{\mathrm{i}\mathrm{i}}& K_{\mathrm{i}\mathrm{b}}\\ K_{\mathrm{b}\mathrm{i}}& K_{\mathrm{b}\mathrm{b}}\end{array}\right]\left\{\begin{array}{c}{u}_{\mathrm{i}}\\ u_{\mathrm{b}}\end{array}\right\}=\left\{\begin{array}{c}0\\ P_{\mathrm{b}}\end{array}\right\},

where M is the mass matrix; C is the damping matrix; Kis the stiffness matrix; ‘i’, ‘e’, and ‘b’ are indices related to the internal, external, and boundary nodes, respectively. Analogously, for {\mathrm{\Omega }}^{+},

(2)\left[\begin{array}{cc}{M}_{\mathrm{b}\mathrm{b}}& M_{\mathrm{b}\mathrm{e}}\\ M_{\mathrm{e}\mathrm{b}}& M_{\mathrm{e}\mathrm{e}}\end{array}\right]\left\{\begin{array}{c}{\ddot{u}}_{\mathrm{b}}\\ {\ddot{u}}_{\mathrm{e}}\end{array}\right\}+\left[\begin{array}{cc}{C}_{\mathrm{b}\mathrm{b}}& C_{\mathrm{b}\mathrm{e}}\\ C_{\mathrm{e}\mathrm{b}}& C_{\mathrm{e}\mathrm{e}}\end{array}\right]\left\{\begin{array}{c}{\dot{u}}_{\mathrm{b}}\\ {\dot{u}}_{\mathrm{e}}\end{array}\right\}+\left[\begin{array}{cc}{K}_{\mathrm{b}\mathrm{b}}& K_{\mathrm{b}\mathrm{e}}\\ K_{\mathrm{e}\mathrm{b}}& K_{\mathrm{e}\mathrm{e}}\end{array}\right]\left\{\begin{array}{c}{u}_{\mathrm{b}}\\ u_{\mathrm{e}}\end{array}\right\}=\left\{\begin{array}{c}-P_{\mathrm{b}}\\ P_{\mathrm{e}}\end{array}\right\}\mathrm{.}

A simple auxiliary problem comprises a half-space without the structure, and site conditions are solved first, thereby yielding the nodal forces P_{\mathrm{b}}^0 of the boundary between the internal and external domains as well as the nodal displacements u_{\mathrm{i}}^0, u_{\mathrm{e}}^0, and u_{\mathrm{b}}^0 of the internal, external, and boundary domains, respectively (see Fig. 4).

Equations of motion for external domain {\mathrm{\Omega }}^{+} in auxiliary problem state

(3)\left[\begin{array}{cc}{M}_{\mathrm{b}\mathrm{b}}& M_{\mathrm{b}\mathrm{e}}\\ M_{\mathrm{e}\mathrm{b}}& M_{\mathrm{e}\mathrm{e}}\end{array}\right]\left\{\begin{array}{c}{\begin{array}{c}{\ddot{u}}_{}\end{array}}_{\mathrm{b}}^0\\ {\begin{array}{c}{\ddot{u}}_{}\end{array}}_{\mathrm{e}}^0\end{array}\right\}+\left[\begin{array}{cc}{C}_{\mathrm{b}\mathrm{b}}& C_{\mathrm{b}\mathrm{e}}\\ C_{\mathrm{e}\mathrm{b}}& C_{\mathrm{e}\mathrm{e}}\end{array}\right]\left\{\begin{array}{c}{\begin{array}{c}{\dot{u}}_{}\end{array}}_{\mathrm{b}}^0\\ {\dot{u}}_{\mathrm{e}}^0\end{array}\right\}+\left[\begin{array}{cc}{K}_{\mathrm{b}\mathrm{b}}& K_{\mathrm{b}\mathrm{e}}\\ K_{\mathrm{e}\mathrm{b}}& K_{\mathrm{e}\mathrm{e}}\end{array}\right]\left\{\begin{array}{c}{\begin{array}{c}{u}_{}\end{array}}_{\mathrm{b}}^0\\ {\begin{array}{c}{u}_{}\end{array}}_{\mathrm{e}}^0\end{array}\right\}=\left\{\begin{array}{c}-{P_{}}_{\mathrm{b}}^0\\ P_{\mathrm{e}}\end{array}\right\}\mathrm{.}

Therefore, the nodal forces P_{\mathrm{e}} in terms of the free-field variables at boundary \mathrm{\Gamma } can be written as

(4)P_{\mathrm{e}}=M_{\mathrm{e}\mathrm{b}}{\ddot{u}}_{\mathrm{b}}^0+M_{\mathrm{e}\mathrm{e}}{\ddot{u}}_{\mathrm{e}}^0+C_{\mathrm{e}\mathrm{b}}{\dot{u}}_{\mathrm{b}}^0+C_{\mathrm{e}\mathrm{e}}{\dot{u}}_{\mathrm{e}}^0+K_{\mathrm{e}\mathrm{b}}{u}_{\mathrm{b}}^0+K_{\mathrm{e}\mathrm{e}}{u}_{\mathrm{e}}^0\mathrm{.}

This formulation does not provide any advantages compared with the traditional finite element approach, as the free-field response in the entire domain {\Omega }^{+} is not known a priori. However, this form can be simplified by rewriting the total displacement u_{\mathrm{e}} as follows:

(5)u_{\mathrm{e}}=u_{\mathrm{e}}^0+w_{\mathrm{e}},

where w_{\mathrm{e}} is the relative displacement with respect to u_{\mathrm{e}}^{\mathrm{0}}. Substituting Eq. (4) into Eq. (2) and adopting Eq. (5), the final form of the system of equations yields

(6)\begin{array}{l}\left[\begin{array}{ccc}{M}_{\mathrm{i}\mathrm{i}}^{\mathrm{\Omega }}& M_{\mathrm{i}\mathrm{b}}^{\mathrm{\Omega }}& 0\\ M_{\mathrm{b}\mathrm{i}}^{\mathrm{\Omega }}& M_{\mathrm{b}\mathrm{b}}^{\mathrm{\Omega }}+M_{\mathrm{b}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}& M_{\mathrm{b}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}\\ 0& M_{\mathrm{e}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}& M_{\mathrm{e}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}\end{array}\right]\left\{\begin{array}{c}{\ddot{u}}_{\mathrm{i}}\\ {\ddot{u}}_{\mathrm{b}}\\ {\ddot{w}}_{\mathrm{e}}\end{array}\right\}+\left[\begin{array}{ccc}{C}_{\mathrm{i}\mathrm{i}}^{\mathrm{\Omega }}& C_{\mathrm{i}\mathrm{b}}^{\mathrm{\Omega }}& 0\\ C_{\mathrm{b}\mathrm{i}}^{\mathrm{\Omega }}& C_{\mathrm{b}\mathrm{b}}^{\mathrm{\Omega }}+C_{\mathrm{b}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}& C_{\mathrm{b}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}\\ 0& C_{\mathrm{e}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}& C_{\mathrm{e}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}\end{array}\right]\left\{\begin{array}{c}{\dot{u}}_{\mathrm{i}}\\ {\dot{u}}_{\mathrm{b}}\\ {\dot{w}}_{\mathrm{e}}\end{array}\right\}+\\ \mathrm{\hspace{1em}\hspace{1em}}\left[\begin{array}{ccc}{K}_{\mathrm{i}\mathrm{i}}^{\mathrm{\Omega }}& K_{\mathrm{i}\mathrm{b}}^{\mathrm{\Omega }}& 0\\ K_{\mathrm{b}\mathrm{i}}^{\mathrm{\Omega }}& K_{\mathrm{b}\mathrm{b}}^{\mathrm{\Omega }}+K_{\mathrm{bb}}^{{\mathrm{\Omega }}^{+}}& K_{\mathrm{b}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}\\ 0& K_{\mathrm{e}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}& K_{\mathrm{e}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}\end{array}\right]\left\{\begin{array}{c}{u}_{\mathrm{i}}\\ u_{\mathrm{b}}\\ w_{\mathrm{e}}\end{array}\right\}=\left\{\begin{array}{c}0\\ -M_{\mathrm{b}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}{\ddot{u}}_e^0-C_{\mathrm{b}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}{\dot{u}}_{\mathrm{e}}^0-K_{\mathrm{b}\mathrm{e}}^{{\mathrm{\Omega }}^{+}}{u}_{\mathrm{e}}^0\\ M_{\mathrm{e}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}{\ddot{u}}_{\mathrm{b}}^0+C_{\mathrm{e}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}{\dot{u}}_{\mathrm{b}}^0+K_{\mathrm{e}\mathrm{b}}^{{\mathrm{\Omega }}^{+}}{u}_{\mathrm{b}}^0\end{array}\right\} ,\end{array}

where superscripts {\mathrm{\Omega }}^{+} and \mathrm{\Omega } denote the external and internal domains, respectively. The resultant system of equations is favorable from the computational perspective as the mass, damping, and stiffness submatrices on the right-hand side of the equation are zeros for all the external domains {\mathrm{\Omega }}^{+}, except for a layer of finite elements at the boundary of the internal and external domains.

A reservoir is considered as a superstructure. Hence, it is only modeled in the second step of the simulation. It is more natural to model the local site such that the entire reservoir is included in it. In addition, the larger the local site, the more realistic is the model. However, the local site should be minimized to reduce the computational cost. Therefore, the near part of the reservoir was included in the simulation, whereas the remainder was truncated at a reasonable distance and replaced by an absorbing boundary to prevent the reflection of outgoing waves.

Algorithm

To implement the formulation, a simple auxiliary problem, similar to the problem presented in Fig. 5, was solved in the first step, and the effective forces were computed at the nodes of a layer of finite elements bounded by {\mathrm{\Gamma }}^{\mathrm{e}} and \mathrm{\Gamma }. Next, the computed effective equivalent forces were utilized in the reduced model (Fig. 6) in the second step, and a problem including the detailed model of the structure and local site was analyzed.

The nodal values of the field variables for the auxiliary problem of the first step can be obtained using deconvolution analysis based on the recorded free-field excitation. Imposing the extracted equivalent seismic forces on the DRM layer and employing a simple viscous boundary for absorbing outgoing waves, the second step can be analyzed to obtain the final result [26]. A flowchart of this procedure is illustrated in Fig. 7.

Verification

ABAQUS commercial software was employed for modeling and performing finite element simulation, whereas an in-house code was developed for other tasks, including deconvolution analysis and providing a correlation between different steps. To verify the model, different cases (defined in Table 1) were introduced and analyzed using both the in-house code and ABAQUS software.

Case I

In the first case, the DRM and its implementation in the finite element program were verified using a wave propagation formulation. In this regard, numerical results from the finite element model were compared with the exact analytical solution of one-dimensional shear wave propagation in an elastic half-space. The material properties of the elastic half-space are listed in Table 2.

A linear two-dimensional domain measuring 240 m × 70 m was discretized using 5 m × 5 m quadrilateral elements to represent the elastic half-space. The incident vertical shear wave was assumed to be a wavelet, and its time history is shown in Fig. 8.

Figure 9 shows a comparison of the numerical response of the free surface with the analytical solution. According to the theory of one-dimensional vertical wave propagation in an elastic half-space, the response obtained at the surface of the half-space is amplified by a factor of two, and a time lag with respect to the incident wave exists. As shown in the figure, the numerical result agreed well with the analytical solution.

Case II

In case II, the numerical results of a vertically propagating shear wave in a viscoelastic half-space were compared with the results obtained using deconvolution. In this regard, the free-field record was assumed to be known at the surface of the viscoelastic half-space, and the response of nodal points was obtained via deconvolution analysis based on one-dimensional vertically propagating shear wave theory, as typically performed in previous studies. The time history of the Koyna earthquake, which is presented in Fig. 10, was employed for this analysis.

A linear two-dimensional viscoelastic half-space domain measuring 240 m × 70 m was discretized using 5 m × 5 m finite elements, and the material properties are listed in Table 3. The time histories of acceleration, velocity, and displacement of the numerical model, which were recorded at the surface of the half-space, were compared with the free-field input of the Koyna earthquake shown in Figs. 11 to 13, respectively. The results indicated good agreement, thereby confirming the robustness of the model in reproducing accurate wave propagation effects in the half-space domain.

Anisotropic Damage–Plastic Model of Concrete

For nonlinear simulations, the anisotropic damage–plastic model of Daneshyar and Ghaemian [61] was employed. The capability of the formulation to reproduce concrete responses and crack trajectories in multidimensional tests has been demonstrated in detail in the original paper. Additionally, the model was successfully employed for the seismic simulation of concrete arch dams in the latter works of the authors (see Ref. [62]). However, for completeness, a brief description of the formulation and some verification examples are provided herein. For detailed information, the reader is referred to the original work of Daneshyar and Ghaemian [61].

Formulation

Based on the concept of effective quantities presented by Kachanov [63], the nominal stress tensor is mapped into the effective space using the following relation:

(7){\overline{\sigma }}_{ij}=M_{ijkl}{\sigma }_{kl},

where {\bar{\sigma }}_{ij} is the effective stress tensor, and M_{ijkl} is the fourth-order mapping tensor such that

(8)M_{ijkl}=M_{ijrs}^{+}{P}_{rskl}^{+}+M_{ijrs}^{-}{P}_{rskl}^{-},

where M_{ijkl}^{+} and M_{ijkl}^{-} are the tensile and compressive fourth-order mapping tensors, respectively; P_{ijkl}^{+} and P_{ijkl}^{-} are projection tensors that project tensile and compressive fourth-order tensors into the principal space, respectively.

The positive part of the projection tensor is expressed as

(9)P_{ijkl}^{+}={\displaystyle \sum _{r=\mathrm{1}}^{\mathrm{3}}H({\sigma }_r)n_i^{(r)}{n}_j^{(r)}{n}_k^{(r)}{n}_l^{(r)}},

where {\sigma }_r is the rth eigenvalue of the stress tensor, n_i^{(r)}is the ith component of the rth eigenvector of the stress tensor, and H(x) is the Heaviside step function, defined as

(10)H(x)=\{\begin{array}{ll}\mathrm{1},& x>\mathrm{0},\\ \mathrm{0},& x\le \mathrm{0}\mathrm{.}\end{array}

Using the tensile projection tensor P_{ijkl}^{+}, the compressive part can be calculated as follows:

(11)P_{ijkl}^{-}=I_{ijkl}-P_{ijkl}^{+},

where I_{ijkl} is the forth-order identity tensor, defined as

(12)I_{ijkl}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\delta }_{ik}{\delta }_{jl}+{\delta }_{jk}{\delta }_{il}),

and {\delta }_{ij} is the Kronecker delta.

The tensile and compressive fourth-order mapping tensors M_{ijkl}^{+} and M_{ijkl}^{-} are defined as

(13)M_{ijkl}^{\pm }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}({\psi }_{ik}^{\pm }{\delta }_{jl}+{\psi }_{jk}^{\pm }{\delta }_{il}),

where {\psi }_{ij}^{+} and {\psi }_{ij}^{-} are the tensile and compressive second-order inverse integrity tensors, respectively, which are obtained by integrating its scalar values over the surface of a unit sphere as follows:

(14){\psi }_{ik}^{\pm }={\displaystyle \frac{\mathrm{3}}{\mathrm{2}\mathrm{\pi }}}{\displaystyle {\int }_{\mathrm{\Omega }}{\psi }^{\pm }{n}_i{n}_j\mathrm{d}\mathrm{\Omega }}\mathrm{.}

The concrete response is rate dependent, which means that its tensile and compressive strengths as well as fracture energy are affected by changes in the strain rate. This behavior was experimentally investigated by Brara and Klepaczko [64] and Bischoff and Perry [65] for different strain rates in tension and compression, respectively. According to their experiments, during seismic loading, in which the strain rate reached 10⁻² 1/s [65], a maximum increase by 1.5 times in concrete strength was observed. This significantly affects the overall seismic response of concrete dams and hence must be considered in comprehensive and realistic models. Therefore, the anisotropic damage–plastic model can be enhanced by the Duvaut–Lions viscous model [66]. Based on the model, the rate of effective visco–plastic strain is defined as

(15){\dot{\overline{\epsilon }}}_{ij}^{\mathrm{vp}}={\displaystyle \frac{\mathrm{1}}{\eta }}{\overline{D}}_{ijkl}^{-\mathrm{1}}({\overline{\sigma }}_{kl}-{\tilde{\sigma }}_{kl}),

where \eta is the viscosity parameter, {\tilde{\sigma }}_{ij}the stress tensor corresponding to the rate-independent behavior of the material, {\bar{\sigma }}_{ij}the current stress tensor, and {\bar{D}}_{ijkl} the stiffness matrix of the intact material, defined as

(16){\overline{D}}_{ijkl}=\overline{\lambda }{\delta }_{ij}{\delta }_{kl}+\overline{\mu }({\delta }_{ik}{\delta }_{jl}+{\delta }_{il}{\delta }_{jk}),

with \bar{\lambda } and \bar{\mu } the Lame constants.

As damage and plasticity are formulated in a coupled and rate-dependent framework, different combinations of input parameters can reproduce similar post-peak responses with similar fracture energies. This means that the selection of parameters is not unique, but it would be correct if the parameters are selected in acceptable ranges. Hence, it is more reasonable to only identify the tensile strength and fracture energy of specimens, instead of reporting a single combination of input parameters. In addition, owing to the coupling of damage and plasticity, a closed-form formulation for relating input parameters to fracture energy cannot be derived; however, the area under the uniaxial tensile stress–strain curve can be interpreted as the fracture energy of the specimens.

Verification

Two-dimensional verification

The experiments of Gálvez et al. [67] on notched concreted beams with different boundary conditions were employed for two-dimensional verification. Two cases were considered, in which specimens with similar sizes and similar concrete mixtures exhibiting the material properties reported in Table 4 were prepared, although different notch depths and different boundary conditions were considered.

First, a notched beam with dimensions and boundary conditions as presented in Fig. 14 was loaded vertically on its upper edge. The notch depth was 37.5 mm in this case. The applied force vs. the crack mouth opening displacement (CMOD) and crack trajectory resulting from the numerical simulation was compared with the experimental results shown in Figs. 15 and 16, respectively. It was evident from the figures that the numerical results complied with the experimental envelopes. The contour of the tensile damage on the deformed configuration of the beam is shown in Fig. 17.

A similar concrete mixture was used in the second test. Therefore, the same properties of concrete were adopted in the second specimen. However, the notch was extended to 45 mm and, as shown in Fig. 18, the vertical displacement of a point on the upper edge of the beam was restricted. A comparison of the numerical force–CMOD curve with the experimental envelope is shown in Fig. 19. The resulting crack profile is shown on the experimental envelope in Fig. 20, excellent agreements were obtained. Finally, the tensile damage distribution and deformed shape of the beam are presented in Fig. 21.

Three-dimensional verification

The pull-out test of Rabczuk et al. [68] was performed to assess the capability of the anisotropic damage–plastic constitutive model in a three-dimensional configuration. The test comprised an anchor embedded in a concrete block. As shown in Fig. 22, a quarter of the specimen was modeled owing to symmetry. The material properties of the specimen are reported in Table 5. Vertical displacement was imposed to mimic anchor movement, and the applied force vs. displacement of the specimen was extracted and compared with the reference result shown in Fig. 23. The response of the specimen based on the anisotropic damage–plastic model agreed well with the reference. The contour of the damaged region and the deformed configuration of the specimen are plotted in Figs. 24 and 25, respectively. The crack trajectory was observable in these figures.

Seismic analysis of concrete gravity dam

The assessment of the model’s capability in reproducing the seismic behavior of a dam-massed foundation system is presented in this section. In this regard, the Koyna concrete gravity dam was analyzed, and the effects of wave propagation in the semi-infinite domain and the radiation damping of the massed foundation on its linear and nonlinear seismic responses were investigated. It is noteworthy that both the DRM and anisotropic damage–plastic model were formulated in general three-dimensional conditions. However, a gravity concrete dam with a two-dimensional configuration was considered in the simulation, which did not violate the generality of the model.

The body of the dam was idealized as an assemblage of planar, four-noded rectangular elements. The reservoir was assumed to be compressible with a depth of 91.7 m, and the foundation was modeled as a semi-infinite, homogeneous, isotropic, viscoelastic half-plane measuring 850 m × 260 m. The finite element model of the system is shown in Fig. 26; it comprised 1480, 5304, and 9152 quadrilateral elements for the body of dam, reservoir, and foundation domain, respectively. Furthermore, the location of the truncated regions in this figure is a schematic representation, and the absorbing boundaries were 390 m from the upstream and downstream faces, and 260 m below the base of the dam. The material properties of both domains and the nonlinear properties of concrete are listed in Tables 6 and 7, respectively. The viscosity parameter was selected such that the tensile strength of concrete was magnified by a factor of 1.2 at a strain rate of 10^–2 (for further details, see Ref. [62]). Additionally, the Koyna earthquake record (see Fig. 10) was selected for these simulations.

Different conditions were considered for the simulations, as summarized in Table 8. For the first four models, the linear behavior of the body was assumed, and the seismic responses of the dam were investigated with massless and massed foundations for empty and full reservoir conditions, whereas similar models were analyzed in the last four simulations, except for the nonlinear response of concrete.

Linear

Empty reservoir

A comparison between the crest relative displacement of models with massless and massed foundation assumptions is shown in Fig. 27. The contour of the maximum principle tensile stress for both simulations is shown in Fig. 28. The maximum values are summarized in Table 9. A maximum relative displacement of 48.7 mm toward the upstream direction was observed in the modeling with the massless foundation assumption, whereas this value decreased to 40.6 mm in the massed foundation model. In both cases, the maximum tensile stress occurred at the location of slope change on the downstream face. The distributions of the maximum principal tensile stress were similar. However, different magnitudes of tensile stress were observed in the two models. As shown in Table 9, the maximum principal tensile stress was 6.40 MPa in the massless foundation model; however, the value decreased to 4.66 MPa when the massed foundation was assumed, i.e., a 27% decrease.

Full reservoir

The full reservoir cases were compared, and the relative crest displacements and contours of the maximum principal tensile stress are shown in Figs. 29 and 30, respectively. The maximum values are reported in Table 10. In contrast to the empty reservoir model, the maximum relative crest displacement of both models occurred in the downstream direction, with a value of 75.6 mm for the massless foundation assumption and 54.4 mm for the model with a massed foundation. The maximum tensile stress was observed at the base of the dam in both models, whereas the value was 12.2 MPa for the massless model and 7.57 MPa for the massed model. Similarly, the stress levels decreased because the foundation mass was considered in the simulation.

Nonlinear

Empty reservoir

The relative crest displacements of the simulations with nonlinear concrete response and empty reservoir assumption for both the massless and massed foundation models were compared, as shown in Fig. 31. The contours of tensile damage for both assumptions at the fourth and tenth second of the analysis were compared as well, as shown in Figs. 32 and 33, respectively. The maximum observed values of tensile damage and relative crest displacement are reported in Table 11. Both simulations show permanent relative displacement of the crest in the upstream direction, which can be attributed to nonlinear deformation resulting from damage and cracking in the body of the dam. Damage growth began at the location of slope change on the downstream face. It propagated within the body and passed through the thickness of the dam in the massless foundation assumption, whereas it halted at the early stages in the massed foundation model. Hence, a more permanent relative displacement was observed in the massless model. The base of the dam remained undamaged in these simulations.

Full reservoir

A comparison between the crest displacement of models with anisotropic damage–plastic behavior of concrete and full reservoir assumption for both the massless and massed foundations are presented in Fig. 34. The tensile damage contours of both models at the fourth and tenth second of the simulation are presented in Figs. 35 and 36, respectively. A summary of the maximum values is shown in Table 12. In the full reservoir model, tensile damage began at the base and downstream face of the body. Stress concentration at both locations triggered damage growth. Propagation continued in the modeling with a massless foundation assumption, whereas it halted at the early stages in the massed foundation model. Therefore, less permanent relative displacement toward the downstream was observed in the modeling with the massed foundation assumption.

A comparison between the duration of analysis for different cases is presented in Table 13. All simulations were performed on an Intel Core i3@2.5 GHz computer with 4 GB of RAM. Slight differences between the massless and massed cases with the linear material response were observed, and the reservoir imposed a negligible effect. Meanwhile, the computational costs in the nonlinear cases depended on the intensity of damage. A higher intensity of damage requires more computation, and vice versa. However, significant discrepancies between the computational costs of different cases was not observed, and all simulations were performed in less than 1 h.

Conclusions

A finite element model for the nonlinear seismic simulation of the dam–foundation–reservoir system was presented. The foundation of the system was regarded as a massed semi-infinite medium, whereas wave propagation effects and radiation damping were considered using the DRM. The nonlinear response of concrete was reproduced using a microplane-based rate-dependent anisotropic damage–plastic model of concrete. Both the DRM and damage–plastic models were first verified and then implemented in a unified finite element model for the seismic simulation of concrete dams. Based on the responses of the Koyna gravity dam under different conditions of concrete behavior, reservoir assumption, and foundation model, the following were observed.

1) Material damping and the far-end of the reservoir were the only energy absorbing mechanisms in the massless foundation assumption. Hence, overestimated results were observed. However, by introducing a foundation mass within the model, wave propagation effects and radiation damping were included, resulting in reduced stress levels and relative displacements.

2) As wave propagation effects were neglected in the massless foundation, even by considering sedimentary materials at the bottom, the far-end of the reservoir was the main absorbing source of hydrodynamic waves, whereas some of them can be transferred to the outside of the domain by propagating within the massed medium. Hence, reduced stress levels was more evident in the full reservoir assumption.

3) The linear material model only provided estimations regarding the location of crack initiation, whereas the nonlinear model delivered detailed crack trajectories and damage profiles.

4) It was observed in the full reservoir models that cracking initiated at the heel of the dam. However, as only the nonlinear response of the body was considered, the crack propagated in a straight line. For a more detailed profile, a more sophisticated model of the dam–foundation interface in addition to a nonlinear description of the foundation response are required.

References

[1] Ghaemian M, Ghobarah A. Staggered solution schemes for dam–reservoir interaction. Journal of Fluids and Structures, 1998, 12(7): 933–948

[2] Ghaemian M, Ghobarah A. Nonlinear seismic response of concrete gravity dams with dam–reservoir interaction. Engineering Structures, 1999, 21(4): 306–315

[3] Chopra K, Chakrabarti P. Earthquake analysis of concrete gravity dams including dam-water-foundation rock interaction. Earthquake Engineering & Structural Dynamics, 1981, 9(4): 363–383

[4] Fenves G, Chopra A K. Earthquake analysis of concrete gravity dams including reservoir bottom absorption and dam-waterfoundation rock interaction. Earthquake Engineering & Structural Dynamics, 1984, 12(5): 663–680

[5] Fenves G, Chopra A K. Simplified earthquake analysis of concrete gravity dams. Journal of Structural Engineering, 1987, 113(8): 1688–1708

[6] Bouaanani N, Lu F Y. Assessment of potential-based fluid finite elements for seismic analysis of dam–reservoir systems. Computers & Structures, 2009, 87(3–4): 206–224

[7] Samii V L, Lotfi V. Application of H–W boundary condition in dam–reservoir interaction problem. Finite Elements in Analysis and Design, 2012, 50: 86–97

[8] Clough R. Non-linear mechanisms in the seismic response of arch dams. In: Proceedings of the International Conference on Earthquake Engineering. Skopje, 1980, 1: 669–684

[9] Tan H, Chopra A K. Earthquake analysis of arch dams including dam-water-foundation rock interaction. Earthquake Engineering & Structural Dynamics, 1995, 24(11): 1453–1474

[10] Tan H, Chopra A K. Dam-foundation rock interaction effects in frequency-response functions of arch dams. Earthquake Engineering & Structural Dynamics, 1995, 24(11): 1475–1489

[11] Chopra K. Earthquake analysis of arch dams: Factors to be considered. Journal of Structural Engineering, 2012, 138(2): 205–214

[12] Lysmer J, Kuhlemeyer R L. Finite dynamic model for infinite media. Journal of the Engineering Mechanics Division, 1969, 95(4): 859–878

[13] Givoli D. High-order local non-reflecting boundary conditions: A review. Wave Motion, 2004, 39(4): 319–326

[14] Katz D, Thiele E, Taffove A. A perfectly matched layer for the absorbing of electromagnetic waves. IEEE Microwave and Guided Wave Letters, 1994, 4: 268270

[15] Basu U, Chopra A K. Perfectly matched layers for time-harmonic elastodynamics of unbounded domains: Theory and finite-element implementation. Computer Methods in Applied Mechanics and Engineering, 2003, 192(11–12): 1337–1375

[16] Basu U, Chopra A K. Perfectly matched layers for transient elastodynamics of unbounded domains. International Journal for Numerical Methods in Engineering, 2004, 59(8): 1039–1074

[17] Ungless R F. Infinite finite element. Dissertation for the Doctoral Degree. Vancouver: University of British Columbia, 1973

[18] Bettess P. Infinite elements. International Journal for Numerical Methods in Engineering, 1977, 11(1): 53–64

[19] Zienkiewicz O, Bicanic N, Shen F. Advances in Computational Nonlinear Mechanics. Vienna: Springer, 1989, 109–138

[20] Wolf J, Hall W. Soil-structure-interaction analysis in time domain. Nuclear Engineering and Design, 1989, 111(3):381–393

[21] Léger P, Boughoufalah M. Earthquake input mechanisms for time-domain analysis of dam-foundation systems. Engineering Structures, 1989, 11(1): 37–46

[22] Wang J, Zhang C, Jin F. Nonlinear earthquake analysis of high arch dam–water–foundation rock systems. Earthquake Engineering and Structural Dynamics, 2012, 41(7):1157–1176

[23] Wang J, Chopra A K. EACD-3D-2008: A Computer Program for Three-Dimensional Earthquake Analysis of Concrete Dams Considering Spatially-Varying Ground Motion. Berkeley, CA: Earthquake Engineering Research Center, University of California, 2008

[24] Bielak J, Christiano P. On the effective seismic input for non-linear soil-structure interaction systems. Earthquake Engineering & Structural Dynamics, 1984, 12(1): 107–119

[25] Bielak J, Loukakis K, Hisada Y, Yoshimura C. Domain reduction method for three-dimensional earthquake modeling in localized regions, part I: Theory. Bulletin of the Seismological Society of America, 2003, 93(2): 817–824

[26] Yoshimura C, Bielak J, Hisada Y, Fernandez A. Domain reduction method for three-dimensional earthquake modeling in localized regions, part II: Verification and applications. Bulletin of the Seismological Society of America, 2003, 93(2): 825–841

[27] Ayari M L, Saouma V E. A fracture mechanics based seismic analysis of concrete gravity dams using discrete cracks. Engineering Fracture Mechanics, 1990, 35(1–3): 587–598

[28] Pekau O, Chuhan Z, Lingmin F. Seismic fracture analysis of concrete gravity dams. Earthquake Engineering & Structural Dynamics, 1991, 20(4): 335–354

[29] Pekau O, Batta V. Seismic crack propagation analysis of concrete structures using boundary elements. International Journal for Numerical Methods in Engineering, 1992, 35(8): 1547–1564

[30] El-Aidi B, Hall J F. Non-linear earthquake response of concrete gravity dams part 1: Modelling. Earthquake Engineering & Structural Dynamics, 1989, 18(6): 837–851

[31] Bhattacharjee S, Leger P. Seismic cracking and energy dissipation in concrete gravity dams. Earthquake Engineering & Structural Dynamics, 1993, 22(11): 991–1007

[32] Guanglun W, Pekau O, Chuhan Z, Shaomin W. Seismic fracture analysis of concrete gravity dams based on nonlinear fracture mechanics. Engineering Fracture Mechanics, 2000, 65(1): 67–87

[33] Hillerborg M, Modéer M, Petersson P E. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. Cement and Concrete Research, 1976, 6(6): 773–781

[34] Bazant Z P, Oh B H. Crack band theory for fracture of concrete. Mat’eriaux et construction, 1983, 16(3): 155–177

[35] Mirzabozorg H, Ghaemian M. Non-linear behavior of mass concrete in three-dimensional problems using a smeared crack approach. Earthquake Engineering & Structural Dynamics, 2005, 34(3): 247–269

[36] Areias P, Rabczuk T. Finite strain fracture of plates and shells with configurational forces and edge rotations. International Journal for Numerical Methods in Engineering, 2013, 94(12): 1099–1122

[37] Areias P, Rabczuk T, Dias-da Costa D. Element-wise fracture algorithm based on rotation of edges. Engineering Fracture Mechanics, 2013, 110: 113–137

[38] Areias P, Rabczuk T. Steiner-point free edge cutting of tetrahedral meshes with applications in fracture. Finite Elements in Analysis and Design, 2017, 132: 27–41

[39] Moës N, Belytschko T. Extended finite element method for cohesive crack growth. Engineering Fracture Mechanics, 2002, 69(7): 813–833

[40] Daneshyar S M, Mohammadi S. Strong tangential discontinuity modeling of shear bands using the extended finite element method. Computational Mechanics, 2013, 52(5): 1023–1038

[41] Afshar A, Daneshyar A, Mohammadi S. XFEM analysis of fiber bridging in mixed-mode crack propagation in composites. Composite Structures, 2015, 125: 314–327

[42] Song J H, Belytschko T. Cracking node method for dynamic fracture with finite elements. International Journal for Numerical Methods in Engineering, 2009, 77(3): 360–385

[43] Ren H, Zhuang X, Rabczuk T. Dual-horizon peridynamics: A stable solution to varying horizons. Computer Methods in Applied Mechanics and Engineering, 2017, 318: 762–782

[44] Sluys L. Wave propagation, localisation and dispersion in softening solids. Dissertation for the Doctoral Degree. Delft: Delft University of Technology, 1992

[45] Bažant Z P, Jir’asek M. Nonlocal integral formulations of plasticity and damage: Survey of progress. Journal of Engineering Mechanics, 2002, 128(11): 1119–1149

[46] Ren H, Zhuang X, Rabczuk T. A nonlocal operator method for solving partial differential equations. Computer Methods in Applied Mechanics and Engineering, 2020, 358: 112621

[47] Løland K. Continuous damage model for load-response estimation of concrete. Cement and Concrete Research, 1980, 10(3): 395–402

[48] Krajcinovic D, Fonseka G. The continuous damage theory of brittle materials. Journal of Applied Mechanics, 1981, 48(4): 809–815

[49] Lee J, Fenves G L. A plastic-damage concrete model for earthquake analysis of dams. Earthquake Engineering & Structural Dynamics, 1998, 27(9): 937–956

[50] Cervera M, Oliver J, Manzoli O. A rate-dependent isotropic damage model for the seismic analysis of concrete dams. Earthquake Engineering & Structural Dynamics, 1996, 25(9): 987–1010

[51] Bourdin G A, Francfort G A, Marigo J J. The variational approach to fracture. Journal of Elasticity, 2008, 91(1–3): 5–148

[52] Miehe F, Welschinger F, Hofacker M. Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations. International Journal for Numerical Methods in Engineering, 2010, 83(10): 1273–1311

[53] Zhou S, Rabczuk T, Zhuang X. Phase field modeling of quasi-static and dynamic crack propagation: Comsol implementation and case studies. Advances in Engineering Software, 2018, 122: 31–49

[54] Zhou S, Zhuang X, Zhu H, Rabczuk T. Phase field modelling of crack propagation, branching and coalescence in rocks. Theoretical and Applied Fracture Mechanics, 2018, 96: 174–192

[55] Zhou S, Zhuang X, Rabczuk T. Phase field modeling of brittle compressive-shear fractures in rock-like materials: A new driving force and a hybrid formulation. Computer Methods in Applied Mechanics and Engineering, 2019, 355: 729–752

[56] Ren H, Zhuang X, Anitescu C, Rabczuk T. An explicit phase field method for brittle dynamic fracture. Computers & Structures, 2019, 217: 45–56

[57] Zhou S, Zhuang X, Rabczuk T. A phase-field modeling approach of fracture propagation in poroelastic media. Engineering Geology, 2018, 240: 189–203

[58] Zhou S, Zhuang X, Rabczuk T. Phase-field modeling of fluid-driven dynamic cracking in porous media. Computer Methods in Applied Mechanics and Engineering, 2019, 350: 169–198

[59] Zhuang X, Zhou S, Sheng M, Li G. On the hydraulic fracturing in naturally-layered porous media using the phase field method. Engineering Geology, 2020, 266: 105306

[60] Goswami S, Anitescu C, Chakraborty S, Rabczuk T. Transfer learning enhanced physics informed neural network for phase-field modeling of fracture. Theoretical and Applied Fracture Mechanics, 2020, 106: 102447

[61] Daneshyar M G, Ghaemian M. Coupling microplane-based damage and continuum plasticity models for analysis of damage-induced anisotropy in plain concrete. International Journal of Plasticity, 2017, 95: 216–250

[62] Daneshyar M G, Ghaemian M. Seismic analysis of arch dams using anisotropic damage-plastic model for concrete with coupled adhesivefrictional joints response. Soil Dynamics and Earthquake Engineering, 2019, 125: 105735

[63] Kachanov L. On creep rupture time. Bulletin of the Academy of Sciences of the USSR, 1958, 8: 26–31

[64] Brara J K, Klepaczko J R. Experimental characterization of concrete in dynamic tension. Mechanics of Materials, 2006, 38(3): 253–267

[65] Bischoff P, Perry S. Compressive behaviour of concrete at high strain rates. Materials and Structures, 1991, 24(6): 425–450

[66] Duvaut G, Lions J. Inequalities in Mechanics and Physics. New York: Springer-Verlag, 1977

[67] Gálvez J, Elices M, Guinea G, Planas J. Mixed mode fracture of concrete under proportional and nonproportional loading. International Journal of Fracture, 1998, 94(3): 267–284

[68] Rabczuk T, Bordas S, Zi G. A three-dimensional meshfree method for continuous multiple-crack initiation, propagation and junction in statics and dynamics. Computational Mechanics, 2007, 40(3): 473–495

RIGHTS & PERMISSIONS

Higher Education Press