Composites Science and Technology 61 (2001) 997±1003
www.elsevier.com/locate/compscitech
An improved dynamic Monte Carlo damage microscopic numerical constitutive model incorporating thermo-mechanical coupling for unidirectional GRP under tensile impact Yuanming Xia, Yang Wang * Department of Mechanics and Mechanical Engineering (5th), University of Science and Technology of China, Anhui, Hefei 230027, PR China Received 13 April 1999; received in revised form 15 March 2000; accepted 19 April 2000
Abstract The phenomenon of thermo-mechanical coupling exists in the adiabatic process of unidirectional glass/epoxy composite (GRP) under tensile impact loading. In the present paper, the frictional work generated by matrix breakage or interfacial debonding is taken into account in order to improve the iteration calculation of thermo-mechanical coupling. An improved dynamic Monte Carlo damage microscopic numerical constitutive model incorporating thermo-mechanical coupling is set up. By means of this improved model, the deformation, damage and failure process of unidirectional GRP under tensile impact is simulated numerically. The apparent tensile stress/strain relationships of GRP at dierent high strain rates are predicted. The numerical results are in good agreement with the experimental data, which proves that the improved model is reasonable and eective. It also reveals that the dominant physical mechanisms of the non-linearity of the stress/strain curve for GRP are the temperature and strain-rate dependence of the mechanical properties of glass ®ber and the thermo-mechanical coupling. The frictional work must be considered during the iteration calculation of the thermo-mechanical coupling and its contribution to the transient temperature rise of GRP under tensile impact cannot be neglected. # 2001 Elsevier Science Ltd. All rights reserved. Keywords: A. Polymer-matrix composite; Tensile impact; Monte Carlo simulation; Thermo-mechanical coupling; Frictional work
1. Introduction The deformation, damage and failure process of unidirectional E-glass-®ber-reinforced epoxy composite (GRP) under tensile impact loading is a very complicated probabilistic process and can also be considered as an adiabatic irreversible working process. The dissipation energy generated by ®ber fracture, matrix breakage, interface debonding and frictional work can produce local transient temperature rises, which lead to non-unique transient temperature ®eld in the composite. These temperature rises can change the mechanical properties of the constituents of GRP, i.e. E-glass ®ber and epoxy (temperature and strain-rate dependence), and aect the microscopic stress ®eld and strain ®eld in the composite. On the other hand, the aected stress and strain ®elds will in¯uence the damage evolution of the GRP and the transient temperature rise. These two processes are *Corresponding author. E-mail address:
[email protected] (Yang Wang).
coupled and we de®ne this phenomenon as thermomechanical coupling. The Monte Carlo simulation technique [1±11] is one of the most eective methods of analyzing such complicated processes, but we must face two problems: (1) the choice of the numerical model, and (2) the determination of the representative volume element. The shear-lag model [12] is a suitable choice because it can re¯ect the natural characteristics of mechanical behavior of a brittle ®ber-reinforced polymermatrix composite under axial tensile loading. The determination of the representative volume element must obey the condition that its mechanical properties can represent those of the practical composite in experiment. From the above consideration, Xia et al. [9,11] proposed a composite element which is composed of ®ne meshes, a sucient number of ®bers, and an adequate composite length to simulate the tensile failure process of GRP. Moreover, they proposed the statistical average iteration method to uncouple the thermo-mechanical coupling [13], and set up a dynamic microscopic numerical constitutive model. The physical foundation
0266-3538/01/$ - see front matter # 2001 Elsevier Science Ltd. All rights reserved. PII: S0266-3538(00)00098-1
998
Y. Xia, Y. Wang / Composites Science and Technology 61 (2001) 997±1003
Fig. 1. Experimental stress/strain curve for loading±unloading tensile impact test on GRP.
where T
t is the magnitude of the temperature rise, t is time, J is Joule's equivalent, Cp is the speci®c heat at constant pressure, is the mass density, T0 is the ambient temperature, and l is the expansion coecient of the materials.
t and "
t are apparent tensile stress and strain of GRP, respectively. The ®rst term in Eq. (1) gives the temperature reduction during adiabatic elastic deformation process and its magnitude is very small. The second term represents the temperature rise arising from the dashed area in Fig. 2. The comparison of the theoretical temperature change curve from Eq. (1) with the experimental measured result is shown in Fig. 3. The solid curve in Fig. 3 is the stress/time curve of the specimen. It must be noted that the above mentioned theoretical results and experimental data are all macroscopic average values for GRP. The theoretical result is in a good agreement with the experimental data, con®rming that: (1) the dissipation work during tensile impact process, i.e. dashed area in Fig. 2, all transforms into heat and the energy analysis of Eq. (1) are valid; (2) the microscopic nonuniform transient temperature ®eld can be described by the macroscopic transient temperature and the interaction between the microscopic transient temperature ®eld and the microscopic stress/strain ®eld can also be described through the interaction between the macroscopic transient temperature and the mechanical properties of the constituents and macroscopic stress/strain curve. The macroscopic stress/strain curve and transient temperature rise are all the statistical average re¯ection of the microscopic non-uniform stress, strain and transient temperature ®eld. The previous simulation results [13] reveals that the strain rate and temperature dependence of the mechanical behavior of glass ®ber and the thermo-mechanical coupling are two dominant factors, which leads to the non-linearity of the stress/strain relation of GRP under tensile impact loading. When the strain rate is less than 300/s for GRP,
Fig. 2. Schematic diagram of stress/strain relation and dissipation work during tensile impact process for GRP.
Fig. 3. Transient temperature change versus time and stress versus time for GRP.
of this method and model is two experimental results of GRP under tensile impact. The ®rst experiment is a loading±unloading tensile impact test [14]. The typical experimental result is shown in Fig. 1. It is obvious that the loading part O±C is non-linear, but the unloading part C±Q is linear and passes through the origin. The second experiment is a transient temperature rise measurement test [15]. From the integral stress/strain curve of GRP (in Fig. 2) and the experimental result of the loading± unloading test, the theoretical transient temperature rise during the tensile impact process can be deduced as follows [16]: l
t T
t T0 1 exp JCp
"
1 1 1
t"
t
td"
t JCp 0 2
Y. Xia, Y. Wang / Composites Science and Technology 61 (2001) 997±1003
the consistency between the simulated results and the experimental data is very good. However, when the strain rate is greater than 800/s and the apparent strain of GRP is greater than 4%, there exists an obvious divergence between the simulated results and experimental results. The reason may be that some mechanism of generating heat during impact loading was not taken into account in the calculation of the thermo-mechanical coupling. The purpose of the present paper, based on the above numerical model, is to improve the iteration calculation of the thermo-mechanical coupling; namely, take the frictional work generated by the matrix breakage or interfacial debonding into account to improve the numerical model and reveal further the dynamic behavior and damage mechanism of unidirectional GRP under tensile impact loading. 2. Improved damage microscopic numerical constitutive model 2.1. Shear-lag model for unidirectional ®ber-reinforced composite The con®guration of the unidirectional GRP used in the present work is shown in Fig. 4. The model consists of N ®bers with length L. In this model, it is assumed that the ®bers sustain the axial force and the matrix transfers through shear the force lost at broken ®bers to the neighboring ®bers. Thus, the force-equilibrium equation for the ®bers is given as follows: di;j A dx i1;j dx
i;j
@2 ui;j Hdx A 2 dx @t
2
where is the density of the composite, H is the thickness of the composite. i;j and "i;j are the normal stress working on the ®ber element from (i; j 1) to (i; j) and the shear stress working on the matrix from (i; j) to (i 1; j), respectively. ui;j is the axial displacement of the nodal point (i; j).
999
The relation between the i;j ; "i;j and displacement ui;j is given by i;j E
ui;j
ui;j x
1
; i;j G i;j G
ui;j
ui D
1;j
3
where E is the modulus of ®ber and G is the shear modulus of matrix. D is the average span between neighboring ®bers. The relation between H and D is simpli®ed as follows: Am Af
1=Vf D=H
1
1; Am HD; Af H2
Vf =Vf
4
where Af is the cross-sectional area of the ®ber, Am is the cross-sectional area of matrix, and Vf is the volume fraction of the ®ber. The left end of the model is ®xed and the right end moves in a constant speed V0 . At the initial time, the model is stable. A ®nite-dierence method of uniform space mesh size x and time step t is used to solve Eq. (2). Thus, the force equilibrium Eq. (2) can be rewritten in the form of displacements of the nodal points. By applying displacement on the right end of composite step by step and by means of the successive over-relaxation algorithm, we can obtain the displacements at all nodal points at any time. The detailed information about the ®nite-dierence method and procedure can be found in Refs. [11,13]. 2.2. The failure criteria of ®ber and matrix The mechanical properties of the E-glass ®ber have temperature and strain-rate dependence. The tensile strength of single glass ®ber follows the double Weibull distribution function [17]. " # 1 2 H
x; 1 exp
5 01 02 x x
Fig. 4. Schematic representation of the model composite.
1000
Y. Xia, Y. Wang / Composites Science and Technology 61 (2001) 997±1003
where H
x; is the failure probability of ®ber element with length x under an applied stress no greater than at a certain temperature and strain rate. 1 , 2 01 02 and x ; x are the shape and scale parameters of double Weibull distribution function, respectively. The scale parameters are dependent on the temperature and strain rate. If a random number is taken from the uniform distribution on the interval (0,1), we can generate the tensile strength of each ®ber element by means of the Newton iteration method. The failure criterion of the ®ber element is, i;j 5Si;j
6
where Si;j is the tensile strength of ®ber element from (i; j 1) to (i; j) at a certain temperature and strain-rate state. The shear strength of all matrix elements are identical. The failure criterion of matrix element is, i;j 5max
7
where max is the shear strength of the matrix at a certain constant strain rate. When the matrix or interface breaks, we assume that a constant friction force, fric , exists between the matrix and ®bers, and we let fric max [9,11,13]. As a result of the limitation of the twodimensional shear-lag model, the matrix breakage and the interface debonding cannot be distinguished. Stated concisely, these two cases are both called matrix breakage through this paper. 2.3. Improved statistical average iteration method Fig. 5 shows the schematic diagram of the thermomechanical coupling iteration calculation of GRP under tensile impact. The curve OABn Dn represents the simulated apparent stress/strain relation of GRP. At time t (point A), there is no damage occurring. At this time the apparent stress and strain of GRP are ! N EVf X c ut M uti M 1 uti 2 uti 1 2Nx i1 i
8 uti ;M V0 t "c L L When the ®ber fracture and the matrix breakage occur at the interval time from t to t t, the dissipation energy of the GRP is the area of OAB0 in Fig. 5. De®ne this dissipation work as,
Wf
"
tt 0
t td"
t t
1
t t"
t t 2
9
Fig. 5. Schematic diagram of the iteration calculation for thermomechanical coupling.
De®ne the work done against friction as, Wfric fric
xH
D
10
where is the change of the shear strain of the matrix element from the maximum level, i.e.
max
max G
fric be the sum of the frictional work for all Let W P fric Wfric . failure matrix elements, then we have W According to Wf and Wfric , the transient temperature rise of GRP at time t t can be obtained as follows, " # fric 1 W Wf T
t t
11 JCp Vc where Vc is the volume for the total composite element. Comparing Eq. (11) with Eq. (1), we can ®nd that there exists two dierences: (1) the term expressed by Kelvin's formula is not included in Eq. (11) because the dynamic elastic modulus of the glass ®ber has been used in the present model; (2) the frictional work term is added in Eq. (11). The reason is that the experimental stress/strain curve of GRP has implicitly included the contribution of frictional work to the transient temperature rise of GRP but, during the numerical simulation process, the frictional work cannot aect the simulated apparent stress/strain curve directly. Therefore, the heat generated by friction work must be considered individually in the numerical simulation.
Y. Xia, Y. Wang / Composites Science and Technology 61 (2001) 997±1003
T01
t t obtained directly from Eq. (11) is de®ned as the initial temperature rise at time t t. Neither the eect of T on the mechanical properties (modulus and strength) of the glass-®ber element and microscopic stress/ strain ®elds and the reaction of the aected stress/strain ®elds on the temperature rise are taken into account in the numerical simulation. Hence, the calculated apparent stress and strain of GRP are not accurate. We must consider such interaction, i.e. redetermine the strength (temperature dependent) of all ®ber elements and recalculate the displacements of all nodes, the stress and strain of all ®ber and matrix elements and the apparent stress and strain of GRP at time t t. Point B1 in Fig. 5 represents the new apparent stress and strain. This calculation process is called the ®rst thermo-mechanical coupling iteration on the interval t; t t. Based on Eq. (11), whether or not there is a new fracture of the ®ber element and a new breakage of the matrix element, the new temperature rise can be calculated. The magnitude of the dissipation work equals the sum of the triangular area OAB1 and the fric with increase of the frictional work, i.e. substitute W
1001
2 1 in Eq. (11), where W 2 and W 1 2 W W tt W fric fric fric fric are new frictional work and ®rst frictional work obtained from the ®rst iteration process. Thus, we can obtain the new nodal displacements and the new apparent stress and strain (corresponding to point B2 ) at time t t. This process is called as the second thermo-mechanical coupling iteration in the interval t; t t. Repeat the above sec fric with ond iteration procedure and substitute W n n n 1 in Eq. (11) until that there is no W W W tt fric fric new breakage of ®ber and matrix element, the slope of ABn n is not less than that of ABn 1 , and W tt tends to zero. Thus, the ®nal apparent stress and strain of GRP at time t t can be determined (corresponding to point Bn), ( N Vf X c
t t c
t E
t t utt utt i;M i;M 1 Nx i1 uti;M uti;M 1 "c
utt V0
t t i;M L L
Fig. 6. Comparison between the numerical results and the experimental data.
12
1002
Y. Xia, Y. Wang / Composites Science and Technology 61 (2001) 997±1003
Fig. 7. Fracture morphology of GRP at strain rate 1100/s.
Similar to the procedure of calculation and iteration from t to t t, the apparent stress and strain at time t 2t can be obtained. The initial temperature rise T02 at t 2 t can be obtained from the area OABnD0 fric and the increase of frictional work, i.e. substitute W t2t tt with Wt2t Wfric Wfric in Eq. (11). Repeat the process of second thermo-mechanical coupling iteration at interval [t t; 2 t] and the ®nal apparent stress and strain at time t 2 t can be determined. Repeat the above procedure and the apparent stress and strain at any time can be obtained until the failure of the composite. This method is the improved statistical average iteration method incorporating the frictional work. 3. Results and discussion The input values of the unidirectional GRP for the present simulation such as Af , Vf , L, M and N are the same as those in Ref. [11]. The mechanical constants for E-glass ®ber and epoxy resin are taken from Refs. [17] and [18], respectively. Fig. 6 shows the comparison between the numerical results incorporating frictional work and the experimental data. The dash curves in Fig. 6 are the numerical results in which the frictional work is not taken into account in the process of thermo-mechanical coupling calculation. It is obvious that there exists apparent divergence between the numerical results and the test data. The main problem is that the calculated maximum stress of GRP is greater and the unstable strain (the strain corresponding to the maximum stress) is smaller compared with the test data. By taking into account the frictional work, we can see that the numerical results are in better agreement with the test ones. This indicates that the improved numerical model incorporating frictional work is reliable and reasonable and it can re¯ect the physical nature of the deformation, damage and failure process of unidirectional GRP under tensile impact loading.
Fig. 7 shows the micro-failure morphology of ®ber and matrix elements when the apparent strain reaches 4%. It is apparent that at this time a great number of broken matrix elements appear. It is unreasonable if the frictional work is not included in the calculation of the dissipation energy. Fig. 8 shows the eect of frictional work on the transient temperature rise of GRP. It can be seen that when GRP is in the elastic deformation state, there is no transient temperature rise in the composite, i.e. there is no heat dissipation, which means that no any microdamage of ®ber fracture and matrix breakage appears. Following the increase of loading, the value of the transient temperature rise increases gradually. This means failure of ®ber and matrix occurs. After the apparent strain exceeds the unstable strain, the transient temperature rises rapidly. At this time, a large number of failure ®bers and matrices have appeared. It is obvious that the frictional work leads to the increase of the transient temperature rise, especially after the apparent strain of GRP is greater than 4%, and the contribution of frictional work to the dissipation energy cannot be neglected.
Fig. 8. The eect of frictional work on the transient temperature rise of GRP.
Y. Xia, Y. Wang / Composites Science and Technology 61 (2001) 997±1003
4. Conclusions 1. An improved dynamic damage microscopic numerical constitutive model incorporating thermomechanical coupling was set up. The contribution of frictional work to the transient temperature rise of GRP is taken into account in the iteration calculation of thermo-mechanical coupling. By means of this improved model, the deformation, damage and failure process of unidirectional GRP under tensile impact loading is simulated and the apparent stress/strain relations of GRP at dierent strain rates, especially the maximum stress and the unstable strain, are predicted numerically. The numerical results have a good consistency with experimental data, which proves that the improved constitutive model is valid and reliable. The consistency also indicates that the dominant physical mechanism of the non-linearity of the stress/strain curve for GRP is temperature and strain-rate dependence of E-glass ®ber and thermo-mechanical coupling. 2. The numerical results reveal that the eect of frictional work on the transient temperature rise cannot be neglected, especially when the apparent strain of GRP exceeds 4%. At this moment, the number of broken matrix elements is large and the existence of the friction work is inevitable. The dissipation energy generated by the frictional work will aect the transient temperature rise, and then aect the apparent stress/strain relation of GRP.
Acknowledgement The present work was supported by National Natural Science Foundation of China.
1003
References [1] Kimpara I, Watanabe I, Okatsu K, Ueda T. The 7th Symp. of Composite Materials. JUSE 1974;169. [2] Kimpara I, Ocaki T. Proc. ISCMS (1st) Beijing, 1986, p. 682. [3] Fukuda H, Kawata K. On the strength distribution of unidirectional composite. Fiber Sci Tech 1977;10:53. [4] Oh KP. A Monte Carlo study of the strength of unidirectional ®ber reinforced composite. J Comp Mater 1979;13:311. [5] Fukuda H, Chou TW. Monte Carlo simulation of the strength of hybrid composite. J Comp Mater 1982;16:371. [6] Manders PW, Bader MG, Chou TW. Monte Carlo simulation of the strength of composite ®ber bundles. Fiber Sci Technol 1982;17:183±204. [7] Lienkamp M, Schwartz P. A Monte Carlo simulation of the failure of a seven ®ber microcomposite. Comp Sci Technol 1993;46:139±46. [8] Goda K, Phoenix SL. Reliability approach to the tensile strength of unidirectional CFRP composite by Monte Carlo simulation in a shear-lag model. Comp Sci Technol 1994;50:457±68. [9] Yuan JM, Xia YM, Yang BC. A note on the Monte Carlo simulation of the tensile deformation and failure process of unidirectional composites. Comp Sci Technol 1994;52:197. [10] Dibenedetto AT, Gurvich MR. Statistical simulation of ®ber fragmentation in a single-®ber composite. Comp Sci Technol 1997;57:543±55. [11] Wang Z, Xia YM, Yuan JM. A dynamic Monte Carlo simulation for the unidirectional composites under tensile impact. Comp Sci Technol 1998;58:487±95. [12] Hedgepeth JM. NASA technical note, D-882, 1961. p. 1. [13] Xia YM, Wang Z. Dynamic damage microscopic numerical constitutive model incorporating thermo-mechanical coupling for unidirectional composite. Comp Sci Technol 1999;59:947±55. [14] Xia YM, Wang X. Consititution equation for unidirectional composites under tensile impact. Comp Sci Technol 1996;56:155. [15] Xia YM, Rao S, Yang B. An infrared transient temperature measuring apparatus and its application to the tensile impact testing. Experimental Mechanics 1990;5:170±6 (in Chinese). [16] Xia YM, Dong L, Ran S, Yang B. Energy transition of ®bre and their unidirectional composite under tensile impact. In: Proc. ICCM-8, Honolulu, 1991, 0±29. p.1±8. [17] Wang Z, Xia YM. Experimental evaluation of the strength of ®bers under high strain rates by bimodal Weibull distribution. Comp Sci Technol 1997;57:1599±607. [18] Xia YM, Li M. Acta Material Composite Sinica 1987;2:59±65 (in Chinese).