Domain information transfer method and its application in quasi-brittle failure analysis

Abstract

This article presents a domain information transfer method for the purpose of quasi-brittle failure analysis. It contains two separated computational levels. On macroscopic level, the material domains are viewed as elementary substances which construct the whole structure. The macroscopic governing equations are established based on this point of view and the boundary displacements of material domains are defined as basic unknowns. Each damaged material domain will be assigned a mesoscale model and the information transfer mechanism from a mesoscale domain to a macroscale domain is established. The damage localization phenomenon is equivalently represented on macroscopic level by boundary deformation of material domains and the macroscale meshes are invariant during analysis. The information exchange element based on nonlocal quadrature scheme is proposed and the efficient upscaling is realized. The proposed method is an efficient and flexible way for quasi-brittle failure analysis.

Introduction

The analysis of quasi-brittle failure is an important topic. The quasi-brittle materials, such as concrete, rocks, and other types of cement-based materials possess complex micro-structure. The damage and fracture of these materials or the structures constructed by them will be affected by this feature significantly. The damage is usually distributed in the materials under relatively low load. With damage evolution, the chaotic, stochastic, and distributed damage pattern will convert to the one with obvious directionality and localized characteristics. The prior-formed damage localization band will depress micro-cracks beside it and annex micro-cracks on its forward path, which leads to the damage localization phenomenon on macroscale and influences the overall performance of structure.

The nature of damage process makes it necessary to consider mesoscopic features of materials.^1–4 Many mesoscopic models have been developed for this purpose, such as lattice model and multi-phase model. The lattice model is a classical model suitable for simulation of fracture phenomenon.^5,6 It has also been applied to the analysis of deterioration of concrete properties caused by alkali-aggregate reaction recently. ⁷ The multi-phase model is another classical numerical model to simulate quasi-brittle fracture, especially for tensile failure behavior. ⁸ In addition, some kind of mesoscopic models introduce interface elements to simulate internal interfaces in materials, such as the model proposed by Caballero and colleagues^9–12 and the model based on mesh fragmentation technique. ¹³ The interface damage law is applied in these models conveniently to predict local fracture. However, micro-cracks can only propagate along the existing internal interfaces, which restrict their distribution in materials. It is clear that mesoscopic modeling provides a perfect way to the analysis of quasi-brittle failure. However, the computational cost is considerable if these models are applied to real engineering specimens or structures. Therefore, the ideal of multiscale computation is necessary.

There are many types of multiscale methods that have been proposed for analysis of heterogeneous materials, such as multiscale finite element method (MsFEM)^14–18 and strong coupling method. ¹⁹ However, concurrent and hierarchical methods are still most popular.

For the concurrent method, the regions with different scale characteristics are combined to be the whole model to balance the accuracy and efficiency.^20–22 The important topic in concurrent method is the connection of different scale regions. ²³ A simple way is to establish the mesh transition zone to realize the variation of mesh density from mesoscopic region to macroscopic region. Another direct method is the enforcement of displacement coordination on the interface between macroscopic and mesoscopic regions. Also, the Lagrange multiplier method is usually applied to the coupling of different scale regions.^24,25 The aim of concurrent method is inserting an elaborate region in the model to ensure the necessary analysis accuracy and to keep the rest of the parts of the model as simple as possible. So, the region with stress concentration or defects, or the region that just needed to be treated carefully should be set to the elaborate region. This strategy guarantees the utmost accuracy because the regions with precise analysis demand can all be restricted to the elaborate area. With damage evolution, some parts of macroscopic region need to be replaced with mesoscopic region, such as the concurrent model for the damage trans-scale evolution analysis in concrete^26,27 and the adaptive concurrent multiscale model with coupling elements. ²⁸ Undoubtedly, the adaptive mechanism is useful in concurrent method. However, it brings some difficulties under concurrent framework. First, the interface between macroscopic and mesoscopic regions continues to change due to the damage evolution, which enlarges the range of mesoscopic region during analysis. It leads to the difficulty to adjust the whole model during analysis and results in the loss of load history, which is a main reason for the demand of repeated computation. ²⁶ Second, if the range of elaborate region keeps enlarging, the degrees of freedom of the whole model will reach a considerable amount which causes the computation to be unsustainable. Therefore, the advantage of concurrent method will vanish on that situation.

The hierarchical framework is a way that injects the information from fine-scale model to coarse-scale model. The most representative method is computational homogenization.^29,30 In this method, the phenomenological macroscopic constitutive law will not be needed. It is extracted from the fine scale, on which a representative volume element (RVE) of the micro-structure is identified.^31–34 For the materials with complex micro-structure, it is an effective way to reflect the influence of microfeatures on the macroscopic mechanical properties. ³⁵ The Hill–Mandel condition is a very important foundation of the method, which ensures the energy consistency between macroscopic and mesoscopic levels.^34,36–39 Many types of mesoscopic models can be applied to construct the RVE, such as Lattice, FEM, and discrete element method (DEM).^40–43 Recently, some adjustments have been applied to the method to make it suitable for damage localization problems,^44–46 including the multiscale aggregating discontinuities method ⁴⁷ and coupled-volume multiscale model.^48,49

The apparent advantage of computational homogenization method is that it transfers the information from mesoscale to macroscale. There is no need to combine different scale regions into a whole model. It brings significant flexibility and the consideration of adjusting the ranges of mesoscopic and macroscopic regions is unnecessary. So, the repeated computation is avoided. This is an efficient method and how to equivalently simulate the damage localization band more accurately and naturally is the hinge to improve the method.

The purpose of this work is to establish a multi-level computational method, the domain information transfer method (DITM), for the analysis of quasi-brittle failure problems. The method contains two different computational levels. On macroscopic level, the material domains are viewed as elementary substances which construct the whole structure. Each damaged material domain will be assigned a mesoscale model with the same size. The mechanical properties of these material domains are obtained by the upscaling method, which transfers information from a mesoscale domain to a macroscale domain. It is an efficient and flexible way which is suitable for materials with complex micro-structure and banded damage pattern.

This article is organized as follows. In section “Domain information transfer method,” the fundamentals of DITM are presented. The discrete governing equations of macroscopic problem are detailed in section “Discretization on macroscopic level.” The downscaling and mesoscopic computational level are detailed in section “Downscaling and mesoscopic computational level.” The information exchange element (IEE) based on nonlocal quadrature scheme will be proposed in section “Quadrature scheme and upscaling,” and the upscaling is detailed. The procedure of the whole computational framework is detailed in section “Damage criterion and two-level computational framework.” The numerical results are discussed in section “Numerical results and discussions,” followed by conclusions in section “Conclusion.”

Domain information transfer method

Here, the material domain, an elementary substance with the ability of independently representing both the deformation and local failure is defined. It will contribute to the analysis from distributed damage to damage localization at macroscopic level.

The data exchange exists between these material domains and corresponding mesoscale models. On macroscopic level, the structure will be considered to be constructed by these material domains. Different from a material point, the intrinsic mechanical properties of material domains are reflected by the relationship between boundary displacement and boundary traction. (Here, the body forces are not under consideration.)

The classical weak form of governing equations for coarse-scale problem or macroscopic problem (without consideration of body force) under quasi-static load are given by

$${\int }_{\Omega }\mathbf{\sigma }:\delta \mathbf{\epsilon }d\Omega ={\int }_{s_{\sigma }}\overline{\mathit{p}}·\delta \mathit{u}\mathrm{ds}$$

(1)

where $\Omega$ denotes the coarse-scale domain or macroscopic domain of the whole problem and $S_{\sigma }$ denotes the external boundary with prescribed traction. $\mathbf{\sigma }$ denotes Cauchy stress and $\overline{\mathit{p}}$ denotes the traction applied to the structure. $\delta \mathbf{\epsilon }$ and $\delta \mathit{u}$ denote virtual strain and virtual displacement, respectively. The essential boundary condition should also be satisfied on $S_u$ as plotted in Figure 1. Based on the view that the structure consists of a number of material domains, the left hand side of equation (1) can be expressed as the sum of items corresponding to each material domain, which gives

$$\sum _{l=1}^{\mathrm{Ne}}{\int }_{{\Omega }_l}\mathbf{\sigma }:\delta \mathbf{\epsilon }d\Omega ={\int }_{s_{\sigma }}\overline{\mathit{p}}·\delta \mathit{u}\mathrm{ds}$$

(2)

where Ne is the amount of material domains in structure and ${\Omega }_l$ denotes the lth material domain. Applying the Gauss divergence theorem and considering the mechanical properties of material domains, which gives

$$\sum _{l=1}^{\mathrm{Ne}}\underset{S_l}{\oint }{\mathit{t}}^{\left[l\right]}·\delta {\mathit{u}}^{\left[l\right]}\mathrm{ds}=\underset{S_{\sigma }}{\int }\overline{\mathit{p}}·\delta \mathit{u}\mathrm{ds}$$

(3)

$${\mathit{u}}^{\left[l\right]}(s)\stackrel{\mathbf{T}}{\to }{\mathit{t}}^{\left[l\right]}(s)$$

(4)

where S_l is the boundary of material domain ${\Omega }_l$ and ${\mathit{t}}^{\left[l\right]}$ , $\delta {\mathit{u}}^{\left[l\right]}$ , and ${\mathit{u}}^{\left[l\right]}$ are traction, virtual displacement, and displacement on the boundary of material domain, ${\Omega }_l$ , respectively. It should be noted that the displacement on common boundary of adjacent material domains is equal. The incremental form of equation (4) should be applied if nonlinear mechanical properties are considered. Here, the material domain is considered to be an elementary substance with independent macroscopic mechanical properties while not just a region. One can note that the concepts of stress and strain vanish on macroscopic level. T is an operator reflecting the intrinsic properties of a material domain, which depends on components and micro-structure of mesoscale. According to equations (3) and (4), a new form of macroscopic governing equations can be expressed as

$$\sum _{l=1}^{\mathrm{Ne}}\underset{S_l}{\oint }\mathbf{T}\left({\mathit{u}}^{\left[l\right]}\left(s\right)\right)·\delta {\mathit{u}}^{\left[l\right]}\mathrm{ds}=\underset{S_{\sigma }}{\int }\overline{\mathit{p}}·\delta \mathit{u}\mathrm{ds}$$

(5)

Here, the basic unknowns convert from the displacement defined at every material point to the displacement defined on boundaries of material domains.

Figure 1.

Domain information transfer method.

Equation (5) is an important fundament of the proposed DITM. It is easy to prove that the combination of equation (5) and essential boundary condition on $S_u$ are the equivalent description of the macroscopic problem. The advantage of this kind of governing equations is the avoidance of representing relatively complex internal status when damage localization phenomenon happens. It uses the boundary deformation to represent damage localization phenomenon which makes the representation more flexible and convenient. The difficulty to solve equation (5) directly is that the explicit expression of operator T is a hard work due to the complex mesoscopic characteristics of quasi-brittle material. However, it provides convenience for the application of information transfer method. A domain information transfer mechanism in DITM is established based on the ideal of upscaling and downscaling to solve the macroscopic problem described by equation (5).

A mesoscopic model with the same size of material domain will be assigned to it to obtain its properties, which replaces the function of operator T. This is a process of upscaling and downscaling, similar to computational homogenization. The difference is that this is an information transfer process from domain to domain. For the macroscopic problem, the boundaries of material domains will be discretized into IEEs which are responsible for data exchange between material domain and mesoscopic model. The macroscopic boundary displacement of each material domain is approximated using the IEEs. By downscaling, this macroscopic boundary displacement will be transferred into mesoscopic level and the solution of a mesoscopic problem will be obtained at that level. Then, the domain forces and domain tangent stiffness matrix of a material domain will be achieved by upscaling. The shape of a material domain could be arbitrary in theory and rectangle is adopted here for the convenience of implementation.

The novelty of the proposed method includes the following aspects:

1. On macroscopic level, material domains are viewed as the elementary substances to construct structure. It is different from the traditional substructure concept because this is the elementary substance with ability to independently reflect both deformation and damage localization phenomenon by itself. The upscaling and downscaling processes are applied to obtain the properties of material domains and corresponding mesoscopic solution is obtained simultaneously. It should be noted that the size of a material domain is not arbitrary. The size should not be too small to reflect macroscopic properties of material. It also should not be too big for the reason that corresponding mesoscopic model needs to be small enough to save computational cost. The governing equations of macroscopic problem are established based on this point of view and the corresponding multi-level domain information transfer mechanism is established to solve the problem.

2. An IEE is proposed based on nonlocal quadrature scheme which realizes the efficient information transfer between material domain and mesoscopic model.

3. A kind of damage criterion is established to determine which material domain enters the damage stage. It relates to where the domain information transfer mechanism is activated. So, the computational efficiency is improved significantly.

4. The boundary displacement of material domain is approximated using IEEs and the damage localization phenomenon is equivalently simulated on macroscopic level, including the phenomenon that the damage localization band gets across boundaries of adjacent material domains.

Discretization on macroscopic level

As basic unknowns on macroscopic level, the boundary displacements of material domains need to be solved. The boundary of each material domain is discretized into several elements, as in macroscopic numerical model depicted in Figure 2. The boundary displacement is interpolated in these elements. Adjacent material domains share the same elements between them. In addition to approximate boundary displacement, the element takes responsibility for information exchange between material domain and mesoscopic model. Thus, it will be called as IEE.

Figure 2.

Macroscopic and mesoscopic computational levels.

The size of IEE is a model parameter, which decides the number of degrees of freedom on the boundary of a material domain. Severe deformation mode is allowed somewhere on the material domain’s boundary by reasonable setting of the size of IEE, which makes the equivalent damage localization band to get across the boundary of adjacent material domains. Compared with the concurrent method, the meshes and degrees of freedom of macroscopic model remain unchanged during the analysis, and the load history can be retained, which significantly improves the flexibility and efficiency of the computational framework. The purpose of this article is to establish a computational framework. It shows some advantages in the description of damage localization phenomenon on macroscopic level, but some mesh dependency appears on the boundary of material domain, as discussed in section “Numerical results and discussions.” It will be the future work to deal with the issue.

Benefit from the same basic displacement mode of each IEE, the damage development near material domain’s boundary will not be differentially induced. Without losing generality, linear interpolation is adopted here. The boundary displacement in IEE k can be expressed as

$${\mathit{u}}^{\left(k\right)}(•)={\mathit{N}}^{\left(k\right)}(•){\mathit{d}}^{\left(k\right)}$$

(6)

where k is number of IEE, ${\mathit{N}}^{\left(k\right)}(•)$ denotes shape functions, ${\mathit{u}}^{\left(k\right)}(•)$ denotes the approximation of material domain’s boundary displacement in IEE k, and ${\mathit{d}}^{\left(k\right)}$ denotes nodal displacement of IEE k. $(•)$ represents the effective global coordinate of a point in IEE. When IEE is parallel to x-axis, $(•)$ is the x-coordinate of point. It is the y-coordinate of point when IEE is parallel to y-axis. In this article, the variables that correspond to IEEs will be indicated with the round brackets superscript and the variables that correspond to material domains will be indicated with the square brackets superscript. Applying equation (6), gives

$$\begin{array}{c}\sum _l{\sum _{k\in \left[l\right]}{\left\{{\overline{\mathit{d}}}^{\left(k\right)}\right\}}^T\underset{S_{l,k}}{\int }[{\mathit{N}}^{\left(k\right)}(•)]}^T{\mathit{t}}^{\left[l\right],\left(k\right)}(•)\mathrm{ds}=\sum _m{\left\{{\overline{\mathit{d}}}^{\left(m\right)}\right\}}^T\underset{S_{\sigma },m}{\int }{[{\mathit{N}}^{\left(m\right)}(•)]}^T\overline{\mathit{p}}(•)\mathrm{ds}\hfill \\ {\mathit{t}}^{\left[l\right]}(•)=\mathbf{T}\left({\mathit{u}}^{\left[l\right]}(•)\right)\hfill \end{array}$$

(7)

Here, ${\overline{\mathit{d}}}^{\left(k\right)}$ denotes virtual nodal displacement of IEE and m is the number of IEE applied with external load. ${\mathit{t}}^{\left[l\right],\left(k\right)}(•)$ is the boundary traction on IEE k of material domain l. It should be noted that ${\mathit{t}}^{\left[l\right],\left(k\right)}(•)$ depends on the boundary displacement, ${\mathit{u}}^{\left[l\right]}(•)$ , of the whole material domain. In other words, it depends on the degrees of freedom, ${\mathit{d}}^{\left[l\right]}$ , of the whole material domain and not on the degrees of freedom of IEE k alone. Let

$${\mathit{F}}^{\left(k\right)}\left({\mathit{d}}^{\left[l\right]}\right)=\underset{S_{l,k}}{\int }{\left[{\mathit{N}}^{\left(k\right)}(•)\right]}^T{\mathit{t}}^{\left[l\right],\left(k\right)}(•)\mathrm{ds}$$

(8)

where ${\mathit{F}}^{\left(k\right)}({\mathit{d}}^{\left[l\right]})$ is the nodal force of IEE k, which belongs to material domain l. The domain force ${\mathit{F}}_i^{\left[l\right]}$ can be written as

$${\mathit{F}}_i^{\left[l\right]}={{\mathit{F}}_i}^{\left(k_i^1\right)}\left({\mathit{d}}^{\left[l\right]}\right)+{{\mathit{F}}_i}^{\left(k_i^2\right)}\left({\mathit{d}}^{\left[l\right]}\right)$$

(9)

where $k_i^1$ and $k_i^2$ are the numbers of two adjacent IEEs connected by common node i and ${{\mathit{F}}_i}^{(k_i^1)}({\mathit{d}}^{\left[l\right]})$ and ${{\mathit{F}}_i}^{(k_i^2)}({\mathit{d}}^{\left[l\right]})$ are element nodal force components on common node i, respectively. Therefore, the discrete governing equations of macroscopic problem can be established as

$$\sum _l\sum _{i\in \left[l\right]}{\left\{{\overline{\mathit{d}}}_i\right\}}^T{\mathit{F}}_i^{\left[l\right]}=\sum _{\mathrm{mp}}{\left\{{\overline{\mathit{d}}}_{\mathrm{mp}}\right\}}^T{\mathit{Q}}_{\mathrm{mp}}$$

(10)

${\overline{\mathit{d}}}_i$ is virtual displacement of domain node and mp is the number of node applied with external nodal force. ${\mathit{Q}}_{\mathrm{mp}}$ and ${\overline{\mathit{d}}}_{\mathrm{mp}}$ are external nodal force and corresponding virtual nodal displacement, respectively. It can be noted that the domain forces of material domain l are determined by ${\mathit{d}}^{\left[l\right]}$ . Considering the arbitrariness of virtual nodal displacement, governing equation (10) can be expressed as

$$\mathit{F}(\mathit{d})=\mathit{Q}$$

(11)

This is the final version of governing equations of macroscopic problem. Here, $\mathit{F}$ denotes the internal nodal force of the whole structure. The nodal displacement on all material domains’ boundaries is denoted by $\mathit{d}$ and $\mathit{Q}$ is the external nodal force.

Downscaling and mesoscopic computational level

As the explicit representation of operator T is quite difficult, the macroscopic mechanical properties of material domain described by T will be equivalently obtained from mesoscopic level. This is a two-way information transfer process. The boundary displacement (boundary displacement increment) of material domain will be converted to the displacement boundary condition of corresponding mesoscopic model by downscaling, as shown in Figure 2. Each IEE determines the nodal displacement of mesoscopic model’s boundary within its range. As plotted in Figure 2, for mesoscopic node p within the range of IEE k, its nodal displacement ${\mathit{d}}_{\mathrm{mesoscopic}}^p$ can be expressed as

$${\mathit{d}}_{\mathrm{mesoscopic}}^p={\mathit{u}}^{\left(k\right)}\left({(•)}_p\right)$$

(12)

Here, $(•{)}_p$ refers to the effective coordinate corresponding to mesoscopic node p.

The solutions of mesoscopic problem will be solved by governing equation (13), which is the governing equation of mesoscopic computational level

$${\mathit{f}}^{\mathrm{int}}\left({\mathit{d}}_{\mathrm{mesoscopic}}\right)={\mathit{f}}^{\mathrm{ext}}$$

(13)

where ${\mathit{d}}_{\mathrm{mesoscopic}}$ , ${\mathit{f}}^{\mathrm{int}}$ , and ${\mathit{f}}^{\mathrm{ext}}$ denote the nodal displacement of mesoscopic model, mesoscopic internal nodal forces, and mesoscopic external nodal forces, respectively. The Newton–Raphson iteration is adopted to solve equation (13).

The tangent stiffness matrix of mesoscopic model is denoted as ${\mathit{K}}_{\mathrm{mesoscopic}}$ . The static condensation is applied to obtain condensed tangent stiffness matrix. From

$$\left[\begin{array}{cc}{\mathit{K}}_{\mathrm{aa}}\hfill & {\mathit{K}}_{\mathrm{ab}}\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{\mathrm{ba}}}& {\mathit{K}}_{\mathrm{bb}}\hfill \end{array}\right]\left[\begin{array}{c}\delta {\mathit{d}}_{\mathrm{mesoscopic}}^a\hfill \\ \multicolumn{1}{c}{\delta {\mathit{d}}_{\mathrm{mesoscopic}}^b}\end{array}\right]=\left[\begin{array}{c}\delta {\mathit{f}}^a\hfill \\ \multicolumn{1}{c}{0}\end{array}\right]$$

(14)

gives

$$({\mathit{K}}_{\mathrm{aa}}-{\mathit{K}}_{\mathrm{ab}}{\mathit{K}}_{\mathrm{bb}}^{-1}{\mathit{K}}_{\mathrm{ba}})\delta {\mathit{d}}_{\mathrm{mesoscopic}}^a=\delta {\mathit{f}}^a$$

(15)

where ${\mathit{K}}_{\mathrm{mesoscopic}}=\left[\begin{array}{cc}{\mathit{K}}_{\mathrm{aa}}\hfill & {\mathit{K}}_{\mathrm{ab}}\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{\mathrm{ba}}}& {\mathit{K}}_{\mathrm{bb}}\hfill \end{array}\right]$ , $\delta {\mathit{f}}^a$ is the variation in boundary nodal forces of mesoscopic model, and $\delta {\mathit{d}}_{\mathrm{mesoscopic}}^a$ and $\delta {\mathit{d}}_{\mathrm{mesoscopic}}^b$ are variations in boundary and internal nodal displacement, respectively. Thus, the condensed tangent stiffness matrix ${\mathit{K}}_T$ can be expressed as

$${\mathit{K}}_T={\mathit{K}}_{\mathrm{aa}}-{\mathit{K}}_{\mathrm{ab}}{\mathit{K}}_{\mathrm{bb}}^{-1}{\mathit{K}}_{\mathrm{ba}}$$

(16)

The solved boundary nodal forces and condensed tangent stiffness matrix of mesoscopic model will be transferred into macroscopic level for the calculation of domain forces and domain tangent stiffness matrix.

Many types of the mesoscopic models could be adopted to describe the details of material, such as lattice model and FEM. The lattice model is adopted here to establish a mesoscopic model which describes a quasi-brittle material containing distributed initial defects. The model is depicted in Figure 2. The defects are randomly inserted in the model with a relatively uniform and distributed pattern by deleting some of the Lattice elements. The type of mesoscopic element is truss element and model parameters are set according to method mentioned in Guo et al., ⁵⁰ which will not be narrated here.

The Mazars damage law is adopted to describe the degradation of mechanical property of material in mesoscopic model and the damage caused by compressive stress is ignored. The damage evolution can be expressed as

$$D=1-\frac{(1-A_t){\epsilon }_{\mathrm{pt}}}{{\epsilon }_t}-\frac{A_t}{e^{B_t({\epsilon }_t-{\epsilon }_{\mathrm{pt}})}}$$

(17)

where ${\epsilon }_{\mathrm{pt}}$ is the strain threshold for damage, ${\epsilon }_t$ is the tensile strain, and $A_t$ and $B_t$ are material parameters, respectively. Let

$$f({\epsilon }_t,\kappa )={\epsilon }_t-\kappa \le 0$$

(18)

The loading criterion can be

$$f=0,\stackrel{·}{D}>0$$

(19)

Equation (18) represents the damage surface and $\kappa$ denotes the variable relating to load history, which represents the maximum that a tensile strain has ever reached during damage evolution. The initial value ${\kappa }_0={\epsilon }_{\mathrm{pt}}$ . The critical value $D_C$ of damage variable is set to 0.5. When $D_C$ is reached in a lattice element, the stiffness of that element will be set to a quite small value to represent the failure of it and the damage variable of the element will be marked with 1.0. The material parameters are shown in Table 1 according to Mazars and Pijaudiercabot ⁵¹ as follows.

Table 1.

Material parameters.

	Young’s modulus E (GPa)	Poisson’s ratio v	${\epsilon }_{\mathrm{pt}}$ (10−6)	At	Bt
Material parameters	30	0.20	100	0.80	20,000

Quadrature scheme and upscaling

In this section, the IEE based on nonlocal quadrature ⁵² is proposed and the domain forces and domain tangent stiffness matrix are derived.

The numerical quadrature scheme is necessary to obtain nodal forces of IEE as expressed in equation (8). If Gauss quadrature is adopted, there may not be any mesoscopic node that happens to correspond to the position of quadrature point and the value on Gauss point should be the one with macroscopic meaning. So, it cannot be the same as that one of corresponding mesoscopic node directly. However, the degrees of freedom are usually considerable for mesoscopic model. Therefore, the calculation efficiency of domain forces and domain tangent stiffness matrix should be improved. Here, nonlocal quadrature is applied for that purpose and the IEE with not only quadrature points but also quadrature neighborhood is established.

The nonlocal quadrature of function $\phi (x)$ can be expressed as

$$I=\sum _L{W}_L\frac{1}{2{\delta \prime }_0}\underset{{\xi }_L-{\delta \prime }_0}{\overset{{\xi }_L+{\delta \prime }_0}{\int }}\phi \left(\xi \right)\mathrm{Jd}\xi$$

(20)

$$I=\sum _L{W}_{L}J\frac{1}{2{\delta }_0}\underset{X_L-{\delta }_0}{\overset{X_L+{\delta }_0}{\int }}\phi \left(x\right)\mathrm{dx}$$

(21)

where equations (20) and (21) are described in parent element and physical element, respectively. $\xi$ denotes the natural coordinate. As depicted in Figure 3, L denotes the number of quadrature point, W_L denotes the nonlocal quadrature weight, and X_L and ${\xi }_L$ are positions of quadrature points in physical and parent elements, respectively. ${\Theta }_L$ and $\Theta {\prime }_L$ denote the quadrature neighborhoods in physical and parent elements. ${\delta }_0$ and $\delta {\prime }_0$ are radius of ${\Theta }_L$ and $\Theta {\prime }_L$ , respectively The Jacobian $J=(b-a)/2$ and $\delta {\prime }_0={\delta }_0/J$ .

Figure 3.

Nonlocal quadrature: (a) physical element and (b) parent element.

Different amount of quadrature points leads to different quadrature parameters. The quadrature parameters of nonlocal quadrature with two quadrature points are

$$\begin{array}{c}{\xi }_{\mathrm{I}}=\sqrt{\frac{1}{3}-\frac{1}{3}{\delta }_0{\prime }^2}\hfill \\ {\xi }_{\mathrm{I}\mathrm{I}}=-\sqrt{\frac{1}{3}-\frac{1}{3}{\delta }_0{\prime }^2}\hfill \\ W_{\mathrm{I}}=1\hfill \\ W_{\mathrm{I}\mathrm{I}}=1\hfill \end{array}$$

(22)

It can be noted that, in a sense, average value of function in the quadrature neighborhood is used as the value on quadrature point, which can smooth local oscillation due to mesoscopic features and obtain a value with macroscopic meaning.

Considering the IEEs shown in Figure 4, each of them contains two quadrature points and corresponding quadrature neighborhoods. The radius $\delta {\prime }_0$ is taken as 0.4 mm. According to the relationship between physical coordinate and natural coordinate, the positions of quadrature points and radiuses of quadrature neighborhoods in every IEE are determined. Applying the nonlocal quadrature, equation (8) can be expressed as

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{F}}^{\left(k\right)}\left({\mathit{d}}^{\left[l\right]}\right)}& =\frac{W_{\mathrm{I}}}{2{\delta }_0}\underset{{(•)}_{\mathrm{I}}-{\delta }_0}{\overset{{(•)}_{\mathrm{I}}+{\delta }_0}{\int }}{\left[{\mathit{N}}^{\left(k\right)}(•)\right]}^T{\mathit{t}}_{\mathrm{mesoscopic}}^{\left[l\right],\left(k\right)}(•)\mathrm{Jds}+\frac{W_{\mathrm{I}\mathrm{I}}}{2{\delta }_0}\underset{{(•)}_{\mathrm{I}\mathrm{I}}-{\delta }_0}{\overset{{(•)}_{\mathrm{I}\mathrm{I}}+{\delta }_0}{\int }}{\left[{\mathit{N}}^{\left(k\right)}(•)\right]}^T{\mathit{t}}_{\mathrm{mesoscopic}}^{\left[l\right],\left(k\right)}(•)\mathrm{Jds}\hfill \\ \multicolumn{1}{c}{}& =\frac{W_{\mathrm{I}}}{2{\delta }_0}J\sum _{p\in {\Theta }_{\mathrm{I}}}{\left[{\mathit{N}}^{\left(k\right)}{(•)}_p\right]}^T{\mathit{t}}_F^p\left({\mathit{d}}^{\left[l\right]}\right)+\frac{W_{\mathrm{I}\mathrm{I}}}{2{\delta }_0}J\sum _{p\in {\Theta }_{\mathrm{I}\mathrm{I}}}{\left[{\mathit{N}}^{\left(k\right)}{(•)}_p\right]}^T{\mathit{t}}_F^p\left({\mathit{d}}^{\left[l\right]}\right)\hfill \end{array}$$

(23)

where ${\mathit{F}}^{(k)}({\mathit{d}}^{\left[l\right]})$ denotes the approximation of nodal forces of IEE k, which belong to material domain l. The mesoscopic traction ${\mathit{t}}_{\mathrm{mesoscopic}}^{\left[l\right],\left(k\right)}(•)$ is located inside the range of IEE k of material domain l. $(•{)}_{\mathrm{I}}$ and $(•{)}_{\mathrm{I}\mathrm{I}}$ denote the coordinates of quadrature points in IEE. ${\mathit{N}}^{\left(k\right)}(•{)}_p$ and ${\mathit{t}}_F^p({\mathit{d}}^{\left[l\right]})$ are shape function of each mesoscopic node in quadrature neighborhood and the mesoscopic traction of that node, respectively. Here, p denotes the number of the mesoscopic node. The quadrature neighborhoods are ${\Theta }_{\mathrm{I}}$ and ${\Theta }_{\mathrm{I}\mathrm{I}}$ in IEE. The mesoscopic nodes are usually distributed densely on the mesoscopic model’s boundary. So, the integral defined on quadrature neighborhoods in equation (23) is simplified to the discrete sum. It should be noted that not all the mesoscopic traction on mesoscopic model’s boundary are utilized to form the nodal forces of IEE. The nonlocal quadrature reduces the necessary information extraction quantity from mesoscopic level. Comparing with the efficiency improvement of calculating the nodal forces of IEE, the efficiency improvement of calculating domain tangent stiffness matrix is more significant for the reason that it simplifies the assembly process, especially when the mesoscopic model is quite fine. One can note that the tractions on mesoscopic model’s boundary are converted to the values with macroscopic meaning on quadrature points in IEEs, which realizes the upscaling from mesoscopic model to material domain.

Figure 4.

IEEs based on nonlocal quadrature.

The quadrature weights are same when two quadrature points are adopted. Let $W=W_{\mathrm{I}}=W_{\mathrm{I}\mathrm{I}}$ , equation (23) can be expressed as

$${\mathit{F}}^{\left(k\right)}({\mathit{d}}^{\left[l\right]})=\frac{W}{2{\delta }_0}J\sum _{p\in \Theta }{\left[{\mathit{N}}^{\left(k\right)}{(•)}_p\right]}^T{\mathit{t}}_F^p({\mathit{d}}^{\left[l\right]})$$

(24)

where $\Theta ={\Theta }_{\mathrm{I}}\cup {\Theta }_{\mathrm{I}\mathrm{I}}$ denotes the whole quadrature neighborhood in IEE.

Inserting equation (24) into equation (9), which gives domain force ${\mathit{F}}_i^{\left[l\right]}$

$$\begin{array}{cc}\multicolumn{1}{c}{{\mathit{F}}_i^{\left[l\right]}=}& \frac{W}{2{\delta }_0}J\sum _{p\in {\Theta }_{k_i^1}}{\left[{\mathit{N}}_i^{\left(k_i^1\right)}{(•)}_p\right]}^T{\mathit{t}}_F^p({\mathit{d}}^{\left[l\right]})\hfill \\ \multicolumn{1}{c}{}& +\frac{W}{2{\delta }_0}J\sum _{p\in {\Theta }_{k_i^2}}{\left[{\mathit{N}}_i^{\left(k_i^2\right)}{(•)}_p\right]}^T{\mathit{t}}_F^p({\mathit{d}}^{\left[l\right]})\hfill \end{array}$$

(25)

As depicted in Figure 4, node i is a common node of two adjacent IEEs, $k_i^1$ and $k_i^2$ , on the boundary of material domain l. The ${\Theta }_{k_i^1}$ and ${\Theta }_{k_i^2}$ are quadrature neighborhoods in IEEs, $k_i^1$ and $k_i^2$ . The nodal forces of IEE cannot be determined by the nodal displacement of itself. So, the real corresponding relationship exists in domain forces and domain nodal displacement.

In addition to domain forces, domain tangent stiffness matrix will be obtained from mesoscopic level. The domain tangent stiffness matrix of material domain l is denoted as ${\mathit{K}}_{\mathrm{domain}}^{\left[l\right]}$ . The variation in domain force ${\mathit{F}}_i^{\left[l\right]}$ due to the variation in domain nodal displacement can be expressed as

$$\delta {\mathit{F}}_i^{\left[l\right]}=\frac{W}{2{\delta }_0}J\sum _{j\in \left[l\right]}\left(\begin{array}{c}{\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^1\right)\hfill & \left(r_j^1\right)\backslash {\eta }_{r_j^1}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)+{\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^1\right)\hfill & \left(r_j^2\right)\backslash {\eta }_{r_j^2}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)\hfill \\ \multicolumn{1}{c}{+{\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^2\right)\hfill & \left(r_j^1\right)\backslash {\eta }_{r_j^1}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)+{\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^2\right)\hfill & \left(r_j^2\right)\backslash {\eta }_{r_j^2}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)}\end{array}\right)\delta {\mathit{d}}_j$$

(26)

Here, ${\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^1\right)\hfill & \left(r_j^1\right)\backslash {\eta }_{r_j^1}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)$ denotes a stiffness assembly symbol and the meaning of which will be detailed in appendix. Define

$${\mathit{K}}_{T(i,j)}^{\left[l\right]}=\frac{W}{2{\delta }_0}J\left(\begin{array}{c}{\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^1\right)\hfill & \left(r_j^1\right)\backslash {\eta }_{r_j^1}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)+{\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^1\right)\hfill & \left(r_j^2\right)\backslash {\eta }_{r_j^2}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)\hfill \\ \multicolumn{1}{c}{+{\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^2\right)\hfill & \left(r_j^1\right)\backslash {\eta }_{r_j^1}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)+{\mathit{K}}_T^{\left[l\right],i}\left(\begin{array}{c}\begin{array}{cc}\left(k_i^2\right)\hfill & \left(r_j^2\right)\backslash {\eta }_{r_j^2}\hfill \end{array}\hfill \\ \multicolumn{1}{c}{j}\end{array}\right)}\end{array}\right)$$

(27)

and note that i and j are global numbers of domain nodes, and the domain tangent stiffness matrix can be expressed as

$${\mathit{K}}_{\mathrm{domain}}^{\left[l\right]}=\left[\begin{array}{cccccccc}{\mathit{K}}_{T({\alpha }_1,{\alpha }_1)}^{\left[l\right]}\hfill & \cdots \hfill & \cdots \hfill & \cdots \hfill & {\mathit{K}}_{T({\alpha }_1,{\alpha }_{\beta })}^{\left[l\right]}\hfill & \cdots \hfill & \cdots \hfill & {\mathit{K}}_{T({\alpha }_1,{\alpha }_{\mathrm{nd}})}^{\left[l\right]}\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \hfill & \hfill & ⋮\hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \hfill & \hfill & ⋮\hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \hfill & \hfill & ⋮\hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{T({\alpha }_{\beta },{\alpha }_1)}^{\left[l\right]}}& \cdots \hfill & \cdots \hfill & \cdots \hfill & {\mathit{K}}_{T({\alpha }_{\beta },{\alpha }_{\beta })}^{\left[l\right]}\hfill & \cdots \hfill & \cdots \hfill & {\mathit{K}}_{T({\alpha }_{\beta },{\alpha }_{\mathrm{nd}})}^{\left[l\right]}\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \hfill & \hfill & ⋮\hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{⋮}& \hfill & \hfill & \hfill & ⋮\hfill & \hfill & \hfill & ⋮\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{T({\alpha }_{\mathrm{nd}},{\alpha }_1)}^{\left[l\right]}}& \cdots \hfill & \cdots \hfill & \cdots \hfill & {\mathit{K}}_{T({\alpha }_{\mathrm{nd}},{\alpha }_{\beta })}^{\left[l\right]}\hfill & \cdots \hfill & \cdots \hfill & {\mathit{K}}_{T({\alpha }_{\mathrm{nd}},{\alpha }_{\mathrm{nd}})}^{\left[l\right]}\hfill \end{array}\right]$$

(28)

where nd denotes the amount of nodes of a material domain. In macroscopic model, ${\alpha }_{\beta }$ denotes the global number of $\beta$ th domain node.

The domain tangent stiffness matrix reflects the current mechanical properties of material domain and is equivalent to the function of operator T. The tangent stiffness matrix, ${\mathit{K}}_{\mathrm{tol}}$ , of the whole macroscopic model is assembled from each domain tangent stiffness matrix

$${\mathit{K}}_{\mathrm{tol}}=\sum _l{\mathit{K}}_{\mathrm{domain}}^{\left[l\right]}$$

(29)

The summation symbol, here, represents the assembly of stiffness matrixes.

Damage criterion and two-level computational framework

The damage criterion of material domain will be described in this section and an efficient two-level computational framework will be established. The mechanical property of intact material domain is assumed to be linear elasticity here.

First, a damage criterion is established based on the deformation of material domain’s boundary. Due to the fact that damage criterion relates to the type of materials and the failure modes of material domain, different damage criterions can be established according to different demands. Here, the failure mode dominated by tensile failure is considered and the boundary strain extremum is defined to establish the damage criterion. For the other failure modes or material types, the damage criterion can also be established according to corresponding characteristics. The boundary strain extremum ${\epsilon }_t^{\left[l\right]}$ is defined as

$${\epsilon }_t^{\left[l\right]}=\underset{\left(k\right)\in \left[l\right]}{max}\left({\epsilon }_t^{\left(k\right),\left[l\right]}\right)$$

(30)

where ${\epsilon }_t^{\left(k\right),\left[l\right]}$ denotes the ratio of tensile deformation of IEE k on the boundary of material domain l with respect to its initial length. If the IEE is compressed, the value of ${\epsilon }_t^{\left(k\right),\left[l\right]}$ is set to zero. The boundary strain extremum represents the maximum tensile deformation degree somewhere on the boundary. Therefore, the damage criterion can be expressed as

$${\epsilon }_t^{\left[l\right]}\le {\epsilon }_{t0}$$

(31)

where ${\epsilon }_{t0}$ denotes a parameter that indicates the maximum permissible value of tensile deformation degree on the boundary of material domain. If this criterion is satisfied, the material domain is considered to be damaged. The ${\epsilon }_{t0}$ is taken as $90\mu \epsilon$ .

For the intact material domain, the mechanical properties are described by domain elastic stiffness matrix, which is assembled by initial condensed stiffness matrix of corresponding mesoscopic model. The domain elastic stiffness matrix is a constant and obtained by the same way of domain tangent stiffness matrix. The work is completed in pre-analysis step before analysis.

The computational framework is depicted in Figure 5. There are two separated computational levels. The macroscopic problem is divided into several macroscopic load steps and each of them is solved by Newton–Raphson iteration. During the solving process, if the material domain is intact, the response of which will be obtained according to its initial elastic property directly. If the material domain reaches the damage criterion, the domain information transfer mechanism will be activated and the response of that material domain will be obtained from the upscaling and downscaling processes.

Figure 5.

Flowchart of domain information transfer method.

It should be noted that the damage criterion is a standard to judge the variation in macroscopic mechanical properties of a material domain, such as the yield surface in plasticity theory. It is not the standard for the substitution of some areas in the model from a macroscopic region to a mesoscopic region. The framework of macroscopic model does not change during the whole analysis. The change exists in the macroscopic properties of material domains. Therefore, different from concurrent method, it is unnecessary to consider the complex variation of the boundary between macroscopic and mesoscopic regions due to damage evolution. However, it realizes the preservation of load history automatically. The load history in the former load steps will not be lost. The variables reflecting load history such as damage variable will be stored in mesoscopic models and be updated after each macroscopic load step. This method saves considerable computational cost because of two main reasons. First, the mesoscopic models will only be assigned to the damaged material domains. Second, load history is preserved automatically so that repeated computation in concurrent method is avoided.

Numerical results and discussions

In this section, results obtained with the DITM are shown and compared with reference solutions obtained by direct numerical simulation (DNS) of the problems. The damage processes of two kinds of specimens consisting of quasi-brittle material are simulated. At mesoscopic level, we consider a quasi-brittle material with initial defects as described in section “Downscaling and mesoscopic computational level.” The DNS model is fully constructed by Lattice elements and its internal structure is periodically repeated by a mesoscopic model (the one assigned to a material domain) with size of 50 mm × 50 mm.

Four-point bending test

For a four-point bending beam subjected to displacement load, the DITM is shown in Figure 6(a). The beam is simply supported. The size of the specimen is 450 mm × 150 mm. The problem is simplified to a plane problem and the result of the sum of reaction forces at loading position is expressed as the force acting on unit thickness. In DITM, the investigated specimen consists of 27 material domains. There are 608 degrees of freedom and 330 IEEs on the boundaries of material domains in the structure. The corresponding DNS for the same problem has 1,36,202 degrees of freedom. The size of each material domain is 50 mm × 50 mm and IEEs with size of 10 mm are applied. Considering the area near mid-span is the main damage region that should be concerned, the damage of domains (16) and (27) which are close to supports is ignored.

Figure 6.

(a) The four-point bending beam modeled by DITM and (b) force displacement curves obtained from DITM and DNS.

The calculated force–displacement curves of DITM and DNS are plotted in Figure 6(b). One can see that DITM result agrees well with the DNS analysis. The damage contours at the end of analysis obtained from DITM and DNS are depicted in Figure 7. The damage pattern obtained from DITM is in relatively good agreement with the one obtained from DNS. It should also be noted that the damage profiles obtained by DITM suggest some spurious damage localization at the horizontal edges of the material domains. It shows some mesh dependency here. One reason is that the size of IEE restricts the range of damage area near it. The deformation along the direction of an IEE is similar; so, the damage is similar there, which results in the range of damage area near IEE close to the size of IEE when damage localization happens. Another reason is that the nonlocal interaction does not extend in the transverse direction of IEE although the nonlocal quadrature is adopted. It is should be pointed out that the damage contours obtained from DITM are results on mesoscopic computational level, assembled by each damaged material domain. The material domains which do not meet damage criterion are also described in Figure 7(a) to construct the whole specimen for the comparison of DITM and DNS results.

Figure 7.

Damage contours obtained from (a) DITM and (b) DNS at the end of analysis.

Figure 8 shows the damage evolution process obtained from DITM. The contours describe the whole process from distributed damage-to-damage localization phenomenon. The damage patterns a–h correspond to reference points a–h depicted in force–displacement curve. The right superscript of each contour indicates the number of damaged material domain. The amount of damaged material domains corresponding to each reference point increases with the continuous increase of displacement load. It should be noted that during the procedure between points a and b, damage distributes around the defects over the bottom of specimen near mid-span areas. The damage localization zone occurs at point c, the peak load. After that, it bifurcates into two main damage localization zones which develop simultaneously. These two damage localization zones coalesce into one at point g and form the obvious damage localization band in the specimen. After point h, the damage localization band penetrates the specimen and results in its failure. It can be seen that damage localization band gets across three material domains and makes the damaged material domains beside these three ones unloaded, which stops the damage evolution there.

Figure 8.

Damage evolution in material domains.

During the DITM analysis, domain information transfer mechanism is activated at six material domains sequentially as plotted in Figure 9. Each subfigure corresponds to a point marked in the force–displacement curve as shown in Figure 9. It shows that up to point a, all of the material domains are intact. With damage evolution, the damaged material domains appear from bottom to top of the specimen. One can see that only the areas near mid-span are damaged while other areas retain intact during analysis.

Figure 9.

Activation procedure of domain information transfer mechanism at material domains.

On macroscopic level, the specimen deformation obtained from DITM at the end of analysis is shown in Figure 10. It can be noted that the deformation of most of material domains is relatively small except for three material domains located on mid-span areas. On the locations where damage localization band gets across the boundaries of adjacent material domains, the deformation is obviously larger than that of other locations on material domains’ boundaries, which indicates the local cracking here. It should be noted that the phenomenon of local failure is described on macroscopic level by local severe deformation of boundaries of material domains. In this method, when and where the damage localization band will generate is solved by DITM automatically. The locations where damage localization band gets across boundaries of adjacent material domains are also obtained automatically. It makes it unnecessary to consider how to insert an equivalent damage localization band (including its length, direction, and pattern) into macroscopic element at proper instant, to represent local failure under hierarchical computational framework. The boundary deformation can represent damage pattern of a material domain containing damage localization band such as domains (2), (5), and (8) in Figure 10. It can be seen that domains (5) and (8) are penetrated by damage localization band while domain (2) is not penetrated totally. According to the deformation pattern of damaged material domains, the region of damage localization band can be located in macroscopic model. One should note that it is a new way to equivalently represent the damage localization band on macroscopic level. Simultaneously, the effect of the failure region on the overall performance of structure is precisely reflected by the mechanical properties of material domains which contain them. If the exact pattern of damage localization band needs to be seen, the result obtained from mesoscopic level in DITM can be assembled to show the pattern, such as shown in Figure 10. (The contours here are fracture contours and the mesoscopic elements with damage value of 1.0 are identified in contours.)

Figure 10.

Specimen deformation on macroscopic level (magnified 200 times) and corresponding mesoscopic fracture pattern obtained from DITM at the end of analysis.

Cantilever beam test

A cantilever beam subjected to a uniform vertical load on its right side is analyzed by DITM and DNS. The material domain is the same as applied in section “Four-point bending test.” The DITM of the specimen with size of 250 mm × 100 mm is depicted in Figure 11(a). The model consists of 10 material domains. The calculated force–displacement curves obtained from DITM and DNS are shown in Figure 11(b). One can see that the DITM result is in relatively good agreement with the DNS analysis result with respect to both the slope of the linear part and the trend of the softening branch. The difference in the peak load is due to coarse-scale discretization of material domain’s boundary. The damage patterns obtained from DITM and DNS at the end of analysis are compared as shown in Figure 12. It can be noted that the DITM result agrees well with DNS analysis result except for some locations in the horizontal edges of the material domains, which shows some mesh dependency as mentioned in section “Four-point bending test.”

Figure 11.

(a) The cantilever beam modeled by DITM and (b) force displacement curves obtained from DITM and DNS.

Figure 12.

Damage contours obtained from (a) DITM and (b) DNS at the end of analysis.

The damage evolution in material domains obtained from DITM is shown in Figure 13. It can be seen that between points a and b, damage is distributed near the upper left end of the specimen, where the tensile stress is relatively high. The damage develops in domains (1) and (2) along lengthwise direction of specimen at first, while damage localization band generates and develops after point c, material domain (2) is unloaded and the damage evolution stops there. At point e, the damage localization band gets across the boundary of domains (1) and (6) and enters into domain (6). It develops until the specimen is almost penetrated as plotted in subfigure g.

Figure 13.

Damage evolution in material domains.

The domain information transfer mechanism is activated at only three material domains during the DITM analysis, as depicted in Figure 14. One can note that up to point a, all the material domains are intact. The activation process and locations of damaged material domains are in consistent with the damage evolution process in Figure 13.

Figure 14.

Activation procedure of domain information transfer mechanism at material domains.

The specimen deformation on macroscopic level obtained from DITM at the end of analysis is shown in Figure 15. It is obvious that domain (1) is penetrated by the damage localization band in the left side of it. The severe deformed IEEs on the two horizontal boundaries of this material domain indicate the locations where the damage localization band gets across material domains’ boundaries. It should be noted that the location of the damage localization band in a material domain can be simulated accurately.

Figure 15.

Specimen deformation on macroscopic level (magnified 200 times) and corresponding mesoscopic fracture pattern obtained from DITM at the end of analysis.

Computational cost comparison

As the initial elastic property is pre-analyzed and assigned to each material domain, the responses of material domains are computed directly according to such property before they meet damage criterion. It saves considerable computational cost. In addition, it is not required to specify where the mechanism should be activated before analysis. In the four-point bending test problem, 6 material domains of the total 27 material domains are activated with the domain information transfer mechanism (≈22%). For the cantilever beam test problem, 30% of the material domains are activated with the domain information transfer mechanism. One should note that the domain information transfer mechanism is activated sequentially, which keeps a relatively lower computational cost before the specimen approaching failure. With the same computational resources, the total computational time of DITM and DNS of the four-point bending test problem are 8 and 60 h, respectively, as plotted in Figure 16. It takes 3 and 26 h of the cantilever beam test problem for DITM and DNS, respectively. Here, models 1 and 2 refer to four-point bending specimen and cantilever beam specimen, respectively. It can be seen that the computational cost saving is significant by DITM.

Figure 16.

Computational cost comparison between DITM and DNS.

Conclusion

The DITM is proposed in this article for the purpose of quasi-brittle failure analysis. It contains two separated computational levels. On macroscopic level, the material domains are viewed as elementary substances which construct the whole structure. The macroscopic governing equations are established based on this point of view and boundary displacement of material domains are defined as basic unknowns. Each damaged material domain will be assigned a mesoscale model and the information transfer mechanism from a mesoscale domain to a macroscale domain is established. The damage localization phenomenon is equivalently represented on macroscopic level by boundary deformation of material domains. Simultaneously, details of damage pattern are obtained from mesoscopic level. The IEE based on nonlocal quadrature scheme is proposed and the efficient upscaling is realized by this type of element.

The numerical results show that the proposed method can properly predict the damage pattern and describe the process form distributed damage-to-damage localization. While, it seems to be some mesh dependency with the size of IEE. Also, if the size of material domain will affect the simulation results or if there is a proper value of the size are still an important issue to be studied. Further research will be conducted to solve the problems mentioned here.