- Email: [email protected]
The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic materials
Journal Pre-proof
The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic materials Junling Hou , Chao Zhang , Qun Li PII: DOI: Reference:
S0167-6636(19)30948-2 https://doi.org/10.1016/j.mechmat.2020.103363 MECMAT 103363
To appear in:
Mechanics of Materials
Received date: Revised date: Accepted date:
31 October 2019 11 February 2020 13 February 2020
Please cite this article as: Junling Hou , Chao Zhang , Qun Li , The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic materials, Mechanics of Materials (2020), doi: https://doi.org/10.1016/j.mechmat.2020.103363
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
Highlights
A new form of M-integral associated with time dependence parameters is presented for viscoelastic materials. Time-dependent M-integral is numerically implemented based on the equivalent domain integral method as an effective and accepted fracture mechanical parameter for damage induced by crack growth in viscoelastic materials. Numerical examples are given to demonstrate the validity and relevance of the time-dependent M-integral in viscoelastic material.
Unmarked The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic materials Junling Hou, Chao Zhang, Qun Li* * Corresponding author, email: [email protected] State Key Laboratory for Strength and Vibration of Mechanical Structures, Shaanxi Key Laboratory of Environment and Control for Flight Vehicle, Shaanxi Engineering Research Center of Nondestructive Testing and Structural Integrity Evaluation, School of Aerospace, Xi’an Jiaotong University, Xi’an 710049, China
Abstract A new form of M-integral associated with time dependence parameters, is presented herein for viscoelastic materials. Based on the equivalent domain integral method, this time-dependent M-integral is numerically implemented as an effective and accepted fracture mechanical parameter for damage induced by crack growth in viscoelastic materials. Based on the linear viscoelastic model defined through Prony series, the conservation of the time-dependent M-integral for viscoelasticity is verified by applying user defined Python scripts. The results show that the newly proposed time-dependent M-integral can be successfully calculated numerically. Furthermore, numerical examples are given to demonstrate the validity and relevance of the time-dependent M-integral in viscoelastic material. In particular, the variations of the time-dependent M-integral for different loading rates are considered, which show that the crack growth behavior over a period of time in viscoelastic material can be evaluated based on the value of the time-dependent M-integral. In addition, the time-dependent M-integral is calculated to assess the damage degree induced by two collinear cracks.
Keywords: Viscoelastic material; M-integral; Time-dependent; Domain integral method; Crack growth.
1. Introduction Owing to the widespread use of time-dependent viscoelastic materials such as solid propellant, Ni-alloy, dispersed fuel, bituminous concrete, and polymer matrix composite etc., it is necessary to have a better understanding of their fracture mechanism for achieving a material design and engineering application. Considering the influence of the manufacturing process, loading condition, and temperature change, it has been found that the crack initiation, evolution and propagation may result in a failure in viscoelastic materials (Yuval et al., 2018; Mohammed et al., 2017; Yi et al., 2014). Thus, finding an effective method becomes critically important for the evaluation of the damage degree and evolution mechanism in viscoelastic materials during the fracture process. For instance, in the framework of linear elastic fracture mechanics, two fracture parameters (stress intensity factor and energy release rate) are employed to investigate the growth and resistance of cracking in viscoelastic materials (Contino et al., 2019). Using precise time-domain expanding algorithm, the viscoelastic cracked problem is solved using converted quasi-elastic variables, such as displacement and stress (Yao et al., 2018). A viscoelastic damage model is proposed for the calculation of crack width in viscoelastic concrete (Sciumè and Benboudjema, 2016). The maximum tensile stress is introduced as a fracture criterion to deal with the crack initiation process in rubber material with viscous-elastic properties (Fukahori et al., 2008). Studies have shown that it is of practical importance to achieve a deeper understanding of the ‘global’ damage state from defect formation, expansion or other factors in viscoelastic materials. That is, for multi-cracked viscoelastic materials, more effective and universal methods are still required. In recent years, the conservation integral method has attracted much attention owing to its advantages in dealing with the complicated damage problems found in different materials (Wang et al., 2016; Yu and Kitamura, 2015; Yu et al., 2009). For the crack problem in both homogeneous and non-homogeneous materials, the domain independent interaction integral method has been successfully used in the calculation of the stress intensity factor (Yu and Kitamura, 2015; Yu et al., 2009). In addition, the energy-conservation integrals, Jk-, M- and L- integrals, have been
proved to be reliable and effective in assessing the fracture properties of materials (Azmi et al., 2017; Li et al., 2017; Motola and Banks-Sills, 2009). In particular, the M-integral, as one of the key concepts in material configurational mechanics, has shown its advantage in quantifying the damage degree for multi-crack problems since it is inherently related to the change of the total potential energy in a damaged system independent of the specific damage details of every defect (Li et al., 2017; Motola and Banks-Sills, 2009). In other words, its value is directly related with all fracture factors such as the material constants, defect characteristic, and mechanical loading. Thus, this approach is much easier in practice and has been widely applied to elastic, elastic-plastic and visco-elastic-plastic problems. Numerous studies have been carried out using the M-integral to evaluate the material damage evolution (Chang and Lin, 2007; Yu et al., 2012; Chang and Wu, 2011; Yu and Li, 2013; Chang et al., 2013; Judt and Ricoeur, 2016). For instance, Chang and Lin (2007) presented a fracture parameter called Mc-integral to investigate the multiple curved crack problem in rubbery materials under large deformation. Using digital image correlation method, the M-integral was measured by Yu et al. (2012) in elastic-plastic samples created using LY12 aluminum. This technique can be extended to evaluate the damage induced by multiple defects in materials. Furthermore, according to a study by Chang and Wu (2011), the energy related M-integral is considered as a damage parameter to describe the irreversible fracture problem in a solid weakened by multi-cracks. Based on the material configurational force theory, Yu and Li (2013) developed the application of M-integral in elastic-plastic material with multi-defects, such as cracks and voids. For a three-dimensional problem, Chang et al. (2013) propose a new energy parameter called Mc-integral for the evaluation of the structural integrity in elastic media with multiple cracks. For engineering applications, the M-integral was introduced by Judt and Ricoeur (2016) to calculate the loading quantities of engineering structures with two cracks. Based on above studies, it is of practical importance to introduce the implementation of M-integral for the assessment of ‘global’ damage state from defect
formation, expansion in materials with a viscous effect. Particularly, the effective and universal methods via the M-integral are necessarily extended for the cracked viscoelastic materials, in which the stress relaxation and the rate correlation as specific material properties of viscoelastic material are taken into account. In contrast with the traditional fracture parameters, for example, stress intensity factor K (only valid in cracking problems of elastic materials) and Jk-integral (can be extended for complicated materials, such as Ct -integral for time-dependent creep behavior, but not feasible for multiple cracks due to their ‘local’ nature associated with a single crack tip). The values of M-integral is directly related with all fracture factors such as material constants, defect characteristic, mechanical loading and etc. The novelty of this paper is to extend the concept of M-integral to deal with the failure of viscoelastic materials. In addition, although the M-integral is a path-independent integral parameter used to measure the damage degree of materials, it is cumbersome to compute in finite element calculation. For numerical implementation, alternative approaches are needed to simplify the numerical computation code. Among them, a numerical technique termed the equivalent domain integral method is recommended. For the equivalent domain method, all necessary variables for numerical computation are easily attainable in a finite element analysis. There is no need to consider the influence of crack tip singularity, or the difficulty in catching the Gauss point on a specific path. Thus, the EDI method appears to be a powerful and attractive tool for transforming the line integral into an equivalent domain integral. For instance, it has been introduced to evaluate the values of conservation line integral in elastic-plastic materials (Li et al., 1985), thermo-elastic materials (Shih et al., 1986), and linear anisotropic elastic materials (Chang and Pu, 1997). The concept of domain integral method is also extended to calculate the value of the M-integral for the measurement of nickel-based superalloys with multi-cracks (Hou et al., 2018). However, in previous study, the definition of M-integral is built in the elastic media. In addition, to evaluate the damage degree of viscoplastic material, the geometry model for finite element calculation is consists of two parts, one of which is the damage zone defined
as viscoplastic media, and the other is a pure elastic zone that surrounds the damage zone. Considering that the M-integral in other fields has been compared with experimental data (Li et al., 2017), numerical results (Hou et al., 2018) in our previous work, it shows that the M-integral is applicable for the validation of the cracked problems in materials. In present study, the work of the time-dependent M-integral is only presented by the theoretical explanations and numerical simulations, without comparison with experimental data. The main purpose of the present paper is to introduce the concept of M-integral into the viscoelastic material owing to its advantage in describing the material fracture behavior under the presence of cracks. In this paper, a computationally efficient and convenient method for the evaluation of the M-integral for defects in two dimensional viscoelastic material is presented. The remainder of this paper is structured as follows. Firstly, the linear viscoelastic model involved in present study is given in Section 2. Then, the derivation of M-integral for viscoelastic material is carried out in Section 3. Next, for numerical implementation, the procedure of converting line M-integral into domain integral is presented in detail. With the newly built viscoelastic domain M-integral, several numerical examples are provided and discussed in Section 5 to illustrate the availability of the approach.
2. Viscoelastic material model The classical viscoelastic model and its application in finite element method are presented in this section. For viscoelastic media, it is known that the stress distribution and material deformation are both time-dependent variables. Specifically, the viscous deformation produces an irreversible flow with the external force work, and all of its energy is converted into internal energy. For linear viscoelastic material, its stress response can be basically described through various combinations of two types of mechanical model. Two classical linear viscoelastic models have been proposed, namely, Maxwell model (related with the stress relaxation) and Kelvin model (related with creep). In the present study, the time domain viscoelasticity defined in ABAQUS is considered under a stress relaxation load, and thus the
Maxwell model is chosen. The basic linear viscoelastic response of the Maxwell model (Gooch, 2011) is described based on the sum of two mechanical models. One is the pure elastic model which obeys Hooke’s law, and the other is a viscous component satisfying Newton’s viscosity law. Under the loading condition (t), the definition formation of the time varying strain can be written as follows:
=
E
(1)
Here, =d / dt represents the strain rate, E is the elastic modulus, and is viscous coefficient. In addition, Eq. (1) is the differential constitutive equation of the Maxwell model. To describe more complicated stress states in viscoelasticity, a generalized Maxwell model (Gooch, 2011) consisting of multiple Maxwell models is proposed. For a two dimensional problem, the mechanical behavior can be illustrated by considering a sudden strain 0 applied to the viscoelastic media and held constant for a long time. That is, considering a varying strain (t) = 0 H(t) applied to the viscoelastic material, the stress response is obtained as follows: n
(t ) 0 Ei e
t Ti
(2)
i 1
Here, H (t) denotes the unit step function, which is defined as follows:
0, t 0 H (t ) 1, t 0
(3)
Here, Ti = i / E corresponds to the relaxation time of i viscous component, while i is the viscosity tensor. And the stress relaxation modulus of the generalized Maxwell model can be written as below (Gooch, 2011) n
E (t ) Ei e
t Ti
(4)
i 1
Eq. (4) is the Prony series form of time-dependent relaxation modulus for viscoelastic material in which the parameters can be directly defined in ABAQUS material.
Otherwise, in finite element analysis, to begin with, an elastic response of the viscoelasticity should be defined. In this paper, the small-strain deformation condition is considered and thus the rate-independent elastic response is defined using a linear elastic material model.
3. Definition of M-integral in viscoelastic material Within the framework of the material configurational force theory (Eshelby, 1975; Li et al., 2013), the physical interpretation of the M-integral is explained as the potential energy change during the process of the self-similar expansion of defects. Thus, the stability of materials with multiple defects, that is, the damage degree of materials can be characterized by the value of the invariant M-integral. The M-integral has recently shown its advantages as an energy fracture parameter in evaluating the degradation of different materials (e. g., elastic-plastic, piezoelectric and general anisotropic materials) with multi-defects (Yu and Li, 2013; Bank-Sill et al., 2008). In the present study, a new form of M-integral for viscoelastic material is defined by considering the crack growth behavior. Considering two-dimensional media with cracks, regardless of the inertia terms and body force in the system, the formula of the M-integral can be written by (Knowles and Sternberg, 1972; Budiansky and Rice, 1973; Herrmann and Herrmann, 1981) M
( w
x jk uk ,i xi ) n j ds
ij i
M
j
n j ds
(i, j , k =1, 2)
(5)
It should be noted that the present study is applicable for two dimensional analyses. For three dimensional conditions, the M3D-integral is written as follows:
1 wij xi jk uk ,i xi jk uk n j da A 2
M3D
(6)
In which, is the total volume surrounded by domain A containing all defects. It is obvious that an additional term 1 jk n j uk must be considered. 2
Based on the concepts of material configurational force theory, the physical meaning of Mj is explained as the material damage extension configuration force,
which is shown as follows
M j w ij xi kj uk ,i xi
(7)
where xi corresponds to the coordinate, ij is the Kronecker symbol, σij refers to the Cauchy stress tensor, uk denotes the displacement tensor, and the subscript prima {},i denotes the corresponding differentiation with respect to the coordinate xi. All cracks are surrounded by the enclosed path with the outer normal vector n j. For viscoelastic material, the strain energy density w(t) can be defined as t
w(t ) ij d ij (t ) ij ij dt ij
0
(8)
Here, εij (t) represents the strain tensor which is related with the load time t. Assuming that the total strain tensor ε could be expressed as follows: (Simo and Hughes, 1998, Chapter 10; Nguyen et al., 2005) ε εe εv
(9)
Here, εe and εv are the elastic and viscous strain tensors, respectively. In the same manner, the following additive split of the stress response tensor is given in the following: σ σ E σ IE
(10)
Here, σE denotes the elastic and equilibrium stress behavior, while σIE corresponds to the time-dependent and inelastic behavior in viscoelasticity. Combining Eqs. (9) and (10), the form of Eq. (8) can be rewritten as follows: w ε σ : ε = σ E σ IE : ε ε v ε v = σ E : ε σ IE : ε ε v σ IE : ε v
(11)
= σ : ε v IE
Here, Ψ (ε, εv) is treated as the free energy density. Then, introducing the equivalent relation in Eq. (11) to Eq. (5), the rate form of the M-integral can be obtained, which is shown as follows u ) x (σ IE : ε v ) x ]n ds M [( ij jk k ,i i ij i j
(12)
Rewriting Eq. (12), the new definition of the time dependent M-integral for viscoelasticity is given as
M (t )
b (t ) x w ij
i
IE
(t ) ij xi n j ds = M Vj n j ds
(13)
where bij (t) denotes the Eshelbian energy momentum for the viscoelastic material, and is defined as follows:
bij (t ) = ij jk uk ,i
(14)
The symbol wIE (t ) is applied to characterize the part of strain energy density induced by the inelastic stress response and viscous strain deformation, that is, wIE (t )=σ IE : ε v . That is, for viscoelastic material, the material damage extension
configuration force M Vj is expressed as follows:
M Vj (t )=bij (t ) xi wIE (t ) ij xi
(15)
Referring to the derivation of this newly built M-integral for viscoelastic material, it is clear that the value of the M-integral is calculated along a path surrounding all the defects but without consideration of the details, such as how many cracks or the shapes of the cracks. Thus, this time-dependent M-integral method can be extended to solve other defects problems such as multi-cracks, interface debonding in viscoelasticity.
4. FEM calculation of M-integral in viscoelastic material
Γs
Viscoelastic medium
dA
Γ
Γ0 Cracks
A
Fig. 1. Description of the equivalent domain in viscoelasticity.
For numerical implementation, the time-dependent M-integral obtained in Eq. (13) can be transformed into a form of a domain integral by introducing the equivalent domain integral method. Initially, a weight function q is introduced for mathematical
manipulation, which represents arbitrary functions that smoothly go through the integral region A. Area A corresponds to the shadow region shown in Fig. 1. Here, the q-function is simply a transition function formed as q (x1, x2) (Becker et al., 2012; Niikishkov and Atluri, 1987), and it has been indicated that the choice of q-function is not sensitive to the results of the M-integral. Thus, the values of q-function are set as 0 and 1 for the outer () and inner (0) contours, respectively. We thus have the following:
M (t ) 1 M Vj n j ds 0 M Vj n j ds 0
V j
M n j qds M Vj n j qds 0
(16)
=
0
M Vj q n j ds
The following identity integral relationship is introduced,
M q A
V j
,j
dA
M A
V j, j
q M Vj q, j dA
(17)
Using the divergence theorem on the left side of Eq. (17) and combining Eq. (16), we obtain the following:
0
M Vj q n j ds = M Vj q dA M Vj , j q M Vj q, j dA A
A
,j
(18)
Considering the free of singularity in zone A, the relation M Vj, j =0 is obtained, and then Eq. (18) can be rewritten as M (t )= M Vj q, j dA A
q = bij (t ) xi wIE (t ) ij xi dA A x
(19)
j
For the convenience of calculation, Eq. (19) can be separated as the superposition of two parts:
M (t ) M I (t ) M II (t )
(20)
with
M I (t ) bij (t ) xi A
q dA x j
(21)
M II (t ) wIE (t ) A
q dA x j
(22)
The form of function q in zone A can be obtained by the shape function interpolation method, which is shown below: 4
q N1q1e N 2 q2e N 3 q3e N 4 q4e Nq e N i ( r , s ) qie
(23)
i 1
where Ni represents the shape function which is only associated with the coordinates of each element, and (r, s) is the position of the corresponding Gauss point. In the present study, four Gauss points are used for the integration in entire equivalent domain. For a single element, rewrite Eqs. (21) and (22) as follows:
M Ie
1
1
I r, s drds
1 1 I
M IIe
1
1
I
1 1 II
r , s drds
(24) (25)
Applying the Gauss integral method, Eq. (24) is approximately equivalent to
M Ie I I r1 , s1 I I r2 , s2 I I r3 , s3 I I r4 , s4
(26)
u u v v q x x y xy x xy x x x y x x x x e I I (r , s ) det( J ) u u v v q y xy x x xy y y y x x y x x y
(27)
with
where det(Je) is the determinant of the Jacobi matrix for the specific element. Similarly, the value of M IIe shown in Eq. (25) equals
M IIe I II r1 , s1 I II r2 , s2 I II r3 , s3 I II r4 , s4
(28)
q q I II (r , s ) xIE xv x1 xyIE xyv x1 xyIE xyv x2 yIE yv x2 det(J e ) x y
(29)
with
In particular, because the integral domain A shown in Fig. 1 is far from the cracked area, the calculation of the M-integral is not sensitive with the crack tip
behavior. Thus, singular mesh is without consideration for crack tip region. For the calculation of Eqs. (27) and (29), the values of ij,
ui q , and can x j x j
be obtained directly from the numerical simulation in ABAQUS postprocessing. Otherwise, the values of ,
IEe and ijv will be calculated from the identity
relation shown in this section. Then, the value of Me in a single element can be obtained by superposing the calculated M Ie and M IIe . Repeat above procedure for every element in zone A, the sum of Me is equal to the corresponding value of the M-integral for chosen integral contour. And with the user defined Python script employed in ABAQUS, the numerical calculation of M-integral can be realized in viscoelastic material.
5. Results and discussions In this section, several numerical examples are presented to show the validity of the newly proposed M-integral for viscoelastic material. For numerical simulation, the crack growth behavior is evaluated using the extended finite element method (XFEM). Differing from the traditional finite element, there is no need to do remeshing for the XFEM. Thus, the mesh continuity can prove the effective transfer of stress and strain. Upon request of the XFEM, the maximum principle stress criterion is set as the damage initiation law to control the behavior of crack growth (max = 1.91 MPa) (Lv, 2013), and the damage evolution law is linear softening stress-strain response whose type is based on displacement, setting u max = 0.25 mm at failure (Lv, 2013). Combining the built-in viscoelastic material model in ABAQUS, numerical calculations of the viscoelastic solid propellant with crack growth could be carried out. The computational cost is between 5~10 minutes for most of the examples simulated using double Intel core i7 CPU (i7-8700 CPU @ 3.2 GHz 3.19GHz). Specifically, for the model containing one single crack described in Section 5.2 and 5.3, a job in ABAQUS could be completed within 5 minutes, and the Python script calculating the visco-elastic integrals can be completed within 5-15 minutes according to the
numbers of the outputting M-integral for different time steps.
5.1 Material properties of viscoelastic solid propellant The viscoelastic material chosen for numerical verification is a common type of solid propellant called hydroxyl-terminated polybutadiene (HTPB), which has been widely used as the fuel in solid rocket engine and shows linear viscoelastic properties under sustained loading. Considering the process of manufacturing, transport or application, it is inescapable that some defects such as cracks will occur in solid propellant owing to the complex load conditions. Research has shown that the damage evolution (e.g. crack growth) in solid propellants will likely lead to a catastrophic accident owing to incomplete combustion. Thus, it is valuable to study the damage mechanism of such viscoelastic solid propellant. Referring to the viscoelastic model presented in Section 2, the Prony series form of time-dependent relaxation modulus is given as follows (Lv, 2013):
E (t ) 4.10 (0.15 0.12e
t
0.04
0.17e
t
0.40
0.25e
t
4
0.32e
t
40
)
(30)
Here, the unit of time t in Eq. (30) is second. In addition, the initial material parameters of viscoelastic solid propellant (HTPB) are shown in Table 1 for completeness.
5.2 Numerical verification of time-dependent M-integral conservation σ
Viscoelastic plate Crack Equivalent a
A1 A2
σ
integral
Fig. 2. Sketch of viscoelastic plate with central crack and equivalent integral domain.
To start with, the validity of a Python script for time-dependent M-integral is examined based on the integral conservation. As shown in Fig. 2, under plane strain condition, the geometry with a central single crack subjected to uniform tension load σ is investigated. The 2D plate has a dimension of 45 mm × 45 mm (length × width) with a central crack (a = 10 mm). The finite element mesh used in this study consists of four-node plane strain elements with full integration in this paper. The applied tensile stress σ is 0.50 MPa over a period of 1 second. In Table 2, the values of time-dependent M-integral for six different areas of equivalent integral domain in viscoelastic material are given. Specifically, the dimensions of six areas vary from 40 to 144 mm2, while the range of length and width are (32-39 mm) and (37-41 mm), respectively. It can be seen from Table 2 that the maximum value of the time-dependent M-integral is 48.33 J/m which is obtained in domain No. 2, whereas the minimum value is 46.47 J/m for domain No. 1. The comparison results show that the error between them (i.e., 4.00%) is controlled at an acceptable level, which means that the results of time-dependent M-integral for different equivalent integral domains show a good conservation nature, which verifies that the numerical implementation of the time-dependent M-integral has been successfully carried out. 5.3 Effects of loading rates on M-integral in viscoelastic material Considering the time dependence of the newly built M-integral in viscoelastic material, the variations of the M-integral under three different loading rates (i.e., 0.90 MPa/s, 0.09 MPa/s, 0.009 MPa/s) are simulated in this section. The model geometry is shown in Fig. 2, that is, a 2D plate with one central crack. The load condition is a purely monotonic increasing loading with a peak value of 0.90 MPa over three time periods (i.e., 1 s, 10 s, 100 s) as shown in Fig. 3, while the time increment is set to the same value as 0.01 s for all simulations. The dimensions of equivalent integral domain are chosen as A1=34.0 and A2=38.0. The curves of the computed time-dependent M-integral under three different
loading rates (Fig. 4) show that the general tendency of time-dependent M-integral increases with increasing load time steps for each loading rate, which confirms that the damage degree of viscoelastic plate increases with the crack growth. For instance, at end of the load time steps, the values of time-dependent M-integral are 113.43 J/m, 110.19 J/m, 110.02 J/m corresponding to the three loading rates as marked in partial enlarged drawing in Fig. 4. (MPa) 0.9
1
0
10
t(s)
100
Fig. 3 Load conditions of different loading rates. 140
M=121.61 J/m
Loading rate = 0.9 MPa/s
M-integral (J/m)
120 100 80
M=99.85 J/m
60 40
M=101.43 J/m
20 0 0.0
0.2
0.4 0.6 0.8 Time step (s) (a)
1.0 1.0
120
M=117.38 J/m
loading rate = 0.09MPa/s
M-integral (J/m)
100 80 M=98.36 J/m 60 40
M=97.93 J/m
20 0 0
2
6 4 Time step (s)
10
8
(b) M=117.02 J/m
Loading rate = 0.009MPa/s
120
M-integral (J/m)
100 80
M=98.36 J/m
60 40
M=98.17 J/m
20 0
00
20 20
40 40
80 60 60 80 Time step (s)
100 100
(c)
Fig. 4. Variation of M-integral against time steps for three loading rates (a) loading rate = 0.90 MPa/s, (b) loading rate = 0.09 MPa/s, and (c) loading rate = 0.009 MPa/s. STATUSXFEM Growth of crack caught by STATUSXFEM
The newly opened crack surface
Original crack length
Original crack length
(a)
(b)
Fig. 5. Crack growth evaluated by STATUSXFEM at time steps (a) t = 8.02 s and (b) t = 8.03 s (deformation scale factor = 10).
It is noteworthy that the curves of the time-dependent M-integral show several
jump points for each loading rate. Referring to the crack growth behavior caught by XFEM, this phenomenon can be explained reasonably. It is because that portion of the damage energy is applied to the building of the new crack surface which may induce a decrease of M-integral associated with the time dependence parameters. However, after the crack growth enters a steady stage, the value of the time-dependent M-integral will raise again owing to the increased damage degree of the viscoelastic plate. For instance, under a loading rate of 0.09 MPa/s, the crack growth behaviors represented by STATUSXFEM (status of XFEM element) in the field output are given in Fig. 5. The value of STATUSXFEM varies from 0 for elements in which damage has not initiated (e.g., no crack growth) to 1 for elements that have completely failed (e.g., no traction across the crack surface). It is clear that from loading time steps 8.02 s to 8.03 s, the new crack surface occurs. At this moment, the value of the time-dependent M-integral has a decrease of 9.47 J/m, that is, it varies from 97.93 J/m to 88.46 J/m. The results indicate that the rewritten time-dependent M-integral can be an effective fracture parameter to measure the damage degree of viscoelastic material induced by crack growth under varied loading rate conditions. 5.4 Propagation of collinear cracks in viscoelastic materials In this section, the proposed time-dependent M-integral of viscoelastic material is extended to the evaluation of the material damage induced by collinear cracks growth. As shown in Fig. 6, the essential information for numerical simulation such as the model geometry, load, and boundary conditions are presented. In detail, the dimensions of viscoelastic plate are 45 45 mm2 in length width with two 5 mm collinear cracks. Four-node plane-strain element meshing is considered. The uniform distributed pressure with a constant value of 0.83 MPa is applied to prove the crack growth over a period of 1 second. Similarly, the integral domain is also chosen as A1 = 34.0 and A2 = 38.0. The distance d shown in Fig. 6 will vary as 2 mm and 3 mm.
σ
Viscoelastic plate Cracks Equivalent
integral
d A1 A2
σ Fig. 6. Finite element mesh and loading conditions for a geometry with two cracks.
See Fig. 7, there are two sharp declines in the line with black symbols, which denotes the variation of time-dependent M-integral for a model with d = 2 mm. One occurs from load time steps 0.94 s to 0.95 s, and at this point the values of the time-dependent M-integral varies from 232.13 J/m to 209.40 J/m. The other is from 0.98 s to 0.99 s, corresponding M value of which is from 199.84 J/m to 198.18 J/m. Combining the variation of the M-integral and crack growth behavior caught by XFEM shown in Fig. 8, it is clear that the decrease of M-integral represents the appearance of the newly built crack surface. This result is consistent with the previous calculation of a model with one central crack. The line with red dots corresponds to the varied M-integral for model with d = 3 mm. The time-dependent M-integral increases along the entire time steps because no new crack surface appears.
250 d=2 mm d=3 mm
M-integral (J/m)
200 150 100 50 0 0.0
0.2
0.4
0.6
0.8
1.0
Time step (s) Fig. 7. Variation of time-dependent M-integral against time steps under = 0.83 MPa. STATUSXFEM
Crack growth
Original cracks
New crack surface
New crack surfaces
Fig. 8. Crack growth evaluated by STATUSXFEM from load time steps of t = 0.94 s to t = 0.95 s (deformation scale factor = 10).
Furthermore, different loading pressures are applied on model with two cracks (d = 2 mm) to see the damage mechanism in viscoelastic material with and without crack growth represented by the time-dependent parameter M-integral. In Fig. 9, the line with black symbols represents the M-integral versus the time steps for a pressure of 0.7 MPa. During this process, no crack growth occurs based on the results of STATUSXFEM. The M-integral has shown a monotonic increasing trend owing to the viscous deformation in material without crack growth. When the pressure increases to 0.8 MPa as shown in the line with red symbols, crack growth occurs but no new surface appears. Comparing the variation of the M-integral for = 0.7 MPa, it is found that the values of the M-integral for = 0.8 MPa are much higher owing to the crack growth behavior. That is, the damage degree induced by crack growth is more
obvious compared with the viscous deformation. For instance, at a time steps of t =1 s, the values of M-integral are 104.98 J/m and 136.84 J/m for = 0.7 MPa and = 0.8 MPa, respectively. However, at this moment, the crack begins to propagate but a new crack surface has not been formed, and thus the M-integral still increases with an increase of time steps.
160 s = 0.7 MPa s = 0.8 MPa
140
M-integral (J/m)
120 100 80 60 40 20 0 -20 0.0
0.2
0.4
0.6
0.8
1.0
Time step (s) Fig. 9. Variation of M-integral versus time steps for two loading pressures.
In summary, it can be concluded that the values of the time-dependent M-integral can effectively describe the damage degree of the viscoelastic material induced by both viscous deformation and crack growth. Moreover, the crack growth behavior has an extremely important influence on the damage of viscoelastic materials based on the variation of the time-dependent M-integral. The contribution of this work is to provide an effective and reliable tool for the damage evaluation of viscoelasticity based on the concept of M-integral method. Herein, the stress relaxation and the rate correlation as specific material properties of viscoelastic material are taken into account. The potential applications of the present method can be used in the damage assessment of viscosity materials.
6. Conclusions In this paper, the application of the M-integral is successfully extended to the
investigation of damage evaluation for viscoelastic material. A new form of time-dependent M-integral, which includes the variables of inelastic stress and viscous strain, has been developed to predict the damage behavior such as crack growth in material. Using Python scripts, several numerical examples are given to evaluate the validity of this newly proposed time-dependent M-integral in tensile fracture. The results show that the time-dependent M-integral can be an effective fracture parameter for the material damage induced by crack growth in viscoelastic material. For instance, the variations of time-dependent M-integral under different loading rates can directly describe the crack growth behavior in viscoelastic material during the whole load time steps. In addition, for a model with two collinear cracks, the tendency of the time-dependent M-integral is also in line with previous calculations, that is, the value of M decreases with the occurrence of a new crack surface, but the overall tendency appears to be upward. Furthermore, the present viscoelastic M-integral method is applicable to evaluating the damage degree of viscoelastic material e.g., solid propellant, Ni-alloy, dispersed fuel with multiple complex defects.
Acknowledge. This work was supported by the National Natural Science Foundation of China (Nos.11772245), the 111 Project (B18040), the Fundamental Research Funds for the Central Universities in China, the Natural Science Basic Research Plan in Shaanxi Province of China (Program No. 2018JC-004), and the Exploration Program-Q of Natural Science Foundation in Zhejiang Province of China (Program No. LQ20A020010). The author Qun Li gratefully acknowledges the support of K.C. Wong Education Foundation.
Author Statement Herewith we submit the revised manuscript “The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic
materials”, which has neither previously nor simultaneously in whole or in part been submitted anywhere else, to you for publication in Mechanics of Materials. We should very much appreciate if you would inform us your decision in your earlier convenience after you review it. You can use our e-mail address for this purpose: [email protected]
Conflict of Interest We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.
References Azmi MSM, Fujii T, Tohgo K, et al. On the Δ J-Integral to characterize elastic-plastic fatigue crack growth. Engineering Fracture Mechanics, 2017, 176:300-307. Banks-Sills L, Motola Y, Shemesh L. The M-integral for calculating intensity factors of an impermeable crack in a piezoelectric material. Engineering Fracture Mechanics, 2008, 75(5):901-925. Becker TH, Mostafavi M, Tait RB, et al. An approach to calculate the J‐integral by digital image correlation displacement field measurement. Fatigue & Fracture of Engineering Materials & Structures, 2012, 35(10): 971-984. Budiansky B, Rice JR. Conservation laws and energy-release rates. Journal of Applied Mechanics, 1973, 40:201. Chang JH, Kang YC, Chung LG. M- and Mc-integrals for multicracked problems in three
dimensions.
Journal
of
Engineering
Mechanics
Asce,
2013,
139(12):1874-1880. Chang JH, Wu WH. Using M-integral for multi-cracked problems subjected to nonconservative and nonuniform crack surface tractions. International Journal of Solids and Structures, 2011, 48(19):2605-2613. Chang JH, Lin JS. Surface energy for creation of multiple curved cracks in rubbery
materials. Journal of Applied Mechanics, 2007, 74(3):488. Chang JH, Pu LP. A finite element approach for J2 calculation in anisotropic materials. Computers & Structures, 1997, 62(4):635-641. Contino M, Andena L, Vincenzo La Valle, et al. A comparison between K and G approaches for a viscoelastic material: the case of environmental stress cracking of HDPE. Mechanics of Time-Dependent Materials, 2019, (18):1-14. Eshelby JD. The elastic energy-momentum tensor. Journal of Elasticity, 1975, 5(3-4):321-335. Fukahori Y, Liang H, Busfield JJC. Criteria for crack initiation during rubber abrasion. Wear, 2008, 265(3):387-395. Gooch JW. Maxwell Model. In: Gooch J.W. (eds) Encyclopedic Dictionary of Polymers. Springer, New York, NY, 2011. Herrmann AG, Herrmann G. On energy-release rates for a plane crack. ASME, Transactions, Journal of Applied Mechanics, 1981, 48:525-528. Hou JL, Zuo H, Li Q, et al. M-integral analysis for cracks in a viscoplastic material with extended finite element method. Engineering Fracture Mechanics, 2018:294-311. Judt PO, Ricoeur A. A new application of M- and L-integrals for the numerical loading analysis of two interacting cracks. ZAMM-Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 2016, 96(1):24-36. Knowles JK, Sternberg E. On a class of conservation laws in linearized and finite elastostatics. Archive for Rational Mechanics and Analysis, 1972, 44(3):187-211. Li Q, Guo YL, Hou JL, et al. The M-integral based failure description on elasto-plastic materials with defects under biaxial loading. Mechanics of Materials, 2017, 112:163-171. Li Q, Hu YF, Chen YH. On the physical interpretation of the M-integral in nonlinear elastic defect mechanics. International Journal of Damage Mechanics, 2013, 22(4):602-613. Li FZ, Shih CF, Needleman A. A comparison of methods for calculating energy
release rates. Engineering Fracture Mechanics, 1985, 21(2):405-421. Lv Z. The analysis of crack growth of solid propellant. PhD thesis, Harbin, China: Harbin Institute of Technology, 2013: 4-16. (In Chinese) El Yaagoubi M, Juhre D, Meier J, et al. Prediction of tearing energy in mode III for filled elastomers. Theoretical and Applied Fracture Mechanics, 2017, 88:31-38. Motola Y, Banks-Sills L. M-Integral for calculating intensity factors of cracked piezoelectric materials using the exact boundary conditions. Journal of Applied Mechanics, 2009, 76(1):011004. Nguyen TD, Govindjee S, Klein PA, et al. A material force method for inelastic fracture mechanics. Journal of the Mechanics and Physics of Solids, 2005, 53(1):91-121. Nikishkov GP, Atluri SN. Calculation of fracture mechanics parameters for an arbitrary three‐dimensional crack, by the ‘equivalent domain integral’ method. International
Journal
for
Numerical
Methods
in
Engineering,
1987,
24(9):1801-1821. Sciumè G, Benboudjema F. A viscoelastic unitary crack-opening strain tensor for crack width assessment in fractured concrete structures. Mechanics of Time-Dependent Materials, 2016, 21:223-243. Shih CF, Moran B, Nakamura T. Energy release rate along a three-dimensional crack front in a thermally stressed body. International Journal of Fracture, 1986, 30(2):79-102. Wang H, Lu W, Barber JR, et al. The roles of cohesive strength and toughness for crack growth in visco-elastic and creeping materials. Engineering Fracture Mechanics, 2016, 160:226-237. Yao W, Li X, Hu X. Viscoelastic crack analysis using symplectic analytical singular element combining with precise time-domain algorithm. International Journal of Fracture, 2018, 214:29-48. Yi J, Shen S, Muhunthan B, et al. Viscoelastic–plastic damage model for porous asphalt mixtures: Application to uniaxial compression and freeze–thaw damage. Mechanics of Materials, 2014, 70(1):67-75.
Yuval M, Giorgio O, Johannes TB, et al. Crack Initiation in Viscoelastic Materials. Physical Review Letters, 2018, 120(26):268002-1-6. Yu HJ, Kitamura T. A new domain-independent interaction integral for solving the stress intensity factors of the materials with complex thermo-mechanical interfaces. European Journal of Mechanics A/Solids, 2015, 49:500–509. Yu HJ, Wu LZ, Guo LC, Du SY, He QL. Investigation of mixed-mode stress intensity factors for nonhomogeneous materials using an interaction integral method. International Journal of Solids and Structures, 2009, 46:3710–3724. Yu NY, Li Q, Chen YH. Measurement of the M-integral for a hole in an aluminium plate or strip. Experimental mechanics, 2012, 52(7):855-863. Yu NY, Li Q. Failure theory via the concept of material configurational forces associated with the M-integral. International Journal of Solids and Structures, 2013, 50(25-26):4320-4332.
List of Table Captions Table 1 Material parameters of HTPB (Lv, 2013) Parameters
Values
Density (kg/mm3)
1.75106
E0 (MPa)
12.32
v
0.50
Table 2 Values of the M-integral for six equivalent domains No.
1
2
3
Dimensions (mm)
A1=39 A2=41
A1=34 A2=37
Area (mm2) M-integral (J/m)
40
53.25
66
46.47
48.33
47.68
4
5
6
A1=34 A2=40
A1=32 A2=40
72
111
144
47.13
47.73
47.23
A1=36 A1=34 A2=39.50 A2=38
The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic materials Junling Hou , Chao Zhang , Qun Li PII: DOI: Reference:
S0167-6636(19)30948-2 https://doi.org/10.1016/j.mechmat.2020.103363 MECMAT 103363
To appear in:
Mechanics of Materials
Received date: Revised date: Accepted date:
31 October 2019 11 February 2020 13 February 2020
Please cite this article as: Junling Hou , Chao Zhang , Qun Li , The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic materials, Mechanics of Materials (2020), doi: https://doi.org/10.1016/j.mechmat.2020.103363
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2020 Published by Elsevier Ltd.
Highlights
A new form of M-integral associated with time dependence parameters is presented for viscoelastic materials. Time-dependent M-integral is numerically implemented based on the equivalent domain integral method as an effective and accepted fracture mechanical parameter for damage induced by crack growth in viscoelastic materials. Numerical examples are given to demonstrate the validity and relevance of the time-dependent M-integral in viscoelastic material.
Unmarked The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic materials Junling Hou, Chao Zhang, Qun Li* * Corresponding author, email: [email protected] State Key Laboratory for Strength and Vibration of Mechanical Structures, Shaanxi Key Laboratory of Environment and Control for Flight Vehicle, Shaanxi Engineering Research Center of Nondestructive Testing and Structural Integrity Evaluation, School of Aerospace, Xi’an Jiaotong University, Xi’an 710049, China
Abstract A new form of M-integral associated with time dependence parameters, is presented herein for viscoelastic materials. Based on the equivalent domain integral method, this time-dependent M-integral is numerically implemented as an effective and accepted fracture mechanical parameter for damage induced by crack growth in viscoelastic materials. Based on the linear viscoelastic model defined through Prony series, the conservation of the time-dependent M-integral for viscoelasticity is verified by applying user defined Python scripts. The results show that the newly proposed time-dependent M-integral can be successfully calculated numerically. Furthermore, numerical examples are given to demonstrate the validity and relevance of the time-dependent M-integral in viscoelastic material. In particular, the variations of the time-dependent M-integral for different loading rates are considered, which show that the crack growth behavior over a period of time in viscoelastic material can be evaluated based on the value of the time-dependent M-integral. In addition, the time-dependent M-integral is calculated to assess the damage degree induced by two collinear cracks.
Keywords: Viscoelastic material; M-integral; Time-dependent; Domain integral method; Crack growth.
1. Introduction Owing to the widespread use of time-dependent viscoelastic materials such as solid propellant, Ni-alloy, dispersed fuel, bituminous concrete, and polymer matrix composite etc., it is necessary to have a better understanding of their fracture mechanism for achieving a material design and engineering application. Considering the influence of the manufacturing process, loading condition, and temperature change, it has been found that the crack initiation, evolution and propagation may result in a failure in viscoelastic materials (Yuval et al., 2018; Mohammed et al., 2017; Yi et al., 2014). Thus, finding an effective method becomes critically important for the evaluation of the damage degree and evolution mechanism in viscoelastic materials during the fracture process. For instance, in the framework of linear elastic fracture mechanics, two fracture parameters (stress intensity factor and energy release rate) are employed to investigate the growth and resistance of cracking in viscoelastic materials (Contino et al., 2019). Using precise time-domain expanding algorithm, the viscoelastic cracked problem is solved using converted quasi-elastic variables, such as displacement and stress (Yao et al., 2018). A viscoelastic damage model is proposed for the calculation of crack width in viscoelastic concrete (Sciumè and Benboudjema, 2016). The maximum tensile stress is introduced as a fracture criterion to deal with the crack initiation process in rubber material with viscous-elastic properties (Fukahori et al., 2008). Studies have shown that it is of practical importance to achieve a deeper understanding of the ‘global’ damage state from defect formation, expansion or other factors in viscoelastic materials. That is, for multi-cracked viscoelastic materials, more effective and universal methods are still required. In recent years, the conservation integral method has attracted much attention owing to its advantages in dealing with the complicated damage problems found in different materials (Wang et al., 2016; Yu and Kitamura, 2015; Yu et al., 2009). For the crack problem in both homogeneous and non-homogeneous materials, the domain independent interaction integral method has been successfully used in the calculation of the stress intensity factor (Yu and Kitamura, 2015; Yu et al., 2009). In addition, the energy-conservation integrals, Jk-, M- and L- integrals, have been
proved to be reliable and effective in assessing the fracture properties of materials (Azmi et al., 2017; Li et al., 2017; Motola and Banks-Sills, 2009). In particular, the M-integral, as one of the key concepts in material configurational mechanics, has shown its advantage in quantifying the damage degree for multi-crack problems since it is inherently related to the change of the total potential energy in a damaged system independent of the specific damage details of every defect (Li et al., 2017; Motola and Banks-Sills, 2009). In other words, its value is directly related with all fracture factors such as the material constants, defect characteristic, and mechanical loading. Thus, this approach is much easier in practice and has been widely applied to elastic, elastic-plastic and visco-elastic-plastic problems. Numerous studies have been carried out using the M-integral to evaluate the material damage evolution (Chang and Lin, 2007; Yu et al., 2012; Chang and Wu, 2011; Yu and Li, 2013; Chang et al., 2013; Judt and Ricoeur, 2016). For instance, Chang and Lin (2007) presented a fracture parameter called Mc-integral to investigate the multiple curved crack problem in rubbery materials under large deformation. Using digital image correlation method, the M-integral was measured by Yu et al. (2012) in elastic-plastic samples created using LY12 aluminum. This technique can be extended to evaluate the damage induced by multiple defects in materials. Furthermore, according to a study by Chang and Wu (2011), the energy related M-integral is considered as a damage parameter to describe the irreversible fracture problem in a solid weakened by multi-cracks. Based on the material configurational force theory, Yu and Li (2013) developed the application of M-integral in elastic-plastic material with multi-defects, such as cracks and voids. For a three-dimensional problem, Chang et al. (2013) propose a new energy parameter called Mc-integral for the evaluation of the structural integrity in elastic media with multiple cracks. For engineering applications, the M-integral was introduced by Judt and Ricoeur (2016) to calculate the loading quantities of engineering structures with two cracks. Based on above studies, it is of practical importance to introduce the implementation of M-integral for the assessment of ‘global’ damage state from defect
formation, expansion in materials with a viscous effect. Particularly, the effective and universal methods via the M-integral are necessarily extended for the cracked viscoelastic materials, in which the stress relaxation and the rate correlation as specific material properties of viscoelastic material are taken into account. In contrast with the traditional fracture parameters, for example, stress intensity factor K (only valid in cracking problems of elastic materials) and Jk-integral (can be extended for complicated materials, such as Ct -integral for time-dependent creep behavior, but not feasible for multiple cracks due to their ‘local’ nature associated with a single crack tip). The values of M-integral is directly related with all fracture factors such as material constants, defect characteristic, mechanical loading and etc. The novelty of this paper is to extend the concept of M-integral to deal with the failure of viscoelastic materials. In addition, although the M-integral is a path-independent integral parameter used to measure the damage degree of materials, it is cumbersome to compute in finite element calculation. For numerical implementation, alternative approaches are needed to simplify the numerical computation code. Among them, a numerical technique termed the equivalent domain integral method is recommended. For the equivalent domain method, all necessary variables for numerical computation are easily attainable in a finite element analysis. There is no need to consider the influence of crack tip singularity, or the difficulty in catching the Gauss point on a specific path. Thus, the EDI method appears to be a powerful and attractive tool for transforming the line integral into an equivalent domain integral. For instance, it has been introduced to evaluate the values of conservation line integral in elastic-plastic materials (Li et al., 1985), thermo-elastic materials (Shih et al., 1986), and linear anisotropic elastic materials (Chang and Pu, 1997). The concept of domain integral method is also extended to calculate the value of the M-integral for the measurement of nickel-based superalloys with multi-cracks (Hou et al., 2018). However, in previous study, the definition of M-integral is built in the elastic media. In addition, to evaluate the damage degree of viscoplastic material, the geometry model for finite element calculation is consists of two parts, one of which is the damage zone defined
as viscoplastic media, and the other is a pure elastic zone that surrounds the damage zone. Considering that the M-integral in other fields has been compared with experimental data (Li et al., 2017), numerical results (Hou et al., 2018) in our previous work, it shows that the M-integral is applicable for the validation of the cracked problems in materials. In present study, the work of the time-dependent M-integral is only presented by the theoretical explanations and numerical simulations, without comparison with experimental data. The main purpose of the present paper is to introduce the concept of M-integral into the viscoelastic material owing to its advantage in describing the material fracture behavior under the presence of cracks. In this paper, a computationally efficient and convenient method for the evaluation of the M-integral for defects in two dimensional viscoelastic material is presented. The remainder of this paper is structured as follows. Firstly, the linear viscoelastic model involved in present study is given in Section 2. Then, the derivation of M-integral for viscoelastic material is carried out in Section 3. Next, for numerical implementation, the procedure of converting line M-integral into domain integral is presented in detail. With the newly built viscoelastic domain M-integral, several numerical examples are provided and discussed in Section 5 to illustrate the availability of the approach.
2. Viscoelastic material model The classical viscoelastic model and its application in finite element method are presented in this section. For viscoelastic media, it is known that the stress distribution and material deformation are both time-dependent variables. Specifically, the viscous deformation produces an irreversible flow with the external force work, and all of its energy is converted into internal energy. For linear viscoelastic material, its stress response can be basically described through various combinations of two types of mechanical model. Two classical linear viscoelastic models have been proposed, namely, Maxwell model (related with the stress relaxation) and Kelvin model (related with creep). In the present study, the time domain viscoelasticity defined in ABAQUS is considered under a stress relaxation load, and thus the
Maxwell model is chosen. The basic linear viscoelastic response of the Maxwell model (Gooch, 2011) is described based on the sum of two mechanical models. One is the pure elastic model which obeys Hooke’s law, and the other is a viscous component satisfying Newton’s viscosity law. Under the loading condition (t), the definition formation of the time varying strain can be written as follows:
=
E
(1)
Here, =d / dt represents the strain rate, E is the elastic modulus, and is viscous coefficient. In addition, Eq. (1) is the differential constitutive equation of the Maxwell model. To describe more complicated stress states in viscoelasticity, a generalized Maxwell model (Gooch, 2011) consisting of multiple Maxwell models is proposed. For a two dimensional problem, the mechanical behavior can be illustrated by considering a sudden strain 0 applied to the viscoelastic media and held constant for a long time. That is, considering a varying strain (t) = 0 H(t) applied to the viscoelastic material, the stress response is obtained as follows: n
(t ) 0 Ei e
t Ti
(2)
i 1
Here, H (t) denotes the unit step function, which is defined as follows:
0, t 0 H (t ) 1, t 0
(3)
Here, Ti = i / E corresponds to the relaxation time of i viscous component, while i is the viscosity tensor. And the stress relaxation modulus of the generalized Maxwell model can be written as below (Gooch, 2011) n
E (t ) Ei e
t Ti
(4)
i 1
Eq. (4) is the Prony series form of time-dependent relaxation modulus for viscoelastic material in which the parameters can be directly defined in ABAQUS material.
Otherwise, in finite element analysis, to begin with, an elastic response of the viscoelasticity should be defined. In this paper, the small-strain deformation condition is considered and thus the rate-independent elastic response is defined using a linear elastic material model.
3. Definition of M-integral in viscoelastic material Within the framework of the material configurational force theory (Eshelby, 1975; Li et al., 2013), the physical interpretation of the M-integral is explained as the potential energy change during the process of the self-similar expansion of defects. Thus, the stability of materials with multiple defects, that is, the damage degree of materials can be characterized by the value of the invariant M-integral. The M-integral has recently shown its advantages as an energy fracture parameter in evaluating the degradation of different materials (e. g., elastic-plastic, piezoelectric and general anisotropic materials) with multi-defects (Yu and Li, 2013; Bank-Sill et al., 2008). In the present study, a new form of M-integral for viscoelastic material is defined by considering the crack growth behavior. Considering two-dimensional media with cracks, regardless of the inertia terms and body force in the system, the formula of the M-integral can be written by (Knowles and Sternberg, 1972; Budiansky and Rice, 1973; Herrmann and Herrmann, 1981) M
( w
x jk uk ,i xi ) n j ds
ij i
M
j
n j ds
(i, j , k =1, 2)
(5)
It should be noted that the present study is applicable for two dimensional analyses. For three dimensional conditions, the M3D-integral is written as follows:
1 wij xi jk uk ,i xi jk uk n j da A 2
M3D
(6)
In which, is the total volume surrounded by domain A containing all defects. It is obvious that an additional term 1 jk n j uk must be considered. 2
Based on the concepts of material configurational force theory, the physical meaning of Mj is explained as the material damage extension configuration force,
which is shown as follows
M j w ij xi kj uk ,i xi
(7)
where xi corresponds to the coordinate, ij is the Kronecker symbol, σij refers to the Cauchy stress tensor, uk denotes the displacement tensor, and the subscript prima {},i denotes the corresponding differentiation with respect to the coordinate xi. All cracks are surrounded by the enclosed path with the outer normal vector n j. For viscoelastic material, the strain energy density w(t) can be defined as t
w(t ) ij d ij (t ) ij ij dt ij
0
(8)
Here, εij (t) represents the strain tensor which is related with the load time t. Assuming that the total strain tensor ε could be expressed as follows: (Simo and Hughes, 1998, Chapter 10; Nguyen et al., 2005) ε εe εv
(9)
Here, εe and εv are the elastic and viscous strain tensors, respectively. In the same manner, the following additive split of the stress response tensor is given in the following: σ σ E σ IE
(10)
Here, σE denotes the elastic and equilibrium stress behavior, while σIE corresponds to the time-dependent and inelastic behavior in viscoelasticity. Combining Eqs. (9) and (10), the form of Eq. (8) can be rewritten as follows: w ε σ : ε = σ E σ IE : ε ε v ε v = σ E : ε σ IE : ε ε v σ IE : ε v
(11)
= σ : ε v IE
Here, Ψ (ε, εv) is treated as the free energy density. Then, introducing the equivalent relation in Eq. (11) to Eq. (5), the rate form of the M-integral can be obtained, which is shown as follows u ) x (σ IE : ε v ) x ]n ds M [( ij jk k ,i i ij i j
(12)
Rewriting Eq. (12), the new definition of the time dependent M-integral for viscoelasticity is given as
M (t )
b (t ) x w ij
i
IE
(t ) ij xi n j ds = M Vj n j ds
(13)
where bij (t) denotes the Eshelbian energy momentum for the viscoelastic material, and is defined as follows:
bij (t ) = ij jk uk ,i
(14)
The symbol wIE (t ) is applied to characterize the part of strain energy density induced by the inelastic stress response and viscous strain deformation, that is, wIE (t )=σ IE : ε v . That is, for viscoelastic material, the material damage extension
configuration force M Vj is expressed as follows:
M Vj (t )=bij (t ) xi wIE (t ) ij xi
(15)
Referring to the derivation of this newly built M-integral for viscoelastic material, it is clear that the value of the M-integral is calculated along a path surrounding all the defects but without consideration of the details, such as how many cracks or the shapes of the cracks. Thus, this time-dependent M-integral method can be extended to solve other defects problems such as multi-cracks, interface debonding in viscoelasticity.
4. FEM calculation of M-integral in viscoelastic material
Γs
Viscoelastic medium
dA
Γ
Γ0 Cracks
A
Fig. 1. Description of the equivalent domain in viscoelasticity.
For numerical implementation, the time-dependent M-integral obtained in Eq. (13) can be transformed into a form of a domain integral by introducing the equivalent domain integral method. Initially, a weight function q is introduced for mathematical
manipulation, which represents arbitrary functions that smoothly go through the integral region A. Area A corresponds to the shadow region shown in Fig. 1. Here, the q-function is simply a transition function formed as q (x1, x2) (Becker et al., 2012; Niikishkov and Atluri, 1987), and it has been indicated that the choice of q-function is not sensitive to the results of the M-integral. Thus, the values of q-function are set as 0 and 1 for the outer () and inner (0) contours, respectively. We thus have the following:
M (t ) 1 M Vj n j ds 0 M Vj n j ds 0
V j
M n j qds M Vj n j qds 0
(16)
=
0
M Vj q n j ds
The following identity integral relationship is introduced,
M q A
V j
,j
dA
M A
V j, j
q M Vj q, j dA
(17)
Using the divergence theorem on the left side of Eq. (17) and combining Eq. (16), we obtain the following:
0
M Vj q n j ds = M Vj q dA M Vj , j q M Vj q, j dA A
A
,j
(18)
Considering the free of singularity in zone A, the relation M Vj, j =0 is obtained, and then Eq. (18) can be rewritten as M (t )= M Vj q, j dA A
q = bij (t ) xi wIE (t ) ij xi dA A x
(19)
j
For the convenience of calculation, Eq. (19) can be separated as the superposition of two parts:
M (t ) M I (t ) M II (t )
(20)
with
M I (t ) bij (t ) xi A
q dA x j
(21)
M II (t ) wIE (t ) A
q dA x j
(22)
The form of function q in zone A can be obtained by the shape function interpolation method, which is shown below: 4
q N1q1e N 2 q2e N 3 q3e N 4 q4e Nq e N i ( r , s ) qie
(23)
i 1
where Ni represents the shape function which is only associated with the coordinates of each element, and (r, s) is the position of the corresponding Gauss point. In the present study, four Gauss points are used for the integration in entire equivalent domain. For a single element, rewrite Eqs. (21) and (22) as follows:
M Ie
1
1
I r, s drds
1 1 I
M IIe
1
1
I
1 1 II
r , s drds
(24) (25)
Applying the Gauss integral method, Eq. (24) is approximately equivalent to
M Ie I I r1 , s1 I I r2 , s2 I I r3 , s3 I I r4 , s4
(26)
u u v v q x x y xy x xy x x x y x x x x e I I (r , s ) det( J ) u u v v q y xy x x xy y y y x x y x x y
(27)
with
where det(Je) is the determinant of the Jacobi matrix for the specific element. Similarly, the value of M IIe shown in Eq. (25) equals
M IIe I II r1 , s1 I II r2 , s2 I II r3 , s3 I II r4 , s4
(28)
q q I II (r , s ) xIE xv x1 xyIE xyv x1 xyIE xyv x2 yIE yv x2 det(J e ) x y
(29)
with
In particular, because the integral domain A shown in Fig. 1 is far from the cracked area, the calculation of the M-integral is not sensitive with the crack tip
behavior. Thus, singular mesh is without consideration for crack tip region. For the calculation of Eqs. (27) and (29), the values of ij,
ui q , and can x j x j
be obtained directly from the numerical simulation in ABAQUS postprocessing. Otherwise, the values of ,
IEe and ijv will be calculated from the identity
relation shown in this section. Then, the value of Me in a single element can be obtained by superposing the calculated M Ie and M IIe . Repeat above procedure for every element in zone A, the sum of Me is equal to the corresponding value of the M-integral for chosen integral contour. And with the user defined Python script employed in ABAQUS, the numerical calculation of M-integral can be realized in viscoelastic material.
5. Results and discussions In this section, several numerical examples are presented to show the validity of the newly proposed M-integral for viscoelastic material. For numerical simulation, the crack growth behavior is evaluated using the extended finite element method (XFEM). Differing from the traditional finite element, there is no need to do remeshing for the XFEM. Thus, the mesh continuity can prove the effective transfer of stress and strain. Upon request of the XFEM, the maximum principle stress criterion is set as the damage initiation law to control the behavior of crack growth (max = 1.91 MPa) (Lv, 2013), and the damage evolution law is linear softening stress-strain response whose type is based on displacement, setting u max = 0.25 mm at failure (Lv, 2013). Combining the built-in viscoelastic material model in ABAQUS, numerical calculations of the viscoelastic solid propellant with crack growth could be carried out. The computational cost is between 5~10 minutes for most of the examples simulated using double Intel core i7 CPU (i7-8700 CPU @ 3.2 GHz 3.19GHz). Specifically, for the model containing one single crack described in Section 5.2 and 5.3, a job in ABAQUS could be completed within 5 minutes, and the Python script calculating the visco-elastic integrals can be completed within 5-15 minutes according to the
numbers of the outputting M-integral for different time steps.
5.1 Material properties of viscoelastic solid propellant The viscoelastic material chosen for numerical verification is a common type of solid propellant called hydroxyl-terminated polybutadiene (HTPB), which has been widely used as the fuel in solid rocket engine and shows linear viscoelastic properties under sustained loading. Considering the process of manufacturing, transport or application, it is inescapable that some defects such as cracks will occur in solid propellant owing to the complex load conditions. Research has shown that the damage evolution (e.g. crack growth) in solid propellants will likely lead to a catastrophic accident owing to incomplete combustion. Thus, it is valuable to study the damage mechanism of such viscoelastic solid propellant. Referring to the viscoelastic model presented in Section 2, the Prony series form of time-dependent relaxation modulus is given as follows (Lv, 2013):
E (t ) 4.10 (0.15 0.12e
t
0.04
0.17e
t
0.40
0.25e
t
4
0.32e
t
40
)
(30)
Here, the unit of time t in Eq. (30) is second. In addition, the initial material parameters of viscoelastic solid propellant (HTPB) are shown in Table 1 for completeness.
5.2 Numerical verification of time-dependent M-integral conservation σ
Viscoelastic plate Crack Equivalent a
A1 A2
σ
integral
Fig. 2. Sketch of viscoelastic plate with central crack and equivalent integral domain.
To start with, the validity of a Python script for time-dependent M-integral is examined based on the integral conservation. As shown in Fig. 2, under plane strain condition, the geometry with a central single crack subjected to uniform tension load σ is investigated. The 2D plate has a dimension of 45 mm × 45 mm (length × width) with a central crack (a = 10 mm). The finite element mesh used in this study consists of four-node plane strain elements with full integration in this paper. The applied tensile stress σ is 0.50 MPa over a period of 1 second. In Table 2, the values of time-dependent M-integral for six different areas of equivalent integral domain in viscoelastic material are given. Specifically, the dimensions of six areas vary from 40 to 144 mm2, while the range of length and width are (32-39 mm) and (37-41 mm), respectively. It can be seen from Table 2 that the maximum value of the time-dependent M-integral is 48.33 J/m which is obtained in domain No. 2, whereas the minimum value is 46.47 J/m for domain No. 1. The comparison results show that the error between them (i.e., 4.00%) is controlled at an acceptable level, which means that the results of time-dependent M-integral for different equivalent integral domains show a good conservation nature, which verifies that the numerical implementation of the time-dependent M-integral has been successfully carried out. 5.3 Effects of loading rates on M-integral in viscoelastic material Considering the time dependence of the newly built M-integral in viscoelastic material, the variations of the M-integral under three different loading rates (i.e., 0.90 MPa/s, 0.09 MPa/s, 0.009 MPa/s) are simulated in this section. The model geometry is shown in Fig. 2, that is, a 2D plate with one central crack. The load condition is a purely monotonic increasing loading with a peak value of 0.90 MPa over three time periods (i.e., 1 s, 10 s, 100 s) as shown in Fig. 3, while the time increment is set to the same value as 0.01 s for all simulations. The dimensions of equivalent integral domain are chosen as A1=34.0 and A2=38.0. The curves of the computed time-dependent M-integral under three different
loading rates (Fig. 4) show that the general tendency of time-dependent M-integral increases with increasing load time steps for each loading rate, which confirms that the damage degree of viscoelastic plate increases with the crack growth. For instance, at end of the load time steps, the values of time-dependent M-integral are 113.43 J/m, 110.19 J/m, 110.02 J/m corresponding to the three loading rates as marked in partial enlarged drawing in Fig. 4. (MPa) 0.9
1
0
10
t(s)
100
Fig. 3 Load conditions of different loading rates. 140
M=121.61 J/m
Loading rate = 0.9 MPa/s
M-integral (J/m)
120 100 80
M=99.85 J/m
60 40
M=101.43 J/m
20 0 0.0
0.2
0.4 0.6 0.8 Time step (s) (a)
1.0 1.0
120
M=117.38 J/m
loading rate = 0.09MPa/s
M-integral (J/m)
100 80 M=98.36 J/m 60 40
M=97.93 J/m
20 0 0
2
6 4 Time step (s)
10
8
(b) M=117.02 J/m
Loading rate = 0.009MPa/s
120
M-integral (J/m)
100 80
M=98.36 J/m
60 40
M=98.17 J/m
20 0
00
20 20
40 40
80 60 60 80 Time step (s)
100 100
(c)
Fig. 4. Variation of M-integral against time steps for three loading rates (a) loading rate = 0.90 MPa/s, (b) loading rate = 0.09 MPa/s, and (c) loading rate = 0.009 MPa/s. STATUSXFEM Growth of crack caught by STATUSXFEM
The newly opened crack surface
Original crack length
Original crack length
(a)
(b)
Fig. 5. Crack growth evaluated by STATUSXFEM at time steps (a) t = 8.02 s and (b) t = 8.03 s (deformation scale factor = 10).
It is noteworthy that the curves of the time-dependent M-integral show several
jump points for each loading rate. Referring to the crack growth behavior caught by XFEM, this phenomenon can be explained reasonably. It is because that portion of the damage energy is applied to the building of the new crack surface which may induce a decrease of M-integral associated with the time dependence parameters. However, after the crack growth enters a steady stage, the value of the time-dependent M-integral will raise again owing to the increased damage degree of the viscoelastic plate. For instance, under a loading rate of 0.09 MPa/s, the crack growth behaviors represented by STATUSXFEM (status of XFEM element) in the field output are given in Fig. 5. The value of STATUSXFEM varies from 0 for elements in which damage has not initiated (e.g., no crack growth) to 1 for elements that have completely failed (e.g., no traction across the crack surface). It is clear that from loading time steps 8.02 s to 8.03 s, the new crack surface occurs. At this moment, the value of the time-dependent M-integral has a decrease of 9.47 J/m, that is, it varies from 97.93 J/m to 88.46 J/m. The results indicate that the rewritten time-dependent M-integral can be an effective fracture parameter to measure the damage degree of viscoelastic material induced by crack growth under varied loading rate conditions. 5.4 Propagation of collinear cracks in viscoelastic materials In this section, the proposed time-dependent M-integral of viscoelastic material is extended to the evaluation of the material damage induced by collinear cracks growth. As shown in Fig. 6, the essential information for numerical simulation such as the model geometry, load, and boundary conditions are presented. In detail, the dimensions of viscoelastic plate are 45 45 mm2 in length width with two 5 mm collinear cracks. Four-node plane-strain element meshing is considered. The uniform distributed pressure with a constant value of 0.83 MPa is applied to prove the crack growth over a period of 1 second. Similarly, the integral domain is also chosen as A1 = 34.0 and A2 = 38.0. The distance d shown in Fig. 6 will vary as 2 mm and 3 mm.
σ
Viscoelastic plate Cracks Equivalent
integral
d A1 A2
σ Fig. 6. Finite element mesh and loading conditions for a geometry with two cracks.
See Fig. 7, there are two sharp declines in the line with black symbols, which denotes the variation of time-dependent M-integral for a model with d = 2 mm. One occurs from load time steps 0.94 s to 0.95 s, and at this point the values of the time-dependent M-integral varies from 232.13 J/m to 209.40 J/m. The other is from 0.98 s to 0.99 s, corresponding M value of which is from 199.84 J/m to 198.18 J/m. Combining the variation of the M-integral and crack growth behavior caught by XFEM shown in Fig. 8, it is clear that the decrease of M-integral represents the appearance of the newly built crack surface. This result is consistent with the previous calculation of a model with one central crack. The line with red dots corresponds to the varied M-integral for model with d = 3 mm. The time-dependent M-integral increases along the entire time steps because no new crack surface appears.
250 d=2 mm d=3 mm
M-integral (J/m)
200 150 100 50 0 0.0
0.2
0.4
0.6
0.8
1.0
Time step (s) Fig. 7. Variation of time-dependent M-integral against time steps under = 0.83 MPa. STATUSXFEM
Crack growth
Original cracks
New crack surface
New crack surfaces
Fig. 8. Crack growth evaluated by STATUSXFEM from load time steps of t = 0.94 s to t = 0.95 s (deformation scale factor = 10).
Furthermore, different loading pressures are applied on model with two cracks (d = 2 mm) to see the damage mechanism in viscoelastic material with and without crack growth represented by the time-dependent parameter M-integral. In Fig. 9, the line with black symbols represents the M-integral versus the time steps for a pressure of 0.7 MPa. During this process, no crack growth occurs based on the results of STATUSXFEM. The M-integral has shown a monotonic increasing trend owing to the viscous deformation in material without crack growth. When the pressure increases to 0.8 MPa as shown in the line with red symbols, crack growth occurs but no new surface appears. Comparing the variation of the M-integral for = 0.7 MPa, it is found that the values of the M-integral for = 0.8 MPa are much higher owing to the crack growth behavior. That is, the damage degree induced by crack growth is more
obvious compared with the viscous deformation. For instance, at a time steps of t =1 s, the values of M-integral are 104.98 J/m and 136.84 J/m for = 0.7 MPa and = 0.8 MPa, respectively. However, at this moment, the crack begins to propagate but a new crack surface has not been formed, and thus the M-integral still increases with an increase of time steps.
160 s = 0.7 MPa s = 0.8 MPa
140
M-integral (J/m)
120 100 80 60 40 20 0 -20 0.0
0.2
0.4
0.6
0.8
1.0
Time step (s) Fig. 9. Variation of M-integral versus time steps for two loading pressures.
In summary, it can be concluded that the values of the time-dependent M-integral can effectively describe the damage degree of the viscoelastic material induced by both viscous deformation and crack growth. Moreover, the crack growth behavior has an extremely important influence on the damage of viscoelastic materials based on the variation of the time-dependent M-integral. The contribution of this work is to provide an effective and reliable tool for the damage evaluation of viscoelasticity based on the concept of M-integral method. Herein, the stress relaxation and the rate correlation as specific material properties of viscoelastic material are taken into account. The potential applications of the present method can be used in the damage assessment of viscosity materials.
6. Conclusions In this paper, the application of the M-integral is successfully extended to the
investigation of damage evaluation for viscoelastic material. A new form of time-dependent M-integral, which includes the variables of inelastic stress and viscous strain, has been developed to predict the damage behavior such as crack growth in material. Using Python scripts, several numerical examples are given to evaluate the validity of this newly proposed time-dependent M-integral in tensile fracture. The results show that the time-dependent M-integral can be an effective fracture parameter for the material damage induced by crack growth in viscoelastic material. For instance, the variations of time-dependent M-integral under different loading rates can directly describe the crack growth behavior in viscoelastic material during the whole load time steps. In addition, for a model with two collinear cracks, the tendency of the time-dependent M-integral is also in line with previous calculations, that is, the value of M decreases with the occurrence of a new crack surface, but the overall tendency appears to be upward. Furthermore, the present viscoelastic M-integral method is applicable to evaluating the damage degree of viscoelastic material e.g., solid propellant, Ni-alloy, dispersed fuel with multiple complex defects.
Acknowledge. This work was supported by the National Natural Science Foundation of China (Nos.11772245), the 111 Project (B18040), the Fundamental Research Funds for the Central Universities in China, the Natural Science Basic Research Plan in Shaanxi Province of China (Program No. 2018JC-004), and the Exploration Program-Q of Natural Science Foundation in Zhejiang Province of China (Program No. LQ20A020010). The author Qun Li gratefully acknowledges the support of K.C. Wong Education Foundation.
Author Statement Herewith we submit the revised manuscript “The concept and numerical evaluation of M-integral based on domain integral method in cracked viscoelastic
materials”, which has neither previously nor simultaneously in whole or in part been submitted anywhere else, to you for publication in Mechanics of Materials. We should very much appreciate if you would inform us your decision in your earlier convenience after you review it. You can use our e-mail address for this purpose: [email protected]
Conflict of Interest We declare that we have no financial and personal relationships with other people or organizations that can inappropriately influence our work, there is no professional or other personal interest of any nature or kind in any product, service and/or company that could be construed as influencing the position presented in, or the review of, the manuscript entitled.
References Azmi MSM, Fujii T, Tohgo K, et al. On the Δ J-Integral to characterize elastic-plastic fatigue crack growth. Engineering Fracture Mechanics, 2017, 176:300-307. Banks-Sills L, Motola Y, Shemesh L. The M-integral for calculating intensity factors of an impermeable crack in a piezoelectric material. Engineering Fracture Mechanics, 2008, 75(5):901-925. Becker TH, Mostafavi M, Tait RB, et al. An approach to calculate the J‐integral by digital image correlation displacement field measurement. Fatigue & Fracture of Engineering Materials & Structures, 2012, 35(10): 971-984. Budiansky B, Rice JR. Conservation laws and energy-release rates. Journal of Applied Mechanics, 1973, 40:201. Chang JH, Kang YC, Chung LG. M- and Mc-integrals for multicracked problems in three
dimensions.
Journal
of
Engineering
Mechanics
Asce,
2013,
139(12):1874-1880. Chang JH, Wu WH. Using M-integral for multi-cracked problems subjected to nonconservative and nonuniform crack surface tractions. International Journal of Solids and Structures, 2011, 48(19):2605-2613. Chang JH, Lin JS. Surface energy for creation of multiple curved cracks in rubbery
materials. Journal of Applied Mechanics, 2007, 74(3):488. Chang JH, Pu LP. A finite element approach for J2 calculation in anisotropic materials. Computers & Structures, 1997, 62(4):635-641. Contino M, Andena L, Vincenzo La Valle, et al. A comparison between K and G approaches for a viscoelastic material: the case of environmental stress cracking of HDPE. Mechanics of Time-Dependent Materials, 2019, (18):1-14. Eshelby JD. The elastic energy-momentum tensor. Journal of Elasticity, 1975, 5(3-4):321-335. Fukahori Y, Liang H, Busfield JJC. Criteria for crack initiation during rubber abrasion. Wear, 2008, 265(3):387-395. Gooch JW. Maxwell Model. In: Gooch J.W. (eds) Encyclopedic Dictionary of Polymers. Springer, New York, NY, 2011. Herrmann AG, Herrmann G. On energy-release rates for a plane crack. ASME, Transactions, Journal of Applied Mechanics, 1981, 48:525-528. Hou JL, Zuo H, Li Q, et al. M-integral analysis for cracks in a viscoplastic material with extended finite element method. Engineering Fracture Mechanics, 2018:294-311. Judt PO, Ricoeur A. A new application of M- and L-integrals for the numerical loading analysis of two interacting cracks. ZAMM-Journal of Applied Mathematics and Mechanics / Zeitschrift für Angewandte Mathematik und Mechanik, 2016, 96(1):24-36. Knowles JK, Sternberg E. On a class of conservation laws in linearized and finite elastostatics. Archive for Rational Mechanics and Analysis, 1972, 44(3):187-211. Li Q, Guo YL, Hou JL, et al. The M-integral based failure description on elasto-plastic materials with defects under biaxial loading. Mechanics of Materials, 2017, 112:163-171. Li Q, Hu YF, Chen YH. On the physical interpretation of the M-integral in nonlinear elastic defect mechanics. International Journal of Damage Mechanics, 2013, 22(4):602-613. Li FZ, Shih CF, Needleman A. A comparison of methods for calculating energy
release rates. Engineering Fracture Mechanics, 1985, 21(2):405-421. Lv Z. The analysis of crack growth of solid propellant. PhD thesis, Harbin, China: Harbin Institute of Technology, 2013: 4-16. (In Chinese) El Yaagoubi M, Juhre D, Meier J, et al. Prediction of tearing energy in mode III for filled elastomers. Theoretical and Applied Fracture Mechanics, 2017, 88:31-38. Motola Y, Banks-Sills L. M-Integral for calculating intensity factors of cracked piezoelectric materials using the exact boundary conditions. Journal of Applied Mechanics, 2009, 76(1):011004. Nguyen TD, Govindjee S, Klein PA, et al. A material force method for inelastic fracture mechanics. Journal of the Mechanics and Physics of Solids, 2005, 53(1):91-121. Nikishkov GP, Atluri SN. Calculation of fracture mechanics parameters for an arbitrary three‐dimensional crack, by the ‘equivalent domain integral’ method. International
Journal
for
Numerical
Methods
in
Engineering,
1987,
24(9):1801-1821. Sciumè G, Benboudjema F. A viscoelastic unitary crack-opening strain tensor for crack width assessment in fractured concrete structures. Mechanics of Time-Dependent Materials, 2016, 21:223-243. Shih CF, Moran B, Nakamura T. Energy release rate along a three-dimensional crack front in a thermally stressed body. International Journal of Fracture, 1986, 30(2):79-102. Wang H, Lu W, Barber JR, et al. The roles of cohesive strength and toughness for crack growth in visco-elastic and creeping materials. Engineering Fracture Mechanics, 2016, 160:226-237. Yao W, Li X, Hu X. Viscoelastic crack analysis using symplectic analytical singular element combining with precise time-domain algorithm. International Journal of Fracture, 2018, 214:29-48. Yi J, Shen S, Muhunthan B, et al. Viscoelastic–plastic damage model for porous asphalt mixtures: Application to uniaxial compression and freeze–thaw damage. Mechanics of Materials, 2014, 70(1):67-75.
Yuval M, Giorgio O, Johannes TB, et al. Crack Initiation in Viscoelastic Materials. Physical Review Letters, 2018, 120(26):268002-1-6. Yu HJ, Kitamura T. A new domain-independent interaction integral for solving the stress intensity factors of the materials with complex thermo-mechanical interfaces. European Journal of Mechanics A/Solids, 2015, 49:500–509. Yu HJ, Wu LZ, Guo LC, Du SY, He QL. Investigation of mixed-mode stress intensity factors for nonhomogeneous materials using an interaction integral method. International Journal of Solids and Structures, 2009, 46:3710–3724. Yu NY, Li Q, Chen YH. Measurement of the M-integral for a hole in an aluminium plate or strip. Experimental mechanics, 2012, 52(7):855-863. Yu NY, Li Q. Failure theory via the concept of material configurational forces associated with the M-integral. International Journal of Solids and Structures, 2013, 50(25-26):4320-4332.
List of Table Captions Table 1 Material parameters of HTPB (Lv, 2013) Parameters
Values
Density (kg/mm3)
1.75106
E0 (MPa)
12.32
v
0.50
Table 2 Values of the M-integral for six equivalent domains No.
1
2
3
Dimensions (mm)
A1=39 A2=41
A1=34 A2=37
Area (mm2) M-integral (J/m)
40
53.25
66
46.47
48.33
47.68
4
5
6
A1=34 A2=40
A1=32 A2=40
72
111
144
47.13
47.73
47.23
A1=36 A1=34 A2=39.50 A2=38