Fracture and Damage Behaviors of Concrete in the Fractal Space

The fracture toughness, the driving force and the fracture energy for an infinite plate with a fractal crack are investigated in the fractal space in this work. The perimeter-area relation is adopted to derive the transformation rule between damage variables in the fractal space and Euclidean space. A plasticity yield criterion is introduced and a damage variable tensor is decomposed into tensile and compressive components to describe the distinct behaviors in tension and compression. A plastic damage constitutive model for concrete in the Euclidean space is developed and generalized to fractal case according to the transformation rule of damage variables. Numerical calculations of the present model with and without fractal are conducted and compared with experimental data to verify the efficiency of this model and show the necessity of considering the fractal effect in the constitutive model of concrete. The structural response and mesh sensitivity of a notched unreinforced concrete beam under 3-point bending test are theoretical studied and show good agreement with the experimental data.


Introduction
Concrete has been widely used in civil engineering for its good in-situ casting and molding abilities.As a quasibrittle material, the fracture behavior of concrete receives much of researchers' concerns.Though classic fracture mechanics which is based on the assumption of smooth cracking in materials can analyze concrete properties and meet the need of structure design in a certain extent, no advance is made to explain the failure mechanism from the change of the concrete internal structure.The difficulty mentioned above has a hindering effect for researchers to improve the mechanical property of concrete.
Fractal geometry is established by Mandelbrot in 1970s [1], which plays an important role in the development of fracture mechanics theory.Researches show that the fracture zone of metal, rock and concrete has fractal characteristics [2][3][4].This leads a widely use of fractal geometry in many fields of material science, for instances, the Sierpinski carpet was adopted by Carpinteri et al. [5] to simulate the composition of concrete cross section, and the fractal effect was also introduced into the cohesive crack model.Another remarkable application of fractal geometry is to describe the roughness of cracks quantitatively.Saouma et al. [6] and Issa et al. [7] investigated the crack profiles of concrete through tests and pointed out that cracks in concrete have an average fractal dimension of 1.1.Meanwhile, Issa et al. [7] analyzed the fracture surface of concrete and found its fractal dimension is about 2.1 to 2.3.
Studies on the damage of concrete point out that the propagation micro defects, i.e. microvoids and microcracks, etc, is the mainly cause of the macro fracture of materials.Based on this fact, a new kind of constitutive model for concrete, called the damage constitutive model, is developed within the framework of continuum damage mechanics.In this model, the choice of damage variable is a key to control the effectiveness and performance.Because of the heterogeneity of concrete, definition of the damage variable still remains at the state that the change of material macro property, such as elastic modulus and stress, are used to reflect the development of damage indirectly, and no direct relation is set with the intrinsic deflects.As the progress in studying fractal phenomena, some researchers try to explore the damage growth by means of fractal geometry.Zhao [8] defined a damage variable as a function of the area of fracture surface, and proposed a new damage constitutive model for rock in which the fractal effect was taken into account; Guarracino [9] gave out a damage variable as the ratio of the porous and REV (or the representative elementary volume) volume fraction of materials and presented a fractal constitutive model for rock.
In this work, the fracture behaviors of a material with fractal cracks are investigated by using fractal geometry.Theoretical expressions of the fracture toughness, the driving force and the fracture energy is derived consequently.The transformation rule of a fractal damage variable in the fractal space and a apparent damage variable in the Euclidean space is obtained by adopting the perimeter-area relation.This rule is introduced into a new plastic damage constitutive model of concrete presented in this research.A notched plain concrete beam under 3-point bending test is simulated to verify the efficiency of the model.

Simplification of Fracture Zone
Figure 1 illustrates an infinite plate with a fractal cut in uniaxial tension.The cut releases the stress in a fracture domain, whose shape can be approximated as an ellipse [10].A standard Koch fractal curve is employed to construct the boundary of the crack, see  , where η is a shape parameter and η=2π when smooth cracking.

Critical Cracking Stress
According to the fractal theory [1], the real length 2a and the apparent length (projected to the axial) 2a 0 of the crack has the following relation: The surface energy of the fracture surface is:  ( ) 4 4 where t is the plate thickness.The perimeter C of the fracture zone is: Thus, one gets the zone area A as follows by adopting the perimeter-area relation where the proportional coefficient m is: Therefore, the strain energy released during the cracking process can be written as: where σ denotes the tensile stress of the plate; E 0 is the elastic modulus.For a plate strain state, one needs to replace E 0 in Equation ( 6) with , where ν 0 is the Poisson's ratio.The Griffith fracture criterion state that: materials cracks when the released elastic energy ΔU equals to the surface energy Π concentrated in the fracture zone for an infinitesimally small increment of the crack length da 0 , i.e., Substituting Equations ( 2) and ( 6) into Equation ( 7), one obtains the critical cracking stress σ c for materials as: For a smooth cut case, D=1.0, η=2π, and we have:

Fracture Toughness
Wunk and Yavari [11] studied the stress field at the fractal crack tip, and presented the component σ y in y direction as shown in Figure 1. where K is the fractal stress intensity factor; α represents the singularity order of the stress field, and can be expressed as: (11) in terms of D for a self-similar fractal crack.
For D=1.0, we have α=1/2 and f I I K K  , and Equation (10) degenerates to the smooth cracking case as: When σ y =σ c , , the crack begins to grow.Referring to Equations ( 8) and (10), one obtains the fractal fracture toughness for materials: dering effect on the cracking of materials.

Driving Force
Classic driving force is defined as the strain energy dissipated to form a unit fracture area [11], and can be expressed as: Submitting Equations ( 1) and (6) into Equation ( 14), one gets: For the mode I cracking case, G f reaches its maximum value max f G when σ increases to the tensile strength σ t of materials: We also have the cracking resistance G IC of materials as: , cracking occurs, and cracks does not grow when max . This is the G-fracture criterion for fractal cracks.

Fracture Energy
The fracture energy G F of materials is defined as the area under the stress vs. the crack open displacement curve in the cohesive law, and represents the energy dissipated on the unitary crack surface [12].G F is usually determined by tests.For a large size specimen, the fracture energy * F G equals to the max driving force max f G expressed in Equation ( 16), approximately.Therefore, fracture energy G F for a normal size specimen has relation with max f G as: where A N and A * are the real areas of the fracture surfaces corresponding to the normal size and large size specimens, respectively, with a transformation as: where the constant b < 1 represents the initially cracked portion of the interface prior to load application, 2a 0 is the initial crack diameter and is assumed to be equal to the maximum aggregate size in concrete.Substituting Equations ( 16) and ( 19) into Equation ( 18), one obtains Here 2a 0 needs is considered to be the final length of the main crack for concrete, and represents a specimen size.Figure 4 shows the behaviors of G F as the changes of a 0 and D. We can find that G F increases with the increasing of these two factors.

Damage Variable
An apparent damage variable is defined as the ratio of the effective bearing area A k and the cross section area A c of materials, and has the following form: For a Euclidean shape, the area A 0 and the perimeter C 0 have the following relation: where m 0 is a shape constant.Substituting Equation (23) into Equation (22), one obtains the expression for the apparent damage variable  as follows: where C k and C c are the perimeters corresponding to A k and A c , respectively.It is reasonable for us to take C k C c as the crack length in one direction and the total length of all cracks in all directions in a unit cell at failure, respectively.By referring to Equation (1), the perimeter C of a fractal damaged surface and the counterpart C 0 in the Euclidean space has the following relation: Substituting Equation (4) into Equation ( 22), concerning Equation (25), and based on the fact that cracks in concrete is statistically self-similar fractal when the yardstick δ ranges from 0.263 to 1 [15], one obtains the following expression for a fractal damage variable in the fractal space: where D k = 1.0 ~ 2.0.From Equation (29), we find that only the fractal dimension is different between the fractal and apparent damage variables, and δ has no influence on this relation.Since A c is the cross section area of materials, which indicates D c = 1.0, thus Equation ( 26) can be rewritten as: For a smooth crack, D k = 1.0, and no fractal effect exists, and we have: From the above discussion, we notice that the apparent damage variable in the Euclidean space is a special case of the fractal damage variable in the fractal space, and the fractal damage variable is the generalization of the apparent damage variable.Figure 5 illustrates the differences between the two kind damage variables in uniaxial compression with an assumption that D k takes the average value of 1.1, and the original evolution data for apparent damage variables is referred from reference [16].

Plastic Damage Constitutive Model for Concrete
Among the various existed constitutive models for concrete, the plastic damage model has a better effect to characterize the stiffness degeneration, the strain softening and the unilateral effect of concrete under various loading conditions.

Decomposition of Effective Stress Tensor
In view of the fact that the typical failure modes of concrete are cracking in tension and crushing in compression, we decompose the effective stress tensor into tensile and compressive parts (denoted by  σ and  σ , respectively) by utilizing spectral decomposition technique [17,18]: : where P + and P -are the fourth-order projection tensors expressed as [19]: where I is the fourth-order identity tensor; represents the Heaviside function calculated for the ith eigenvalue i  of σ ; P ij is the second-order tensor and is defined as: where n i is the ith normalized eigenvector corresponding to i  .

Plasticity
We adopt a plasticity yield function f and a plastic potential function F p as: where is the algebraically maximum effective principal stress.α, β and c are parameters with the following forms [20]: where 0 b f  and 0 f  are the initial equibiaxial and uni- axial compressive yield stresses, respectively.
lies between 1.10 and 1.20 from experiments, therefore, α varies from 0.08 to 0.12.

 
c  κ represent the inner cohesion, and

 
f  κ are the evolution stresses (positive values are used here in compression) in the effective stress space due to plastic hardening or softening under uniaxial tension and compression, respectively.1 I is the first invariant of the effective stress tensor, 2 J is the second invariant of the effective deviatoric stress tensor.α p ≥ 0 is a dilation parameter with 0.2 ≤ α p ≤ 0.3 for concrete.
According to the flow rule, the rate of the effective plastic strain p   can be written as: where p   is a plastic consistency factor, and can be determined by the consistency condition for the yield surface f, which can be expressed in the Kuhn-Tucker form as: In this study, the linear isotropic hardening rules are introduced to describe the change of the yield surfaces in the effective stress space, and can be expressed in simple forms as [21]: where y f  are the effective yield strengths in uniaxial tension and compression, and have approximate values as , respectively.0 f  are the uniaxial yield strengths of concrete corresponding to tension and compression, respectively.p E  are the effective plastic hardening modulus in uniaxial case, and have relation with the elastoplastic tangent modulus ep E  as [22]: According to the plasticity consistency condition, we have: and obtain the rate form of the constitutive equation as follows: where C 0,ijkl is the fourth-order undamaged elastic stiff-ness tensor.
Accounting for the coupling of tension and compression, κ  can be written as [22]: ˆî where ˆi  are the principle stresses.
Substituting Equations ( 42) and (43) into Equation (41), we obtain: where p h  can be expressed in the following forms: where max  and min  are the maximum and minimum effective principle stresses, respectively.Therefore, we can rewrite the rate form of the relation Equation ( 42) for the effective stress and the strain as: where C ijkl is the elasto-plastic tangent stiffness tensor and has the following form: 0, 0, 0,

Free Energy
A damage constitutive model of a material is based on the second law of thermodynamics which states that all the selected internal variables must satisfy the Clausius-Duhem inequality for any irreversible process under an isothermal condition, and has a simple form as: : 0 where σ and ε are the stress and strain tensors,  is the total Helmholtz free energy (HFE) which can be considered as the sum of the elastic part e  and the plastic part p  , that is: where Φ is the damage variable tensor.
We decompose e  into tensile and compressive parts as: The incremental form of Equation ( 55) can be expressed as: Equations ( 56) and (48) form the final plastic damage constitutive equations for the plain concrete.
Similarly, we can rewrite p  as: Referring to the fact that the contribution to the plastic HFE from plastic strains of concrete in tension is much smaller comparing to the one in compression, we assume that 0 p    .Thus, we have: Substituting Equations ( 51), ( 52) and (58) into Equation (50), we get: : Since the above inequality must be satisfied for any elastic strain ε e , we have: where Y is the damage energy release rate and can be expressed as: According to the above discussions about the total HFE, we can rewrite Y as: We define the damage criteria in tension and compression for concrete, respectively, as: where r  are the current damage thresholds whose initial values is denoted by 0 r  .Equation (65) indicates that damage is initiated when Y  exceed the corresponding damage thresholds r  .Initial strain energy of materials can be written as: where 0 e   is the initial strain energy of materials in compression, and can be expressed as: where parameters Assume that concrete is in biaxial compression, which indicates σ 3 ≡0, and we have: For an uniaxial compression case, we denote the uniaxial compressive ultimate strength as 0 f  , and have . We can derive the initial damage threshold 0 r  as: And for an equibiaxial compression case, we use 0 b f  to denote the biaxial compressive ultimate strength.One gets σ 1 =σ 2 = 0 b f  , and derives 0 r  as: Noticing that the initial damage threshold is unique for a material, we obtain the expression for Ω from Equations (73) and (74): For concrete, one has , and finally have Y  ≥0.Accounting for the fact that damage is irreversible, we get   ,   ≥ 0. Therefore, the total HFE defined in Equation (51) satisfies the thermodynamic consistency, or satisfies the inequality of Equation (61).
Equation (75) reaches the limit state Ω→+∞ when , and the undamaged area of a material in compression can be characterized by the following inequality: Figure 6 illustrates the changes of the undamaged areas with α p .We find that the increase of α p will lead an expansion of the undamaged area.

Evolution Laws
The evolution laws proposed by Faria et al. [18] are adopted in this research.Damage variables (in the Euclidean space) and their rate forms are expressed as:   where A -is a material parameter, and can be determined by the uniaxial compression test; B ± can be expressed as: where F G  is the fracture energy of a material and can be determined by tests or from Equation (21).
Replacing   in the above model by   in Equation (27), we can generalize the Euclidean constitutive model for concrete to the fractal space.

Comparison of Constitutive Models
In this section, numerical simulations of the present model considering fractal effect are performed for concrete under different loading conditions, and comparisons of the results are done with some experimental data, i.e. the unaxial loading test by Karsan and Jirsa [23] and the biaxial loading one by Kupfer et al. [24].Material parameters for concrete are listed in Table 1.The values of E 0 , ν 0 , α p and α are obtained from the study of Lee and Fenves [18]; Average values of l ch,N , * ch l , η, a 0 , b and D k are used here because of their narrow range intervals.
Therefore, we can obtain the fracture energy for concrete from Equation ( 21) as: F G  =45.3N/m in tension and F G  =1497N/m in compression.

Uniaxial Tension
Both the predictions of the present model with and without the fractal effect and the test data obtained by Karsan and Jirsa [23] for concrete under uniaxial tension are illustrated in Figure 7(a).We can find that the two kind results of the present model well with the test data in the stress hardening stage.Comparing with the case of no fractal, the prediction considering the fractal is more coincident with the test.In the last large deformation stage, the three behaviors are close with each other.[23].Comparing with the uniaxial tension case, both of the two kind results of the model are more efficient to simulate the concrete compressive behaviors.The predictions considering fractal is slightly lower than the two other data.

Biaxial Tension
In this simulation, equibiaxial tensile condition is adopted, and material parameters are taken as the same as those listed in Table 1. Figure 8(a) gives the behaviors obtained from the present model and the test in biaxial tension [24].We find that the two kind theoretical results are coincident with the test data; in the initial strain softening stage, the result concerning no fractal agrees with the test better, but there is obviously a difference between them comparing with the data ing the fractal effect.All these show the superiority of the proposed model.

Biaxial Compression
In equibiaxial compression, concrete shows a good plastic deformation ability which can be found both in the proposed model and the test [24]; see Figure 8(b).The three curves are well close with each other.Comparing with the three other loading cases discussed above, the peak stress increases obviously, accompanied by the slowest decreasing softening stage in biaxial compression, which indicates a good compressive capacity of concrete under the confining pressure condition.Specially, model with fractal damage variables is more accurately than the one with apparent damage variables.

Structural Analysis
An unreinforced notched concrete beam under 3-point bending is simulated to verify the efficiency of the present concrete damaged plasticity model.This problem has been studied extensively both experimentally by Petersson [25] and analytically by Meyer et al. [26], among others.This beam is simply supported at both ends with concentrated force acting at the center.Its sketch is illustrated in Figure 9 (unit: m).
The present constitutive model considering the fractal effect of concrete is adopted for the theoretical analysis.The model parameters are taken as: F G  = 138 N/m, 0 f  = f t = 3.33 MPa, 0 f  = 30 MPa, and other properties are same as that listed in Table 1.The beam is under the plane stress condition.Accounting for the symmetry of both the structure and the load, only one half of the beam is modeled.The model is meshed with 280 4-node bilinear, reduced integration plane elements.The beam is loaded by prescribing the vertical displacement at the center of the beam until it reaches a value of 0.0015 m.The Riks method is used to solve this problem.
Figure 10 illustrates the variation curves of the concentrated force and the center displacement of the beam calculated from the theoretical analysis and the Petersson's test [25].We can note that the theoretical result is coincidence with the test at the loading stage and slightly higher then the test at the unloading branch.Figure 11 shows the distribution of the principal tensile stress of the structure as well.
Mesh sensitivity is investigated in this study by meshing the structure into coarse and fine grids, respectively, with 70 and 1120 same type elements with the medium mesh case.Resolving the above problem at the same condition, we get the relation between the load and displacement in center.Figure 12 represents the relation curves corresponding to the three kind meshes.We find that: comparing with the coarse mesh case, the structural responses agree well for the other two meshes.The structure is not sensitive to the mesh size in general.

Conclusions
In this research, the fracture toughness, the driving force and the fracture energy of a material with fractal cracks are investigated and their theoretical expresses in the fractal space are derived based on fracture mechanics and fractal geometry.The surface energy and the strain energy in the fractal fracture zone are theoretical expressed in the fractal space.The transformation rule of damage variables in the fractal space and the Euclidean space is obtained which indicates that the apparent damage variable in the Euclidean space is a special case of the fractal one in the fractal space with the fractal dimension of cracks equals to 1.We introduce a plastic yield function and decompose the damage variable tensor into tensile and compressive parts to establish a plastic damage constitutive model for concrete in the Euclidean space.Generalization of this model to the fractal space is done by utilizing the damage variable transformation rule.Comparisons of the results obtained from the present model and tests for concrete under different loading conditions are done to verify the efficiency of this model and show the necessity of considering the fractal effect in the constitutive model of concrete.The present model considering the fractal effect is used to analyze a notched plain concrete beam under 3-point bending.Mesh sensitivity is also concerned.The numerical results show the efficiency and validation of the present model for structural analysis.

Figure 2 . 0 a
Figure1illustrates an infinite plate with a fractal cut in uniaxial tension.The cut releases the stress in a fracture domain, whose shape can be approximated as an ellipse[10].A standard Koch fractal curve is employed to construct the boundary of the crack, see Figure 2. n denotes the construction step.Keep the area of frac ture zone as a constant of 2 0 a  , then the fractal dimension D of the crack is independent on the yardstick 0 3 n a 

Figure 1 .
Figure 1.A fractal crack in the infinite plate.

Figure 2 .
Figure 2. Construction of the crack boundary with a standard Koch fractal curve.

Figure 3 Figure 3 .
Figure 3 illustrates the relation among f IC K , D and a 0 .It can be noted from Figure 3 that f IC K decreases with the increasing of the crack length; If a 0 =0, f IC K tends to be infinite, and no crack exists in materials, which verifies the concept that the growth of initial deflects leads to the final failure of materials.The greater the D value, the more bifurcated the cracks are, and the larger f IC K is, which indicates that roughness has a hin- lengthes corresponding to A * and A N , respectively.Usually, * ch l is taken to be the minimum threshold of characteristic lengths.For concrete, * ch l = 0.15 mm when 0 f  = 100 MPa and * ch l = 0.15 mm when 0 f  = 200 MPa [13], and * ch l can be obtained by linear interpolation for other strengths.

Figure 4 .
Figure 4. Influences of D and a 0 on G F.

Figure 5 .
Figure 5. Evolutions of damage variables for concrete in uniaxial compression cases.
materials;   are the tensile and com- pressive components of Φ .Therefore, we derive:


is the second invariant of the compressive effective deviatoric stress tensor  s ; 1 I  and 1 I  are the first invariant of  σ and  σ  has the following form in compression case: 16 and p  = 0.2~0.3.Assuming ν 0 =0.2, and substituting the above three parameters into Equation (74), we find that Ω ≥ 0 is always hold.Thus we have 0 p   ≥0 (seeEquation (69)

Figure 6 .
Figure 6.Influences of α p on undamaged domain.

Figure 7 (Figure 7 .
Figure 7.Comparison of the constitutive curves of the present model with test results (Karsan and Jirsa 1969) in the uniaxial loading conditions.(a) Uniaxial tension; (b) Uniaxial compression.

Figure 10 .
Figure 10.Comparison of the theoretical and experimental data for the load vs. displacement curve.

Figure 11 .
Figure 11.Distribution of the principle tensile stress.

Figure 12 .
Figure 12.Influences of mesh size on structural response. ,