# Numerical simulation of ductile fiber-reinforced cement-based composite

## Numerical simulation of ductile fiber-reinforced cement-based composite

Journal of Computational and Applied Mathematics 270 (2014) 433–442 Contents lists available at ScienceDirect Journal of Computational and Applied M...

Journal of Computational and Applied Mathematics 270 (2014) 433–442

Contents lists available at ScienceDirect

Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam

Numerical simulation of ductile fiber-reinforced cement-based composite Jan Vorel a,∗ , W.P. Boshoff b a

Department of Mechanics, Faculty of Civil Engineering, Czech Technical University in Prague, Thákurova 7, 166 29 Praha 6, Czech Republic b

Department of Civil Engineering, Stellenbosch University, Private Bag X1 Matieland 7602, South Africa

article

info

abstract Strain Hardening Cement-based Composite (SHCC) is a type of High Performance Concrete (HPC) that was developed to overcome the brittleness of conventional concrete. Even though there is no significant compressive strength increase compared to conventional concrete, it exhibits superior behavior in tension. The primary objective of the presented research is to develop a constitutive model that can be used to simulate structural components with SHCC under different types of loading conditions. In particular, the constitutive model must be efficient and robust for large-scale simulations. The proposed model, based on previous research Vorel and Boshoff (2012), for plane stress is outlined and the further focus of this paper is on the mesh objectivity of the model. It is shown to be mesh size independent. © 2013 Elsevier B.V. All rights reserved.

1. Introduction At the beginning of the 21st century, civil engineers more than ever faced the often-contradictory demands for designing larger, safer and more durable structures at lower cost and shorter time. Concrete has been used for many centuries as a safe and durable building material. Two of the main advantages of concrete are its high compressive strength and that it can be cast on the construction site into a variety of shapes and sizes. The most prominent disadvantages of concrete and other cementitious materials are their brittle failure behavior in tension and the low tensile strength. The low tensile strength is usually compensated for with steel reinforcement, but wide cracks leading to the corrosion of the steel reinforcement still occur during the normal use of concrete . These cracks lead to durability problems and cause structural degradation to occur more rapidly. Fiber Reinforced cement-based composites (FRC or FRCC) is a large group of composites with variety of properties. The reason for adding fibers is to overcome the brittleness of the concrete by improving the post-cracking behavior and enhancing ductility. The paper deals with the group of Strain Hardening Cement based Composite (SHCC), a type of High Performance Concrete (HPC), that exhibits excellent mechanical behavior showing tensile strain hardening and multiple fine cracks [2,3]. It has been shown to reach a tensile strain capacity of more than 4% during a strain hardening phase [2,3] caused by the fine, closely spaced multiple cracks with crack widths normally not exceeding 100 µm [2,3]. These fine cracks, compared to large (greater than 100 µm) localized cracks found in conventional concrete, have the advantage of increased durability. For a further discussion of the mechanical properties of SHCC, the reader is referred to [4,5]. Several scholars have simulated SHCC mechanical behavior with the Finite Element Method (FEM) [6–9], Extended Finite element Method (XFEM)  or lattice model [11,12]. Kabele  formulated a model to simulate the mechanical behavior

Corresponding author. Tel.: +420 224354495. E-mail addresses: jan.[email protected] (J. Vorel), [email protected] (W.P. Boshoff).

434

J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442

Fig. 1. Coordinates and transformation angle.

J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442

435

ε = {ε11 , ε22 , γ12 = 2ε12 } T and a principal strain e = {e1 , e2 , 0} T or a global stress σ = {σ11 , σ22 , τ12 } T and principal stress s = {s1 , s2 , 0} T , respectively s = [T]σ σ.

e = [T]ε ε,

(1)

Using the standard transformation rule the tensors are

n211 [T]ε =  n221 2n11 n21

n212 n222 2n12 n22

n11 n12 , n21 n22 n11 n22 + n12 n21

cos α − sin α

 n=

sin α , cos α

(2)

with the relations between [T]ε and [T]σ given as [T]σ

[T]ε T = [T]σ −1 .

= [T]ε −1 and

T

(3)

The rotation angle α can be obtained by means of a standard relation

α=

arctan [γ12 / (ε11 − ε22 )] .

1 2

(4)

Preserving coaxiality between principal stress and strain the tangential shear modulus in the linearized incremental stress–strain law (in the crack orientation) reads  s1 − s2

G12 =

2 (e1 − e2 )

.

(5)

The derivation of whole tangential stiffness matrix in the rotating crack framework can be found in . 2.1. Poisson’s ratio effect and equivalent principal stresses The rotating crack approach does not automatically include the effect of Poisson’s ratio as the stress is evaluated on the basis of individual principal strains. In [21,8] the concept of equivalent strain is utilized to take this effect into account. Such an approach is reliable when a model formulation does not permit residual deformations by cyclic loading, i.e. by changing state (tension to compression and vice versa). However, permanent (residual) deformations are allowed in the presented  model. Therefore, a new approach is employed to treat the effect of Poisson’s ratio. The effective principal strain eˆ is used



to determine the equivalent stress sˆ from the simplified uniaxial stress–strain diagram (see Sections 2.2 and 2.3). The effective principal strain eˆ is based on the principal strain (e) which is free of inelastic deformations associated with the stress state change. The final principal stresses are consequently evaluated as

 

  s1 s2

=

1 1 − ν12 ν21

ν12 = ν0 E1 /E0 ,

1

ν21

ν12

  sˆ1 , 1 sˆ2

ν21 = ν0 E2 /E0 ,

(6) (7)

where E0 and ν0 stand for Young’s modulus and Poisson’s ratio of the undamaged material respectively. The parameters E1 , E2 , ν12 and ν21 represent the characteristics of the damaged material in a given direction and are defined in Sections 2.2 and 2.3. The isotropic elastic material is represented in the state without cracks (E1 = E2 = E0 , ν12 = ν21 = ν0 ) and the orthotropic when the crushing or cracking starts sˆ1 , sˆ2

T

  = E1 eel1 , E2 eel2 T .

(8)

The stiffness matrix introduced by this approach satisfies the condition of symmetry for orthotropic materials. Combining Eqs. (6) and (8) further gives

  s1 s2

=

1 1 − ν12 ν21

E1 ν21 E1

ν12 E2 E2

  el 

e1 , eel 2

(9)

where ν12 E2 = ν21 E1 and superscript ·el represents the elastic part. 2.2. Tension The equivalent stress state in principal direction is determined by the stress function sˆt (c ) as a function of the effective principal strain and associated history parameters. The stress function is based on the uniaxial stress–strain diagrams in tension (compression). The experimental data are idealized (approximated) to obtain a suitable mathematical representation of this constitutive model as shown in Figs. 2 and 4. In particular, the true picture is given of the strain hardening and softening in tension as well as softening in compression.

436

J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442

a

b

c

Fig. 3. Evolution of inelastic strain in tension.

The material response for virgin loading in tension (Fig. 2(a)) is defined by means of total effective strains for the ascending branches and as a cohesive crack for the softening part to ensure the proper energy dissipation and mesh independency. The elastic part is assumed to be linear, the hardening branch is the Hermite polynomial function and the softening is described by the power function. The power function is a modification of the linear tension softening diagram as defined in  for ordinary concrete. Note that this choice was made to capture the proper shape of the softening branch, see Fig. 7(a), and the area under the stress–displacement diagram is equal to the fracture energy needed when the element (specimen) is fully saturated by micro-cracks. The complete loading path in tension then reads

E eˆ 0     2  3       eˆ − εt0 eˆ − εt0     σt0 + σtp − σt0 3 εtp − εt0 − 2 εtp − εt0 sˆt eˆ ≥ εtmax =     w at    σtp 1 −   wt ,cr  0

0 ≤ eˆ ≤ εt0

εt0 < eˆ ≤ εtp (10)

εtp < eˆ ≤ εtu εtu < eˆ

where w stands for the crack opening, wt ,cr is a parameter with the dimension of length, which controls the ductility of the material, and at ≥ 1 is the material parameter. In this study the degree of polynomial softening is assumed to be the same as the degree of unloading (Eq. (11)). The ultimate tensile strain, characterized by the opening of one crack at complete

J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442

437

  = bt εtp − εt0 + wt ,cr /h, where bt is the material characteristic

failure and closing of the remaining cracks, reads εtu driving the evolution of damage and cracking strains, see Eq. (13) and Fig. 3. The remaining model parameters are depicted in Fig. 2(a)–(b). In the context of the crack-band approach, the crack opening corresponds to the inelastic (cracking) strain multiplied by the effective thickness h of the crack band. In this paper the effective thickness is estimated by projecting the finite element into the direction normal to the crack at its initiation. Note that for the given strain the evaluation of corresponding stress is no longer explicit and the Newton iteration is utilized to take into account the nonlinear unloading as well. The unloading and reloading scheme (e˙ < 0 and e˙ ≥ 0, where e˙ is the effective principal strain increment) shown in Fig. 2(c) is based on the experiments presented in  and is written as

 E0 eˆ     at   eˆ − εtul ∗  σtmax  ∗ sˆt eˆ < εtmax = εtmax − εtul   ∗   eˆ − εtul   ∗ ∗  σtul + σtmax − σtul ∗ εtmax − εtul

0 ≤ εtmax ≤ εt0

εt0 < εtmax , e˙ < 0

(11)

εt0 < εtmax , e˙ ≥ 0.

  ∗ εtmax = min εtmax , εtprl ,   ∗ εtul = max εtul , εtpul ,

(12)

where σtmax , σtul are associated stresses and εtmax is the maximum strain experienced in previous steps with stress σtmax . The evolution of inelastic strain εtul is assumed to be linearly dependent on εtmax . This simplification correlates well with cyclic tensile results done at Stellenbosch University (Fig. 3) and results presented for ordinary concrete in . ∗

0

εtul = max [previously reached εtul ; bt (εtmax − εt0 )] max. reached εtul

0 ≤ εtmax ≤ εt0

εt0 < εtmax ≤ εtu εtu < εtmax .

(13)

By considering the standard definition of the damage parameter ω Et = (1 − ωt )E0 ,

(14)

where Et denotes the actual elastic modulus, the damage variable can be determined by introducing Eq. (13) into Eq. (14) as

ωt = 1 −

Et E0

=1−

σtmax . (εtmax − εtul ) E0

(15)

The damage is assumed to be irreversible and the largest reached value in tension and compression is assumed (i.e., ω = max (ωt ; ωc )). The transverse strain ratio in Eq. (7) can then be evaluated as

νij = ν0 (1 − ω) .

(16)

This definition ensures the decreasing influence of Poisson’s ratio while the material cracks or crushes. Such behavior is very important for the SHCC which undergoes large strain in tension. 2.3. Compression In the present formulation, the so called ‘‘overlapping crack model’’ is used for the analysis of the post-peak response of SHCC in compression to ensure the mesh insensitivity. This approach was introduced in [19,24] for the analysis of reinforced concrete beams and is similar to the well-known cohesive crack model used to analyze the behavior of quasi-brittle materials in tension. The function similar to a stress–strain diagram presented in  is utilized to describe the overlapping. The prepeak stage is characterized by the initial linear part followed by the nonlinear hardening. The virgin compression loading response is shown in Fig. 4(a) and is defined mathematically as

 E0 eˆ     εcp −εc0      εcp − eˆ E0 σcp −σc0     σcp − σcp − σc0 sˆc eˆ ≤ εmin = εcp − εc0      β 1 + w/w  c c ,cr    β σcp βc − 1 + 1 + w/wc ,cr c

0 > eˆ ≥ εc0

εc0 > eˆ ≥ εcp εcp > eˆ

(17)

438

J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442

where wc ,cr and βc are the material parameters characterizing the compression crushing and therefore have to be determined experimentally or by fitting to match the prescribed post-peak energy. The unloading and reloading scheme (e˙ > 0 and e˙ ≤ 0) is depicted in Fig. 4(b) and is based on similar assumptions as for tension 0 > εcmin ≥ εc0

 E0 eˆ    a c    eˆ − εcul ∗  σcmin  ∗ sˆc eˆ > εcmin = εcmin − εcul   ∗   eˆ − εcul   ∗ ∗  + σcmin − σcul σcul ∗ εcmin − εcul

εc0 > εcmin , e˙ > 0

(18)

εc0 > εcmin , e˙ ≤ 0

where ac ≥ 1 is the material parameter controlling the unloading and

  ∗ εcmin = max εcmin , εcprl ,   ∗ εcul = min εcul , εcpul

(19)

where σcmin , σcul are associated stresses and εcmin is the minimum strain reached in previous steps with stress σcmin . The evolution of inelastic strain (crushing) is again assumed to be linearly dependent on εcmin . ∗

εcul =

0 min [previously reached εcul ; bc (εcmin − εc0 )]

0 > εcmin ≥ εc0 εc0 > εcmin

(20)

where bc is the material parameter characterizing the speed of inelastic strains evolution and for that reason have to be determined from experimental tests. The other possible nonlinear definition is, e.g., presented in , where the evolution of compressive irrecoverable strain is based on exponential function. The damage parameter is determined in a similar fashion as for tension (Eq. (15)) and reads

ωc = 1 −

Ec E0

=1−

σcmin . (εcmin − εcul ) E0

(21)

2.4. Biaxial behavior To demonstrate the complex behavior of the proposed approach the failure envelope in the space of the principal stresses is shown in Fig. 5(a) for the material properties summarized in Table 1. The biaxial interaction is influenced by the damage evolution in both directions and agrees reasonably well with the shape of the experimental results presented in [26,27]. Nevertheless, the failure criterion for better capturing of the biaxial behavior will be included at a later stage as this is currently under experimental investigation. 2.5. Cyclic loading The above described model is adjusted for cyclic loading when the orientation of principle stresses changes (Fig. 6). The residual deformations during the stress state change are assumed to be dependent on the inelastic strain. Therefore, a simple linear definition is employed and the permissible closing (opening) strain is evaluated as

εtcl(c ) = bclt(c ) εtul(cul) ,

(22)

cl where bcl t and bc are material parameters and can therefore be calculated from reverse cyclic loading tests. The trajectories of reloading after the stress state change are in good agreement with experimental results presented in [29,30,28]. The ability of the model to capture this loading is depicted in Fig. 5(b) where the model response is compared with the experimental data for cyclic uniaxial loading presented in [28, Figs. 3–39] (material model parameters are pretty similar to those presented in Table 1 and are not presented herein). To limit the paper length, we have restricted the definition to the tension behavior after stress state change only which is defined as ∗∗

 ∗∗  ∗ ∗ σ (˙e ≥ 0) = σtul + σtmax − σtul ∗ σ (˙e < 0) = σtmax

∗∗ eˆ − εtul ∗ σtmax /Et

at

,

∗ eˆ − εtul ∗∗ ∗ εtmax − εtul

−εtul Et εtmax ∗∗ ∗

σtmax −σtul

,

(23)

(24)

J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442

439

model experiment

Fig. 5. (a) Numerical model failure envelope in the normalized principal stress space compared with experimental data presented in , (b) cyclic loading—comparison of the model response with the experimental data presented in .

Fig. 7. (a) Tensile tests and model approximation, (b) scheme and dimensions of three-point bending test.

where the driving parameter eˆ is again shifted to correspond with the diagram for a virgin loading, i.e., to eliminate the ∗∗ ∗∗ inelastic deformations associated with the stress state change, and εtmax = max (εt0 , εtmax ) with associated stress σtmax , i.e., ∗∗ ∗∗ the stress experienced for given εtmax . The inelastic strain εtul is assumed to be equal to

 ∗  ∗∗ ∗ εtul = min εtmax − σtmax /Et , bt (εtmax − εt0 ) .

(25)

The stress evolution in compression after stress state change is derived in the same fashion. The corresponding relations can be obtained by substitution of tensile driving parameters for compressive variables and replacement of the maximum (max) with the minimum (min) values and vice versa in Eqs. (23)–(25). Note that during the loading after stress state change the old cracks are reopened (closed) and the tangent modulus increases to reach the value of the previously experienced modulus Et (c ) (Eqs. (15) and (21)).

440

J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442 Table 1 Material model parameters. General

Tension

Param.

Value

E

9200a 0.23a

ν

MPa

Compression

Param.

Value

σt0 εtp σtp wt ,cr

2.0a 0.043a (0.060a ) 2.4a (3.0a ) 2.9a 3.0 0.8a 0.9

at bt bcl t

MPa MPa mm

Param.

Value

|σc0 | |εcp | |σcp | |wc ,cr |

12.25 0.0042 17.5a 2.1 3.0 0.8 0.98 5.1

ac bc bcl c

βc

MPa MPa mm

(·) parameters of the surface layer with aligned fibers. a

Parameters experimentally obtained or presented in .

To demonstrate the model response, a loading change from tension to compression to tension (A–H) is shown in Fig. 6(a):

• • • • • • •

Bear in mind that if the loading follows the stress state change, the loading path has a tangent equal to the actual modulus (Eq. (23)), intervals F–K and L–G in Fig. 6(b). The unloading at this stage is defined in Eq. (24), see interval K–L in Fig. 6(b). 3. Numerical simulations As mentioned in the previous section, the constitutive model is implemented in the open source finite element code OOFEM . Isoparametric four-node quadrilateral plane-stress finite elements with four integration points (PlaneStress2d) are employed for the discretization of the numerical models presented in this paper. Moreover, the unstructured mesh based on the triangular three-node constant strain plane-stress finite elements with one integration point (trplanestress2d) is utilized to study the effects caused by the element type. Due to the lack of a reverse cyclic loading data some parameters are set up using the engineering judgment of the authors as this does not have a significant influence on the presented results. The available tensile tests for the same mixture as beams are used to set up the parameters describing tension (Table 1 and Fig. 7(a)). The compressive strength is determined from the available experiments. The rest of the parameters are set up according the general semi-empirical equation available in literature for ordinary concretes , where the elastic part is n σcp limited as σc0 = 0.7σcp . According to  the strain at the peak under uniaxial compression is estimated as εcp = n− , 1 E |σ |

cp where n = 0.8 + 17 . The presented simulations are focused on the mesh objectivity study of the proposed model intended to be used for larger structures. To examine the feasibility of the presented approach, the aforementioned numerical model is used to obtain the force–deflection diagrams of three-point bending test. To trigger the non-symmetric localization, a small irregularity is introduced into the structure (node shift equal to 1% of the beam height). Different numerical studies of beam in Fig. 7(b) are presented in this paper. First, the mesh-dependency of the model is investigated by varying the element size (1.5 mm × 1.5 mm, 5 mm × 1.5 mm and 5 × 3 mm) and element type (Fig. 8). The smaller dimension of elements is chosen in accordance with the theory introduced in  to deal with the crack spacing and therein presented mesh dependency. Second, the edge-effect caused by the aligned fibers along the bottom surface of beams is studied. The tensile strength and the strain hardening capacity in the region of the aligned fibers increase, see Table 1 and , and so the 3 mm thick layer with enhanced material properties is introduced at the bottom of the finite element (FE) models. The obtained results are compared with experimental data.

4. Results and discussions The presented crack rotating model (RCM) is utilized to obtain force–deflection diagrams of the three-point bending test for different finite element meshes (Fig. 8). The results are plotted in Fig. 9 together with the experimental data. As can be seen, the influence of the extra layer representing the aligned fibers along the surface of the beam is quite strongly pronounced and have to be taken into account for relatively thin structures. However, it can be seen from the presented results that the models with enhanced layer slightly underpredict the deformation capacity. This is attributed to the aligned fibers which probably occupy wider layer than is assumed in the study. Such characteristics still have to be confirmed by the CT scans of the specimens. Moreover, the model is set up according the tensile test where the specimen is damaged

J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442

441

Fig. 8. Three-point bending test numerical model: (a) mesh 1.5 mm × 1.5 mm, (b) mesh 5 mm × 1.5 mm, (c) mesh 5 mm × 3 mm (d) triangular mesh (1.5 mm).

Fig. 9. Loading path for three-point bending test: (a) without enhanced layer, (b) with enhanced layer.