- Email: [email protected]
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 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
Article history: Received 7 October 2013 Received in revised form 29 November 2013 Keywords: Strain Hardening Cement-based Composite (SHCC) Rotating crack model Damage Cyclic loading Nonlinear unloading
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 [1]. 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) [10] or lattice model [11,12]. Kabele [9] formulated a model to simulate the mechanical behavior
∗
Corresponding author. Tel.: +420 224354495. E-mail addresses: [email protected] (J. Vorel), [email protected] (W.P. Boshoff).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.12.021
434
J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442
Fig. 1. Coordinates and transformation angle.
of SHCC using a smeared crack approach. Despite acknowledging that a discrete cracking model would be best for the final localizing crack, Kabele decided to use a smeared cracking approach for the localization. This is due to the uncertainty of the position of the final localizing crack. Another model was proposed by [8]. This model was created to simulate the behavior of SHCC under cyclic loading to test the improvement in structural response if SHCC elements are used to dissipate energy during earth-quake loadings. Computational modeling of SHCC was also performed by [13] who used an embedded discontinuity approach for the final material softening. This method would have the same kinematic characterization as one obtained with interface elements for discrete cracking, but does not require remeshing procedures. Their conclusion was that it did not simulate the experimental results of SHCC satisfactorily due to the simplicity of the model. Boshoff [6] created a simple damage mechanics based model for the tensile behavior of SHCC. This was implemented numerically using the FEM. Even though numerous shortcomings still exist, the model showed relatively good results. Remaining issues include an unresolved mesh dependence and the under prediction of the deformation when analyzing a structure with a strain gradient. Radtke et al. [10] employed discrete fiber distributions taken into account by employing the partition of unity property of finite element shape functions without explicitly meshing them. In this approach the definition of the constitutive behavior of the matrix material, the fiber material and the fiber-matrix bond is assumed. Although this method showed satisfactory results and versatility, the simulated specimens must be limited in size as the method would be too computationally demanding for real world structures. 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 [14,15], for plane stress is outlined and the further focus of this paper is on the mesh objectivity of the model. The paper is organized as follows. The general features and definitions of the model are described in Section 2. Sections 3 and 4 deals with the finite element simulations and results of studied experiments. Finally, the concluding remarks and discussion are provided in Section 5. In the following sections the tension and compression related parameters and variables are denoted by subscripts starting with t and c, respectively. The superscript ·∗ is utilized to characterize the driving parameters for unloading and reloading when no stress state change occurs, while ·∗∗ denotes the driving parameters when stress state change takes place. 2. Numerical model definition In this section the main features of the utilized numerical model are described. To model the specific behavior of SHCC in tension, the application of classical constitutive material models used for quasi-brittle materials is not straightforward. The proposed numerical model is based on a rotating crack assumption to capture the strain hardening and softening, the multiple cracking, the crack localization and multiple orthogonal crack patterns [16]. A schematic representation of orthogonal cracking using the rotating crack model is shown using global and local axes in Fig. 1. In heterogeneous materials where micro-cracking occurs prior to the formation of a macro-crack, the rotating crack model may be more realistic than the fixed crack model. Micro-cracks are formed orthogonal to the major principal stress when the tensile strength is first violated. However, upon rotation of the principal stress axes new micro-cracks arise in the ‘‘rotated’’ direction and it is most likely that upon termination of the stress rotation, the latter micro-cracks will grow into macro-cracks. This justifies the choice of a rotating crack model from a physical perspective. A complete description of the rotating crack model can be found, e.g., in [17]. The presented model is implemented in the open source finite element code OOFEM [18] for plane stress elements using a coaxial rotating crack method (RCM) with two orthogonal cracks as described in [8]. This numerical approach is classified as the smeared crack model with the softening defined by means of the cohesive crack and overlapping crack model [19]. The rotating crack model evaluates a given strain state and generates the inelastic strain in the principal directions of the strain. Therefore, it is inevitably required to introduce a transformation tensor ([T]ε , [T]σ ) interconnecting a global strain
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 [17] s1 − s2
G12 =
2 (e1 − e2 )
.
(5)
The derivation of whole tangential stiffness matrix in the rotating crack framework can be found in [20]. 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. 2. Tensile response: (a) virgin loading—elastic and hardening part, (b) virgin loading—softening, (c) loading/unloading.
Fig. 3. Evolution of inelastic strain in tension.
Fig. 4. Compressive response: (a) virgin loading—elastic and hardening part, (b) virgin loading—softening, (c) loading/unloading.
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 [22] 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 [23] 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.
The unloading curve is based on the polynomial function of degree at and the reloading is assumed to be linear. The partial unloading and reloading is incorporated using the stress–strain values experienced during the unloading
∗ ε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 [21]. ∗
∗
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 [25] 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 [21], 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 [26], (b) cyclic loading—comparison of the model response with the experimental data presented in [28].
Fig. 6. Schematic cyclic behavior: (a) model response, (b) loading/unloading after stress state change.
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 [6].
To demonstrate the model response, a loading change from tension to compression to tension (A–H) is shown in Fig. 6(a):
• • • • • • •
A–B: initial virgin loading, Eq. (10), B–C: unloading, Eq. (11), C–D: cracks closing and compressive loading, D–E: virgin loading, Eq. (17), E–F: unloading, Eq. (18), F–G: cracks reopening and tensile loading, Eq. (23), G–H: virgin loading, Eq. (10).
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 [18]. 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 [31], where the elastic part is n σcp limited as σc0 = 0.7σcp . According to [32] 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 [6] 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 [6], 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.
in the weakest point of stretched zone, but for the three-point bending test the weakest point does not need to be found in the mid-span where the largest stresses are experienced. Therefore, it would be more accurate to use the four-point bending test where the constant moment is spread over the larger zone to avoid this effect. Nevertheless, the numerical models demonstrate good agreement with experimental data in the elastic as well as the hardening part. The discrepancy is detected for the softening part if the wall effect is not considered. The mesh objectivity is fulfilled for all models except the ones with triangular three-node constant strain finite elements as expected (even quite coarse mesh 5 × 3 mm provides sufficient results, especially if one looks at the peak load value). However, the triangular elements still provide reasonable predictions which differ only by less then 6% from the results obtained by means of the quadrilateral mesh. This is the promising feature of the model, which predestinates it for the complex simulation of larger structures. 5. Conclusions A two-dimensional numerical model for Strain Hardening Cement-based Composites was introduced and can be easily extended for three-dimensional problems. The presented model takes into account: (i) strain hardening and softening in tension as well as in compression; (ii) nonlinear unloading; (iii) nonlinear loading after stress state change including crack closing, reopening; (iv) the effect of a decreasing Poisson’s ratio with cracking or crushing. To capture the material behavior in different loading stages, 17 parameters need to be specified. However, most of them can be effortlessly determined from simple tensile and compressive tests or used as their default values, i.e., using linear laws for softening and unloading (at (c ) = 1), and full crack closing and hence full crack reopening for a stress state change (bcl t (c ) = 1), the elastic part in compression can be limited by σc0 = 0.7σcp [28,31] and according to [32] the strain at the peak under uniaxial compression σ
|σ |
cp cp n is estimated as εcp = n− , where n = 0.8 + 17 . 1 E The presented results show the mesh objectivity examined using the simulations of a three-point bending test. It was found that the model is mesh size independent and can therefore be used for large scale analysis. Further work is required to incorporate the true biaxial behavior.
Acknowledgment The financial support provided by the GAČR grant No. 105/12/P353 is gratefully acknowledged.
442
J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442
References [1] M. Otieno, M. Alexander, H.-D. Beushausen, Corrosion in cracked and uncracked concrete—influence of crack width, concrete quality and crack reopening, Mag. Concr. Res. 62 (6) (2010) 393–404. [2] V. Li, S. Wang, Tensile strain-hardening behavior of PVA-ECC, ACI Mater. J. 98 (6) (2001) 483–492. [3] W. Boshoff, G. van Zijl, Time-dependant response of ECC: characterisation of creep and rate dependence, Cement Concr. Res. 37 (2007) 725–734. [4] W. Boshoff, V. Mechtcherine, G. van Zijl, Characterising the time-dependant behaviour on the single fibre level of SHCC: part 1: mechanism of fibre pull-out creep, Cement Concr. Res. 39 (2009) 779–786. [5] W. Boshoff, V. Mechtcherine, G. van Zijl, Characterising the time-dependant behaviour on the single fibre level of SHCC: part 2: the rate effects on fibre pull-out tests, Cement Concr. Res. 39 (2009) 787–797. [6] W. Boshoff, Time-dependant behaviour of engineered cement-based composites, Ph.D. Thesis, University of Stellenbosch, 2007. [7] W. Boshoff, G. van Zijl, A computational model for strain-hardening fibre-reinforced cement-based composites, J. South Afr. Inst. Civ. Eng. 49 (2) (2007) 24–31. [8] T.-S. Han, P. Feenstra, S. Billington, Simulation of highly ductile fiber-reinforced cement-based composite components under cyclic loading, ACI Struct. J. 100 (6) (2003) 749–757. [9] P. Kabele, Assessment of structural performance of engineered cemetitious composites by computer simulation, Ph.D. Thesis, Czech Technical University in Prague, a Habilitation Thesis, 2000. [10] F. Radtke, A. Simone, L. Sluys, A partition of unity finite element method for simulating non-linear debonding and matrix failure in thin fibre composites, Internat. J. Numer. Methods Engrg. 86 (4–5) (2011) 453–476. [11] A. Spagnoli, A micromechanical lattice model to describe the fracture behaviour of engineered cementitious composites, Comput. Mater. Sci. 46 (1) (2009) 7–14. [12] Z. Qian, E. Schlangen, 3D modeling of fracture in cement-based materials, J. Multiscale Model. 1 (2) (2009) 245–261. [13] A. Simone, L. Sluys, P. Kabele, Combined continuous/discontinuous failure of cementitious composites, in: Proceedings of EURO-C 2003: Computational Modelling of Concrete Structures, 2003, pp. 133–137. [14] J. Vorel, W. Boshoff, Numerical modelling of engineered cement-based composites, in: Proceedings for Engineering Mechanics 2012, 2012, pp. 1555–1564. [15] W. Boshoff, G. van Zijl, Mesh-objectivity of crack modelling in SHCC, in: Proceedings of the International Conference Nonlocal Modelling of Material’s Failure, 2007, pp. 57–70. [16] B. Suryanto, K. Nagai, K. Maekawa, Influence of damage on cracking behavior of ductile fibre-reinforced cementitious composite, in: Proceedings of 8th International Conference on Creep, Shrinkage and Durability of Concrete and Concrete Structures, 2008, pp. 495–500. [17] J. Rots, Computational modeling of concrete fracture, Ph.D. Thesis, Delft University of Technology, 1998. [18] B. Patzák, Z. Bittnar, Design of object oriented finite element code, Adv. Eng. Softw. 32 (10–11) (2001) 759–767. [19] A. Carpinteri, M. Corrado, M. Paggi, G. Mancini, Cohesive versus overlapping crack model for a size effect analysis of RC elements in bending, in: Proceedings of 6th International FraMCoS, 2007, pp. 655–663. [20] M. Jirásek, Z. Zimmermann, Analysis of rotating crack model, J. Engrg. Math. 124 (8) (1998) 842–851. [21] W. He, Y.-F. Wu, K. Liew, A fracture energy based constitutive model for the analysis of reinforced concrete structures under cyclic loading, Comput. Methods Appl. Mech. Engrg. 197 (51–52) (2008) 4745–4762. [22] H. Reinhardt, Fracture mechanics of an elastic softening material like concrete, Heron 29 (2) (2010) 1–42. [23] V. Mechtcherine, P. Jůn, Stress–strain behaviour of strain-hardening cement-based composites (SHCC) under repeated tensile loading, in: Fracture Mechanics of Concrete Structures, 2007, pp. 1441–1448. [24] A. Carpinteri, M. Corrado, M. Paggi, An analytical model based on strain localisation for the study of size-scale and slenderness effects in uniaxial compression tests, Strain 47 (4) (2011) 351–362. [25] A. Ezeldin, Normal- and high-strength fiber-reinforced concrete under compression, J. Mater. Civ. Eng. 4 (4) (1992) 415–429. [26] K. Molapo, The behaviour of strain-hardening cement composites under biaxial compression, Master’s Thesis, Stellenbosch University, 2010. [27] S. Swaddiwudhipong, P.E.C. Seow, Modelling of steel fiber-reinforced concrete under multi-axial loads, Cement Concr. Res. 36 (7) (2006) 1354–1361. [28] K. Kesner, S. Billington, Tension, Compression and Cyclic Testing of Engineered Cementitious Composite Materials, Tech. Rep., 2004. [29] S. Billington, Damage-tolerant cement-based materials for performance-based earthquake engineering design: research needs, in: Fracture Mechanics of Concrete Structures, 2004, pp. 53–60. [30] G. Fischer, V. Li, Deformation behavior of fiber-reinforced polymer reinforced engineered cementitious composite (ECC) flexural members under reversed cyclic loading conditions, ACI Struct. J. 100 (1) (2003) 25–35. [31] A. Carpinteri, M. Corrado, M. Paggi, An integrated cohesive/overlapping crack model for the analysis of flexural cracking and crushing in RC beams, Int. J. Fract. 161 (2) (2010) 161–173. [32] R. Selby, F. Vecchio, A constitutive model for analysis of reinforced concrete solids, Can. J. Civil Eng. 24 (3) (1997) 460–470.
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
Article history: Received 7 October 2013 Received in revised form 29 November 2013 Keywords: Strain Hardening Cement-based Composite (SHCC) Rotating crack model Damage Cyclic loading Nonlinear unloading
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 [1]. 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) [10] or lattice model [11,12]. Kabele [9] formulated a model to simulate the mechanical behavior
∗
Corresponding author. Tel.: +420 224354495. E-mail addresses: [email protected] (J. Vorel), [email protected] (W.P. Boshoff).
0377-0427/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cam.2013.12.021
434
J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442
Fig. 1. Coordinates and transformation angle.
of SHCC using a smeared crack approach. Despite acknowledging that a discrete cracking model would be best for the final localizing crack, Kabele decided to use a smeared cracking approach for the localization. This is due to the uncertainty of the position of the final localizing crack. Another model was proposed by [8]. This model was created to simulate the behavior of SHCC under cyclic loading to test the improvement in structural response if SHCC elements are used to dissipate energy during earth-quake loadings. Computational modeling of SHCC was also performed by [13] who used an embedded discontinuity approach for the final material softening. This method would have the same kinematic characterization as one obtained with interface elements for discrete cracking, but does not require remeshing procedures. Their conclusion was that it did not simulate the experimental results of SHCC satisfactorily due to the simplicity of the model. Boshoff [6] created a simple damage mechanics based model for the tensile behavior of SHCC. This was implemented numerically using the FEM. Even though numerous shortcomings still exist, the model showed relatively good results. Remaining issues include an unresolved mesh dependence and the under prediction of the deformation when analyzing a structure with a strain gradient. Radtke et al. [10] employed discrete fiber distributions taken into account by employing the partition of unity property of finite element shape functions without explicitly meshing them. In this approach the definition of the constitutive behavior of the matrix material, the fiber material and the fiber-matrix bond is assumed. Although this method showed satisfactory results and versatility, the simulated specimens must be limited in size as the method would be too computationally demanding for real world structures. 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 [14,15], for plane stress is outlined and the further focus of this paper is on the mesh objectivity of the model. The paper is organized as follows. The general features and definitions of the model are described in Section 2. Sections 3 and 4 deals with the finite element simulations and results of studied experiments. Finally, the concluding remarks and discussion are provided in Section 5. In the following sections the tension and compression related parameters and variables are denoted by subscripts starting with t and c, respectively. The superscript ·∗ is utilized to characterize the driving parameters for unloading and reloading when no stress state change occurs, while ·∗∗ denotes the driving parameters when stress state change takes place. 2. Numerical model definition In this section the main features of the utilized numerical model are described. To model the specific behavior of SHCC in tension, the application of classical constitutive material models used for quasi-brittle materials is not straightforward. The proposed numerical model is based on a rotating crack assumption to capture the strain hardening and softening, the multiple cracking, the crack localization and multiple orthogonal crack patterns [16]. A schematic representation of orthogonal cracking using the rotating crack model is shown using global and local axes in Fig. 1. In heterogeneous materials where micro-cracking occurs prior to the formation of a macro-crack, the rotating crack model may be more realistic than the fixed crack model. Micro-cracks are formed orthogonal to the major principal stress when the tensile strength is first violated. However, upon rotation of the principal stress axes new micro-cracks arise in the ‘‘rotated’’ direction and it is most likely that upon termination of the stress rotation, the latter micro-cracks will grow into macro-cracks. This justifies the choice of a rotating crack model from a physical perspective. A complete description of the rotating crack model can be found, e.g., in [17]. The presented model is implemented in the open source finite element code OOFEM [18] for plane stress elements using a coaxial rotating crack method (RCM) with two orthogonal cracks as described in [8]. This numerical approach is classified as the smeared crack model with the softening defined by means of the cohesive crack and overlapping crack model [19]. The rotating crack model evaluates a given strain state and generates the inelastic strain in the principal directions of the strain. Therefore, it is inevitably required to introduce a transformation tensor ([T]ε , [T]σ ) interconnecting a global strain
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 [17] s1 − s2
G12 =
2 (e1 − e2 )
.
(5)
The derivation of whole tangential stiffness matrix in the rotating crack framework can be found in [20]. 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. 2. Tensile response: (a) virgin loading—elastic and hardening part, (b) virgin loading—softening, (c) loading/unloading.
Fig. 3. Evolution of inelastic strain in tension.
Fig. 4. Compressive response: (a) virgin loading—elastic and hardening part, (b) virgin loading—softening, (c) loading/unloading.
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 [22] 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 [23] 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.
The unloading curve is based on the polynomial function of degree at and the reloading is assumed to be linear. The partial unloading and reloading is incorporated using the stress–strain values experienced during the unloading
∗ ε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 [21]. ∗
∗
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 [25] 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 [21], 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 [26], (b) cyclic loading—comparison of the model response with the experimental data presented in [28].
Fig. 6. Schematic cyclic behavior: (a) model response, (b) loading/unloading after stress state change.
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 [6].
To demonstrate the model response, a loading change from tension to compression to tension (A–H) is shown in Fig. 6(a):
• • • • • • •
A–B: initial virgin loading, Eq. (10), B–C: unloading, Eq. (11), C–D: cracks closing and compressive loading, D–E: virgin loading, Eq. (17), E–F: unloading, Eq. (18), F–G: cracks reopening and tensile loading, Eq. (23), G–H: virgin loading, Eq. (10).
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 [18]. 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 [31], where the elastic part is n σcp limited as σc0 = 0.7σcp . According to [32] 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 [6] 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 [6], 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.
in the weakest point of stretched zone, but for the three-point bending test the weakest point does not need to be found in the mid-span where the largest stresses are experienced. Therefore, it would be more accurate to use the four-point bending test where the constant moment is spread over the larger zone to avoid this effect. Nevertheless, the numerical models demonstrate good agreement with experimental data in the elastic as well as the hardening part. The discrepancy is detected for the softening part if the wall effect is not considered. The mesh objectivity is fulfilled for all models except the ones with triangular three-node constant strain finite elements as expected (even quite coarse mesh 5 × 3 mm provides sufficient results, especially if one looks at the peak load value). However, the triangular elements still provide reasonable predictions which differ only by less then 6% from the results obtained by means of the quadrilateral mesh. This is the promising feature of the model, which predestinates it for the complex simulation of larger structures. 5. Conclusions A two-dimensional numerical model for Strain Hardening Cement-based Composites was introduced and can be easily extended for three-dimensional problems. The presented model takes into account: (i) strain hardening and softening in tension as well as in compression; (ii) nonlinear unloading; (iii) nonlinear loading after stress state change including crack closing, reopening; (iv) the effect of a decreasing Poisson’s ratio with cracking or crushing. To capture the material behavior in different loading stages, 17 parameters need to be specified. However, most of them can be effortlessly determined from simple tensile and compressive tests or used as their default values, i.e., using linear laws for softening and unloading (at (c ) = 1), and full crack closing and hence full crack reopening for a stress state change (bcl t (c ) = 1), the elastic part in compression can be limited by σc0 = 0.7σcp [28,31] and according to [32] the strain at the peak under uniaxial compression σ
|σ |
cp cp n is estimated as εcp = n− , where n = 0.8 + 17 . 1 E The presented results show the mesh objectivity examined using the simulations of a three-point bending test. It was found that the model is mesh size independent and can therefore be used for large scale analysis. Further work is required to incorporate the true biaxial behavior.
Acknowledgment The financial support provided by the GAČR grant No. 105/12/P353 is gratefully acknowledged.
442
J. Vorel, W.P. Boshoff / Journal of Computational and Applied Mathematics 270 (2014) 433–442
References [1] M. Otieno, M. Alexander, H.-D. Beushausen, Corrosion in cracked and uncracked concrete—influence of crack width, concrete quality and crack reopening, Mag. Concr. Res. 62 (6) (2010) 393–404. [2] V. Li, S. Wang, Tensile strain-hardening behavior of PVA-ECC, ACI Mater. J. 98 (6) (2001) 483–492. [3] W. Boshoff, G. van Zijl, Time-dependant response of ECC: characterisation of creep and rate dependence, Cement Concr. Res. 37 (2007) 725–734. [4] W. Boshoff, V. Mechtcherine, G. van Zijl, Characterising the time-dependant behaviour on the single fibre level of SHCC: part 1: mechanism of fibre pull-out creep, Cement Concr. Res. 39 (2009) 779–786. [5] W. Boshoff, V. Mechtcherine, G. van Zijl, Characterising the time-dependant behaviour on the single fibre level of SHCC: part 2: the rate effects on fibre pull-out tests, Cement Concr. Res. 39 (2009) 787–797. [6] W. Boshoff, Time-dependant behaviour of engineered cement-based composites, Ph.D. Thesis, University of Stellenbosch, 2007. [7] W. Boshoff, G. van Zijl, A computational model for strain-hardening fibre-reinforced cement-based composites, J. South Afr. Inst. Civ. Eng. 49 (2) (2007) 24–31. [8] T.-S. Han, P. Feenstra, S. Billington, Simulation of highly ductile fiber-reinforced cement-based composite components under cyclic loading, ACI Struct. J. 100 (6) (2003) 749–757. [9] P. Kabele, Assessment of structural performance of engineered cemetitious composites by computer simulation, Ph.D. Thesis, Czech Technical University in Prague, a Habilitation Thesis, 2000. [10] F. Radtke, A. Simone, L. Sluys, A partition of unity finite element method for simulating non-linear debonding and matrix failure in thin fibre composites, Internat. J. Numer. Methods Engrg. 86 (4–5) (2011) 453–476. [11] A. Spagnoli, A micromechanical lattice model to describe the fracture behaviour of engineered cementitious composites, Comput. Mater. Sci. 46 (1) (2009) 7–14. [12] Z. Qian, E. Schlangen, 3D modeling of fracture in cement-based materials, J. Multiscale Model. 1 (2) (2009) 245–261. [13] A. Simone, L. Sluys, P. Kabele, Combined continuous/discontinuous failure of cementitious composites, in: Proceedings of EURO-C 2003: Computational Modelling of Concrete Structures, 2003, pp. 133–137. [14] J. Vorel, W. Boshoff, Numerical modelling of engineered cement-based composites, in: Proceedings for Engineering Mechanics 2012, 2012, pp. 1555–1564. [15] W. Boshoff, G. van Zijl, Mesh-objectivity of crack modelling in SHCC, in: Proceedings of the International Conference Nonlocal Modelling of Material’s Failure, 2007, pp. 57–70. [16] B. Suryanto, K. Nagai, K. Maekawa, Influence of damage on cracking behavior of ductile fibre-reinforced cementitious composite, in: Proceedings of 8th International Conference on Creep, Shrinkage and Durability of Concrete and Concrete Structures, 2008, pp. 495–500. [17] J. Rots, Computational modeling of concrete fracture, Ph.D. Thesis, Delft University of Technology, 1998. [18] B. Patzák, Z. Bittnar, Design of object oriented finite element code, Adv. Eng. Softw. 32 (10–11) (2001) 759–767. [19] A. Carpinteri, M. Corrado, M. Paggi, G. Mancini, Cohesive versus overlapping crack model for a size effect analysis of RC elements in bending, in: Proceedings of 6th International FraMCoS, 2007, pp. 655–663. [20] M. Jirásek, Z. Zimmermann, Analysis of rotating crack model, J. Engrg. Math. 124 (8) (1998) 842–851. [21] W. He, Y.-F. Wu, K. Liew, A fracture energy based constitutive model for the analysis of reinforced concrete structures under cyclic loading, Comput. Methods Appl. Mech. Engrg. 197 (51–52) (2008) 4745–4762. [22] H. Reinhardt, Fracture mechanics of an elastic softening material like concrete, Heron 29 (2) (2010) 1–42. [23] V. Mechtcherine, P. Jůn, Stress–strain behaviour of strain-hardening cement-based composites (SHCC) under repeated tensile loading, in: Fracture Mechanics of Concrete Structures, 2007, pp. 1441–1448. [24] A. Carpinteri, M. Corrado, M. Paggi, An analytical model based on strain localisation for the study of size-scale and slenderness effects in uniaxial compression tests, Strain 47 (4) (2011) 351–362. [25] A. Ezeldin, Normal- and high-strength fiber-reinforced concrete under compression, J. Mater. Civ. Eng. 4 (4) (1992) 415–429. [26] K. Molapo, The behaviour of strain-hardening cement composites under biaxial compression, Master’s Thesis, Stellenbosch University, 2010. [27] S. Swaddiwudhipong, P.E.C. Seow, Modelling of steel fiber-reinforced concrete under multi-axial loads, Cement Concr. Res. 36 (7) (2006) 1354–1361. [28] K. Kesner, S. Billington, Tension, Compression and Cyclic Testing of Engineered Cementitious Composite Materials, Tech. Rep., 2004. [29] S. Billington, Damage-tolerant cement-based materials for performance-based earthquake engineering design: research needs, in: Fracture Mechanics of Concrete Structures, 2004, pp. 53–60. [30] G. Fischer, V. Li, Deformation behavior of fiber-reinforced polymer reinforced engineered cementitious composite (ECC) flexural members under reversed cyclic loading conditions, ACI Struct. J. 100 (1) (2003) 25–35. [31] A. Carpinteri, M. Corrado, M. Paggi, An integrated cohesive/overlapping crack model for the analysis of flexural cracking and crushing in RC beams, Int. J. Fract. 161 (2) (2010) 161–173. [32] R. Selby, F. Vecchio, A constitutive model for analysis of reinforced concrete solids, Can. J. Civil Eng. 24 (3) (1997) 460–470.