A Static Damage Constitutive Model of Concrete Based on Microscopic Damage Mechanism

Abstract

In this article, a microscopic constitutive model is established that includes friction, plastic, and spring elements and has clear physical meaning. The friction unit reflects the mutual friction between crack surfaces, the plastic unit reflects the development of concrete plasticity, and the fracture of the spring unit reflects the formation and expansion of interior cracks in concrete. In addition, the integration of the random field theory into this model uncovers the physical underpinnings of the relationship between concrete’s nonlinearity and randomness. The multi-scale modeling of the concrete static damage constitutive model is then realized once the parameters of the random field are discovered using the macro test results. In order to apply the model’s applicability in finite element programs, a subroutine was ultimately constructed. The experimental data and the anticipated values from the numerical simulation are in good agreement, supporting the model’s realism.

1. Introduction

In recent decades, academics have concentrated on constitutive modeling and nonlinear analysis of concrete materials. There has been significant advancement in the aforementioned problems thanks to the efforts of researchers. However, given the complexity of concrete materials, this subject still presents many difficulties [1,2,3,4]. The primary characteristic of concrete materials is the non-uniformity of its internal structure. The concrete material’s stress–strain curve exhibits a clear nonlinear connection when it is subjected to external load [5,6,7]. The cause is the nucleation and propagation of microcracks and microdefects in concrete materials due to external loads, which causes localized stress concentration and degradation of the stress–strain curve. Therefore, due to the widespread use of continuum damage mechanics (CDM) to describe and model the issues with concrete materials over the past 20 years, numerous significant advancements have been made [8,9,10,11,12]. It can accurately realize the building of concrete material constitutive relationships since it is based on the developed irreversible thermodynamics theory system. The special properties of concrete materials, such as strain hardening, strain softening, residual stress, etc., are effectively reflected in the theoretical results based on continuum damage mechanics, and the engineering application is also quite successful. However, the majority of continuum damage mechanics-based theoretical models can only employ the empirical damage development equation to derive the damage variables, which results in the lack of parameters with a clear physical meaning.

Therefore, because of the precise physical description of their properties, fiber bundle models (FBMs) have received increased attention within the context of damage mechanics. This model was created using random properties and micromechanics of materials. It contends that throughout the creation process of concrete materials, numerous distributed microcracks and pores gradually occur. These nucleated microcracks and pores will cause the formation of random damage when subjected to load. The FBM first appeared in the paper published by Peirce in 1926 [13], and then Daniels [14] made a detailed static analysis of it. Later, some researchers also further developed the FBMs [15,16]. However, in the above research, only the elastic damage (stiffness degradation) caused by fiber fracture in concrete is considered, while the plastic damage is ignored. In 1982, Krajcinovic and Silva [17] used the joint probability density method to couple the elastic-plastic damage in concrete materials and proposed a model represented by a series of parallel micro units with the same stiffness. The elastic damage variable of the model is defined as the fracture probability of microscopic units, while plasticity is characterized by the units in the toughness model proposed by Iwan [18].

After that, researchers have also studied and developed FBMs, considering plastic damage [19,20,21,22,23,24,25]. For example, reference [19] proposed a kind of plastic fiber bundle model (PFBM) that can still bear a part of the force less than its breaking strength after the spring chain is broken. Le [20] et al. proposed a stochastic microscopic damage model and predicted the static strength and crack development of concrete materials under load. Based on the fiber chain element proposed by Li and Ren [12], Feng et al. [21] connected the plastic slider in series to characterize the hysteretic properties of concrete materials. Yu and Shan [22] introduced an irrecoverable plastic element to characterize the elastic damage and irrecoverable damage characteristics. Chen and Li [23] conducted dimensionality reduction parameter identification based on constrained optimization problems for the model proposed in reference [12]. Guo and Kuang [24] proposed a uniaxial stochastic micromechanical damage model to describe the elastoplastic response of quasi-brittle materials subjected to monotonic and cyclic loading. Voyiadjis et al. [25] developed an anisotropic elasto-plastic damage model for quasi-brittle materials within the thermodynamics framework.

However, there are some defects in the above models. The model in the literature [19] lacks clear physical meaning for the characterization of plasticity in concrete materials. In the model proposed by Feng et al. [22], the plastic slider is only used to characterize the macroscopic plastic strain and stress, and the corresponding result is based on the empirical formula, which leads to the low accuracy of the calculated result. In the model proposed by Shan and Yu [23], the damage evolution equation is determined by assuming a probability density cumulative distribution function with multiple coefficients. Due to the complexity of their establishment methods, other FBMs considering plastic damage are difficult to determine in numerical simulation. The method proposed in reference [24] requires a large number of uniaxial tension–compression complete curves of concrete with different strengths. The model proposed in reference [25] cannot represent the hysteresis loading with minor cyclic loops and is limited to one dimension.

Therefore, based on the above conclusions, a coupled microscopic model of elastic-plastic damage has been established, which can consider the plastic development of concrete and the generation, development, and interaction of internal cracks in concrete. At the same time, random field theory was introduced into the microscopic model to reveal the physical essence of the coupling between nonlinearity and randomness in concrete. This model can be used to describe the physical properties of concrete materials under uniaxial and biaxial loading and can be better used in the mechanical analysis of concrete structures.

2. Constitutive Model of Static Damage of Concrete

2.1. Proposal of Meso Physical Model

Using damage mechanics and meso statistical mechanics, a new meso physical model is created to explain the nonlinearity and unpredictability of concrete materials. The model consists of three elements: a spring element, friction element, and plastic element, as shown in Figure 1. The spring element is compressed and deformed by an external tension when under the influence of an external load. When the spring element reaches the fracture threshold and the spring snaps, the friction element begins to bear tension and generate sliding friction. The introduction of plastic components is a sign of the expansion of material plasticity. As the spring element draws closer to the yield threshold, the plasticity starts to slide. In addition, this work also introduces the random field theory to reflect the unpredictable nature of microfracture and element production. It is assumed that the yield threshold of each plastic element and the fracture threshold of each spring element both follow a random distribution. In other words, the way a micro-spring quits working and the way plastic yields are both random events.

The model has three switches that regulate the operation of various elements. When switch 1 is closed and switches 2 and 3 are open, only the spring element operates in the entire model to represent the linear elastic component of concrete. When switch 2 is closed and switches 1 and 3 are open, the series element made up of a spring element and a plastic element acts in the model. At this time, the plastic element slides, which indicates the plastic development of concrete materials. When switch 3 is closed and switches 1 and 2 are opened, only the friction element acts in the whole model to express the residual stress caused by the friction between crack surfaces after concrete cracking. So, when the concrete material is in the ideal initial lossless state, the model is in the state of switch 1 closed and switch 2 and 3 opened. Then, with the increase in external load, switches 2 and 3 are triggered, respectively. In the whole force process, only one of the three switches of the model is closed, and the other two are opened.

It is assumed that the threshold of spring fracture in the model is ${\Delta }_r$, and the yield threshold of plastic sliding chain is ${\Delta }_y$. When element strain $\epsilon \le {\Delta }_r,\epsilon \le {\Delta }_y$, the spring element does not break, and the plastic element does not yield. At this time, only the spring element is in the working state, that is, switch 1 is closed and switches 2 and 3 are opened. So, the control function of switch 1 can be expressed as

$$S1=H({\Delta }_r-\epsilon )H({\Delta }_y-\epsilon )$$

(1)

where $H(\cdot )$ denotes Heaviside function:

$$H(x)=\left\{\begin{array}{l}0 (x\le 0)\\ 1 (x>0) \end{array}\right.$$

(2)

So, the Formula (1) can be expressed as

$$S1=H({\Delta }_r-\epsilon )H({\Delta }_y-\epsilon )=\left\{\begin{array}{l}1 \mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{c}\mathrm{h} \mathrm{o}\mathrm{n}\\ 0 \mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{c}\mathrm{h} \mathrm{o}\mathrm{f}\mathrm{f}\end{array}\right.$$

(3)

When $\epsilon \ge {\Delta }_y,\epsilon \le {\Delta }_r$, the spring element does not break, and the plastic element begins to yield. At this time, the series element composed of spring element and plastic element in the model acts, that is, switch 2 is closed, and switches 1 and 3 are opened. So, the control function of switch 2 can be expressed as

$$S2=H(\epsilon -{\Delta }_y)H({\Delta }_r-\epsilon )=\left\{\begin{array}{l}1 \mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{c}\mathrm{h} \mathrm{o}\mathrm{n}\\ 0 \mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{c}\mathrm{h} \mathrm{o}\mathrm{f}\mathrm{f}\end{array}\right.$$

(4)

When $\epsilon \ge {\Delta }_r,\epsilon \le {\Delta }_y$, the spring element is broken, and the friction element is sliding. At this time, only the friction element in the model is in the working state, that is, switch 3 is closed and switches 1 and 2 are opened. So, the control function of switch 3 can be expressed as

$$S3=H(\epsilon -{\Delta }_r)H({\Delta }_y-\epsilon )=\left\{\begin{array}{l}1 \mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{c}\mathrm{h} \mathrm{o}\mathrm{n}\\ 0 \mathrm{s}\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{c}\mathrm{h} \mathrm{o}\mathrm{f}\mathrm{f}\end{array}\right.$$

(5)

2.2. Mechanical Characteristics of Single Element

According to the above definition, the stress–strain relationship of each element in the model can be expressed as follows:

(1)

The spring element:

$${\sigma }_e=H({\Delta }_r-\epsilon )H({\Delta }_y-\epsilon )E\epsilon$$

(6)

where E denotes the elastic modulus.

(2)

The plastic element

It has been pointed out in the literature [26] that the bearing capacity of concrete material will gradually decrease with the increase in plastic strain after crossing the peak point of the stress–strain curve. The law can be expressed as

$$\frac{d\sigma }{d{u}_p}=-K$$

(7)

where $u_p$ denotes the plastic deformation, and K is a constant. The relationship between $u_p$ and plastic strain can be expressed as

$${\epsilon }_p=\frac{u_p}{L}$$

(8)

where $L$ denotes the length of the specimen. So, the Formula (7) can be expressed as

$$\frac{d\sigma }{d{\epsilon }_p}=-KL=-Y$$

(9)

In order to describe the fact that the bearing capacity of concrete gradually decreases with the development of plasticity, the model assumes that when $\epsilon \ge {\Delta }_y$, the spring element shrinks with the sliding of plasticity. At this time, the stress of the series element composed of the spring element and plastic element decreases with the increase in strain. For simplicity, it is assumed that the model stress decreases linearly with the increase in strain, and the stress characteristics of the element are shown in Figure 2.

Therefore, the stress–strain relationship of the plastic element can be expressed as

$${\sigma }_p=[E{\Delta }_y+Y(\epsilon -{\Delta }_y\left)\right]H[E{\Delta }_y+Y(\epsilon -{\Delta }_y\left)\right]H(\epsilon -{\Delta }_y)H({\Delta }_r-\epsilon )$$

(10)

(3)

The friction element

According to the literature [27], the relationship between the residual stress ${\sigma }_r$ and the fracture stress ${\sigma }_c$ of concrete can be expressed as follows:

$${\sigma }_r={\eta }_s{\sigma }_c$$

(11)

where ${\eta }_s$ denotes the shear lag coefficient. So, from Equations (6) and (11), it can be concluded that

$${\sigma }_r={\eta }_s{\sigma }_c={\eta }_{s}E{\Delta }_{r}H(\epsilon -{\Delta }_r)H({\Delta }_y-\epsilon )$$

(12)

Therefore, the stress characteristics of the friction block element can be shown in Figure 3.

For the sake of simplification, the following assumptions can be made:

$$\begin{array}{l}{\phi }_1=H\left({\Delta }_r-\epsilon \right)H\left({\Delta }_y-\epsilon \right)\\ {\phi }_2=H\left({\Delta }_r-\epsilon \right)H\left(\epsilon -{\Delta }_y\right)H\left[E{∆}_y+Y\left(\epsilon -{\epsilon }_y\right)\right]\\ {\phi }_3=H(\epsilon -{\Delta }_r)H({\Delta }_y-\epsilon )\end{array}$$

(13)

Therefore, the total stress–strain relationship of a single meso physical model can be expressed as

$$\sigma ={\phi }_{1}E\epsilon +{\phi }_2[E{\Delta }_y+Y(\epsilon -{\Delta }_y\left)\right]+{\phi }_3{\eta }_{s}E{\Delta }_r$$

(14)

2.3. Parallel Bar System

On the basis of the previous section, N meso units are combined into a parallel system to simulate the stress of concrete, as shown in Figure 4.

Therefore, the stress of the kth meso element ${\sigma }^k$ can be expressed as

$${\sigma }^k=\sigma \left(x_k\right),k=\mathrm{1,2},3\dots \dots N$$

(15)

where $x_k$ denotes the space coordinates of the kth meso element, and N denotes the total number of meso units.

Assuming that each meso element has the same section form, the stress of N parallel meso elements can be expressed as

$${\sigma }_e=\frac{1}{N}{\sum }_{k=1}^{N}H\left[{\Delta }_r\right(x_k)-\epsilon ]H\left[{\Delta }_y\right(x_k)-\epsilon ]E\epsilon$$

(16)

$${\sigma }_p=\frac{1}{N}{\sum }_{k=1}^N\left\{\begin{array}{l}\left[E{\Delta }_y\right(x_k)+Y\epsilon -Y{\Delta }_y(x_k\left)\right]\times H\left[E{\Delta }_y\right(x_k)+Y\epsilon -Y{\Delta }_y(x_k\left)\right]\\ \times H[\epsilon -{\Delta }_y(x_k\left)\right]\times H\left[{\Delta }_r\right(x_k)-\epsilon ]\end{array}\right\}$$

(17)

$${\sigma }_r=\frac{1}{N}{\sum }_{k=1}^N{\eta }_{s}E{\Delta }_r\left(x_k\right)H[\epsilon -{\Delta }_r(x_k\left)\right]H\left[{\Delta }_y\right(x_k)-\epsilon ]$$

(18)

When $N\to \infty$, the above formula can be expressed as follows:

$${\sigma }_e={\int }_0^{1}H({\Delta }_r-\epsilon )H({\Delta }_y-\epsilon )E\epsilon dx$$

(19)

$${\sigma }_p={\int }_0^1[E{\Delta }_y+Y(\epsilon -{\Delta }_y\left)\right]H[E{\Delta }_y+Y(\epsilon -{\Delta }_y\left)\right]H(\epsilon -{\Delta }_y)H({\Delta }_r-\epsilon )dx$$

(20)

$${\sigma }_r={\int }_0^1{\eta }_{s}E{\Delta }_{r}H(\epsilon -{\Delta }_r)H({\Delta }_y-\epsilon )dx$$

(21)

From a single meso element to a discrete system with N meso elements, and then to a continuous system with infinite meso elements, the mechanical characteristics of the whole process model are shown in Figure 5.

In addition, the concept of randomness is introduced to reflect the uncontrollable nature of micro-fracture and yield. In this model, it is considered that the fracture threshold ${\Delta }_r$ of spring element satisfies $0\le {\Delta }_r\le {\Delta }_r^{\max }$, which probability density function is $f_r\left({\Delta }_r\right)$ and distribution function is $F_r\left({\Delta }_r\right)={\int }_0^{+\infty }{f}_r\left({\Delta }_r\right)d{\Delta }_r$. It is considered that the yield threshold ${\Delta }_y$ of plastic element satisfies $0\le {\Delta }_y\le {\Delta }_y^{\max }$, which probability density function is $f_y\left({\Delta }_y\right)$ and distribution function is $F_y\left({\Delta }_y\right)={\int }_0^{+\infty }{f}_y\left({\Delta }_y\right)d{\Delta }_y$. Additionally, $f\left({\Delta }_r,{\Delta }_y,\gamma \right)$ is the two-dimensional joint probability density function of ${\Delta }_r$ and ${\Delta }_y$. $f_r\left({\Delta }_r\right)$ and $f_y\left({\Delta }_y\right)$ are two independent probability distributions, as shown in Figure 6.

The expected values of ${\sigma }_e$ are calculated as follows:

$$\begin{gathered} \mu \left({\sigma }_e\right)=\mu \left[{\int }_0^{1}H({\Delta }_r-\epsilon )H({\Delta }_y-\epsilon )E\epsilon d\epsilon \right] \\ ={\int }_0^1\mu \left[H\right({\Delta }_r-\epsilon \left)H\right({\Delta }_y-\epsilon )E\epsilon ]d\epsilon \end{gathered}$$

(22)

From the commutativity of the expectation operator and the integral operator, it can be seen that

$$\begin{array}{l}\mu \left[H\left({\Delta }_r-\epsilon \right)H\left({\Delta }_y-\epsilon \right)E\epsilon \right]\\ =1\times P\left({\phi }_1=1\right)\times E\epsilon +0\times P\left({\phi }_1=0\right)\times E\epsilon \\ =P\left({\phi }_1=1\right)\times E\epsilon \\ =P\left(\epsilon \le {\Delta }_r\cap \epsilon \le {\Delta }_y\right)\times E\epsilon \\ =E\epsilon {\int }_{\epsilon }^{+\infty }{\int }_{\epsilon }^{+\infty } f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{r}d{\Delta }_y\end{array}$$

(23)

It can be assumed that

$$H_1={\int }_{\epsilon }^{+\infty }{\int }_{\epsilon }^{+\infty } f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{r}d{\Delta }_y$$

(24)

So,

$$\mu ({\sigma }_e)={\int }_0^1\left[{\int }_{\epsilon }^{+\infty }{\int }_{\epsilon }^{+\infty } f({\Delta }_r,{\Delta }_y;\gamma )E\epsilon d{\Delta }_{r}d{\Delta }_y\right] dx=H_{1}E\epsilon$$

(25)

Similarly, the expected values of ${\sigma }_r$ and ${\sigma }_y$ can be expressed as

$$\begin{array}{ll}\mu ({\sigma }_r)& =\mu \left[{\int }_0^1{\eta }_{s}E{\Delta }_{r}H(\epsilon -{\Delta }_r)H({\Delta }_y-\epsilon )dx\right]\\ & ={\int }_0^1\mu \left[{\eta }_{s}E{\Delta }_{r}H\right(\epsilon -{\Delta }_r\left)H\right({\Delta }_y-\epsilon \left)\right]dx\\ & ={\eta }_{s}E{\int }_0^1\left[{\int }_{\epsilon }^{+\infty }{\int }_0^{\epsilon }{\Delta }_{r}f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{r}d{\Delta }_y\right] dx\end{array}$$

(26)

It can be assumed that

$$H_2 ={\int }_{\epsilon }^{+\infty }{\int }_0^{\epsilon }{\Delta }_{r}f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{r}d{\Delta }_y$$

(27)

So,

$$\mu ({\sigma }_r)={\eta }_{s}E{\int }_0^1 H_{2}dx={\eta }_{s}E{H}_2$$

(28)

$$\begin{array}{ll}\mu ({\sigma }_p)& =\mu \left\{ {\int }_0^1\left[E{\Delta }_y+Y(\epsilon -{\Delta }_y)\right]H\left[E{\Delta }_y+Y(\epsilon -{\Delta }_y)\right]H(\epsilon -{\Delta }_y)H({\Delta }_r-\epsilon ) \right\} dx\\ & ={\int }_0^1\mu \left\{ \left[E{\Delta }_y+Y(\epsilon -{\Delta }_y)\right]H\left[E{\Delta }_y+Y(\epsilon -{\Delta }_y)\right]H(\epsilon -{\Delta }_y)H({\Delta }_r-\epsilon ) \right\} dx\\ & ={\int }_0^1\left\{ {\int }_{\epsilon }^{+\infty }{\int }_0^{\epsilon }\left[E{\Delta }_y+Y(\epsilon -{\Delta }_y)\right]H\left[E{\Delta }_y+Y(\epsilon -{\Delta }_y)\right] f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{y}d{\Delta }_r\right\} dx\end{array}$$

(29)

It can be assumed that

$$H_3={\int }_{\epsilon }^{+\infty }{\int }_0^{\epsilon }\left[{\Delta }_y+Y(\epsilon -{\Delta }_y)/E\right]H\left[{\Delta }_y+Y(\epsilon -{\Delta }_y)/E\right] f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{y}d{\Delta }_r$$

(30)

Then, Equation (21) can be simplified as follows:

$$\mu ({\sigma }_p) =E{\int }_0^1 H_3 dx=H_{3}E$$

(31)

Then, the expected value of total stress $\sigma$ can be expressed as follows:

$$\begin{array}{ll}\mu \left(\sigma \right)& =\mu \left({\sigma }_e\right)+\mu \left({\sigma }_p\right)+\mu \left({\sigma }_r\right)\\ & =H_{1}E\epsilon +H_{2}E{\eta }_s+H_{3}E\end{array}$$

(32)

2.4. Damage Definition

Each meso element in this parallel meso model has the same spring stiffness, and the plastic yield and fracture thresholds of the spring are random distribution field functions. Concrete is considered to have experienced elastic damage when a spring element fractures. Concrete is thought to have suffered plastic degradation from plastic elements sliding. Then, the total damage of concrete is defined as the ratio of the damaged area to the total area, which can be expressed as follows:

$$d=\frac{A_d}{A}$$

(33)

$$d_e(\epsilon )=\frac{1}{N}{\sum }_{k=1}^{N}H[\epsilon -{\Delta }_r(x_k)]\cdot H[{\Delta }_y(x_k)-\epsilon ]$$

(34)

$$d_p(\epsilon )=\frac{1}{N}{\sum }_{k=1}^{N}H[\epsilon -{\Delta }_y(x_k)]\cdot H[{\Delta }_r(x_k)-\epsilon ]$$

(35)

When $N\to \infty$, the above formulas can be expressed as

$$d_e(\epsilon )={\int }_0^{1}H[\epsilon -{\Delta }_r(x)]\cdot H[{\Delta }_y(x)-\epsilon ] dx$$

(36)

$$d_p(\epsilon )={\int }_0^{1}H[\epsilon -{\Delta }_y(x)]\cdot H[{\Delta }_r(x)-\epsilon ] dx$$

(37)

$$d(\epsilon )=d_e(\epsilon )+d_p(\epsilon )$$

(38)

After considering the damage, the constitutive model can be expressed as

$$\sigma =(1-d)E\epsilon +{\int }_0^1{\phi }_2[E{\Delta }_y+Y(\epsilon -{\Delta }_y)] dx+{\int }_0^1{\phi }_3{\eta }_{s}E{\Delta }_{r}dx$$

(39)

The mean values of damage variables $d_e$, $d_p$ and $d$ can be calculated as follows, respectively:

$$\begin{array}{ll}\mu (d_e)& =\mu \left\{{\int }_0^{1}H[\epsilon -{\Delta }_r(x)]\cdot H[{\Delta }_y(x)-\epsilon ] dx\right\}\\ & ={\int }_{\epsilon }^{+\infty }{\int }_0^{\epsilon }f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{r}d{\Delta }_y\end{array}$$

(40)

$$\begin{array}{ll}\mu (d_p)& =\mu \left\{{\int }_0^{1}H[\epsilon -{\Delta }_y(x)]\cdot H[{\Delta }_r(x)-\epsilon ] dx\right\}\\ & ={\int }_{\epsilon }^{+\infty }{\int }_0^{\epsilon }f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{y}d{\Delta }_r\end{array}$$

(41)

$$\begin{array}{ll}\mu (d)& =\mu \left\{{\int }_0^{1}H[{\Delta }_y(x)-\epsilon ]\cdot H[{\Delta }_r(x)-\epsilon ] dx\right\}\\ & =1-{\int }_{\epsilon }^{+\infty }{\int }_{\epsilon }^{+\infty }f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{r}d{\Delta }_y\end{array}$$

(42)

Therefore, $\mu \left(\sigma \right)$ can be expressed as follows:

$$\mu (\sigma )=[1-\mu (d)]E\epsilon +H_{2}E{\eta }_s+H_{3}E$$

(43)

When $\mu (d)=0$, $\mu (\sigma )=E\epsilon$, which means it is in a non-destructive state. When $\mu (d)=1$, the stress of all plastic elements in the model is 0, and $\mu (\sigma )=H_{3}E$, which means it is in a completely damaged state.

2.5. Biaxial Static Damage Constitutive Model of Concrete

Under the condition of two-dimensional plane stress (${\sigma }_3=0$), the full variable form of the damage constitutive model of concrete can be expressed as follows:

$$\mu \left(\mathit{\sigma }\right)=[\mathit{I}-\mu (\mathit{d}\left)\right]:\mathit{E}:\mathit{\epsilon }+{\eta }_s{H}_2\left(\mathit{\epsilon }\right):\mathit{E}+H_3\left(\mathit{\epsilon }\right):\mathit{E}$$

(44)

where E denotes the initial material stiffness tensor of concrete, which can be expressed as follows:

$$\left[\mathit{E}\right]=\frac{E}{1-{{\nu }_0}^2}\left[\begin{array}{ccc}1& {\nu }_0& 0\\ {\nu }_0& 1& 0\\ 0& 0& (1-{{\nu }_0}^2)G/E\end{array}\right]$$

(45)

where ${\nu }_0$ and $G$ denotes Poisson’s ratio and the initial shear modulus of concrete.

Then, for different stress states, the damage constitutive relation of concrete can be expressed as follows:

(1)

Biaxial tensile stress state (${\sigma }_1>0,{\sigma }_2>0$):

$$\mu \left(\mathit{\sigma }\right)=[\mathit{I}-\mu ({\mathit{d}}^{+}\left)\right]:\mathit{E}:\mathit{\epsilon }+{\eta }_s{H}_2^{+}\left(\mathit{\epsilon }\right):\mathit{E}+H_3^{+}\left(\mathit{\epsilon }\right):\mathit{E}$$

(46)

(2)

Biaxial compression stress state (${\sigma }_1<0,{\sigma }_2<0$):

$$\mu \left(\mathit{\sigma }\right)=[\mathit{I}-\mu ({\mathit{d}}^{-}\left)\right]:\mathit{E}:\mathit{\epsilon }+{\eta }_s{H}_2^{-}\left(\mathit{\epsilon }\right):\mathit{E}+H_3^{-}\left(\mathit{\epsilon }\right):\mathit{E}$$

(47)

(3)

Biaxial tension and compression stress state (${\sigma }_1>0,{\sigma }_2<0$ or ${\sigma }_1<0,{\sigma }_2>0$):

$$\mu \left(\mathit{\sigma }\right)=\sqrt{[\mathit{I}-\mu ({\mathit{d}}^{+}\left)\right]:[\mathit{I}-\mu ({\mathit{d}}^{-}\left)\right]}:\mathit{E}:\mathit{\epsilon }+{\eta }_s\sqrt{H_2^{+}\left(\mathit{\epsilon }\right):H_2^{-}\left(\mathit{\epsilon }\right)}:\mathit{E}+\sqrt{H_3^{+}\left(\mathit{\epsilon }\right):H_3^{-}\left(\mathit{\epsilon }\right)}:\mathit{E}$$

(48)

3. Identification of Random Field Parameters Based on Global Genetic Algorithm

Now, the constitutive curve of concrete is established by the method of meso random damage of concrete. The nonlinear development process of concrete macro specimens involves a complex process of loading–random damage–stress redistribution–damage evolution, as can be shown from the examination of the meso physical mechanism. Therefore, in order to model the link between meso stress and strain, the concept of multi-scale analysis is developed. This system has three fundamental scales: the micro-, meso-, and macro-scales. The stress–strain connection is determined at the micro-scale by random fracture and yield processes. Then, utilizing meso element modeling, the macroscopic mean stress–strain relationship is provided. System identification (inverse problem solving) can be used to identify the parameters of micro-fracture strain and yield strain random field from the test data of macro specimens, and the probability description of meso damage variables can then be provided. As a result, the probability description of the constitutive relation at the meso level is provided.

Using the modeling method based on the test sample set, the algorithm is divided into two processes:

First, the data of the uniaxial test of a batch of concrete specimens are interval discretized with respect to strain ${\epsilon }_i$, and the average stress value $({\epsilon }_i;\mu ({\sigma }_i^o\left)\right)$ of the sample set at the corresponding point is obtained based on the midpoint of the interval.

Second, the random field parameters of micro-fracture strain and micro yield strain are assumed. Then, the above theoretical model is used to calculate the mean stress $\mu \left({\sigma }_i\right)$ when the strain is ${\epsilon }_i$. Finally, the unknown parameters of the random field are obtained by inverse solution with the minimum square sum of the difference between the corresponding theoretical value and the experimental value as the optimization criterion, which can be shown as

$${\sum }_{i=1}^N\left[\mu \right({\sigma }_i)-\mu ({\sigma }_i^o){]}^2\to min$$

(49)

Because there are many random field parameters to be identified in the model and the calculation is complex, in order to obtain an ideal optimization solution, the global genetic algorithm is used to identify the random field parameters. The calculation flow of random field parameter identification is shown in Figure 7.

4. Model Validation

It is assumed that the fracture threshold ${\Delta }_r\left(x\right)$ of the spring element and the yield threshold ${\Delta }_y\left(x\right)$ of the plastic element obey the lognormal distribution. Then, $Z_r\left(x\right)=\mathit{ln}{\Delta }_r\left(x\right)$ and $Z_y\left(x\right)=\mathit{ln}{\Delta }_y\left(x\right)$ obey normal distribution, which can be shown as $Z_r\left(x\right):N({\mu }_r,{\sigma }_r)$ and $Z_y\left(x\right):N({\mu }_y,{\sigma }_y)$.

Then, standardize $Z_r\left(x\right)$ and $Z_y\left(x\right)$, respectively. Let $\alpha =\left(\mathit{ln}{\mathsf{\Delta }}_r-{\lambda }_r\right)/{\zeta }_r$ and $\beta =\left(\mathit{ln}{\mathsf{\Delta }}_y-{\lambda }_y\right)/{\zeta }_y$, so $\alpha$ and $\beta$ all follow the standard normal distribution, which can be shown as $Z\left(\alpha \right):N\left(\mathrm{0,1}\right)$ and $Z\left(\beta \right):N\left(\mathrm{0,1}\right)$. So, the distribution functions and probability density functions of ${\Delta }_r\left(x\right)$ and ${\Delta }_y\left(x\right)$ ca n be shown as follows:

$$F_r=\phi \left(\frac{\mathit{ln}{\Delta }_r-{\lambda }_r}{{\zeta }_r}\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\frac{\mathit{ln}{\Delta }_r-{\lambda }_r}{{\zeta }_r}}{e}^{-t^2/2} dt$$

(50)

$$F_y=\phi \left(\frac{\mathit{ln}{\Delta }_y-{\lambda }_y}{{\zeta }_y}\right)=\frac{1}{\sqrt{2\pi }}{\int }_{-\infty }^{\frac{\mathit{ln}{\Delta }_y-{\lambda }_y}{{\zeta }_y}}{e}^{-t^2/2} dt$$

(51)

$$f_r({\Delta }_r)=F_r^{\prime }=\frac{1}{\sqrt{2\pi }{\Delta }_r{\zeta }_r}\exp \left[-1/2(\frac{\mathit{ln}{\mathsf{\Delta }}_r-{\lambda }_r}{{\zeta }_r}{)}^2\right]$$

(52)

$$f_y({\Delta }_y)=F_y^{\prime }=\frac{1}{\sqrt{2\pi }{\Delta }_y{\zeta }_y}\exp \left[-1/2(\frac{\mathit{ln}{\mathsf{\Delta }}_y-{\lambda }_y}{{\zeta }_y}{)}^2\right]$$

(53)

According to the above, ${\Delta }_r\left(x\right)$ and ${\Delta }_y\left(x\right)$, respectively, obey one-dimensional lognormal distribution (as shown in Figure 8), and $f({\Delta }_r,{\Delta }_y;\gamma )$ is a two-dimensional lognormal distribution function. When ${\Delta }_r\left(x\right)$ and ${\Delta }_y\left(x\right)$ are independent of each other, it can be determined that

$$f({\Delta }_r,{\Delta }_y;\gamma ) =f_r\left({\Delta }_r\right)f_y\left({\Delta }_y\right)$$

(54)

Then, $\mu \left(d\right)$ can be simplified as follows:

$$\begin{array}{l}\mu (d)=1-{\int }_{\epsilon }^{+\infty }{\int }_{\epsilon }^{+\infty }f({\Delta }_r,{\Delta }_y;\gamma ) d{\Delta }_{r}d{\Delta }_y\\ =1-{\int }_{\epsilon }^{+\infty }{\int }_{\epsilon }^{+\infty }{f}_r\left({\Delta }_r\right) f_y\left({\Delta }_y\right) d{\Delta }_{r}d{\Delta }_y\\ =1-{\int }_{\epsilon }^{+\infty }{\int }_{\epsilon }^{+\infty }\frac{1}{2\pi {\Delta }_r{\Delta }_y}\mathit{exp}\left\{-\frac{1}{2}\left[(\frac{\mathit{ln}{\Delta }_r-{\lambda }_r}{{\zeta }_r}{)}^2+(\frac{\mathit{ln}{\Delta }_y-{\lambda }_y}{{\zeta }_y}{)}^2\right]\right\} d{\Delta }_{r}d{\Delta }_y\\ =1-{\int }_{\epsilon }^{+\infty }\left\{{\int }_{\epsilon }^{+\infty }\frac{1}{\sqrt{2\pi }{\Delta }_r}\mathit{exp} \left[-\frac{1}{2}(\frac{\mathit{ln}{\Delta }_r-{\lambda }_r}{{\zeta }_r}{)}^2\right] d{\Delta }_r\right\}\frac{1}{\sqrt{2\pi }{\Delta }_y}\mathit{exp} \left[-\frac{1}{2}(\frac{\mathit{ln}{\Delta }_y-{\lambda }_y}{{\zeta }_y}{)}^2\right] d{\Delta }_y\end{array}$$

(55)

4.1. Concrete Uniaxial Tension and Compression

In this paper, the random field parameter identification is carried out by using the uniaxial compression and tension test data of self-compacting concrete in reference [28]. The elastic modulus of concrete measured in the test is $E_0=23.22 \mathrm{GPa}$. The parameters to be identified in the model are ${\lambda }_r, {\zeta }_r, {\lambda }_y, {\zeta }_y, Y, {\eta }_s$. According to the algorithm proposed in Section 3, the parameter identification results are obtained, as shown in Figure 9 and Figure 10 and

Table 1. Random field parameters identification results.

Parameters	${\mathit{\lambda }}_{\mathit{r}}$	${\mathit{\zeta }}_{\mathit{r}}$	${\mathit{\lambda }}_{\mathit{y}}$	${\mathit{\zeta }}_{\mathit{y}}$	$\mathit{Y}$	${\mathit{\eta }}_{\mathit{s}}$	R
Compression constitutive	7.653	0.247	7.5	0.18	−0.316	0.077	1.019
Tensile constitutive	4.968	0.496	3.817	0.799	−0.08	0.06	0.613

.

From the above recognition results, it can be seen that the R-value of the objective function of the genetic algorithm is very small, close to 0. Therefore, the model proposed in this paper is in good agreement with the test results, so the rationality of this model can be verified.

Further model validation was carried out based on the experimental results of uniaxial cyclic compression [29] and uniaxial cyclic tension [30], as shown in Figure 11. Compared with the experimental results, it can be seen that the model values are in good agreement with the test values.

4.2. Concrete under Biaxial Compression

By substituting the identified random field parameters into the concrete biaxial constitutive equation in the above summary, the concrete biaxial strength envelope predicted by this model can be obtained. Figure 12 is a comparison diagram of model prediction results and test results [12] (Kupfer et al., 1969). It can be seen from the figure that the two are in good agreement. Therefore, it can be preliminarily verified that the model proposed in this paper can be used for the response and damage analysis of structures under biaxial loads.

Further, according to the concrete biaxial compression tests (${\sigma }_3=0,{\sigma }_1/{\sigma }_2=-1/-1$) conducted by Kupfer et al., the material parameters measured in the tests are elastic modulus E = 31,000 MPa and Poisson’s ratio ν = 0.2. Then, Figure 13 shows the comparison between the predicted biaxial stress–strain response of this model and the experimental results. The comparison results show that the predicted values of the model are in good agreement with the experimental results, which verifies the rationality of the biaxial damage constitutive model of concrete in this paper.

4.3. Analysis of Plain Concrete Structures

The test of simply supported beams with notches in plain concrete under concentrated loads at the mid-span provides a good validation standard for the constitutive model of concrete materials [31]. The size and loading form of the test model are shown in Figure 14, and

Table 2. Material parameters of the specimen.

E (MPa)	ft (MPa)	fc (MPa)	ν	GF (N/m)
27,413	2.886	34.0	0.18	40.29

shows the material parameters of the specimen.

Then, establish a finite element model (Figure 15) using the constitutive model proposed in this article and calculate the load–deflection curve of a simply supported beam at mid-span, as shown in Figure 16. It can be seen that the experimental values are in good agreement with the model values, which verifies the feasibility of applying the model proposed in this paper to the calculation of concrete structures.

5. Conclusions

This article establishes a microscopic constitutive model with clear physical significance, which is composed of a spring element, plastic element, and friction element. The fracture of the spring unit represents the generation and expansion of internal cracks in concrete, the plastic unit represents the development of concrete plasticity, and the friction unit represents the mutual friction between crack surfaces. In addition, the random field theory is introduced into this model, which reveals the physical essence of the coupling of nonlinear and randomness of concrete. Then, the parameters of the random field are identified using the macro test data, and the multi-scale modeling of the concrete static damage constitutive model is realized. Finally, the model was developed as a subroutine to implement its application in finite element programs. The predicted values obtained from numerical simulation are in good agreement with the experimental values, which verify the rationality of the model.