ICF13 Beijing (China) 2013 Vol. C
13th International Conference on Fracture June 16–21, 2013, Beijing, China -1- The interphase elasto-plastic damaging model Giuseppe Giambanco1,*, Giuseppe Fileccia Scimemi1, Antonino Spada1 1 Department of Civil, Environmental, Aerospace and Materials’ Engineering, University of Palermo, 90128, Italy * Corresponding author: giuseppe.giambanco@unipa.it Abstract Heterogeneous materials present a mechanical response strongly dependent on the static and kinematic phenomena occurring in the constituents and at their joints. In order to analyze this kind of materials it is a common practice to distinguish a macroscopic length scale of interest from a mesoscopic one, where the mesoscopic length scale is of the order of the typical dimensions of the constituents. At the mesoscopic level the interaction between the units is simulated by mean of apposite mechanical devices. Among these devices is popular the zero thickness interface model where contact tractions and displacement discontinuities are the primary static and kinematic variables respectively. However, in heterogeneous materials the response also depends on joint internal stresses as much as on contact stresses. The introduction of internal stresses brings to the interphase model or an enhancement of the classical zero-thickness interface. With the term ‘interphase’ we shall mean a layer separated by two physical interfaces from the bulk material or a multilayer structure with varying properties and several interfaces. Different failure conditions can be introduced for the physical interfaces and for the joint material. The interphase model has been implemented in an open-source research-oriented finite element analysis program for 2D applications. Numerical simulations are provided to show the main features of the model. Keywords Interphase element, Damage, Elastoplasticity, Finite element 1. Introduction The mechanical response of all those structures that are constituted by heterogeneous materials is dependent by different static and kinematic phenomena occurring in each constituent and at their joints. Material degradation due to nucleation, growth and coalescence of microvoids and microcracks is usually accompanied by plastic deformations as decohesion and sliding that cause strain softening and induced anisotropy. The mesoscopic approach is by now the most diffused technique to understand this kind of materials, because it overcomes the problems associated with the strong simplifications that have to be introduced when the macroscopic approach is applied. In particular, with the mesoscopic approach all the material constituents are modelled individually and their interactions are regulated by using appropriate devices able to reproduce the inelastic phenomena that usually occur at the physical interfaces. In literature, these mechanical devices are generally called contact elements, normally distinguished between link elements, thin layer elements and zero-thickness interface elements (ZTI). In the last decades interface elements have been applied in several engineering applications due to their simple formulation and to their easiness to be implemented in finite element codes [1-10]. The interface constitutive laws are expressed in terms of contact tractions and displacement discontinuities which are considered as generalized joints strains. In order to model the nonlinear behaviour caused by plastic deformations and damage evolution the constitutive laws of the interface elements are formulated making use of concepts borrowed by theory of Plasticity and Continuum Damage Mechanics. In many cases the structural response depends also on internal stresses and strains within the joint. It is sufficient to think to the fracture that appears in the middle of masonry blocks caused by the horizontal tangential contact stresses between the mortar and the block when a masonry assembly is subjected to a pure compressive load. These tangential stresses cannot be captured by the classical ZTI model. Therefore, the usual assumption used in zero-thickness interface elements, where the response is governed by contact stress components, may require a correction by introducing the
13th International Conference on Fracture June 16–21, 2013, Beijing, China -2- effect of the internal stresses into the analysis. This enhancement of the ZTI is known as interphase model, for the first time proposed by Giambanco and Mròz [11]. The interphase element has been formulated by authors as a new contact element and introduced in a scientific oriented finite element code. Patch tests have been carried out in elasticity to investigate the numerical performance and convergence of the element. All the results are shown in the paper written by Giambanco el al. [12]. In particular, in that paper is shown how strategies such as the Reduced Selective Integration or the Enhanced Assumed Strain methods are necessary to avoid shear locking effects of the element. In this work the interphase element is implemented for nonlinear applications by introducing separate limit conditions for the joint bulk material and for the physical interfaces. In particular the damage mechanics theory is used to simulate the formation and propagation of fractures in the bulk material. The elastoplastic limit condition of Mohr-Coulomb type with a tensile cut-off is adopted to describe the decohesion process of the interfaces. The contact strains are subdivided in an elastic and a plastic part. The overall model is thermodynamically consistent and the flow rules are derived by applying the Lagrangian method. With the aim to show the effectiveness of the model the interphase constitutive laws have been implemented in an open-source research-oriented finite element analysis program for 2D applications and by using the Selective Reduced Integration. The paper is organized as follows. In Section 2 the general assumptions of the model are reported and the expression of the Helmholtz free energy is furnished. In Section 3 the state equations and the flow rules are determined on the base of the thermodynamically consistent theory. Section 4 is finally dedicated to numerical applications in order to show the effectiveness of the proposed model. 2. General assumptions and Thermodynamics. Let us consider, in the Euclidean space 3ℜ referred to the orthonormal frame ( ) , , , O 1 2 3 i i i , a structure formed by two adherents +Ω , −Ω connected by a third material Ω in contact with the two bodies by means of the two physical interfaces +Σ and −Σ respectively, as in Fig. 1. Ω+ Ω− Σ+ Σ− Ωj Γ+ u Γ− u h Γ+ t Γ− t i1 i2 i3 Ω+ Ω− Σ Γ+ u Γ− u Γ+ t Γ− t i1 i2 i3 e 1 e2 e3 (a) (b) Figure 1. (a) Mechanical scheme of a third body interposed between two adherents; (b) Interphase mechanical scheme. It is assumed that the thickness h of the joint is small if compared with the characteristic dimensions of the bonded assembly. The boundary of the two adherents is divided in the two parts u ±Γ and t ±Γ , where kinematic and loading conditions are specified respectively. The joint interacts with the two adherents through the following traction components: 1 1 2 2 3 3 t t t ± ± ± ± = + + t e e e (1) which can be considered as the external surface loads for the joint.
13th International Conference on Fracture June 16–21, 2013, Beijing, China -3- In Eq. (1) 1e , 2e and 3e are the unit vectors of the local reference system, with 3e oriented along the normal to the middle surface Σ and directed towards the adherent +Ω . The joint can be regarded as an interphase model. It is assumed that the fibers inside the interphase and directed along 3e are maintained rectilinear during the deformation process. In view of this hypothesis the interphase displacement field u can be easily obtained from the displacement +u and −u of the interfaces +Σ and −Σ , thus ( ) ( ) 3 3 1 2 3 1 2 1 2 1 1 ( , , ) , , 2 2 x x x x x x x x x h h + − = + + − u u u (2) with 1x , 2x and 3x the Cartesian coordinates in the orthonormal frame ( ) 1 2 3 , , , Oe e e . Since the thickness of the joint is generally small if compared to the characteristic dimensions of the adherents, we can assume the strain state ε uniform along the 3e direction and given by: ( ) 2 1 2 1 2 3 3 2 1 ( , ) , , h s h x x x x x dx h − = ∇ ∫ ε u (3) Substituting the Eq. (2) we have: [ ] [ ] ( ) ( ) 1 2 1 1 ( , ) 2 2 s x x h + − = ⊗ + ⊗ + ∇ + ε u n n u u u (4) where [ ] + − = − u u u , n is the unit normal vector to the interphase plane and s∇ is the symmetric gradient operator defined as ( ) 1 2 s T ∇ = ∇+∇ . Let us note that in the interphase model the joint curvatures generated by displacement field (2) and the related flexural effect are neglected. Equilibrium equations are derived by applying the principle of virtual displacements (PVD) that asserts that the external work produced by the contact tractions equals the internal work developed in the joint. According to the hypothesis of a constant strain state, by applying the divergence theorem and assuming that + − Σ=Σ =Σ , the PVD leads to the following local equilibrium relations of the interphase model: ; , 2 2 h h div div on + − − ⋅ + = + ⋅ + = Σ t σ n σ 0 t σ n σ 0 (5) on ⋅ = Γ m σ 0 . (6) The basic kinematical hypotheses are the additive decomposition of total strain in the internal (i) and contact (c) parts and, for the contact strain only, a further decomposition in elastic (e) and inelastic (p) parts: i c = + ε ε ε (7) c ce cp = + ε ε ε (8) with i = ε AεA (9) being = − ⊗ A I n n the unit second order tensor. In order to comply with thermodynamic requirements, the interphase Helmholtz free energy is introduced in the following form: ( ) ( ) ( ) ( ) , , , , , , , , , , , , , i c cp i i c c cp i c i c cp d p d p ξ ξ ξ ξ Ψ =Ψ +Ψ +Ψ ε ε ε ω ε ω ε ε ε ε ε ω ; (10)
13th International Conference on Fracture June 16–21, 2013, Beijing, China -4- where iΨ and cΨ represent the free energies related to the internal and contact parts of the strain state respectively and ic Ψ is the mixed term of the free energy which takes into account the co-presence of the contact and internal strains. ξd and ξp are the damage and plastic internal variables, respectively. The principle considered for developing the constitutive laws is that damage occurs in the bulk material, therefore the damage tensor ω appears in the two terms of the total free energy that are functions of the internal strains also. In this way the constitutive model takes into account the onset of microvoids and fractures along the thickness of the joint. On the other hand, debonding of the joint from the adherents, sliding and fractures developing on surfaces parallel to the middle plane of the interphase are modelled using elastoplasticity and the inelastic contact strains cp ε are the related internal variables. In this work a single scalar damage variable ω governs the loss of stiffness of the bulk material. It ranges from 0 to 1, with the inferior and superior limits having the meaning of a pristine and a fully damaged bulk material, respectively. The explicit expression of the components of the free energy is given below: ( ) ( ) ( ) 2 1 1 tr 2 : ln 1 2 i i i i d d d h ω λ μ ξ ξ Ψ = − + − + − ε ε ε (11) ( ) ( ) ( ) 2 2 1 1 tr 2 : 2 2 c c cp c cp c cp p p h λ μ ξ Ψ = − + − − + ε ε ε ε ε ε (12) ( ) ( ) ( ) , i c cp 1 tr tr i c ω λ Ψ = − − ε ε ε (13) where λ and μ are the Lamè’s constants, hd is a material parameter which governs the softening response associated to the damage onset and growth, and hp is a material parameter specifying isotropic hardening/softening interface response. 3. State equations and flow rules. In order to derive the interphase constitutive equations, the second principle of thermodynamics, taking into account also the balance equation (first principle) can be applied in the form of the Clausius-Duhem inequality. This inequality for isothermal purely mechanical evolutive process reads as : 0 D= −Ψ≥ σ εɺ ɺ (14) where D is the interphase dissipation (for unit surface) or net entropy production. From the assumed form of the Helmholtz free energy (Eq. 10) its general rate has the following expression: , , , , : : : : i i c c i c c i c i c cp i i c c cp cp i i c i c d p d p ω ξ ξ ω ω ξ ξ ∂Ψ ∂Ψ ∂Ψ ∂Ψ ∂Ψ ∂Ψ Ψ= + + + + + + ∂ ∂ ∂ ∂ ∂ ∂ ∂Ψ ∂Ψ ∂Ψ ∂Ψ + + + + ∂ ∂ ∂ ∂ ε ε ε ε ε ε ε ε ε ɺ ɺ ɺ ɺ ɺ ɺ ɺ (15) Particularizing Eq. (14) for a purely elastic incremental deformation process ( , 0 cp d p ω ξ ξ = = = = ε 0 ɺ ɺ ɺ ɺ ), assuming the decomposition of the stress state similar to that used for the strain state
13th International Conference on Fracture June 16–21, 2013, Beijing, China -5- , being i c i = + = σ σ σ σ AσA (16) and considering the adopted expressions of the free energy parts (Eq. 11-13), the elastic stress-strain equations can be derived, thus ( ) ( ) ( ) { } 1 tr tr 2 i i c cp i σ ω λ μ = − + − + ε ε ε A ε (17) ( ) ( ) ( ) ( ) ( ) 1- tr tr 2 c i c cp c cp σ λ ω μ = + − − + − ε ε ε I A ε ε (18) p p p h χ ξ = (19) 1 d d d d h ξ χ ξ = − (20) ( ) ( ) ( ) 2 1 tr 2 : tr tr 2 i i i i c cp ζ λ μ λ = + + − ε ε ε ε ε ε (21) where χp and χd are the static variables conjugate of the internal variables ξp and ξd respectively, and ζ the thermodynamic force conjugate of the damage variable ω. Making use of the elastic strain-stress equation and of the previous positions, the final expression of the instantaneous dissipation is obtained: : 0 c cp p p d d D χ ξ χ ξ ζω = − − + ≥ σ ε ɺ ɺ ɺ ɺ . (22) In order to regulate the activation of each dissipative mechanism, two different yield functions are defined in the space of the proper static variables, namely: ( ) ( ) , 0, , 0 c p p d d χ ζ χ Φ ≤ Φ ≤ σ (23) where Φp is the classical plastic yield function specifying the elastic contact domain assumed convex and Φd is the damage activation function also assumed convex. The activation of each dissipation mechanism can be effectively described by a variational formulation which is represented by the generalized principle of maximum intrinsic dissipation: ( ) , , , max : c p d c cp p p d d χ χ ζ χ ξ χ ξ ζω − − + σ σ ε ɺ ɺ ɺ ɺ (24) subject to the constraints (Eq. 23). The Kuhn-Tucker conditions of the maximum constrained problem provide the plastic and damage evolution laws of the interphase: ; p cp d p d c λ ω λ ζ ∂Φ ∂Φ = = ∂ ∂ ε σ ɺ ɺ ɺ ɺ (25) ; p d p p d d p d ξ λ ξ λ χ χ ∂Φ ∂Φ − = − = ∂ ∂ ɺ ɺ ɺ ɺ (26) ( ) ( ) , 0, 0, , 0 c c p p p p p p χ λ λ χ Φ ≤ ≥ Φ = σ σ ɺ ɺ (27) ( ) ( ) , 0, 0, , 0 d d d d d d ζ χ λ λ ζ χ Φ ≤ ≥ Φ = ɺ ɺ (28) being pλɺ and dλɺ the plastic and damage multiplier, respectively.
13th International Conference on Fracture June 16–21, 2013, Beijing, China -6- In the present study the elasto-plastic convex domain is defined by the intersection of the classical Mohr-Coulomb bilinear function with a tension cut-off: ( ) ( ) 1 0 , tan 1 c c c p p p c χ σ φ χ Φ = + − − σ τ (29) ( ) ( ) 2 0 , 1 c c p p p χ σ σ χ Φ = − − σ (30) where c τ and cσ are the tangential stress vector and the normal stress component of the contact stresses, φ is the friction angle, c0 and σ0 the cohesion and tensile strength of the virgin interfaces. The two yield functions are depicted in Fig. 2. The following four zones can be distinguished: I elastic zone: 1 2 0, 0 p p Φ < Φ < II plastic activation in shear: 1 2 0, 0 p p Φ = Φ < III plastic activation in tension: 1 2 0, 0 p p Φ < Φ = IV plastic activation in tension and shear (corner): 1 2 0, 0 p p Φ = Φ = . The damage activation function is linear and the first activation occurs when the thermodynamic force reaches the relative threshold value ζ0: ( ) ( ) 0 , 1 d d d ζ χ ζ ζ χ Φ = − − (31) Figure 2. Plastic yield conditions represented in the stress space. 4. Numerical applications. The interphase model presented in Sections 2 and 3 has been implemented in an open-source research-oriented finite element analysis program for 2D applications. With the aim to run a step-by-step integration, flow rules given in rate form were rewritten as discrete laws. The implicit backward-Euler difference method was applied to obtain results within the time step
13th International Conference on Fracture June 16–21, 2013, Beijing, China -7- [ ] [ ] 1 , 0, n n t t T + ⊂ . In particular the nonlinear solution at time 1 nt + has been calculated by means of an elastic prediction – plastic and/or damaging correction procedure. The interphase element has four nodes and zero-thickness. The integration of the stiffness matrix has been solved by applying the Reduced Selective Integration method, that is two Gauss are used for the integration in the direction normal to the interphase plane while one Gauss point is used in the tangential one. The numerical applications presented in this work regard uniaxial compression and diagonal compression tests on masonry specimens. All numerical examples have been carried out under the hypothesis of plane stress state. 4.1. Uniaxial compression tests on masonry. Uniaxial compression tests have been carried out to assess the performance of the interphase element. A brick-mortar-brick system uniformly compressed (Fig. 3) has been considered with the aim to show the formation of damage in the bulk material when the stiffness of the bricks is lower than that one of mortar. In this case, in fact, the bricks tend to expand more than the mortar in the horizontal direction because of the Poisson effect, but the higher stiffness of the mortar opposes to this displacement. The bricks are therefore subjected to compressive stresses while the mortar to tensile ones along the x-axis. Figure 3. Uniaxial compression test on a masonry block. The parameters used for tests are reported in Table 1, while the results are shown in Fig. 4. The numerical test has been performed under displacement control. In this paper the results obtained with a load multiplier equal to 3 and for a FE model with 80 interphase elements are reported in terms of stresses along the mortar layer. Table 1. Parameters used for compression test. Ebrick 500 MPa vbrick 0.3 Emortar 15000 MPa vmortar 0.2 c0 4.5 MPa φ 30° σ0 1.0 MPa hp 1 ζ0 0.001 MPa hd 0.0025 MPa
13th International Conference on Fracture June 16–21, 2013, Beijing, China -8- The diagrams depicted in Fig. 4 clearly show the points were the damage nucleation takes place. The first fracture forms in the middle of the layer and grows during the following steps. As a consequence, a redistribution of stresses takes place in the two parts on sides of the fracture. When the thermodynamic force in another Gauss point or group of Gauss points overcomes the threshold value a new fracture appears symmetrically with respect to the middle of the mortar joint. On the other hand, while σx stresses suddenly fall down and tend to zero, the σz values are little influenced by fracture. This is due to the fact that damage affects only internal and mixed terms of the free energy. 0 4 8 12 16 20 d [cm] 0 10 20 30 40 50 60 σx [MPa] x 10-1 0 4 8 12 16 20 d [cm] -20 -40 -60 -80 -100 -120 σz [MPa] 0 4 8 12 16 20 d [cm] -0.6 -0.4 -0.2 0 0.2 0.4 0.6 τ [MPa] x 10-1 x 10-1 Figure 4. Uniaxial compression test on a masonry block: results at a load multiplier equal to 3. (a) normal stresses along x-axis; (b) vertical stresses; (c) tangential stresses. 4.2. Diagonal compression tests on masonry. Diagonal compression numerical tests have been carried out on a masonry panel and compared with the experimental results obtained in laboratory. With reference to Fig. 5, the specimen is made up of four courses of sandstone blocks with calcium-cement mortar. It has a squared shape with a length of 67 cm for each side. A single block is 33 cm long and 16 cm high. The mortar layer has a thickness of 1 cm. A total number of 256 plane stress 2D solid elements and 72 interphase elements has been used to create the finite element model. hinge hydraulic jack Load cell Figure 5. Diagonal compression test on a masonry panel. Experimental setup and finite element model. The mechanical variables for blocks and mortar and the parameters used for the finite element model are reported in Table 2. (a) (b) (c)
13th International Conference on Fracture June 16–21, 2013, Beijing, China -9- Table 2. Parameters used for diagonal compression test. Ebrick 11141.67 MPa vbrick 0.25 Emortar 37000 MPa vmortar 0.14 c0 5.0 MPa φ 35° σ0 1.0 MPa hp 0.04 ζ0 1·10 -5 MPa h d 0.001 MPa The test has been run with 400 steps under displacement control. At each step the values of the total vertical load F and the displacement discontinuities h h h δ δ δ + − = + and v v v δ δ δ + − = + , in horizontal and vertical direction respectively, have been evaluated and reported in the load-displacement curves of Fig. 6. A good agreement has been obtained with experimental results. -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5 displacement [mm] 0 4000 8000 12000 16000 Load [daN] exper V1 exper V2 numeric δv ++δ v - exper H1 exper H2 numeric δh -+δ h + Figure 6. Load-displacement curves: comparison between numerical and experimental results. Fig. 7 finally shows a comparison between the collapse experimental and numerical configurations. A very good agreement can be observed, even if the numerical model is not able to simulate the fracture inside the blocks because an elastic behavior was chosen for bricks. Figure 7. Comparison between numerical and experimental collapse configurations. 5. Conclusions. The present paper deals with the mesomodelling of heterogeneous structures by means of interphase elements, that can be considered as an enhancement of the common interface elements. The possibility to distinguish internal and external contact stresses and strains permits to introduce separate failure conditions for the bulk material and for contact tractions. In particular an isotropic
13th International Conference on Fracture June 16–21, 2013, Beijing, China -10- damage model has been used to model the nonlinear response of the bulk material, while an elastoplastic bilinear domain governs the evolution of plasticity for contact tractions. The interphase model has been written in the framework of a thermodynamically consistent theory. State equations and flow rules have been derived and rewritten in a discrete form to be suitable to be used for finite element implementation. Two numerical applications on masonry structural elements have been conducted. In particular a compression test and a diagonal test have been carried out. Ongoing and future efforts are devoted to the introduction of plastic activation functions on the physical interfaces between mortar and block, and to the possibility to introduce a new damage model to catch horizontal fractures also. Acknowledgements The authors acknowledge the financial support given by the Italian Ministry of Education, University and Research (MIUR) under the PRIN09 project 2009XWLFKW_005 , “A multiscale approach for the analysis of decohesion processes and fracture propagation”. References [1] G. Giambanco, S. Rizzo, R. Spallino, Numerical analysis of masonry structures via interface models. Comput. Methods Appl. Mech. Engrg., 190 (2001) 6493-6511. [2] G. Alfano, E. Sacco, Combining interface damage and friction in a cohesive-zone model. Int. J. Numeric. Methods Engrg., 68 (2006) 542-582. [3] I. Einav, G. T. Houlsby, G. D. Nguyen, Coupled damage and plasticity models derived from energy and dissipation potentials. Int. J. Solids Struct., 44 (2007) 2487-2508. [4] M. R. Salari, S. Saeb, K. J. William, S. J. Patchet, R. C. Carrasco, A coupled elastoplastic damage model for geomaterials. Comput. Methods Appl. Mech. Engrg., 193 (2004) 2625-2643. [5] P. B. Lourenco, J. Rots, Multisurface interface model for analysis of masonry structures. Journal Engrg. Mech., 123 (7) (1997) 660-668. [6] A. Spada, G. Giambanco, P. Rizzo, Damage and plasticity at the interfaces in composite materials and structures. Comput. Methods Appl. Mech. Engrg., 198 (2009) 3884-3901. [7] G. Alfano, M. A. Crisfield, Finite element interface models for the delamination analysis of laminated composite: mechanical and computational issues. Int. J. Numer. Meth. Engrg., 50 (2001) 1701-1736. [8] K. William, I. Rhee, B. Shing, Interface damage model for thermomechanical degradation of heterogeneous materials. Comput. Methods Appl. Mech. Engrg., 193 (2004) 3327-3350. [9] H. R. Lofti, P. B. Shing, Interface model applied to fracture of masonry structures. J. Struct. Eng. (ASCE), 120 (1) (1994) 63-80. [10] F. Parrinello, B. Failla, G. Borino, Cohesive-frictional interface costitutive model. Int. J. Solids and Structures, 46 (2009) 2680-2692. [11] G. Giambanco, Z. Mroz, The interphase model for the analysis of joints in rock masses and masonry structures. Mechanica, 36 (1) (2011) 111-130. [12] G. Giambanco, G. Fileccia Scimemi, A. Spada, The interphase element. Comp. Mech., 50 (3) (2012) 353-366. [13] E. Sacco, F. Lebon, A damage-friction interface model derived from micromechanical approach. Int. J. Solids and Structures, 49 (26) (2012) 3666-3680.
13th International Conference on Fracture June 16–21, 2013, Beijing, China -1- The damage evolution of 2D-Woven-C/SiC composite under tension loading Min-ge Duan1,2, Fei Xu1,*, Zhong-bin Tang1 1 School of Aeronautics, Northwestern Polytechnical University, Xi’an, 710072, China 2 State Key Laboratory of Explosion Science and Technology, Beijing Institute of Technology, Beijing, 100081, China * Corresponding author: fay.xu@nwpu.edu.cn Abstract A meso-numerical model is set up to investigate the damage evolution of 2D woven C/SiC based on RVE(representative volume element) because the woven composite has a periodic in-plane structure. Periodic boundary conditions are imposed to the model. Three material models for the fibers, the matrix, and the yarn/matrix interface are considered. Basic damage modes and their evolutions are observed in the simulation. Firstly, the matrix is damaged at the global applied stress around 50Mpa, which agrees well with the tests. The stress/strain relation obtained from the meso-numerical model can represent the behavior of the 2D woven C/SiC under tension loading. Secondly, the void in the matrix plays an important role. The bundle/matrix interface is essential to investigate the interface debonding. Keywords 2D-C/SiC composites, RVE, damage mechanisms 1. Introduction C/SiC, as one of the most promising high temperature structural materials for its performance, such as low coefficient of thermal expansion, excellent thermal shock resistance, abrasion resistance, and high specific strength/modulus, has been applied in structures such as air-breathing engines, hot-gas valves and aerospace thermal structures. Consequently, the mechanical properties of the C/SiC has been an active research field for several decades. Both experimental studies and macro- or micromechanical methods have been used to obtain the mechanical properties of the C/SiC[1-13,20]. In experimental studies, the mechanical characterization of the C/SiC were studied by classical strain measurements, ultrasonic method, AE(acoustic emission) technique , infrared thermography and microstructural observations under the load of uniaxial tension[1-6], compression[1,5], in-plane shear[7], tension fatigue[8,9], thermal fatigue[10,11] and low-velocity impact[12]. Based on the expensive experiments, the damage modes found by microstructural observations are: matrix cracking, transverse bundle cracking, bundle/matrix and inter-bundle debonding, fiber cracking, ply delamination, bundle splitting and matrix wear. The damage evolution or development mention in the former studies are all deduced from the possible damage modes, which can only reflect the overall damage event and the detailed mechanical behavior of the constituents identification is beyond the detection. Fortunately, it is possible to identify detailed mechanical behavior of the constituents in FEM. However, there is not so much study on C/SiC with simulation method as experiment method. The prediction of the stiffness tensor has been carried out with FEM[13]. While, little work was done on the prediction of damage evolution and the strength[4]. Recently, simulation of plain woven fabrics based on photomicrograph measurement and idealized sinusoidal representation of the weave structure has been successfully applied to determine their mechanical properties by taking the advantages of their inherent periodicity of woven fabric architecture[14-19] at mesoscale in-plane. Both elastic moduli [14,15], damage[16] and fracture[18] behavior of the woven fabrics have been successfully predicted for the case of two-phase heterogeneous materials without considering the fiber-matrix interface using FEM by RVE (representative volume elements). While, it has been confirmed that interface obviously has great
13th International Conference on Fracture June 16–21, 2013, Beijing, China -2- influence on the macroscopic behavior, the toughness and strength of CMCs[20]. Interfaces are very narrow regions between matrix and fibers and are responsible for a variety of key properties, because of its significance role of stress transfer between fibers and matrix[20-23]. Recently, CZM(cohesive zone model) is being increasingly used in describing interfacial debonding and separation[23]. It was first proposed to analyze brittle fracture, and developed to model debonding of inclusions, crack growth processes in elastic-plastic solids, delamination process in composites subjected to low velocity impacts and so on. It is not so easy to simulate the bundle/matrix interface, because the mechanical properties of the interface is very difficult to obtain and the geometry of the interface in woven fabrics is curved surface. The aim of the present contribution is to introduce the interface into the RVE model in terms of CZM, and to analyze the damage evolution in 2D woven C/SiC with the damage behavior of different constituents identified clearly. This paper is organized as follows: firstly, the geometry model and boundary conditions are presented. Next, the meso-model was outlined with the description of material models of matrix, fiber yarns and interface-matrix interface. Then, the damage of different constituents was predicted based on elastic model. In the last but one section, the damage initiation and evolution are described with different damage modes identified. Finally, conclusions are summarized in the final sections. 2. RVE and periodic boundary conditions 2D woven C/SiC is plain weave fabrics, which are formed by weaving of yarns, as shown in Fig. 1. The black yarn are in X(longitudinal) direction and the grey yarn are in Z(transverse) direction. Figure 1. Sketch of 2D plain weave single lamina(left) and its RVE without matrix(right) In order to modeling the fabric-reinforced laminates using FEM, the representative volume elements (RVE) of the respective configurations are chosen. The RVE is the repeated element that can represent the whole composite fabric structure by a periodical array. Then characteristics of the macrostructure composites can be represented by the RVE with a certain proper boundary conditions. And the data points are fit to a sinusoidal function of the form[15]: where, a, the amplitude of the yarn curve is 0.125mm; b=2π/1.8≈3.49, and 1.8 is the pitch of the yarn path curve in millimeters; d, the offset depicted is ±0.125mm; c is the phase adjusting factor. Figure 2. The photomicrograph of the 2D plain weave C/SiC
13th International Conference on Fracture June 16–21, 2013, Beijing, China -3- Fig. 2. is the photomicrograph of the 2D plain weave C/SiC by SEM. According to the micrograph the geometry of the fiber yarns is determined, AB=0.25mm, AE=AD=1.8mm, with the woven hole in the square of 0.2mm×0.2mm. The periodic boundary conditions meets the requirement of displacement periodicity and continuity of the proper boundary conditions, which implies that no separation or overlap will be found between neighboring of RVEs. Each RVE in the composites has the same deformation mode. As stated by Z.Xia[18], the periodic conditions on the boundary V is where, ik are the average strains, i u is the periodic part of the displacement components on the boundary surfaces and it is related to the applied global loads. Aboudi has developed a unified micromechanical theory based on the study of interacting periodic cells, and its boundary conditions, plane-remains-plane, were applied to the RVE models in the normal traction loading conditions. Now based on these the proper boundary conditions are the planes ABFE and CDGH keeps plane during the loading in x direction on EFHG plane. The ABCD plane keeps zero displacement in X direction. And to avoid the rigid motion, the displacement components of the center point of the ABCD plane are assumed to be zero. 3. The meso-model In the model, RVE that used to investigate the 2D woven C/SiC under tension loading, is shown in Fig. 3. The matrix, fiber and the yarn/matrix interface are taken into account with their specific damage behavior. (a) Matrix (b) Interface (c) Fiber bundles Figure 3. The RVE model and its constituents 3.1. Matrix The elastic constants for the matrix material SiC: Young's modulus E, Poisson's ratio v are found to be E=430Gpa, v=0.3. The brittle cracking material constitutive in ABAQUS is chosen to describe the damage and failure behavior of the isotropic brittle matrix. The failure stress and displacement are 200Mpa and 0.02mm, respectively. In order to consider void in matrix, the elastic constants is assumed to reduced to 10% of the original matrix, and the model with void is shown in Fig. 3(a), the small squares in the RVE is the void region. 3.2. Fiber yarns The fiber bundles in C/SiC is regarded as transverse isotropic materials. And the properties are listed in Table 1. Table 1. The mechanical properties of the fiber yarns E1(Gpa) E2(Gpa) v12 v23 G12(Gpa) G23(Gpa) f(Mpa) f 220 13.8 0.2 0.18 9 4.8 800 1.5%
13th International Conference on Fracture June 16–21, 2013, Beijing, China -4- Figure 4. The flow chart for the VUMAT To study the damage process, the damage model is introduced to the fiber yarns with VUMAT by Fortran, and the flow chart is shown in Fig. 4. Here the maximum stress criterion and linear stiffness degradation are chosen as the damage initiation criterion and damage evolution law of the fiber yarns. 3.3. Yarn/matrix interface The interface zone plays a key role in the mechanical behavior of composite materials. Its properties can be influenced by different fiber coatings such as PyC, BN or SiC and its thickness. And the most important mechanism for improving the toughness of CMCs is the crack deflection along the interface after the initiation of the matrix cracking. In this meso-level FEM model, the inner-yarn interface is not considered, and fiber bundles is simplified as transversely isotropy. In this work, bilinear models both for normal and tangential separation is used to simulate the response of the yarn/matrix interfacial zone. The curve of normal and tangential tractions with respect to n andt are shown in Fig.5, where max and max are the interface normal and tangential strength, are 200Mpa and 150Mpa respectively;max is the interface characteristic length parameter; n andt denote the non-dimensional normal and tangential displacement respectively. Figure 5. The curve of normal(left) and tangential(right) tractions with respect to n and t In the FE model, CZM was realized by either cohesive element or cohesive behavior and maximum stress criterion was used to predict the damage initiation in Abaqus/explicit, and detailed constitutive model for the CZM can be found in the help documentation[24]. 4. The result analysis
13th International Conference on Fracture June 16–21, 2013, Beijing, China -5- 4.1. The mesh sensibility and validation Firstly, in order to get accurate results, the mesh sensibility was studied with mesh size 10um, 30um and 50um in the quarter model. The result of the overall force-displacement response and the mises contour are shown in Fig.6. The overall force-displacement response shows very small difference. The mises contour pictures corresponding to the loading of 1.12×10-3mm for 10um, 30um and 50um are listed from left to right(white means 0Mpa, and black grey means 200Mpa, and black means exceeding 200Mpa), the maximum stress in matrix is between 203.3Mpa and 208.7Mpa, while it is 250.8Mpa and 273.1Mpa for the yarns, which shows little difference. So, 50um was accepted for its time-efficiency. (a) Overall force-displacement (b) The mises contour Figure 6. The effect of mesh size The prediction of the elastic modulus is about 75Gpa. And the elastic modulus from experiments is almost between 70Gpa[10,11] and 81.2Gpa[3]. Furthermore, the damage in RVE initiates at the stress about 53Mpa in the matrix, which agrees well with the initial proportional limit of the composites[3], around 50Mpa. 4.2. The damage prediction White-Black 0-200Mpa 0-150Mpa 0-800Mpa Figure 7. The mises contour of matrix(5th step)/interface/fibers(37th step) From the elastic RVE model, the matrix is probably damaged first at the thinnest region that close to the yarns, the damage of the yarns will initiate at the woven point, and the interface is easy to get damaged at the region of woven edge, shown in black in Fig. 7. When the stress in matrix reaches the strength, the average stress in the RVE is about 53Mpa. When the stress in the yarns reaches the strength, the average stress in the RVE is about 330Mpa,
13th International Conference on Fracture June 16–21, 2013, Beijing, China -6- corresponding to strain about 0.005 in Fig. 8 of the grey curve, higher than the failure stress of the composites, which indicates the limitation of the elastic model. So the damage models are introduced to the constituents of the C/SiC, and the corresponding stress-strain curve compared with no damage model is shown in Fig .8. The detailed damage process of the three constituents are shown in subsequent sections. Figure 8. The stress-strain curve 4.3. The matrix damage evolution The brittle cracking was introduced to matrix for its brittleness in fracture process. Fig. 9 shows the damage initiation and propagation process of the matrix in terms of STATUS, which represents the state of the element (1.0 means the element is active, shown in white; and 0.0 means the element is inactive, shown in black). It is found that the damage first occurs at the thinnest region close to the yarns in loading directions as predicted in the elastic model. In the upper(Y+) matrix, the damage initiates at the edge center of the RVE, while in the lower(Y-) matrix, it initiates at the center of the RVE. Then, the damage propagates in Z direction, vertical to loading. And the damage first runs through the Z direction in the lower matrix at the strain of 2200, while it happens to the upper matrix at the strain of 2600. The damage has already run through the cross section of the RVE at the strain of 2600. 800 900 1100 2200 2600 3500 10000 Figure 9. The damage process of the matrix 4.4. The fiber damage evolution The simulation of fiber damage was introduced by VUMAT in Abaqus, which aims at investigating the possible damage process of fiber bundle damage. Fig. 10 shows the damage initiation and propagation process of the matrix in terms of STATUS. It is found that the damage first occurs at the weaving center of the yarns in loading directions as predicted in the elastic model. Then, the damage propagates in Z direction, vertical to loading. And the damage first runs through the Z direction in longitudinal yarns at the strain of 5800, while it happens to the center yarns in
13th International Conference on Fracture June 16–21, 2013, Beijing, China -7- longitudinal at the strain of 6000. 5700 5800 5900 6000 6500 10000 y+ y- Figure 10. The damage process of the fiber bundle 4.4. The yarn/matrix interface damage evolution The CZM was introduced to the RVE by cohesive element in Abaqus, which aims at investigating the possible damage process of yarn/matrix interface. 500 1000 2000 4000 6000 8000 10000 Figure 11. The damage process of the interface Fig. 11 shows the damage initiation and propagation process of the interface in terms of MAXSCRT(0.0 means no damage occurs, shown in white grey; and 1.0 corresponding to failure, shown in black grey). It is found that the damage initiates at the weaving edge of the interface around the yarns in longitudinal direction as shown in the MAXSCRT contour of 6000, then it spread to the weaving point of the RVE in longitudinal(X) direction, marked with circles and arrows. The interface around yarns in transverse(Z) direction get damaged at the strain of 10000, later than it in longitudinal direction. 4.6. The effect of the void in matrix The porosity in C/SiC may have great influence on the damage initiation, evolution. So, the RVE with void in matrix is investigated. From the photomicrograph of the 2D plain weave C/SiC, it can be found that the void in matrix is mainly around the weaving gap between yarns, so the mechanical properties of the matrix in this region, showing as four small rectangles in the matrix (the volume is about 6% of the whole volume of the RVE) is reduced to 10%. The result of the matrix damage evolution is shown below in Fig. 12. The damage of the matrix initiates around the void, which differs from the damage evolution in matrix without void compared with Fig. 8. And it is also found that the through damage in the model with void comes earlier than that without void. So the void in
13th International Conference on Fracture June 16–21, 2013, Beijing, China -8- matrix plays an important role in the damage evolution. 1200 1400 2000 2500 5000 5900 10000 Figure 12. The damage process of the matrix with void 5. Conclusions This paper focuses on the 2D plain woven C/SiC, and the conclusions are as follows: 1) The geometry model of RVE is established based on the photomicrograph of the 2D plain weave C/SiC by SEM with the yarn/matrix interface successfully introduced. And the damage models for the matrix, yarn/matrix interface and yarns are considered. 2) The established meso-level FEM is able to describe mechanical behavior of different constituents in the plain-woven C/SiC. Especially, the damage initiates in the matrix at the overall stress about 50Mpa, which agrees well with the experiments. 3) The damage of the matrix initiates near the void and spreads in transverse directions, and spreads through the section of the RVE; the damage of the fiber begins at the very region of the weaving point in longitudinal yarns and spreads in transverse direction; the damage of the yarn/matrix occurs near the weaving edge and spreads both in transverse and longitudinal directions. 4) The void in the matrix has great influence on the damage initiation and evolution of the matrix in the RVE model. It is found that the damage of the matrix without void initiates near the yarn in Z direction at the thin region of the matrix, while the damage of the matrix with void initiates around the void. Acknowledgements This work is supported by the Open Foundation provided by the State Key Laboratory of Explosion Science and Technology (KFJJ12-14M), the Natural Science Foundation of China (90916027), and the Natural Science Fund for Young Scholars of China (11102168). References [1] C. Gerald, G. Laurent, B. Stephane, Development of damage in a 2D woven C/SiC composite under mechanical loading: I. mechanical characterization. Composites Science and Technology, 56 (1996) 1363–1372. [2] E.B. Rachid, B. Stephane, C. Gerald, Development of damage in a 2D woven C/SiC composite under mechanical loading: II. ultrasonic characterization. Composites Science and Technology, 56 (1996) 1373–1382. [3] M. Hui, L.F. Cheng, L.T. Zhang, Y.D. Xu, Z.X. Meng, C.D. Liu, Damage evolution and microstructure characterization of a cross-woven C/SiC composite under tensile loading. Journal of the Chinese ceramic society, 35 (2007) 137–143. (in Chinese) [4] C.P. Yang, G.Q. Jiao, B. Wang, Tensile damage behavior of 2D-C-SiC composite. Journal of Aeronautical Materials, 30 (2010) 87–92. (in Chinese) [5] G.Y. Guan, G.Q. Jiao, Z.G. Zhang, Uniaxial macro-mechanical property and failure mode of a 2D-woven C/SiC composite. Acta Material Composite Sinica, 22 (2005) 81–86. (in Chinese) [6] M.D. Wang, C. Laird, Characterization of microstructure and tensile behavior of a cross-woven C-SiC composite. Acta mater, 44 (1996) 1371–1387.
13th International Conference on Fracture June 16–21, 2013, Beijing, China -9- [7] G.Y. Guan, G.Q. Jiao, Z.G. Zhang, In-plane shear fracture characteristics of plain-woven C/SiC composite. Mechanical Science and Technology, 24 (2005) 515–521.(in Chinese) [8] M.D. Wang, C. Laird, Tension-tension fatigue of a cross-woven C/SiC composite. Material Science and Engineering, A230(1997) 171–182. [9] C.D. Liu, L.F. Cheng, X.G. Luan, B. Li, J. Zhou, Damage evolution and real-time non-destructive evaluation of 2D carbon-fiber/SiC-matrix composites under fatigue loading. Material Letters, 62 (2008) 3922–3924. [10] M. Hui, L.F. Cheng, L.T. Zhang, X.C. Luan, J. Zhang, Behavior of two-dimensional C/SiC composites subjected to thermal cycling in controlled environments. Carbon, 44 (2006) 121–127. [11] M. Hui, L.F. Cheng, Damage analysis of 2D C/SiC composites subjected to thermal cycling in oxidizing environments by mechanical and electrical characterization. Material Letters, 59 (2005) 3246–3251. [12] L.J. Yao, Z.S. Li, Y.Q. Cheng, X.Y. Tong, Damage behavior of 2D C/SiC composites under low velocity impact. Journal of Inorganic Materials, 25(2010)311-314. (in Chinese) [13] Y.J. Chang, Study on damage constitutive model and mechanical properties of C/SiC composite. PhD thesis, Northwestern Polytechnical University, Xi'an, 2008. (in Chinese) [14] E.J. Barbero, J. Trovillion, J.A. Mayugo, K.K. Sikkil, Finite element modeling of plain weave fabrics from photomicrograph measurement. Composites Structures, 73 (2006)41–52. [15] E.J. Barbero, T.M. Damiani, J. Trovillion, Micromechanics of fabric reinforced composites with periodic microstructure. International Journal of Solids and Structures, 42 (2005) 2487–2504. [16] E.J. Barbero, P. Lonetti, K.K. Sikkil, Finite element continuum damage modeling of plain weave reinforced composites. Composites: Part B, 37 (2006)137–147. [17] H. Ismar, F. Schroter, F. Streicher, Modeling and numerical simulation of the mechanical behavior of woven SiC/SiC regarding a three-dimensional unit cell. Computational Materials Science, 19 (2000) 320–328. [18] S.I. Haan, P.G. Charalambides, M. Suri, A specialized finite element for the study of woven composites. Computational Mechanics, 27 (2001) 445–462. [19] Z. Xia, Y. Zhang, F. Ellyin, A unified periodic boundary conditions for representative volume elements of composites and applications. International Journal of Solids and Structures, 40 (2003) 1907–1921. [20] J.F. Despres, M. Monthioux, Mechanical properties of C/SiC composites as explained from their interfacial features. Journal of the European Ceramic Society, 15 (1995) 209–224. [21] R.R. Naslain, The design of the fiber-matrix interfacial zone in ceramic matrix composites. Composites: Part A, 29A (1998) 1145–1155. [22] P.C. Yang, G.Q. Jiao, Effects of interface on tensile properties of fiber reinforced ceramic matrix composite. Acta Material Composite Sinica, 27 (2010) 116–121. (in Chinese) [23] N. Chandra, H. Li, C. Shet, H. Ghonem, Some issues in the application of cohesive zone models for metal-ceramic interfaces. International Journal of Solids and Structures, 39 (2002) 2827–2855. [24] ABAQUS documentation.
13th International Conference on Fracture June 16–21, 2013, Beijing, China -1- Structure of micro-crack population and damage evolution in concrete Andrey P Jivkov1,* 1 School of Mechanical, Aerospace and Civil Engineering, The University of Manchester, Manchester M13 9PL, UK * Corresponding author: andrey.jivkov@manchester.ac.uk Abstract Tensile behaviour of concrete is controlled by the generation and growth of micro-cracks. A 3D lattice model is used in this work for generating micro-crack populations. In the model, lattice sites signify solid-phase grains and lattice bonds transmit forces and moments between adjacent sites. The meso-scale features generating micro-cracks are pores located at the interfaces between solid-phase grains. In the model these are allocated to the lattice bonds with sizes dictated by an experimentally determined pore size distribution. Micro-cracks are generated by removal of bonds when a criterion based on local forces and pore size is met. The growing population of micro-cracks results in a non-linear stress-strain response, which can be characterised by a standard damage parameter. This population is analysed using a graph-theoretical approach, where graph nodes represent failed bonds and graph edges connect neighbouring failed bonds, i.e. coalesced micro-cracks. The evolving structure of the graph components is presented and linked to the emergent non-linear behaviour and damage. The results provide new insights into the relation between the topological structure of the population of micro-cracks and the macroscopic response of concrete. They are applicable to a range of quasi-brittle materials with similar dominant damage mechanisms. Keywords Concrete porosity; Lattice model; Cracking graphs; Macroscopic damage 1. Introduction The mechanical behaviour of quasi-brittle materials, such as concrete, graphite, ceramics, or rock, emerges from underlying microstructure changes. At the engineering length scale it can be described with continuum constitutive laws of increasing complexity combining damage, plasticity and time-dependent effects [1-4]. In these phenomenological approaches the damage represents reduction of the material elastic constants. From the microstructure length scale perspective damage is introduced by the nucleation and evolution of micro-cracks. While the population of micro-cracks formed under loading could be sufficiently well captured by various continuum damage models, the latter cannot help to understand the effects of the population on other important physical properties of the material. In many applications the quasi-brittle materials have additional functions as barriers to fluid transport via convection/advection and/or diffusion. It is therefore important to take a mechanistic view on the development of damage by modelling the evolution of micro-crack population, which can inform us about changes in the transport properties. Such a mechanistic approach needs to account for the material microstructure in a way corresponding to the mechanism of micro-crack formation [5]. Micro-cracks typically emerge from pores in the interfacial transition zone between cement paste and aggregate in cement-based materials [6]. Discrete lattice representation of the material microstructure seems to offer the most appropriate modelling strategy for analysis of micro-crack populations. This is a meso-scale approach, where the material is appropriately subdivided into cells and lattice sites are placed at the centres of the cells. Discrete lattices allow for studies of distributed damage without constitutive assumptions about crack paths and coalescences that would be needed in a continuum finite element modelling. The deformation of the represented continuum arises from the interactions between the lattice sites. These involve forces resisting relative displacements and moments resisting relative rotations between sites. Two conceptually similar approaches have been proposed to link local interactions to continuum response. In the first one, the local forces are related to the stresses in the continuum cell, e.g. [7, 8]. In the second one, the interactions are represented by structural beam elements, the stiffness coefficients of which are determined by equating the strain energy in the discrete and the
www.icfweb.orgRkJQdWJsaXNoZXIy MjM0NDE=