Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 UDK - UDC 620.172.21:678.7 Izvirni znanstveni clanek - Original scientific paper (1.01) Predvidevanje nelinearnega lezenja izdelkov Prediction of the Nonlinear Creep Deformation of Plastic Products Jan Spoormaker1 - Ihor Skrypnyk2 - Anton Heidweiller3 (1Spoormaker Consultancy, The Netherlands; 2Goodyear Corporate Research, USA; 3Delft University of Technology, The Netherlands) Prispevek na primeru nelinearnega lezenja materiala vstopnika zraka prikazuje sodobne zmožnosti analiz z MKE, zapletenih 3D viskoelasticnih deformacij izdelkov iz plastike. Predstavljena je pomembnost tovrstnih zmožnosti za konstruiranje zapletenih komponent iz plastike. Model nelinearne viskoelasticnosti je vkljucen kot podprogram v tržne programske rešitve, ker le-te še ne dajejo tovrstne uporabe. V prispevku so podane posebnosti modela in njegova numericna izvedba v tržnem programskem paketu MSC.MARC. © 2007 Strojniški vestnik. Vse pravice pridržane. (Kljucne besede: nelinearno lezenje, sprostitve napetosti, analize koncnih elementov, polimerni izdelki) Based on an example of the non-linear creep deformations of an air inlet, this paper demonstrates modern capabilities in the FEA modeling of complex 3D visco-elastic deformations in relation to the design of plastic products. The importance of such capabilities for designing complex plastic components is discussed. Because commercial FEA packages do not yet render these capabilities "off the shelf, the nonlinear visco-elasticity model is incorporated through a user subroutine. The specifics of the constitutive model and its numerical implemen-tation are outlined for the case of implementation in the commercial FEA package MSC.MARC. © 2007 Journal of Mechanical Engineering. All rights reserved. (Keywords: non-linear creep, stress relaxation, FEA modeling, plastic product design) 0 INTRODUCTION The everyday design practice for plastic components relies on relatively simple but effective tools. Just a decade ago these tools would typically include linear elastic finite-element analysis (FEA). In the case when a long-term prediction was required, these FEA calculations would make use of the isochronous curves. Due to the technical limitations of computer hardware even linear visco-elastic FEA calculations were not usually possible. It is well known, however, that most engineering plastics feature pronounced non-linear visco-elastic behaviour. Accounting for such long-term behaviour is important for designing sustainable, reliable plastic products. The advances in computer hardware and FEA techniques over the years brought the capabilities of visco-elasticity modelling within the reach of designers of plastic products. The main obstacle that remains here, however, is the lack of well-established models for the time-dependent non-linear behaviour that would be readily available in commercial FEA software. This paper considers a relatively simple example of the product design of a plastic air inlet (Fig. 1). One of the main design constraints in this case is related to the long-term creep of the product. Based on this design case the capabilities and the importance of non-linear visco-elastic FEA modelling for plastic design are demonstrated. The prediction of the creep behaviour of the air inlet was a demonstration project for a company that designs and manufactures plastic products for the automotive industry. The upper beam of the air inlet is loaded with 60 N for a period of 10 hours. The latter task itself consists of two challenges: experimental measurement of the time-dependent behaviour of plastic material and a FEA 834 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 Fig. 1. Photograph of the air inlet simulation using a non-linear visco-elastic constitutive model. In order to answer these types of questions, it is necessary to be capable of predicting non-linear visco-elastic deformations of complex three-dimensional solids. The creep behaviour of the applied plastic, PP with talc, had been measured in an inaccurate way and resulted then in a discrepancy between the predicted behaviour and the behaviour determined in the experimental verification. After accurately measuring the creep behaviour the discrepancy between the prediction and the experimental verification was sufficiently small to demonstrate the possibilities of predicting non-linear creep deformation. 1 NON-LINEAR VISCO-ELASTICITY MODEL OF THE RLAXATION TYPE The traditional constitutive model for visco-elasticity (Leaderman (1943), Schapery (1969)) is a creep-type dependency, where strains are given as functions of stress, time and, possibly, temperature. This convention comes from the comparative ease and cost efficiency of creep experiments (vs. their relaxation counterpart). It is necessary to mention, however, that the commercial FEA software packages are displacement based, meaning that a constitutive dependence should be provided in the form of stresses as functions of strains, time and temperature (MARC, Vol. D (1998)). In an incremental form the dependence will look as follows: Ds = DDe +G(qi ,Dt,K) (1). Predvidevanje nelinearnega lezenja - Prediction of the Nonlinear Creep Deformation 835 This means that a creep-type dependency should be inverted at each increment of the solution process. The inversion of the non-linear tensorial equation will easily lead to non-convergence (Skrypnyk, Spoormaker & Smit (2000)), making the task of implementing the creep-type models into the FEA very complicated. It is much easier (Skrypnyk & Spoormaker (2000)) to handle the implementation of the relaxation-type models, which originate from the Maxwell relaxation model (Fig. 2). These models have a structure similar to (1), and thus do not require inversion. Earlier authors proposed a generalization to the Schapery model of the relaxation type: a(t) = E0g0 (Ł) + Xg 1 i (e)/exp[-Ai (Ci -0^^di i 0 (2). J ai(e) J ai(e) Here Ç(t) and f(|) denote strain-reduced functions. Although the functions in Equation (2) might be chosen from among the wide range of function types (polynomial, exponential, etc.), the following forms seem to be preferable for a description of the experimental data for many plastics: E0-g0{e ) = [Ee + D0e^\ (3) g1 i[e] = exp[/?ie] (4) g2i [e] = D1 iŁai (5) a[e\= exp[/ie] (6). Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 E 1 I_______ E n V2 ri 'n Fig. 2. Dashpot-spring schematics for the generalized Maxwell model 2 MODEL CALIBRATION The relaxation-type models, like Equation (2), usually require material characterization based on relaxation tests, which is more complex and expensive than creep experiments. This is because Equation (2) seamlessly downgrades to an analytical approximation of a relaxation curve: s[ek, t ] = E0g0[ek] + +^exp(-l i t ai[ek]) gi 1[ek]gi2[ek] (7), i when the strain history is kept constant: e (t) = const. In order to overcome the impediment of relaxation testing, a special calibration procedure has been developed (Skrypnyk & Spoormaker (2000)) that makes it possible to estimate the parameters in Equation (2) based on creep tests. The procedure is based on the minimization of the relative deviation between the experimental data (creep-recovery tests) and the corresponding model prediction. This relative deviation can be considered as a functional that depends on the material functions (3) to (6). In other words, it may be determined as a function of the model parameters I{Di, a, ß, y} and, therefore, the parameter identification procedure involves the problem of multi-variables function minimization. In order to accomplish this, however, the simulations of creep-recovery tests are necessary. Unlike with the creep-type model, it cannot be accomplished in analytic form, but must be carried out numerically. Since the numerical scheme for calculations of an arbitrary loading history of visco-elastic materials using hereditary integrals, the so-called Henriksen procedure (1984), is already well established and features a high rate of calculations, ę[ffk,t ] -Ł[cyk,t; { Di,rç,A,yJ di (8) it can be used within the parameter-identification procedure (for multiple evaluations). The above procedure itself will be described in the next paragraph, concerning the model implementation into the FEM program. Here, let us simply assume that e [ak, t;{Di, ai, ßi, y}] are the results of calculations of the creep-relaxation curve for the {Di, a, ß, 7} set of model parameters and the stress level ak. The correspondent test data would be denoted as e [ok, t]. Then, the functional of the total error can be written as follows: I{Di ,oi,ß,yi}= Niend ka0 Ł]Pk, t N is the number of loading levels o<7k. The minimization algorithm developed is based on Powell’s method, described in Numerical Recipes (1989) for multi-variables’ function minimization. For the evaluation of the integral in (8) the Romberg integration procedure described in Numerical Recipes (1989) is deployed. This procedure usually implies an equal distribution of the integration points t over the integration period [0, tend] with a continuously increasing number of integration points, until the preset accuracy of the integration is reached. However, when the normal time scale is used, this integration procedure puts 90% of the integration points within the last time decade. Thus, the minimization procedure would produce a resulting set that describes the material behaviour very accurately only for the last time decade, but may be far from reality for most of the other decades. In order to overcome this obstacle, the integration was performed over the logarithmic time scale: log(t) t end log(t end ) (9). 836 Spoormaker J. - Skrypnyk I. - Heidweiller A. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 In this case, the Romberg procedure distributes the integration points equally over all time decades. The choice of fixed X. -parameters (inverted relaxation times) gives additional advantages for the mini-mization procedure developed: 1. It makes the procedure unconditionally stable (the program runs, until a minimum is found). 2. It decreases the number of parameters (to be estimated) and as a result makes the procedure faster. 3. If the {X = 10- i} set is chosen, it covers the entire creep-recovery period and yields a fair approximation for each time decade. Another point that significantly influences the efficiency of the minimization process is a proper choice of the “starting point” for minimization. This was accomplished in three steps. Firstly, available creep data were approximated by the following relations: ę[ffk,t ] =A 0, k ]+ZA i,k] ( 1~exp i=1 -t /10i (10). [i,k] i (sk ) A combination of the principle of locally separated variables and the Prony series form for time functions provides a good approximation for an arbitrary creep-data set. Next, the representation (10) is used to find a variation of stresses in time for constant strain levels (which were chosen from the range [0, 0.3e ], max where e - is the maximum creep strain measured max in the test): B0(p(t)f0 ¦^^^Bii^vn ( 1—exp (—t/10ijj= const i=1 (11). The non-linear Equation (11) was solved for different time moments t using Brent’s method, Numerical Recipes (1989). The last stage consisted of an approximation of the obtained relaxation dependencies g (t) by a generalized representation similar to (7): a[Łk,t] = C[0, k ] + ]TC[i , k]exp(-t/10i) C[i,k]=Di(ŁkT The obtained parameter set {Di, a} was chosen as the initial for the minimization procedure. The starting va-lues for the parameters {ß, 7} were set to zero. The above procedure makes it possible to start the minimization from 20-60% of the total deviation for an arbitrary set of experimental data. Creep tests on plastics cannot be reproduced to (12). within 3 to 5 % and the errors in the FEA simulation are usually of the order of 5 to 10%. There is, therefore, no reason to strive for a fit better than 3 to 5% of total deviation between the data and the model prediction. 3 IMPLEMENTATION OF THE RELAXATION MODEL INTO FEA MSC.MARC 3.1 Extension of the model to the 3-D formulation Below, the 3-D representation of the developed model is briefly presented. The assumptions are as follows: • the strains are small enough to accept a conventional definition of stress and strain tensors; • the behaviour of the visco-elastic material is, to a certain extent, similar to the elastic behaviour; • the material is compressible and originally isotropic; • the deviatoric and hydrostatic part of the deformation process are completely uncoupled; • the rate of viscous flow is proportional to the deviatoric part of the deformation. As a result, the following 3-D representation is proposed: cr = M \v M Eg1 i [l]/exp[-Ai (C-C') ö[g2i[e]e 9( (13) di C{t^I^Ü ( t ). (14). Here, the functions in Equation (13) are modified in such a way that for the case of uniaxial loading: fel=g 0 fęle g 1 i [e] = g 1 i [e]e (15). The following vector and matrix notations are used above: xy Tyz TzxT (16) xy ^lyz 7zxT (17) ° = { Vxx Vyy ^ - [ ^xx ^ yy ^zz M\v\ = \ 1-2n 2n if i = j and i, j < 3 if i ^ j and i, j < 3 1-2*/ (18). 2, if i = jand i,j> 3 0, if i ^ j and i, j > 3 n 0 0 Predvidevanje nelinearnega lezenja - Prediction of the Nonlinear Creep Deformation 837 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 3.2 Main elements of the numerical algorithm for model implementation into the FEA code The MSC.MARC package was chosen for the implementation of the relaxation model (13)-(14), because of its open structure (which enables easy access to the variables, like stresses, strains, time increment, etc. via user subroutines), and extensive features to handle the geometrically and physically non-linear problems. The constitutive relations in MSC.MARC can be modelled using the capabilities of the user subroutine HYPELA2. As mentioned above, the model should be presented in incremental form as follows: Aâ = DAë + G(e,,At, (19). Here, D - is Jacobian D = 9<7ij/oŁi ; G denotes the vector of stress increment due to the variation of the internal parameters 9 i . Therefore, the relation (13) has to be rewritten in an incremental form. For sufficiently regular functions g2i(e), such that d2\g2i\Ł)Ł\ldt2 «1, the hereditary integral with the exponential kernel functions: 9i{t) = Jexp ( -\i{Ç-t) ) ^ d (20) 0 can be calculated recurrently, following the Henriksen scheme (Henriksen (1984)): 9i(t) = exp(_ AiAt)0i(t - At)- A [gi(Ł jejT^At) where: r(A0 1-exp(-zlAt) lAt (21), (22). To write the constitutive Equation (13) in the incremental form (19), the total differential of the right-hand side of Equation (13) has to be derived. Zhang (1995) published similar derivations for the Schapery model (Schapery (1969)). As a result, a quite cumbersome expression with low convergence was obtained. At the same time, certain assumptions can considerably simplify the resulting relations with only a small loss of accuracy (Lai (1995)). If we assume that the effective strain is constant within the time increment: de ---- de 0 (23) Aa = M[^}g0(ę)Ae + n + M h ] ]T [g1 i (ę) (exp (-Ai At) - 1)9i {t-At) (24) + A[g2i(ę)Aë]ri(At)] and expressions for the Jacobian D and the vector G can be derived: n D = M[I/0]g0(ę) + M[I/1]J)g2i(ę)ri(At) i=1 G = M[^]^g1 i(ę)(exp(-\At)-1)ei(t-At) (25). i=1 If necessary, the reduced time increment can be estimated as follows: AC= At a(e(t-At) (26). the stress increment can be represented as follows: The above presented numerical scheme (24) to (26) is recurrent. To estimate the stress field at the end of the current time increment, only the data for the strain field e and the internal parameters 9 from the previous step are needed. Depending on whether the effective strain in relations (24) to (26) is estimated for the beginning or for the end of an increment, the numerical scheme is explicit or implicit. The STATE VAR option is used to store an array of 9 data. In this case the FEA code handles the 9 array itself, which enables restarting and adaptive meshing technology to be used together with the non-linear visco-elasticity material model. 4 EXPERIMENTAL STUDIES The experimental part of the study consisted of: • measuring the creep of the air-inlet material: polypropylene (PP) with 40% talc; • designing a test fixture and measuring the deflection of the middle part of the upper beam of the air inlet. 4.1 Results from the creep experiments Two sets of creep experiments were performed. The first set was carried out on “dog-bone” specimens (5 cm in length), deploying a large tensile-test set-up. The MTS extensometer (Fig. 3) has a range of 10% and this is a rather high range for the lowest strain levels tested. The creep data obtained from the first test set was fitted to model (2) to (6) using the i=1 838 Spoormaker J. - Skrypnyk I. - Heidweiller A. Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 23 C 26 MPal 20 MPa 23 MPa . ¦ 1 16.5 MPa ^^ V /18.9 MPa^ 10 ^rrS^ ^-*"""^^ ^^^^^^ 10 MPa ^-"~"*^ 5 MPa ^^¦H^ _^-»^v_^ ^•^^^^^ ¦ - experiment ------- model prediction 10 Fig. 3. Creep/recovery measurements set-up with “dog-bone” specimen and MTS extensometer calibration procedure described above (see paragraph 3). The test results and the model predictions are depicted in Fig. 4. The agreement between the experimental and modeling results for this test set was rather poor. It seemed that the reason for such a disagreement is erroneous test data obtained for the stress level of 20 MPa and, possibly, some errors in the testing methodology at the initial stages of creep. This initial conclusion, however, turned out to be not correct. One should keep in mind that the calibration procedures like the one described in paragraph 3 are based on the optimization of the total error (Equation (8)), i.e., the deviation between the data and the prediction at all times and load levels. If the test data are rather irregular, the obtained prediction results might not necessarily deviate in the areas where the data are erroneous, but where the data are correct. The latter happens because the model parameters (3) to (6) are time, s Fig. 4. Creep-test data and model predictions for PP with 40% talc estimated for all sets of test data and, therefore, influence the modeling results at all levels. As it proved to be in this case, the test data were unreliable at low stress levels, not at 20 MPa. This becomes obvious when plotting (Fig. 5) the creep moduli for the three lowest stress levels (5, 10 and 16.5 MPa). If we assume that at this stress range the material behaves as linear visco-elastic, the creep moduli should be almost the same. The creep moduli for 16.5MPa are neither similar to those for 5MPa and 10MPa nor are they completely different. Nevertheless, the initial experimental study (Fig. 4) made it possible to obtain approximate model parameters for the model (2) to (6). Using this set of parameters an initial FEA simulation of the air inlet was performed. The FEA mesh and the applied boundary conditions are depicted in Fig. 14. A sketch of the typical results is given in Fig. 6. From the FEA results it became clear that the air inlet does not experience stress levels higher Fig. 5. Creep moduli for PP with 40% talc for three loading levels Fig. 6. Distribution of Von Misses stresses of the air inlet loaded at the middle of the upper beam 100 10 10 10 10 10 10 3000 2000 1000 100 101 102 103 104 105 time (s) Predvidevanje nelinearnega lezenja - Prediction of the Nonlinear Creep Deformation 839 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 Creep modulus ^ ^5 "--- -_.^ "^-stS^-fc ^^ -^ ^-* 100 101 102 103 104 105 Fig. 7. Creep-test data for PP with 40% talc for the loading levels of up to 15 MPa, measured using long specimens than 15 MPa. At the same time, an initial set of creep tests were performed for the levels up to 26 MPa and the data for lower stress levels (which are of interest) appeared to be unreliable. Therefore, the creep experiments were performed again for the three lowest levels (5MPa, 10MPa and 15MPa) using relatively large specimens. The resulting creep curves are shown in Fig. 7 and the creep moduli are given in Fig. 8. The latter plot proves that the material remains in the area of linear visco-elasticity for the stress levels of 5MPa and 10MPa, but becomes non-linear visco-elastic above that range of loading. The creep-test results from the second test were deployed to calibrate the model (2) to (6) using the procedure described in paragraph 3. The results of the model calibration are given in Fig. 9. The resulting set of model parameters for the model (2) to (6) are given in Table 1. Because of the small number of creep curves (as well as the shortened time of the creep tests), Fig. 8. Creep moduli recalculated from Figure 7 1000 time, s Fig. 9. Creep-test data and model predictions for PP with 40% talc for the loading ranges of up to 15 MPa Table 1. Set of model parameters for PP with 40% talc (at 23 °C) Term, i E 0 0.27095427104 D i a i ßi Ji Ai 1 0.31352655103 0.88982058 2.5333602 -0.66960937102 IO2 2 0.3819852M02 0.86928058 0.76228033103 -0.15235806104 IO3 3 0.26457767103 0.78604436 -0.51668278102 -0.57895746103 IO4 4 0.68070737104 1.3374238 -0.27413947103 -0.330784-IO3 IO5 840 Spoormaker J. - Skrypnyk I. - Heidweiller A. 0.010 0.008 0.006 0.004 0.002 0.000 101 time (s) 10 100 10000 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 Fig. 10. Experimental setup with the air inlet the fitting results do deviate from experimental data, especially for longer testing times. The integral deviation between the fitting results and the test data, however, is only 6.6%, which can be considered satisfactory for engineering applications. Further improvement of the fitting results was considered to be an overfitting, since in the process of achieving it some parameters of the model got too high values, which we consider as “non-physical” and such that would disturb the numerical robustness of the model. 4.2 Measuring the creep deflection of the air inlet The air inlet was fixed to a testing setup, as depicted in Fig. 10. A special fixture was developed in order to position the upper plane in the horizontal plane, as depicted. The upper beam was loaded in the middle of the upper beam with a constant load of 60N for 10 hours (Fig. 11). The upper beam in the middle is not completely flat. Therefore, a cylinder was used to exert the load and allow some twisting of this beam. The vertical deflection of the beam at the point of loading was recorded (to be compared with the FEA predictions). 4.3 FEA verification of the developed approach The accuracy of the created approach to the prediction of the long-term behavior of the air inlet was verified using MSC.MARC 2000. In order to achieve a high accuracy for the FEA prediction the complete air inlet was meshed with shell elements (13,000 elements with 7 layers), see Fig. 14. The Fig. 11. Close-up photo of loaded upper beam of the air inlet out-of-core solver procedure was applied (for this MSC.MARC option ELSTO was deployed). The accuracy requested for the FEA calculation was 5%. A higher accuracy would significantly increase the calculation time, but was considered not necessary, since the model has been fitted with only a 6.6% accuracy. The results of the FEA predictions for the middle of the upper beam were compared to the experimental data (Fig.12). The deviation between the test data and the prediction was less than 5.5% for the entire modeling period. This should be considered as satisfactory, since the accuracy of the constitutive model was around 6.6% and the accuracy of the FEA analysis was 5%. S**1 HlAAAi - experimental data - FEA simulation 0 5000 10000 15000 20000 25000 30000 35000 40000 Time, s Fig. 12. Long-term deformation of the upper beam of the inlet at the middle. Comparison of experimental results vs. FEA prediction. -0.5 Predvidevanje nelinearnega lezenja - Prediction of the Nonlinear Creep Deformation 841 Strojniški vestnik - Journal of Mechanical Engineering 53(2007)12, 834-843 Fig. 13. Finite-element model 5 CONCLUSIONS 1. The presented example of a design problem related to a long-term constraint demonstrated that the long-term deformation of complex 3-D plastic products can be successfully predicted by using commercially available FEA software and modern constitutive models for non-linear visco-elasticity. 2. For successful FEA modeling it is critical to have reliable experimental data for a range of parameters (stresses, strains, time period and temperature), in which the designed product is supposed to operate. For non-linear plastic materials (which are rather typical) at least two iterations between the tests and the FEA simulations are usually required to reach an accurate predictive capability. 6 NOMENCLATURE a, ßi , y, D coefficient in material functions gi(e) e i strain (one-dimensional formulation) Ł x. V % a, a(t) G e. i A ] [i,k B experimental creep data for stress level