calcium carbonate composites under uniaxial tension and three-point bending

calcium carbonate composites under uniaxial tension and three-point bending

Composite Structures 171 (2017) 370–381 Contents lists available at ScienceDirect Composite Structures journal homepage: www.elsevier.com/locate/com...

3MB Sizes 1 Downloads 24 Views

Composite Structures 171 (2017) 370–381

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

An analytical study of mechanical behavior of polypropylene/calcium carbonate composites under uniaxial tension and three-point bending Zhongqiang Xiong a, Yuqi Li a,⇑, Lulu Pan a, Jinhong Yu b, Shaorong Lu a,⇑ a Key Laboratory of New Processing Technology for Nonferrous Metals and Materials, Ministry of Education, College of Materials Science and Engineering, Guilin University of Technology, Guilin 541004, China b Key Laboratory of Marine New Materials and Related Technology, Zhejiang Key Laboratory of Marine Materials and Protection Technology, Ningbo Institute of Material Technology & Engineering, Chinese Academy of Sciences, Ningbo 315201, China

a r t i c l e

i n f o

Article history: Received 18 November 2016 Revised 13 March 2017 Accepted 15 March 2017 Available online 18 March 2017

a b s t r a c t Analytic solutions of test curve are beneficial to analyze mechanical properties and to understand the reason behind mechanical behaviors of materials. Herein, two relationships about both average true strain and average true stress on necking were obtained. Further, the homogenizing constitutive equation was obtained according to the linear viscoelastic constitutive equation becoming the foundation of solving analytic solutions of material behavior. The applicability of analytic solutions in the entire test range was discussed and improved by experimental results of both uniaxial tension and three-point bending. Analytic solutions of test curve in this work conform well to the experiments. The moduli and strengths of polypropylene/calcium carbonate composites were predicted according to this model. This mathematical model indicates that mechanical behavior of linear viscoelastic material was the outcome of competition between stress-increasing of extension and stress-decreasing of relaxation. It is hoped that the empirical formulae of two routine test curves can bring convenience for an analytical study of mechanical behavior of linear viscoelastic materials. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction Polymeric materials are typical viscoelastic materials. In particular, the particle-reinforced composites are macroscopically isotropic and simple structure, which are linear viscoelastic materials characterized in terms of the integral representation of the Boltzmann superposition principle as well as the differential representation based on the generalized Maxwell or Kelvin model. Two kinds of constitutive equations are equivalent. The linear viscoelastic constitutive equation of tridimensional anisotropy is also developed to characterize different components based on the correspondence principle to the theory of elasticity [1]. Furthermore, the fractional viscoelastic constitutive equation is employed [2–5], which fewer parameters are required to represent material viscoelastic behavior than traditional integer-order models [4]. Mechanical behaviors of materials are understood by the macromechanics theory that has a set of complete equations [1,6]. Accordingly, the mechanical problems of different materials are discussed in detail by macromechanics theory, so does viscoelastic theory, for instance, dynamic problems [7–10], anisotropic problems [11–14] and non-linear problems [5,12,15–17]. ⇑ Corresponding authors. E-mail addresses: [email protected] (Y. Li), [email protected] (S. Lu). http://dx.doi.org/10.1016/j.compstruct.2017.03.054 0263-8223/Ó 2017 Elsevier Ltd. All rights reserved.

Meanwhile, mathematical operations are usually difficult due to the complex boundary conditions or constructions in engineering, but some powerful numerical methods are developed to simulate the mechanical behaviors of materials under different situations, such as the finite difference method [18–21] or the finite element (FE) method [5,9,22–27]. Most engineering problems have been solved. The most common ways of mechanical performances test of materials, the uniaxial tension and three-point bending, are a testing standard to obtain the material constants with stationary boundary conditions. If there are some analytical solutions of mechanical behaviors of material, it will be beneficial to execute test in detail. This is meaningful for analyzing what the role of parameters is understood based on various functional relations rather than numerical calculations. In this paper, two relationships on necking were obtained from viscoelastic theory, which are the average true strain hex ðtÞi and average true stress hrx ðtÞi related to the nominal strain e0 and nominal stress r0 , respectively. Further, the homogenizing constitutive equation between hex ðtÞi and hrx ðtÞi was obtained and became the foundation of solving geometrically nonlinear mechanical behavior question. The analytic solutions of test curve both uniaxial tension and three-point bending were solved by the differential equations. And then the applicability of analytic solutions

371

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

was discussed emphatically via experimental results. The model was corrected by experiments and empirical formulae were built in two routine test to suit large deformed flexible polymeric materials. The moduli and strengths of polypropylene/calcium carbonate composites were predicted according to this model. The impact of different test rates to materials performance was discussed, which was convenient to convert among different test conditions. And the changing tendency of mechanical behavior in the entire test range influenced by parameters was shown. The establishment of the model will be beneficial to quantitative research the relationship among a series of ingredients of linear composites, such as guide formulation design of linear composites when the relationships among model parameters and the dosage of a series of raw material was known. 2. The homogenizing constitutive equation on necking 2.1. The definition of average true strain Viscoelastic materials are mostly high ductile materials, which are produced large deformation when it suffers the sustained stress. The large deformation is characterized by geometric equation. Corresponding to the condition of uniaxial tension (Detailed process are shown in the Appendix A), the elongation Dl of a specimen is equal to

Z

Dl ¼

Dl

Z

l0

du ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 2ex ðx; y; z; tÞ þ 1  f ðx; y; z; tÞ  1 dx

ð1Þ

0

0

where f ¼ ð@ v =@xÞ2 þ ð@w=@xÞ2 , ex is the true strain. There is and ex ðx; y; z; tÞ ¼ ex ðtÞ before necking generated. Hence the Eq. (1) reduces as 2

v ¼w¼0

1 2

ex ðtÞ ¼ e0 ðtÞ þ e20 where e0 ðtÞ ¼ Dl=l0 is the nominal strain. As is shown in part A of Fig. 1, however, v ; w – 0 and ex is related to the coordinates in the place where necking occur. (Fig. 1 is the finite element numerical simulation) There is also a formula resemblance to the relation before necking generated, which was proven in Appendix A.

hex ðtÞi ¼ e0 ðtÞ þ

1 2 e ðtÞ 2 0

where hex ðtÞi is called average true strain and hex ðtÞi 2 ½emin ; emax . emax is true strain of the smallest columnar shape in the specimen and emin is the biggest. Hence hex ðtÞ ¼ emin ðtÞi ¼ ex ðtÞ ¼ emax ðtÞ when there is nonexistent necking. 2.2. The definition of average true strain The volume variation ratio g is changing with the process of stretching as a function of the nominal strain

V0 ¼ gðe0 Þ V

ð2Þ

Hence, the true stress is obtained based on Eq. (2) when the condition is uniaxial tension before necking generated

F A

rx ðtÞ ¼ ¼

Fðl0 þ DlÞgðe0 Þ ¼ r0 ð1 þ e0 Þgðe0 Þ A0 l 0

where A0 , l0 are initial size of cross-sectional area and gauge length respectively. r0 ¼ F=A0 is the nominal stress. After necking, there is also a formula resemblance to the relation before necking, which was proven in Appendix B.

hrx ðtÞi ¼ r0 ð1 þ e0 Þgðe0 Þ where hrx ðtÞi is called average true stress and hrx ðtÞi 2 ½rmin ; rmax . rmax is the stress of the smallest columnar shape in the specimen and rmin is the biggest. Hence hrx ðtÞi ¼ rmin ðtÞ ¼ rx ðtÞ ¼ rmax ðtÞ when there is nonexistent necking. Some detailed investigations [28,29] that is volume variation process of materials during tensile are carried out. Assuming that the function gðe0 Þ is enough to form the Taylor’s series, the higher order in Oðe20 Þ can be omitted if the volume changed is not obvious in the finite tensile range. We say the linear variation, gðe0 Þ ¼ 1 þ ae0 . Therefore

hrx i ¼ r0 ð1 þ e0 Þð1 þ ae0 Þ 2.3. The constitutive equation related to the average true strain and average true stress The constitutive equation of linear viscoelastic materials does relate to ex ðtÞ and rx ðtÞ rather than hex ðtÞi and hrx ðtÞi. In this paper, the homogenizing constitutive equation on necking was proven according to the definition of the average true strain and average true stress based on the linear viscoelastic constitutive equation in Appendix C

hrx ðtÞi þ

i i m n X X d d pi i hrx ðtÞi ¼ qi i hex ðtÞi dt dt i¼1 i¼0

3. A uniaxial tension model and three-point bending model 3.1. The relationship between the nominal stress and nominal strain It was expedient to consider the standard linear solid constitutive equation (he_ x i denotes dhex i=dt)

hrx i þ p1 hr_ x i ¼ q0 hex i þ q1 he_ x i

ð3Þ

that included first-order derivative and only one relaxation time but a satisfactory conclusion was still given out, which was illustrated by the experiments in this paper. For uniaxial tension, e0 ðtÞ ¼ e_ 0 t. e_ 0 is the strain rate of extension, e_ 0 ¼ v 0 =l0 . According to Section 2.1, it is easy to obtain

hex i ¼ e_ 0 t þ

1 2 2 e_ 0 t ; 2

he_ x i ¼ e_ 0 þ e_ 0 2 t

ð4Þ

Substituting Eq. (4) into the (3) and solving it (Details of solving process are shown in the Appendix D), we get

   e0 þ C t e20 þ Dt e0 hrx i ¼ At 1  exp  Bt where At ¼ ð1  p1 e_ 0 Þðq1  p1 q0 Þe_ 0 ; Bt ¼ p1 e_ 0 ; C t ¼ q0 =2; Fig. 1. The varying stress in the necking A otherwise well-distributed stress respectively (Draw in ABAQUS 6.14-2).

Dt ¼ 2C t þ At =ð1  Bt Þ:

ð5Þ

372

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

The parameter Bt has a physical meaning that is relaxation strain corresponding to relaxation time p1 . Especially, there is p1 ; q1 ! 0þ for the linear elastic materials. The true stress is rx ¼ q0 e0 þ q0 e20 =2 Combining to the discussion of Section 2.2, we directly obtain the conclusion between the nominal stress and the nominal strain that is the most common curve in the uniaxial tension test

h

r0 ¼

 i At 1  exp  eB0t þ C t e20 þ Dt e0 ð1 þ e0 Þð1 þ ae0 Þ

0

e0 ¼0

¼

ð6Þ

dF Ab ¼ þ Cb dw w¼0 Bb

ð14Þ

the flexural modulus

At þ Dt Bt

ð7Þ

when the Eq. (5) is allowed, the (7) can simplify into Et ¼ q1 =p1 . the yield strength

ðA =B Þexpðe0 =Bt Þ þ 2C t e0 þ Dt ry ðAt ; Bt ; C t ; Dt ; aÞ ¼ t t 1 þ a þ 2ae0 where

3

of simply supported elastic beam in the midpoint is w ¼ Fl =ð48q0 Iz Þ. Based on the formula (12), further, it is convenient to deduce the performance parameters of materialthe elastic ratio

K b ðAb ; Bb ; C b Þ ¼

Based on the formula (6), further, it is convenient to deduce the performance parameters of material:the Young’s modulus

dr0 Et ðAt ; Bt ; Dt Þ ¼ de

The parameter Bb has a physical meaning that is relaxation deflection corresponding to relaxation time p1 . Especially, there is p1 ; q1 ! 0þ for the linear elastic materials. The deflection equation

ð8Þ

3

Eb ðAb ; Bb ; C b Þ ¼

e is the yield strain.the tensile strength

rt ðAt ; Bt ; C t ; Dt ; aÞ ¼ maxðr0 je0 ¼bi ; r0 je0 ¼eend Þ

ð15Þ

the flexural strength (Appendix F shows the details)

rb ðAb ; Bb ; C b Þ ¼

   M max l Ab ¼ Ab þ Bb C b 1 þ ln  4W z Wz Bb C b

ð16Þ

where W z is the flexural factor of section. When the Eq. (13) is allowed, the Eq. (15) can simplify into

Eb ¼

 0

l K b ðAb ; Bb ; C b Þ 48Iz

q1 ¼ Et p1

ð17Þ

and the Eq. (16) can simplify into

where ðdr0 =de0 Þje0 ¼bi ¼ 0, i ¼ 1; 2; . . . n. n is the number of extremums.

rb ¼ kv 0

ð18Þ

where k is a proportional constant just related to the nature and geometry of materials.

the fracture toughness

xt ðAt ; Bt ; C t ; Dt ; aÞ ¼

R eend 0

3.2.1. The error estimation According to the Eq. (11), the form of deflection equation is

r 0 d e0

2

ð4x2  3l Þx hðtÞ 48Iz

3.2. The relationship between the force and deformation

wðx; tÞ ¼

For pure bending, the flexural behavior of viscoelastic materials is characterized by the approximate deflection equation corre_ 00 denotes @ 3 w=ð@x2 @tÞ) sponding to the Eq. (3) (w

which ought to be limited by condition of w02 max  1 of Eq. (9). In we consider order to obtain the value of w0max , w00 ðx; tÞ ¼ hðtÞx=ð2Iz Þ ¼ 0 with hðtÞ–0 for t–0 and combine with wðx; tÞjx¼l=2 above. Hence

_z _ 00 ¼ Mz þ p1 M q0 Iz w00 þ q1 Iz w

ð9Þ

which is limited by condition < tan2 ðhmax Þ ¼ w02 max  1. I z is moment of inertia. Three-point bending of beam is the transverseforce bend that contain shear stress, which dissatisfy the condition of Eq. (9). Long beam, however, the maximum size of cross-section is so less than the span’s, which the bending shear stress is neglected. Hence, the equation is still used in the three-point bending test [30]. The bending moment of three-point bending is h2max

1 M z ðx; tÞ ¼ FðtÞx 2

ð10Þ

Substituting (10) into (9), using the boundary conditions and the initial conditions (Appendix E shows the details), an relationship is derived

wðx; tÞ ¼

      2 ð4x2  3l Þx p1 1 q0 p1 t  FðtÞ FðtÞ þ  2 exp  48Iz q1 q1 =q0 q1 q1

ð11Þ There is wðx; tÞjx¼l=2 ¼ v 0 t in the three-point bending test leading us to get a relationship between the force and deflection. The minus of the relationship was omitted, which just express the opposite direction with the y-axis. Its form is changed into

   w þ Cbw F ¼ Ab 1  exp  Bb

ð12Þ

w0max ðtÞ

2 12x2  3l ¼ hðtÞ 48Iz

¼ x¼0

3 wðx; tÞjx¼l=2 l

Thus, the deflection of midpoint should be limited by wðxÞjx¼l=2 ¼ w0max l=3. The formation of error  in the Eq. (9) is owing

to neglect the w0 in the curvature formula. The value of  is known based on a given deflection w

max ¼

w00  w00 ð1 þ w02 max Þ

3=2

3=2

w00 ð1 þ w02 max Þ

"

¼ 1þ



2 #3=2 3w 1 l

or calculate allowed deflection w on a given

l w¼ 3

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð1þmax Þ2=3  1

max ð19Þ

Therefore, a specific estimation to the condition w02 max  1 is got based on the model. This is independent with flexural rate. The domain of formula (12) is directly given under different error max by Eq. (19) as long as to know the span of beam such as l ¼ 60:00 mm. w 2 ½0; 3:64 mm of formula (12) when max ¼ 5%; w 2 ½0; 5:12 mm of formula (12) when max ¼ 10%. 4. Experiments 4.1. Materials and preparation

where Ab ¼ 48ðq1  q0 p1 Þv 0 Iz =l ¼ C b ðq1 =q0  p1 Þv 0 ; Bb ¼ p1 v 0 ; 3

3

C b ¼ 48Iz q0 =l :

ð13Þ

All materials used in this paper were commercially available and dried in a thermotank with 60 °Call night.

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

The polypropylene/calcium carbonate composites were prepared to indicate the results of preceding discussion by the following method: 57.55 wt% Polypropylene (PPH-T03, Standard: Q/ SHPRD253-2009, with a melt flow rate (MFR) of 1.31 g/10 min (at 200 °C, 10 kg) and a density of 0.88 g/cm3), 35.56 wt% CaCO3 (Light calcium carbonate of edible grade, Standard: GB18982007, with 1500 mesh and a density of 2.47 g/cm3), 6.60 wt% PPg-MAH (Grafting-PP LSS308; Full Name: Maleic anhydride grafted Polypropylene) and 0.29 wt% Antioxidant 1010 (SK1010) were put into twin-screw extruder (TLE16-4) homogeneous mixing with 200 °C, which the screw speed was maintained at 65 rpm. The product was pelleted and dried. Then specimens were shaped by the injection molding machine (TW-25V-1S) with 200 °C. Meanwhile, the mold temperature was maintained at room temperature.

373

As we can see from Fig. 2, the mathematical model is useful for uniaxial tension of polypropylene/calcium carbonate composites according to the experimental results. The phenomenon of mechanical behavior is comprehensive result by multiple factors such as the geometric nonlinearity considered in Section 2.1, the cross-section variation considered in Section 2.2, the viscoelastic attributes and the extension rates considered in Section 3.1. To explain the effect of different test rates for extension behavior of materials, the parameters of formula (6) how to change with the extension rate were described. In fact, there are still some simple relations in routine test shown in Fig. 3(b) (c) (d) (e) and (f) by experiment. The value of parameters is influenced by the fitting range, which is the characteristic of fitting method. For guaranteeing comparability, there is the same fitting range in different test rate for uniaxial tension as well as three-point bending.

4.2. Test procedure All tests of specimens were executed at room temperature. The uniaxial tension and three-point bending properties were measured by an electronic universal testing machine (SUNS UTM14483) according to the National Standard of China. The tensile test standard was GB/T 1040.2-2006 and the bending test standard was GB/T 9341-2008. A series of tension speed 2, 4, 7, 10, 20, 30, 40 and 50 mm/min was adopted in tensile test. A series of bending speed was 2, 4, 6, 8, 12 and 14 mm/min. All stress-strain curves and material constant data were obtained by the software related to the testing machine. For each group of specimen, all measurements were repeated to insure that there were 5 valid data. Then all experimental data were processed in the Origin 9.0.0 in this paper. The mean value and the standard deviation of the mechanical properties were reported. 5. Comparing model results and experimental data 5.1. Uniaxial tension 5.1.1. The stress-strain curve According to the conclusion of Section 3.1, the formula was h  i e

written as

r0 ð1 þ e0 Þ ¼

At 1exp  B0 t

þC t e20 þDt e0

ð1þae0 Þ

to fit the experimen-

tal data with different extension rates. (Incidentally, the unit of right side e0 is ‘‘%” while the left is ‘‘one”). It was beneficial to fit when Dt was regarded as an independent parameter, which Dt should have been calculated by At Bt and C t according to the Eq. (5). Herein, the uniaxial tension model was improved by the experiment directly and Eq. (5) became the most special situations. Fitting results are shown in Fig. 2(a). Converted into the relationship between nominal stress and nominal strain as known to all, the results are shown in Fig. 2(b).

5.1.2. The Young’s modulus and yield strength The Young’s modulus was calculated by the formula (7) according to the experimental relationships of parameters related to extension rates, results of Fig. 3(b) (c) (e). There is the same tendency of model results related to the extension rate. However, the value of model results is far away from the experimental as shown in Fig. 4(a). Checked carefully, the possible reasons are that there is the initial small orientation of clamp or an accelerated process from the quiescence to the setting test rate in test machine. Every single specimen is in test state of instability at initial measurement as we can see from Fig. 4(b). The Young’s modulus was tested according to the GB/T 1040.2-2006 within small strain 0.0025 precisely. On the other hand, there was a hypothesis e0 ðtÞ ¼ e_ 0 t without the accelerated process as well as the Young’s modulus was obtained by the first-order derivative of formula (6) at e0 ¼ 0 in model shown in Section 3.1. Hence, this is the main source of discrepancy. The yield strength was calculated according to the experimental relationships of parameters related to extension rates, results of Fig. 3(b) (c) (d) (e) and (f). It can be estimated by the relation of little changed yield strain e0 according to the formula (8) and also can be solved accurately via the first-extremums of formula (6). As we can see from the Fig. 5, there is about 1 MPa deviation than experiment. Indeed, the model is not good at explaining the yield behavior as shown in Fig. 3(a). The model results higher and lower than experiment are probably led by the change of microstructures [29] with strain shown in Fig. 4(b). This is the aspect neglected by the homogenizing approach. 5.2. Three-point bending 5.2.1. The force-deflection curve According to the conclusion of Section 3.2, the formula (12) was used to fit the experimental data with different flexural rates. It

Fig. 2. (a) The exhibition of fitting results; (b) The comparison of model results and raw data Note: It is taken on unloading that the Photo Report marked in the stress-strain curve shown in (b).

374

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

Fig. 3. (a) The details of fitting results about different extension rates in the part ② to ③ of Fig. 2(b); (b) (c) (d) (e) (f) The relationship between extension rates and each parameter of formula (6): At ¼ 0:0016v 20 þ 0:042v 0 þ 34:16; Bt ¼ 0:0079v 0 þ 3:35; C t ¼ 2:16v 0 þ 247:35, Dt ¼ 0:78v 0 þ 55:30; a ¼ 1:79v 0 þ 490:74. Note: For guaranteeing comparability, there is the same fitting range e0 2 ½0; 40% in different rate.

Fig. 4. (a) The Young’s modulus related to different extension rates; (b) The reason of deviation of model results.

was commendable to fit on entire range as shown in Fig. 6(a), which the formula (12) should have been used in small deflection with max according to the error estimation in Section 3.2.1. Herein, the three-point bending model was extended by the experiment

directly and Eq. (13) became the most special situations. It is worth noting that parameters lost original meaning such as C b turned into a negative by extending. The error estimation was not suitable anymore owing to the shaky Eq. (13).

375

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

previously. The answer is explicit according to the experimental data in this section. For rectangle, parallel axis theorem indicates

Iz ¼ Izc þ g2 A ¼

3

bh þ g2 bh 12

ð20Þ

where g is the distance apart from the centroidal axis. Here some typical size of bending specimen: width b ¼ 10:02 mm, thickness h ¼ 4:12 mm and the span of beam l ¼ 60:00 mm. Experimental results of elastic ratio K b and flexural modulus Eb are shown in Appendix G. According to the Eq. (15), the

c ¼ K b =Eb ¼ 48Iz =l3 does not change with test rate. The value of c was calculated by experimental data of Appendix G and shown in Table 3. The mean value of c, calculated by Table 3, is c ¼ 0:01325 mm. Combining with the Eq. (20), we obtain (Minus means opposite direction)

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 cl3 h g¼  ¼ 0:17 mm 48bh 12

Fig. 5. The comparison of yield strength between model accurate results and experimental data.

To explain the effect of different test rates for flexural behavior of materials, the parameters of formula (12) how to change with the flexural rate were described. In fact, there are still some linear relations in routine test shown in Fig. 6(b) (c) (d) by experiment. 5.2.2. The location of inertia axis It only had vague notions about the moment of inertia Iz where to find its inertia axis that the three-point bending was discussed

The flexural strength was about 28 MPa in experiment, so negative solution was given up. The location of inertia axis on the cross section is shown in Fig. 7.

Table 3 The value of c in different extension rates.

v0 (mm/min)

2

4

6

8

10

12

14

c (102 mm)

1.330

1.327

1.318

1.320

1.329

1.323

1.325

Fig. 6. (a) The fitting results of formula (12) on entire range and different flexural rates; (b) (c) (d) The relationship between flexural rates and each parameter of formula (12): Ab ¼ 0:64v 0 þ 71:02; Bb ¼ 0:02v 0 þ 3:75; C b ¼ 0:01v 0  1:81. Note: It is taken on unloading that the Photo Report marked in the force-deflection curve shown in (a) and there is the same fitting range w 2 ½0; 23 mm.

376

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

deformed in the big deflection w in the three-point bending test. Hence, the value of W z was changing with the process of bending no longer a rectangle’s. The mean value of W z was discussed

Wz ¼

Iz ¼ 28:02 mm3 ymax

ð22Þ

The predicted result based on Eq. (22) is shown in Fig. 8(b), which was executed the same operational process as rectangle. 6. The effect of parameters Fig. 7. The location of inertia axis on the cross section Note: Point O is the symcenter and the distance g is disproportionate to h for expressing more clearly.

5.2.3. The flexural modulus and flexural strength The flexural modulus was calculated by the formula (15) according to the experimental relationships of parameters related to flexural rates, results of Fig. 6(b) (c) (d). The result is shown in Fig. 8(a). According to the location of inertia axis on the cross section, the flexural factor of section was calculated

bh =12 þ g2 bh ¼ 26:72 mm3 h=2 þ g 3

Wz ¼

ð21Þ

The model was established based on the viscoelastic theory to suit two common tests, further, corrected and improved by experiments. The effect of parameters to the model was briefly discussed in this section. Here just mentioned the influence of each parameter and the others were constant at the level shown in Tables 1 and 2. The variation tendencies were shown in Figs 9 and 10. In fact, parameters are synergistic effect on model and different situations are needed to fit. The parameters changing tendency is beneficial to us to understand this model better. For instance, the results of Fig. 6(b) (c) (d) are parametric equations that are a space curve in space of Fig. 10(d). Hence, the flexural strength is known as shown in Fig. 8(b). 7. Discussions

Combining Eqs. (16) and (21) with fitting results of Fig. 6 (b) (c) (d), the flexural strength was predicted based on the rectangular cross-section implied in Eq. (21). The comparison of flexural strength between model predicted results and experimental data is shown in Fig. 8(b). The prediction was about 1.5 MPa higher than the experimental data based on the rectangular cross-section, but the tendency was satisfactory. We noted that the cross-section

The mathematical model deduced from viscoelastic theory is phenomenological. This model was shown the mechanical behavior of linear viscoelastic material, which was the outcome of competition between the stress-increasing of extension and stress-decreasing of relaxation. The concept of yield may be not needed anymore to explain the reason why the nominal stress or

Fig. 8. (a) Predicted flexural modulus by formula (15); (b) The comparison of flexural strength between model predicted results and experimental data. Note: Dash curve result based on the rectangular cross-section implied in Eq. (21) and solid curve result based on the deformed cross-section implied in Eq. (22).

Table 1 The value of fitting parameters of formula (6) (10 mm/min). At (MPa)

Bt (%)

Ct (104MPa)

Dt (102 MPa)

a (102)

Adj. R2

F Value

Prob > F

34.7683

3.6118

0.0192

1.5341

0.0765

0.9994

5.5405E7

0.0000

Table 2 The value of fitting parameters of formula (12) (10 mm/min). Ab (N)

Bb (mm)

Cb (N/mm)

Adj. R2

F Value

Prob > F

72.9466

3.4503

1.8633

0.9996

2.74E7

0.0000

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

377

Fig. 9. The effect of each parameter on the formula (6) in uniaxial tension.

force would decrease to a certain degree. And the reason why Young’s modulus is not equal to flexural modulus tested by experiment (Shown in Appendix G) instead of Eq. (17) was cleared by the experimental results of Figs. 3 and 6. The fitting parameters A; B; C; D were concerned in this paper, which was convenient to drive the performance parameters of materials and satisfactory to grasp the mechanical behavior of linear viscoelastic materials. It will be beneficial to guide formulation design of linear composites when the relationships among model parameters and the dosage of a series of raw material were known. Of course, there are some deficiencies of this model. For instance, the place where necking and fracture generate cannot be provided by this model. It is not equality that material constants p1 ; q0 ; q1 come from two kinds of tests based on the same standard linear viscoelastic solid constitutive equation. There are some reasons to explain: The different assumptions and approximations

were used for different test; besides, the model was corrected by experiments. It is more adaptable when the value of fitting parameters A; B; C; D was considered than p1 ; q0 ; q1 . Otherwise, the highorder including fractional order constitutive equation or viscoelastic spectrum and its traits [31,32] are needed to note as well as numerical calculations are necessary.

8. Conclusions The polypropylene/calcium carbonate composites mechanical behavior of entire range under both two kinds of tests was simulated. The analytic solutions were given out, in addition, corrected and improved by experiments. Hence, the empirical formulae were built for an analytical study of mechanical behavior of linear viscoelastic materials:

378

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

Fig. 10. (a) (b) (c) The effect of each parameter on the formula (12) in three-point bending; (d) The effect of parameters to the flexural strength plotted by formula (16) in MATLAB 2012b.

The relationship (formula (6)) of stress-strain curve of linear viscoelastic materials in the entire cold-drawing range is h  i e

r0 ð e0 Þ ¼

At 1exp  B0 t

þC t e20 þDt e0

ð1þe0 Þð1þae0 Þ

The relationship (formula 12) of force-deflection curve of linear viscoelastic materials in the entire flexural range is h  i FðwÞ ¼ Ab 1  exp  Bw þ C b w b

Model constants (the uniaxial tension has five and the threepoint bending has three) needed to fit with experimental data ultimately. These solutions were performed well by experiment and then the moduli and strengths were illustrated by the model. Consequently, it was convenient to convert among different test conditions. This mathematical model put the mechanical behavior of materials coming down to a few parameters which contained most of mechanical information such as the Young’s modulus, the tensile strength, the flexural modulus, the flexural strength, etc. It is hoped that the empirical formulae of two routine tests can bring convenience for an analytical study of mechanical behavior of linear viscoelastic materials. Acknowledgements The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China (51463007 and 51605109), the Natural Science Foundation of Guangxi Province, China (2014GXNSFDA118006), Guangxi Ministry-Province Jointly-Constructed Cultivation Base for State Key Laboratory of Processing for Non-ferrous Metal and Featured Materials (15AA, 15KF-9), Young teachers based capacity building project of

Guangxi (KY2016YB198), Science and Technology Development Program of Nanning city (20151267). Appendix A The

large

deformation

is

characterized

by

equation

eij ¼ ðrj;i þ ri;j þ rk;i rk;j Þ=2, r ¼ ðu; v ; wÞ (ri;j denotes index convention) [33] and the conditions correspond to uniaxial tension. The true strain is

"  2  2  2 # 1 @u @u @v @w 2 þ þ þ 2 @x @x @x @x  2 @u 1 @u 1 2 þ f ðx; y; z; tÞ ¼ þ @x 2 @x 2

ex ¼

2

where f implies the measurement of necking. @u=@x is solved out and give up the negative solution (It means the compression of materials)

@u ¼ @x

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2ex þ 1  f  1

Thus, the elongation Dl of the specimen is equal to

Z

Dl

Dl ¼

Z

du ¼

0

l0

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð 2ex þ 1  f  1Þdx

0

ex is related to the time ðtÞ and the coordinates ðx; y; zÞ, but we can get the following conclusion

rðx; y; z; tÞ ¼ r

ðA:1Þ

379

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

where l0 ¼

Pn

i¼0 ðliþ1

 li Þand mark

Xn pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2exi þ 1ðliþ1  li Þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i¼0 Xn ¼ 2ex ðtÞ þ 1 ðl  li Þ i¼0 iþ1

Fig. 11. The phenomenon of necking.

The formula (A.2) is reduced as from the rij;j þ f i ¼ q€r i [33] which disregard body force f i and inertial force q€ri ¼ 0 because it is uniform velocity process, e0 ðtÞ ¼ e_ 0 t, in the uniaxial tension condition. It is easy to obtain the strain eðx; y; z; tÞ ¼ eðtÞ via the constitutive equation of linear viscoelastic materials limited by Eq. (A.1) and there is v ¼ w ¼ 0 before necking. 2

Hence f ðx; y; z; tÞ ¼ 0.

hex ðtÞi ¼ e0 þ

1 2 e 2 0

where hex ðtÞi is called average true strain and hex ðtÞi 2 ½emin ; emax . emax is the strain of the smallest columnar shape in the specimen, and emin is the biggest. Hence ex ðtÞ ¼ emin ðtÞ ¼ hex ðtÞi ¼ emax ðtÞ when

ex ¼ ð1 þ Þ  ¼ e0 þ e20

there is nonexistent necking. Equivalently, in a simple way, we can obtain the answer directly from the mean value theorem for integrals

where the nominal strain is written as



Dl l0

1 2

2

1 2

1 2

Z

2

varying with time and coordinates as well as f ðx; y; z; tÞ–0. Supposing that there is n-segment necking (such as Fig. 11), Marking li express i th place arising necking in the original specimen when i is odd number, and liþ1 express the end of i th necking (Here marking in the deformed specimen for convenience to illustrate). So ðli ; liþ1 Þ is normal columnar shape when i is even number. Thus, the elongation Dl of the specimen is equal to

¼

liþ1

X Z

i is odd

þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð 2ex ðx; y; z; tÞ þ 1  f ðx; y; z; tÞ  1Þdx

li

i¼0

liþ1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 ð 2ex ðx; y; z; tÞ þ 1  f ðx; y; z; tÞ  1Þdx

X Z

liþ1

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð 2ex þ 1  1Þdx

li

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 2ex ðx; y; z; tÞ þ 1  f ðx; y; z; tÞ  1 dx

i is ev en

i is odd

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ð 2ex þ 1  1Þ

We can omit the value of f ðni ; tÞ when f ðni ; tÞ  2ex ðni ; tÞ þ 1. It implies small necking. Regard ex ðni ; tÞ and ex as exi when i is odd and even respectively. Hence 2

 Xn pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2exi þ 1  1 ðliþ1  li Þ i¼0 Xn Xn pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðl  li Þ 2exi þ 1ðliþ1  li Þ  ¼ i¼0 i¼0 iþ1 Xn pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2exi þ 1ðliþ1  li Þ i¼0 Xn  1 ¼ e0 ðl  li Þ i¼0 iþ1

Correspondingly, average true stress is driven from the compressible condition implied in Eq. (2). The force F is the same in the whole specimen despite necking generated.

V gðe0 Þ ¼ rx Ax V0

V ¼F l

l V F gðe0 Þ ¼ r0 ð1 þ e0 Þgðe0 Þ V V0

where hrx ðtÞi is called average true stress and hrx ðtÞi 2 ½rmin ; rmax . rmax is the stress of the smallest columnar shape in the specimen, and rmin is the biggest. Hence hrx ðtÞi ¼ rmin ðtÞ ¼ rx ðtÞ ¼ rmax ðtÞ when there is nonexistent necking. It is worth noting that the definition of hrx ðtÞi come from the general mean value theorem for integrals

R l0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2 ð 2ex þ 1  f Þdx l 1 ¼ rx Ax hrx ðtÞi ¼ rx Ax ¼ rx Ax R 0 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi V Af l0 2 Ax ð 2ex þ 1  f Þdx 0 where f 2 X0 . Appendix C

Dl ¼

Dl ¼ l0

Appendix B

hrx ðtÞi ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi X 2 ðliþ1  li Þð 2ex ðni ; tÞ þ 1  f ðni ; tÞ  1Þ þ ðliþ1  li Þ

2

ex ðtÞ ¼ en

Hence

where ni ¼ ðxi ; yi ; zi Þ 2 Xi . Xi is the body of specimen between li and liþ1 .Hence

X

Define

hrx ihAx i ¼ hrx i

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 2ex ðni ; tÞ þ 1  f ðni ; tÞ  1 ¼ ðliþ1  li Þ

Dl ¼

1 2

where rx ; Ax are the true stress and present size of each crosssection, respectively. Defined

li

From the mean value theorem for integrals, we can get liþ1

2

en ðtÞ ¼ e0 þ e20

F¼F

li

i is ev en

Z

We omit the value of f ðn; tÞ when f ðn; tÞ  2en þ 1, n 2 X0 . It implies small necking. 2

In fact, we can also discuss the Dl after necking. The specimen is no longer columnar shape in the place of necking. Therefore, r is

Dl ¼

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 ð 2ex þ 1  f Þdx ¼ l0 2en þ 1  f ðn; tÞ

0

DlðtÞ e0 ðtÞ ¼ l0

n Z X

l0

For stretching, the form of constitutive equation of linear viscoelastic materials is given by [34]

ðA:2Þ

rx ðtÞ þ

Xm

p i¼1 i

i

d

dt

i

rx ðtÞ ¼

Xn

q i¼0 i

i

d

dt

i

ex

ðC:1Þ

380

Z. Xiong et al. / Composite Structures 171 (2017) 370–381

The Eq. (C.1) does relate to ex ðtÞ and rx ðtÞ rather than hex ðtÞi and hrx ðtÞi. Processing the Laplace transformation to Eq. (C.1) with the initial conditions ex ð0Þ ¼ rx ð0Þ ¼ 0, we get

1þ ~ex ðsÞ ¼

m X

pi si

i¼1

n X qi si

r~ x ðsÞ ¼ ~JðsÞr~ x

Appendix E The bending moment of three-point bending is Eq. (10). Due to the symmetry, we just need to consider half scope of beam, x 2 ½0; l=2Þ. The Eq. (10) was substituted into Eq. (9). We can derived

_ 00 ¼ q0 Iz w00 þ q1 Iz w

i¼0

Hence, the inverse Laplace transformation of ~ex ðsÞ is

ex ðtÞ ¼ JðtÞ  rx

ðC:2Þ

where JðtÞ ¼ L1 f~JðsÞg and is independent with coordinates. Combining the definition of rx ðtÞ shown in Appendix B with the Eq. (C.2), we get

Af hex iðtÞ ¼ JðtÞ  hrx ðtÞi Ax

ðC:3Þ

Combining the definition of hex ðtÞi shown in Appendix A with the Eq. (C.3) as well as noting that the JðtÞ; Af ; rx ðtÞ are independent with coordinates and convolution is an operation relate to time, we obtain

Af hex ðtÞi ¼ en ðtÞ ¼ JðtÞ  hrx ðtÞi An

hex ðtÞi ¼ JðtÞ  hrx ðtÞi ¼

Z

t

This is an equation whose the variables have been separated. Carrying out two integrations with respect to x and using the _ 0 ðl=2Þ ¼ 0 and wð0Þ ¼ wð0Þ _ ¼0 boundary conditions, w0 ðl=2Þ ¼ w and using the initial conditions, Fð0Þ ¼ wð0Þ ¼ 0, proceeding the Laplace transform:

~ ¼ FðsÞ

48Iz q0 þ q1 s ~ w 2 ð4x2  3l Þx 1 þ p1 s

~ We find the inverse transformation of FðsÞ:

FðtÞ ¼

     q1 q0 q1 t exp   wðtÞ wðtÞ þ  2 p1 p1 p21 ð4x2  3l Þx p1

Jðt  sÞhrx ðsÞids

48Iz

      2 ð4x2  3l Þx p1 1 q0 p1 t  FðtÞ FðtÞ þ  2 exp  48Iz q1 q1 =q0 q1 q1

For three-point bending test, wðx; tÞjx¼l=2 ¼ v 0 t, we can get

FðtÞ ¼ 

48Iz

0

that is similar to Eq. (C.2), which is an integral constitutive equation with JðtÞ ¼ hrx ðtÞ ¼ 0i for t 6 0. Using the Laplace transformation

Pn

P i i and LfJðtÞg ¼ 1 þ m i¼1 pi s = i¼0 qi s We obtain

  Xn  Xm ~x ¼ 1þ p si r q si h~ex i i¼1 i i¼0 i

Xm

p i¼1 i

i

d

dt

hrx ðtÞi ¼ i

Xn

i

q i¼0 i

d

48Iz l

3

   w q0 w þ ðq1  q0 p1 Þv 0 1  exp  p1 v 0

dt

   w  Cbw F ¼ Ab 1  exp  Bb Ab ¼ 48ðq1  q0 p1 Þv 0 Iz =l ¼ C b ðq1 =q0  p1 Þv 0 , 3

Bb ¼ p1 v 0 ,

C b ¼ 48Iz q0 =l . Appendix F

1 2 1 e ¼ e_ 0 t þ e_ 20 t2 ; 2 0 2

he_ x i ¼ e_ 0 þ e_ 20 t

ðD:1Þ

where e_ 0 ¼ v 0 =l0 . v 0 is the rate of extension. Substituting Eq. (D.1) into the Eq. (3) and using the initial conditions, rx ð0Þ ¼ 0, the Laplace transform method is employed to find the true stress rx

rx ðtÞ ¼ ½q0 þ ðq1  q0 p1 Þe_ 0 e_ 0 t þ

q0 2 2 e_ 0 t 2

þ ð1  p1 e_ 0 Þðq1  q0 p1 Þe_ 0 

FðwÞ ¼ 

   t q0 v 0 t þ ðq1  q0 p1 Þv 0 1  exp  p1

3

According to the relationship between the true strain and nominal strain with the nominal strain e0 ¼ e_ 0 t for uniaxial tension, it directly drives

Simplified into

3

l

and

where

hex i i

Appendix D

hex i ¼ e0 þ



Simplified into

Eventually, we prove the differential constitutive equation does relate to hex ðtÞi and hrx ðtÞi

hrx ðtÞi þ



We can also derive

wðx; tÞ ¼

where Af ; An 2 ½Amin ; Amax . Hence, Af =An  1 when there is small necking just as Appendix A. Therefore

1 1 _ FðtÞx þ p1 FðtÞx 2 2



rx ðe0 Þ ¼ At 1  exp 

e0 Bt



rb ¼

Mz;max F max l ¼ 4W z Wz

where the F max is the maximum (the extremum) of Eq. (12). Hence

  dF Ab w þ Cb ¼ 0 exp  ¼ dw Bb Bb

   Ab F max ¼ Ab þ Bb C b 1 þ ln  Bb C b thereby

  t 1  exp  p1

rb ¼

þ C t e20 þ Dt e0

where At ¼ ð1  p1 e_ 0 Þðq1  p1 q0 Þe_ 0 , Dt ¼ 2C t þ At =ð1  Bt Þ

According to the Eq. (10) and the definition of the flexural strength

   l Ab Ab þ Bb C b 1 þ ln  4W z Bb C b

Appendix G Bt ¼ p1 e_ 0 ,

C t ¼ q0 =2,

According to the Eqs. (7) and (15), the material constant change with test rate. The experimental results were shown in Tables 4–6,

381

Z. Xiong et al. / Composite Structures 171 (2017) 370–381 Table 4 The Young’s modulus and yield strength related to different extension rates.

v0 (mm/min)

2

4

7

10

20

30

40

50

Et (MPa) STDEV (MPa) ry (MPa) STDEV (MPa)

823.90 10.22 25.72 0.29

842.92 16.44 26.02 0.25

855.57 11.94 26.37 0.22

875.83 16.56 26.59 0.37

915.25 14.14 27.01 0.36

941.88 18.09 27.21 0.36

964.36 8.30 27.40 0.25

970.42 18.69 27.50 0.30

Table 5 The flexural modulus and flexural strength related to different flexural rates.

v0 (mm/min)

2

4

6

8

10

12

14

Eb (MPa) STDEV (MPa) rb (MPa) STDEV (MPa)

1354.82 34.85 26.88 0.49

1392.37 19.95 27.91 0.78

1428.42 31.99 28.08 0.78

1476.28 9.65 28.77 0.63

1514.44 22.14 29.60 0.91

1541.76 23.81 29.37 0.64

1600.33 6.10 30.06 0.52

Table 6 The elastic ratio related to different flexural rates.

v0/(mm/min)

2

4

6

8

10

12

14

Kb/(N/mm) STDEV/(N/mm)

18.02 0.42

18.48 0.61

18.82 0.62

19.48 0.39

20.12 0.33

20.40 0.60

21.20 0.39

which was the mean value of five sets of valid data under each of test rate. The Young’s modulus is not equal to flexural modulus shown by experiment. (STDEV denotes standard deviation) References [1] El Kouri M, Bakkali A, Azrar L. Mathematical modeling of the overall timedependent behavior of non-ageing viscoelastic reinforced composites. Appl Math Model 2016;40:4302–22. [2] Heymans N, Podlubny I. Physical interpretation of initial conditions for fractional differential equations with Riemann-Liouville fractional derivatives. Rheol Acta 2006;45:765–71. [3] Koeller RC. Application of fractional calculus to the theory of viscoelasticity. J Appl Mech 1984;51:299–307. [4] Liu JG, Xu MY. Higher-order fractional constitutive equations of viscoelastic materials involving three different parameters and their relaxation and creep functions. Mech Time-Depend Mater 2006;10:263–79. [5] Muller S, Kastner M, Brummund J, Ulbricht V. On the numerical handling of fractional viscoelastic material models in a FE analysis. Comput Mech 2013;51:999–1012. [6] Christensen R. Theory of viscoelasticity. 2nd ed. Academic Press; 1982. [7] Vaz MA, Ariza AJ. Quasi-static response of linear viscoelastic cantilever beams subject to a concentrated harmonic end load. Int J Nonlinear Mech 2013;54:43–54. [8] Graziani A, Cardone F, Virgili A. Characterization of the three-dimensional linear viscoelastic behavior of asphalt concrete mixtures. Constr Build Mater 2016;105:356–64. [9] Eratli N, Argeso H, Calim FF, Temel B, Omurtag MH. Dynamic analysis of linear viscoelastic cylindrical and conical helicoidal rods using the mixed FEM. J Sound Vibr 2014;333:3671–90. [10] El Hachemi M, Koutsawa Y, Nasser H, Giunta G, Daouadji A, Daya EM, et al. An intuitive computational multi-scale methodology and tool for the dynamic modelling of viscoelastic composites and structures. Compos Struct 2016;144:131–7. [11] Rasolofosaon PNJ. Generalized anisotropy parameters and approximations of attenuations and velocities in viscoelastic media of arbitrary anisotropy type theoretical and experimental aspects⁄. Geophys Prospect 2010;58:637–55. [12] Ryou H, Chung K, Yu W-R. Constitutive modeling of woven composites considering asymmetric/anisotropic, rate dependent, and nonlinear behavior. Compos Part A: Compos Part A-Appl S 2007;38:2500–10. [13] Chung K, Ryou H. Development of viscoelastic/rate-sensitive-plastic constitutive law for fiber-reinforced composites and its applications. Part I: theory and material characterization. Compos Sci Technol 2009;69:284–91. [14] Rebouah M, Chagnon G. Extension of classical viscoelastic models in large deformation to anisotropy and stress softening. Int J Nonlinear Mech 2014;61:54–64. [15] Zhu HL, Muliana A, Rajagopal KR. On the nonlinear viscoelastic deformations of composites with prestressed inclusions. Compos Struct 2016;149:279–91.

[16] Abdelrahman AA, El-Shafei AG, Mahmoud FF. Analysis of steady-state rolling contact problems in nonlinear viscoelastic materials. J Tribol-Trans ASME 2015;137:10. [17] Trivaudey F, Placet V, Guicheret-Retel V, Boubakar ML. Nonlinear tensile behaviour of elementary hemp fibres. Part II: modelling using an anisotropic viscoelastic constitutive law in a material rotating frame. Compos Pt A-Appl Sci Manuf 2015;68. [18] Zhao WJ, Chen LQ, Zu JW. Finite difference method for simulating transverse vibrations of an axially moving viscoelastic string. Appl Math Mech-Engl Ed 2006;27:23–8. [19] Asvadurov S, Knizhnerman L, Pabon J. Finite-difference modeling of viscoelastic materials with quality factors of arbitrary magnitude. Geophysics 2004;69:817–24. [20] Kristek J, Moczo P. Seismic-wave propagation in viscoelastic media with material discontinuities: a 3D fourth-order staggered-grid finite-difference modeling. Bull Seismol Soc Amer 2003;93:2273–80. [21] Schmalholz SM, Podladchikov YY, Schmid DW. A spectral/finite difference method for simulating large deformations of heterogeneous, viscoelastic materials. Geophys J Int 2001;145:199–208. [22] Zobeiry N, Malek S, Vaziri R, Poursartip A. A differential approach to finite element modelling of isotropic and transversely isotropic viscoelastic materials. Mech Mater 2016;97:76–91. [23] Hartmann S, Hamkar AW. Rosenbrock-type methods applied to finite element computations within finite strain viscoelasticity. Comput Methods Appl Mech Eng 2010;199:1455–70. [24] Zhang HH, Li LX. Modeling inclusion problems in viscoelastic materials with the extended finite element method. Finite Elem Anal Des 2009;45:721–9. [25] Fritzen F, Bohlke T. Reduced basis homogenization of viscoelastic composites. Compos Sci Technol 2013;76:84–91. [26] Fliegener S, Hohe J, Gumbsch P. The creep behavior of long fiber reinforced thermoplastics examined by microstructural simulations. Compos Sci Technol 2016;131:1–11. [27] Tan JF, Jia YJ, Li LX. A series–parallel mixture model to predict the overall property of particle reinforced composites. Compos Struct 2016;150:219–25. [28] Addiego F, Dahoun A, Sell CG, Hiver JM. Volume variation process of highdensity polyethylene during tensile and creep tests. Oil Gas Sci Technol 2006;61:715–24. [29] Parenteau T, Ausias G, Grohens Y, Pilvin P. Structure, mechanical properties and modelling of polypropylene for different degrees of crystallinity. Polymer 2012;53:5873–84. [30] Gere JM, Goodno BJ. Mechanics of materials. 8th ed. Cengage Learning; 2012. [31] McDougall I, Orbey N, Dealy JM. Inferring meaningful relaxation spectra from experimental data. J Rheol 2014;58:779–97. [32] Anderssen RS, Davies AR, de Hoog FR. On the sensitivity of interconversion between relaxation and creep. Rheol Acta 2008;47:159–67. [33] Banks HT, Hu SH, Kenz ZR. A brief review of elasticity and viscoelasticity for solids. Adv Appl Math Mech 2011;3:1–51. [34] Findley WN, Lai JS, Onaran K. Creep and relaxation of nonlinear viscoelastic materials. North-Holland Publishing Company; 1976.