domain: phys sde

text:
Abdul-Hameed Hemin

Roguet Eléonore

Brusselle-Dupend Nadège

Boulharts Habiba

Cangémi Laurent

Explicit integration methods for constitutive equations of a mean-stress dependent elastoviscoplastic model: Impact on structural finite element analyses

Keywords: Viscoplasticity, Polymer, Deformation behavior, Structure analysis, Constitutive model, Integration method

The strong dependent behavior of semi-crystalline polymers can lead to the use of simplified material laws in Finite Element structural calculations for reasons of robustness to the detriment of the quantitative response of the models. This work focuses on numerical integration methods as a solution to overcome the possible convergence and robustness limitations of mean-stress dependent elastoviscoplastic material laws, typical of the semi-crystalline polymers mechanical behavior.

What is proposed here is a rational application of three explicit integration methods (fourth and second order Rung-Kutta method, a hybrid schema between Runge-Kutta and Euler method) in engineering structural calculations, which provide a reliable solution for constitutive models of semi-crystalline polymer. These methods are examined for structure creep test and tensile test, in comparison with experimental data. The investigations have been done in terms of the stability toward convergence, the accuracy of results, the plastic consistency, and CPU time efficiency. This work, proposes an easy implementation of integration methods in any computational Finite Element code. It also provides a flexible modular implementation which is applicable to any different constitutive equations. An integration step sub-division technique is recommended. It is a powerful technique to improve the convergence of solution and accuracy of result by damping oscillation around stress Gauss point integration solution. The results obtained illustrate the effect of numerical integration schemas on structural analysis and provide an insight to select suitable method.

Introduction

Finite Element Analysis engineers appear as the end users of the material laws developed by the solid mechanics community. Easy calibration and robustness of material laws are required for designing structural applications and predicting their service life from Finite Element calculations. But the mechanical behavior of semicrystalline polymers (SCPs) is so strongly dependent (time, mean stress, temperature) that simplified material laws are generally favored regarding the convergence of the structural calculations to the detriment of their response. This work focuses on numerical integration methods as a solution to overcome the possible convergence and robustness limitations of mean-stress dependent elastoviscoplastic material laws, typical of the semi-crystalline polymers mechanical behavior.

The contributions of semi-crystalline polymers (SCPs) in complex industrial structure, need a reliable constitutive model since their mechanical behavior is time-dependent and strongly nonlinear once undergo large deformation. An accurate computing of mechanical response is necessary to predict lifetime of structures. The sophistication of the constitutive model depends on how to interpret phenomena occuring during deformation (including elasticity, plasticity, relaxation, and creep) with respect to internal (structure of polymer at different scales) and external (temperature, rate of deformation) factors, besides, numerical approaches and mathematical framework employed are important to attain accurate mechanical response.

Over decades, many continuum-level constitutive models have emerged to predict the mechanical behavior of SCPs, taking into account nonlinearities and rate-dependence, as reviewed in many works . However, the work of Brusselle and Cangemi [START_REF] Brusselle-Dupend | A two-phase model for the mechanical behaviour of semicrystalline polymers. Part I: Large strains multiaxial validation on HDPE[END_REF] has been considered in with work, for its capacity of covering main features of SCPs such as time and mean stress dependence, volume change, cavitation development, strain hardening under large deformation.

A complex mechanical behavior of SCPs requires an elaborate constitutive model incorporating all aspects of mechanical behavior. The choice of the numerical integration method is thus important.

Importance of integration methods for the constitutive model:

The time-dependent elasto-viscoplastic equations can be considered as a system of discontinues, nonlinear Differential Algebraic Equation (DAE), due to the yield condition and evolution of hardening variables with plastic strain [START_REF] Büttner | Runge-Kutta methods in elastoplasticity[END_REF]. Therefore, the mathematical structure of the constitutive law becomes very complex and mathematically stiff which causes a key problem for numerical integration methods to update the stress Gauss point and evolution of internal state variables.

The difficulties associated with a reliable implementation of integration schemas have limited the employment of complex model [START_REF] Hiley | Explicit Runge-Kutta methods for the integration of rate-type constitutive equations[END_REF][START_REF] Liu | A scheme of automatic stress-updating on yield surfaces for a class of elastoplastic models[END_REF][START_REF] Scalet | Computational Methods for Elastoplasticity: An Overview of Conventional and Less-Conventional Approaches[END_REF].

Modern constitutive laws aim to improve the quality of engineering structure calculations by involving complex models to capture the mechanical behavior of materials. Providing a reliable tool does not only focus on a reliable constitutive model based on a physical description, it also enhances the numerical integration tool for improving the constitutive models in structural computation [START_REF] Besson | Continuum Models of Ductile Fracture: A Review[END_REF]. Accurate result of complex models in simulations, conversely, depends on the robustness of solution of integration methods for the stress Gauss point integration. Therefore, constitutive models need to be provided with a robust numerical integration method to deliver an accurate solution.

In structural computation, the stress updating at each Gauss point for specific strain increment in a nonlinear Finite Element is a key solution. This update occurs through integrating constitutive model equations. The constitutive model solution and internal state variables update need to reevaluate several thousand times during the structure calculation, depending on the complexity of structure design. Thus, the computational cost changes considerably with the complexity of the constitutive model and with the integration method schemas employed.

In consequent, the accuracy and efficiency of variable update in each load step increment affect the entire structure solution [START_REF] Rezaiee-Pajand | Accurate and approximate integrations of Drucker-Prager plasticity with linear isotropic and kinematic hardening[END_REF].

Criteria for choosing an integration schema: The choice of numerical integration methods in practice is important, however, an efficient integration method relies on: I) the numerical stability of the method; II) the accuracy of solution; III) the incremental plastic consistency; IV) the adaptability to implement in constitutive model equations; V) a low CPU time cost [START_REF] Hiley | Explicit Runge-Kutta methods for the integration of rate-type constitutive equations[END_REF][START_REF] Scalet | Computational Methods for Elastoplasticity: An Overview of Conventional and Less-Conventional Approaches[END_REF]. The first two criteria are vital to satisfy the convergence of the numerical solution. To limit the accumulative error produced as a consequence of increment stepsize, a step sub-division procedure has been proposed which allows the increment stepsize to be divided into sub-increments and to reduce the error. This technique was proposed by Sloan [START_REF] Sloan | Substepping schemes for the numerical integration of elastoplastic stress-strain relations[END_REF]. The third criterion satisfies the keeping stress updated value at a time (t n+1 ) on the subsequent yield surface during incremental plastic deformation. The fourth condition is a practical limitation: it does not change the results but represents major criterion for engineering use. Finally, computing time cost is also a significant factor in structure calculation within a Finite Element package.

A lot of works have been achieved to apply integration methods to solve Initial Value Problem for DAE.

Meanwhile, few works have been done to implement numerical integration methods to update stress and its corresponding internal variables to capture the evolution of elasto-viscoplastic behavior of SCPs materials under large deformation within a Finite Element framework. In this contribution, the modular implementations of different integrations methods are proposed for an elastic-viscoplastic constitutive model of SCPs, this manner of implementation of integration schemas can be applied to any constitutive model. The varieties of integration methods implemented in this work provide an insight into a suitable integration tool for a different structure calculation loading conditions.

In this work, explicit numerical integration schemas are employed. The explicit methods are straightforward, not need to build a matrix of algebraic equations system or iterative loop to obtain integrated variables, solutions are divided in time (t), the solution in current step time (t n+1 ) is computed directly from previous step time (t n ).

The explicit integration schemas that have been employed in this study are a second-order explicit Runge-Kutta, a fourth-order explicit Runge-Kutta method, and a Hybrid method (between second-order Runge-Kutta and explicit forward Euler method). These methods have been chosen due to the possibility to apply to elastoplasticity constitutive equations and their enhanced property comparing with other methods in terms of criteria explained above [START_REF] Büttner | Runge-Kutta methods in elastoplasticity[END_REF][START_REF] Scalet | Computational Methods for Elastoplasticity: An Overview of Conventional and Less-Conventional Approaches[END_REF][START_REF] Ashour | A new very simple explicit method for the integration of midly stiff ordinary differential equations[END_REF][START_REF] Hairer | Solving ordinary differential equations I: Nonstiff problems[END_REF]. The very high order Runge-Kutta method is not recommended due to computing efficiency vs. obtained accuracy and stability. This work is organized as follows. Section 2 provides details of the two-phase elastic-viscoplastic constitutive model considered for SCPs. Then section 3 details the modular implementation of explicit numerical integration methods employed in this work into constitutive model within a Finite Element code. Section 4 presents experimental database and comparison with structural calculations. Finally, the work concluded in Section 5.

Two-phase elasto-viscoplastic constitutive model for semi-crystalline polymers

In this work, an elasto-viscoplastic constitutive model based on the thermodynamic framework of porous media is used. It has been originally proposed in [START_REF] Cangémi | A two-phase model for the mechanical behaviour of semicrystalline polymers[END_REF] and extended to large strains in [START_REF] Brusselle-Dupend | A two-phase model for the mechanical behaviour of semicrystalline polymers. Part I: Large strains multiaxial validation on HDPE[END_REF]. The proposed two-phase model intends to predict the large strains tridimensional behavior of SCPs, by also considering the other features of the mechanical behavior met during the life of a structure with plastic components (mean stress sensitivity, volume change evolution, short-term and the long-term viscous behavior). Considering their microstructure, semi-crystalline polymers can be represented as two-phase materials at mesoscopic scale. The first phase is a solid network corresponding to the crystalline phase and tied amorphous zones of the polymer, and the second one is a "fluid-like" phase corresponding to the very soft amorphous phase (free amorphous). The solid network can be compared to an open "porous" skeleton which is saturated by free amorphous phase. Moreover, a plastic yield surface inspired by porous-like materials (clays for instance) has been selected to describe the sensitivity of the mechanical behavior of SCPs to mean stress revealed for several decades in these materials. A nonassociated viscoplastic law also permits to catch both the plastic volume change and the viscous response of SCPs.

Kinematical fra mework

The solid network deformation of continuum media, described by Lee [START_REF] Lee | Elastic-Plastic Deformation at Finite Strains[END_REF], introduces the unstressed configuration of the representative elementary volume  d , represented by R d , which is obtained by unloading from the deformed configuration t d . This concept leads to the following decomposition of the deformation gradient:

(1)

(vp) (e) F F F .  , F is considered the deformation gradient transforming R d into t d . (e)
F is the elastic part of F , and (vp) F

the viscoplastic part of F transforming 0  d (initial configuration of  d ) into R d .
For large strain, the strain rate tensor can be expanded as below:

(2)

s s s ) ( ) ( ) ( 1 1 1 1             (e) (vp) (vp) (e) (e) (e) F F F F F F F F D   
. However, we can consider for SCPs that the elastic part of the deformation remains small at large strains, which leads to the following approximation:

  s ) ( 1       (vp) (vp) (e) vp (e) F F ε D D D   .

2.2

Constitutive equations of the two -phase model

Mass balance and coupling equations

A porous medium saturated by a fluid-like phase has been considered by the authors to introduce a natural coupling between the amorphous phase and the solid skeleton. This mobile and soft phase is supposed to be free and can move through the skeleton. Then, a mechanical loading imposes a differential movement of the saturating fluid with the skeleton, and the volume of matter  d , which follows the movement of its skeleton, is able to exchange fluid matter with the exterior. The balance of the fluid mass contained in the volume  d thus holds:

(5)

        0 ) ( , ,     d m m d t t m a a t a X X X x
where a m ~ is fluid mass variation between the reference and the actual state. In the porous media, the fluid mass per unit volume depends on the volumic mass a  and on the connected porosity a  , can be formulated as: It allows formulating the plastic porosity as following:

(9) 0 (vp) (vp) ~a a m   
, where 0 a  is the volumic mass at the initial state. Furthermore, the viscoplastic porosity (vp)  can be related to the total viscoplastic deformation of the two-phase material by assuming a linear dependence with β . The second order tensor β introduces an anisotropy effect linked to a non-isotropic initial microstructure. This tensor reduces to a scalar parameter once the material is isotropic.

The previous considerations about fluid movements inside the porous crystal skeleton lay the foundation for a general framework which could be able to deal with mass exchanges of gas or liquids in the polymer. However, the exact formulation of this phenomenon should introduce the notion of molecular diffusion through the free amorphous phase. An attempt to construct a more precise approach can be found in [START_REF] Baudet | A damage model for the blistering of polyvinylidene fluoride subjected to carbon dioxide decompression[END_REF]. In the context of this work, there are no exchanges supposed between the constituent sub-volumes  d and the external environment for SCPs, as a result, it considered a closed system, therefore 0 ~ a m . Then, the mass variation of the plastic part per unit volume is results of deformation of the soft phase (free amorphous macromolecules) imposed by the deformation of the solid network. From this point of view, it is possible to consider the intrinsic evolution of volumic mass of the soft phase during plastic deformation of the skeleton.

For infinitesimal strains, classical thermodynamics extended to porous continuous media and a quadratic choice of the free energy lead to the following state laws and dissipation inequality: The Eqs. [START_REF] Hadj Hamouda | Viscoplastic behaviour of a medium density polyethylene (MDPE): Constitutive equations based on double nonlinear deformation model[END_REF] and [START_REF] Brusselle-Dupend | A two-phase model for the mechanical behaviour of semicrystalline polymers. Part I: Large strains multiaxial validation on HDPE[END_REF] describe the interactions between the solid network and the embedded soft phase as two springs acting in parallel, so that B, M, and β appear as coupling parameters, where M is linked to the compressibility of the free amorphous phase within the solid network. The tensors B and β allow to represent the evolution of the intensity of the macromolecular interactions between the hard phase and the soft phase through the tied amorphous zone.

The elastic tensor of the SCPs system can be deduced from

B B C C s    M
. In Eq. ( 12), the terms J A are the affinities of internal variables J  .

Through Eq. ( 12) and free enthalpy equivalent to the pressure acting on free amorphous, an effective viscoplastic stress (vp) s σ of the two-phase material associated with the viscoplastic strain can be presented as:

(

) a P β σ σ (vp) s   . 14
A direct extension of the equations ( 11), ( 13) and ( 14) can be done from the previous kinematic framework at finite strains:

(15)

a a p G J 0 - B σ σ s   so that     (vp) s D D : C σ s     , a p a g J G 
, and

F det  p J .    s
σ is the so-called Jaumann derivative, which can be interpreted as the time derivative of the Cauchy stress tensor σ with respect to a rotating frame. D  and

(vp) D 
are respectively the strain rate tensor and its viscoplastic part.

The second coupling equation takes the following form:

(16)       (vp) vp E : B β D : B β E : B D : B              ) ( 0 M M M a a a G G  ,
In Eq.( 16),

V E ln

 where V is the left stretch tensor, and the viscoplastic strain

   t 0 d (vp) (vp) D E where s ) ( 1    (vp) (vp) (vp) F F D  .
Finally, the viscoplastic effective stress becomes:

(17) a a p G J 0  β σ σ (vp) s  
, and can also be expressed as follows:

(18) a a G 0 ) (  B β σ σ s s (vp)   

Viscoplastic modeling of the solid phase

The sensitivity to hydrostatic pressure and the volume change evolution generally observed in SCPs has led Cangémi and Meimon [START_REF] Cangémi | A two-phase model for the mechanical behaviour of semicrystalline polymers[END_REF] to consider SCPs as porous-like materials similarly to clay materials. In these materials, the volume change is associated with hydrostatic pressure. SCPs also show this mechanical behavior feature, and [START_REF] Cangémi | A two-phase model for the mechanical behaviour of semicrystalline polymers[END_REF] admitted to result from the difference in compressibility between the solid network and the free amorphous inside. A generalized "cam-clay" surface has thus naturally been considered as an appropriate plastic yield surface (f ) intended to the solid network of SCPs: At s b = 0, the Eq. ( 19) is reducing to a Drucker-Prager like criterion. The evolution of the yield surface also depends on the functions r and s p . They correspond to dual affinities J A associated with the following choice of internal state variables J  , which are linked to hardening and to the strain energy stored in the material:

(19)   s s eq p g f s     with   r p p p g c s s                  
(20)     (vp) vol (vp) d J J E E , 2 , 1    with  d E t (vp) d ) ( I 0 2   (vp) D and  d E t (vp) vol ) ( I 0 1   (vp)

D

Note that for simplification of writing, the following notations

D

thereafter.

(vp) d E and (vp)

vol E are the cumulative viscoplastic strain according to Von Mises, and to the volume viscoplastic strain respectively.

(vp) d E and r introduce a classical isotropic hardening effect as a consequence of cumulative plastic strain, whereas (vp) vol E and c p introduce a volume-dependent effect with plastic volume strain.

The stored energy contribution of non-reversible phenomena occurring in the material outcomes of the following equations for r and c p :

(21) 0 ( 0 r d vp) d d (vp) d E a E       (22)   ) s 0 c exp p vp vol c E p (    0 d  is the initial size of the yield surface,
a is a parameter controlling the rate of isotropic hardening and affects the transition domain between zero stress and perfect yielding, before the hardening stage at large strains, 0 c p refers to residual confining stresses from process conditions related to the morphology of the material microstructure. It can evolve with texture structure evolution. Finally, s  is a control parameter.

A viscoplastic flow rule, especially a non-associated flow rule has been naturally proposed by [START_REF] Cangémi | A two-phase model for the mechanical behaviour of semicrystalline polymers[END_REF] because of the viscous effects and the volume strain evolution accompanying plastic yielding in SCPs in general:

(23)     (vp) s (vp) σ Λ D f    1
where  acts as a long relaxation time and the function,  

           s k f f Φ sinh (25)     1 σ dev σ Λ (vp) s (vp) s               s s s c 3 2 3 s s eq s eq p  
In Eq. ( 24), s k is a parameter that controls the rate sensitivity at short times, the sinh function is based on the Eyring's molecular approach to the viscosity phenomena within the framework of thermal activation theory [41].

In Eq. ( 25), the parameter s  controls the intensity of the volume strain. The sign of the volume strain depends on the evolution of the yield surface in the (  s eq , s p ) plane.

Most of the parameters of the material law concern the plastic yield surface making it possible to grasp the dependence of the mechanical behavior on mean stress in order to cover a wide range of multiaxial stress states possibly encountered in structural calculations. In addition, parameters relating to the viscous behavior and to the non-associated nature of the viscoplastic flow rule are involved.

Features at large strains

Whatever the external loading, the progressive orientation of the crystalline and amorphous parts during large strains accompanies the progressive destruction of the initial isotropic spherulitic mesostructure into an oriented crystalline structure. This can lead to two competing phenomena at large strains: macroscopic stress hardening by the evolution of the crystalline texture of the SCP and cavitation mechanisms in the amorphous phase to account for the macroscopic loading, as suggested in [START_REF] Brusselle-Dupend | A two-phase model for the mechanical behaviour of semicrystalline polymers. Part I: Large strains multiaxial validation on HDPE[END_REF].

The model is based on the interpretation of the microstructure evolution at a mesoscopic level considering the evolution of a solid network containing a softer phase, like in soil-like materials. An interaction relation between the two phases is obtained from coupling equations between the solid network and the free amorphous inside through B (Eq. 10). A specific yield surface is moreover associated with the solid network, which is supposed to be at the origin of the viscoplastic evolution of the material. The specific hardening behavior at large strains is supposed to result from the evolution at a mesoscopic scale of both the intensity of the interactions between the two phases, supposed isotropic initially (Eq. 26), and the residual confining pressure (Eq. 27). The reader can refer to [START_REF] Brusselle-Dupend | A two-phase model for the mechanical behaviour of semicrystalline polymers. Part I: Large strains multiaxial validation on HDPE[END_REF] for details. Finally, a global energy criterion has been proposed to catch a possible process of cavities nucleation in the free amorphous phase of SCPs. The consequence of the cavitation process is traduced in the model by a sudden decrease of the compressibility modulus M which involves an increase of the volume plastic strain and a global stress softening.

(26)       3 , 1 , j i j i ij e e b B with n 1 0 k tanh L b                     Ε (vp) ij ij b for j i  , and 0  ij b for j i  .
In the present work, B and 0 c p have been considered constant and the energy cavitation criterion not activated for the sake of simplicity, also because it wasn't the scope of the present work. However, applications of the models for large strains have been achieved in [START_REF] Brusselle-Dupend | A two-phase model for the mechanical behaviour of semicrystalline polymers. Part I: Large strains multiaxial validation on HDPE[END_REF] and more recently in [START_REF]Diffuso-mechanical modelling of semicrystalline polymer under water pressure[END_REF][START_REF]Fragilisation mécanique du polyamide 11 en condition d'hydrolyse : mécanismes et modélisation, French Thesis IFPEN[END_REF] to deal with some mechanical specificities encountered with SCPs.

A summary of the whole parameters involved in the constitutive model is finally presented in Table 1.

Table 1. Parameters of the material law.

Elasticity parameters of the solid network

Implementation of numerical integration schemas

The constitutive elastic-viscoplastic model presented in section 2 has been implemented into an UMAT (User MATerial subroutine in the Abaqus Software®). In a numerical calculation step, the mechanical subroutine is called at the beginning of each increment. Then, the stress tensor through the Eq. ( 18) and their corresponding internal variables are calculated according to a given numerical integration process. For each element and each increment, the updated stress tensor and its corresponding internal variables are stored at the end of the increment, and is recalled at the beginning of the next one.

The mechanical constitutive behavior established as a set of rate equations, has the general form of an initial value problem. In this constitutive model, the Eq. ( 18) is the main equation that required to be updated, indeed with its state variables. Each time increment ( ) is divided into J increments to become a sub-increment ( ) or ( ), where ( ∑ ⁄ ) and ( ⁄ . The J parameter is an integer value equal or greater than one.

The plastic consistency condition depends on stepsize increments: enlarging the stepsize will increase the error during the entire integration process. Thus in with work, the step sub-division technique is employed, which divides each increment stepsize into sub-increments. The number of sub-increments (J) is controlled manually (as an input value), in order to avoid extra calculation and increase the CPU time.

Explicit numerical integration methods

An explicit integration method is a direct calculation of dependent variables at a time step (t i+1 ) that can be computed in term of known quantities at the previous time step (t i ). In this section, the employed integration methods are presented briefly and the details of their implantations are provided.

Forward Euler and Runge-Kutta methods implementation

For the initial value problem of a general nonlinear system of first-order ordinary differential equations, such as ( ) at ( ), the solution can be obtained from: ( ) ∫ ( ) . The forward Euler method, thus can be formulated as:

(28) ( )
The Runge method has modified the Euler method by adding further steps, in order to provide better accuracy.

Afterward, the method was generalized by Kutta. The general form of Runge-Kutta method of s-stage can be written as follow:

The Butcher tables for a s-stage Runge-Kutta order ( ) are shown in Table2, where is stage evaluation, is a step size, and are real coefficients, can be determined as ∑ . The stability function ( ) for explicit Runge-Kutta method is given by [START_REF] Hairer | Solving ordinary differential equations I: Nonstiff problems[END_REF]:

(30) ( )
where z , and is corresponding eigenvalues. The stability of the method is in positive correlation with value of ( ) . Runge-Kutta has a structure of one-step methods, offering a possibility of multistage per step; this led to better accuracy and superior stability in general. The increasing number of the stages per step will increase evaluations of functional at each step, which is costly from a CPU time viewpoint [START_REF] Press | Numerical Recipes in C: The Art of Scientific Computing (Fortran version)[END_REF]. The rigorous mathematical description of modifying the Euler method, expanding stages to upgrade to a higher order, generalizing of the methods, and evaluating the stability functions are out of the scope of this work. Readers are invited to refer to some books [START_REF] Hairer | Solving ordinary differential equations I: Nonstiff problems[END_REF][START_REF] Hairer | Solving ordinary differential equations II: Stiff and Differential-Algebraic Problems[END_REF] for more details.

The explicit Runge-Kutta general form of one stage of evaluation is equivalent to Forward Euler method. The numerical implementation of stress update with corresponding internal variables is shown in Fig. 1. For nonelastic deformation, the time interval is divided into sub-increments (J sub-increments are here considered). The corresponding variables of the constitutive model are recalling at the beginning of each increment, then estimation processes start and initialize the temporary memory of local arrays. In this case, The processes indicated in the box with the red dashed line are the core numerical integration processes, for next numerical integration methods, the only processes in red dashed line box will be modified and repeated (as the other parts are similar for all methods). The equations presented in the red dashed line box are main equations (in this constitutive model) which need to be updated at each Gauss point and at each time-step. This flexible way of numerical integration implementation allows applying the integration method to any constitutive model in future works.

Fig.1.

Flowchart of stress-strain updates implementation in a finite element, according to the forward Euler integration method.

j j j J       s s s σ σ σ  ( t ( ( 1 _ j j j J ) ( t ) ( ( (vp) (vp) (vp) _1 D E E     j j ) ( ) ( 1 _ s s σ σ   Computing Eq. (18): a a j j G 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j vp vol j vp vol j vp vol D J E E ) ( t ) ( ) ( 1 _ ) ( ) ( ) (    j vp d j vp d j vp d D J E E ) ( t ) ( ) ( 1 _ ) ( ) ( ) (   

Second order Runge-Kutta

The Runge-Kutta second order can be obtained by taking the initial derivative at each step to obtain midpoint across the interval. According to the generalized Butcher table for second-order ( ), Runge-Kutta can be reduced as shown in Table 3, and thus, the Eq. ( 29) can be derived as:

(31) ( ) ( )
Table 3. Butcher table for second-order Runge-Kutta method.

There are two evaluation stages in the Runge-Kutta second order method: first evaluation ( ) at the initial point and second evaluation ( ) at the middle of the interval. The process of stress and their corresponding internal variables update for both evaluation stages is carried out as detailed in Fig. 2. The Fig. 2 presents the integration process to evaluate the at the beginning of the interval . Hence the estimation of occurs at the mid-interval, the time interval /J is divided by two : ( )/2J. To obtain variables update at (t i+1 ) the Butcher table (Table 3)

and Eq. ( 31) are used. The variables

Fourth-order Runge-Kutta

For each step in fourth-order Runge-Kutta, the derivative is evaluated four times, one time at the initial point, twice at midpoints, and once at the endpoint. As a result, the fourth order Runge-Kutta method is reduce as shown in below equation [START_REF] Scalet | Computational Methods for Elastoplasticity: An Overview of Conventional and Less-Conventional Approaches[END_REF], and the corresponding Butcher table is depicted in Table 4.

where k 1 is the slope at the beginning of the step, k 2 is the slope at the middle of the step using k 1 , k 3 is the slope at the middle of the step using k 2 , and k 4 is the slope at the end of the step using k 3 .

First estimation ( )

-Stress and strain integration:

j j j J       1 _ 1 _ ( t ( ( s s s σ σ σ  j j j J Δt ) (D ) (E (E (vp) _1 (vp) (vp) _1    j j ) ( ) ( 1 _ s s σ σ  
Computing Eq. ( 18):

a a j j G 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j (vp) vol_ j (vp) vol j (vp) vol_ D J E E ) ( t ) ( ) ( 1 1    j (vp) d_ j (vp) d j (vp) d_ D J E E ) ( t ) ( ) ( 1 1

  

Second estimation ( ) -Stress and strain integration:

j j j J       2 _ 2 _ ( 2 t ( ( s s s σ σ σ  j j j J Δt ) (D ) (E (E (vp) _2 (vp) (vp) _2 2    j j ) ( ) ( 2 _ s s σ σ  
Computing Eq. ( 18):

a a j j G 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j (vp) vol_ j (vp) vol j (vp) vol_ D J E E ) ( 2 t ) ( ) ( 2 2    j (vp) d_ j (vp) d j (vp) d_ D J E E ) ( 2 t ) ( ) ( 2 2

  

Recapitulation of estimated stages according to Eq. ( 31): The Fig. 3 depicts the estimation procedures for all four stages and final estimation according to the Eq. ( 32). The indexes 1, 2, 3 and 4 refer to the first, second, third and fourth estimation respectively.

j j j J         2 _ 1 _ 1 ( 2 t ( ( s s s s σ σ σ σ   ; j j j J ) ( t ) ( ( ) ( ) ( 1 _ ) ( 1 ) ( vp _2 vp vp vp D D E E       a G a j j 0 ) ( ) ( ) ( 1 1  B β σ σ s (vp) s      ; j (vp) vol_ (vp) vol_ j (vp) vol j (vp) vol D D J E E ) ( t ) ( ) ( 2 1 1      j (vp) d_ (vp) d_ j (vp) d j (vp) d D D J E E ) ( t ) ( ) ( 2 1 1     

Fig.3.

Flowchart showing only stress -strain integration processes, according to Runge-Kutta fourth-order.

Hybrid integration method (Forward Eul er and Runge-Kutta second-order)

This Hybrid method is an explicit one, based on both forward Euler and Runge-Kutta second-order methods.

This combination has been proposed by [START_REF] Ashour | A new very simple explicit method for the integration of midly stiff ordinary differential equations[END_REF], in order to improve the stability property. Explicit forward Euler is First estimation ( ) -Stress and strain integration:

j j j J       1 _ 1 _ ( t ( ( s s s σ σ σ  j j j J ) ( t ) ( ( (vp) _1 (vp) (vp) _1 D E E     j j ) ( ) ( 1 _ s s σ σ  
Computing Eq. ( 18):

a a j j G 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j (vp) vol_ j (vp) vol j (vp) vol_ D J E E ) ( t ) ( ) 1 1 (    j (vp) d_ j (vp) d j (vp) d_ D J E E ) ( t ) ( ) 1 1 (   

Fourth estimation ( ) -Stress and strain integration:

j j j J       4 _ 4 _ ( t ( ( s s s σ σ σ  j _ j j J ) ( t ) ( ( ) ( 4 ) ( ) ( vp vp vp _4 D E E     j j ) ( ) ( 4 _ s s σ σ  
Computing Eq. ( 18):

a a j j G 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j (vp) vol_ j (vp) vol j (vp) vol_ D J E E ) ( t ) ( ) 4 4 (    j (vp) d_ j (vp) d j (vp) d_ D J E E ) ( t ) ( ) 4 4 (   
Recapitulation of estimated stages according to Eq. ( 32): ( )

j s j j J           4 _ 3 _ 2 _ 1 _ 1 2 2 ( 6 t ( ( σ σ σ σ σ σ s s s s s     ; j j j J ) 2 2 ( 6 t ) ( ( ) ( 4 _ ) ( ) ( ) ( ) ( 1
) ( vp vp _3 vp _2 vp _1 vp vp D D D D E E         ; a G a j j 0 ) ( ) ( ) ( 1  B β σ σ s (vp) s     ; j (vp) vol_ (vp) vol_ (vp) vol_ (vp) vol_ j (vp) vol j (vp) vol D D D D J E E ) 2 2 ( 6 t ) ( ) ( 4 3 2 1 1        ; j (vp) d_ (vp) d_ (vp) d_ (vp) d_ j (vp) d j (vp) d D D D D J E E ) 2 2 ( 6 t ) ( ) ( 4 3 2 1 1       

Second estimation ( ) -Stress and strain integration:

j j j J       2 _ 2 _ ( 2 t ( ( s s s σ σ σ  j _ j j _ J ) ( 2 t ) ( ( ) ( 2 ) ( ) ( 2 vp vp vp
D E E     j j ) ( ) ( 2 _ s s σ σ  
Computing Eq. ( 18):

a a j j G 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j (vp) vol_ j (vp) vol j (vp) vol_ D J E E ) ( 2 t ) ( ) 2 2 (    j (vp) d_ j (vp) d j (vp) d_ D J E E ) ( 2 t ) ( ) 2 2 (   

Third estimation ( ) -Stress and strain integration:

j j j J       3 _ 3 _ ( 2 t ( ( s s s σ σ σ  j j j J ) ( 2 t ) ( ( (vp) _3 (vp) (vp) _3 D E Ε     j j ) ( ) ( 3 _ s s σ σ  
Computing Eq. ( 18):

a a j j G 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j (vp) vol_ j (vp) vol j (vp) vol_ D J E E ) ( 2 t ) ( ) 3 3 (    j (vp) d_ j (vp) d j (vp) d_ D J E E ) ( t ) ( ) 3 3 (   
not recommended to be employed alone due to the accuracy of the method. This Hybrid method offers an advantage over forward Euler method by being more accurate, and over Runge-Kutta second-order method by having a wider stable zone without adding extra evaluation stages. Thus, the computation time-cost remains the same as the second order Runge-Kutta but with a higher stability. This stability zone depends on a ( ) parameter that manages the importance of the forward Euler or the Runge-Kutta method in this Hybrid method, as shown in Fig. 4. The Eqs. ( 28) and ( 31) can be reformulated as follows;

(33) ( ) ( )= (34) ( ) ( ) { [ ] [ ( ) ]}
The hybrid method thus can be obtained using Eqs. ( 33) and ( 34) as follows:

(35) ( ) ( )( )
( ) {( ) [ ] ( ) [ ( ) ]}
Fig. 4. Stability region of (Forward Euler/Runge-Kutta) Hybrid method [START_REF] Ashour | A new very simple explicit method for the integration of midly stiff ordinary differential equations[END_REF].

where z ,with , the increment stepsize and the corresponding eigenvalue.

Two stages of estimations are required in this Hybrid method; at the first one, the explicit forward Euler method is computed, the result of this first evaluation becomes the input for the second evaluation using Runge-Kutta method. These evaluations are done as explained in the Fig. 5 with respecting the time interval for each evaluation. The ( ) parameter is a percentage ratio between forward Euler and the second-order Runge-Kutta methods. It is set equal to 0.74 in this work for stability reason(Fig. 4).

Fig.5.

Flowchart showing only stress-strain integration processes, according to Hybrid method.

Experimental verification and evaluation of the integration schemas performance

The constitutive model and numerical integration methods presented respectively in the Section 2 and the Section 3 have been coded via FORTRAN 90 and implemented in a commercial finite elements package ABAQUS TM /Standard (2014) through a user subroutine UMAT [START_REF]ABAQUS/Standard, Version 6.14[END_REF].

Two mechanical loading cases have been simulating with the mechanical behavior law calibrated on a polyethylene industrial grade. The results obtained successively with the three explicit integration methods have been analyzing and comparing in order to evaluate their performance. The integration methods have been examined according to the stability of numerical integration method, the accuracy of the obtained solution, the plastic consistency criterion, and the CPU time efficiency.

The first case of simulation deals with a classic uniaxial tensile test. The second one, a creep test, concerns a configuration related to a polymer structure in an engineering context.

Uniaxial tensile test

The typical results of experimental tensile stress tests performed on a polyethylene industrial grade at room temperature 23°C are depicted in Fig. 6. The tensile force F was measured and converted to true axial stress according to the following relationship:

First estimation ( ) -Stress and strain integration:

j j j J       1 _ 1 _ ( t ( ( s s s σ σ σ  j j j J ) ( t ) ( ( (vp) _1 (vp) (vp) _1 D E E     j j ) ( ) ( 1 _ s s σ σ   Computing Eq. (18): a G a j j 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j (vp) vol_ j (vp) vol j (vp) vol_ D J E E ) ( t ) ( ) ( 1 1    j (vp) d_ j (vp) d j (vp) d_ D J E E ) ( t ) ( ) ( 1 1

  

Recapitulation of estimated stages according to Eq. ( 35):

( ) {( ) [ ] ( ) [ ( ) ]} j j j J           2 _ 1 _ 1 ) 1 ( ) 1 (( 2 t ( ( s s s s σ σ σ σ     ; j j j J ) ) 1 ( ) 1 (( 2 t ) ( ( ) ( ) ( ) ( 1 vp _2 vp _1 vp (vp) D D E E           ; a G a j j 0 ) ( ) ( ) (  B β σ σ s (vp) s    ; j (vp) vol_ (vp) vol_ j (vp) vol j (vp) vol D D J E E ) ) 1 ( ) 1 (( 2 t ) ( ) 2 1 1 (          j (vp) d_ (vp) d_ j (vp) d j (vp) d D D J E E ) ) 1 ( ) 1 (( 2 t ) ( ) ( 2 1 1         

Second estimation ( ) -Stress and strain integration:

j j j J       2 _ 2 _ ( 2 t ( ( s s s σ σ σ  j j j J ) ( 2 t ) ( ( (vp) _2 (vp) (vp) _2 D E E     j j ) ( ) ( 2 _ s s σ σ   Computing Eq. (18): a G a j j 0 ) ( ) ( ) (  B β σ σ s (vp) s    -Internal variables integration: j (vp) vol_ j (vp) vol j (vp) vol_ D J E E ) ( 2 t ) ( ) ( 2 2    j (vp) d_ j (vp) d j (vp) d_ D J E E ) ( 2 t ) ( ) ( 2 2    (36) ) 1 ( 2     o o S S o S o F   
where S is current sample cross-section, S o is the initial cross section,  o is the nominal stress,  o is the nominal strain and  is the Poisson's ratio chosen to be equal to 0.4. Following a classical analytical deterministic method, a set of parameters have been determined for the two-phase elasto-viscoplastic constitutive model. Calculations are performed on one linear cubic element C3D8. Three numerical simulation scenarios are conducted at 23°C with strain rate equal to 0.083 s -1 (equal to the one used in the experiment). The first scenario is conducted without sub-division technique. The second scenario is conducted with relatively coarse subdivision, and a third scenario uses a finer sub-division. The maximum number of automatic cutbacks allowed for an increment and the increasing factor when two consecutive increments converge are set to 5 and 1.5 respectively.

It should be noted that the goal of this work is not to have a perfect material behavior (material parameters should be optimized to this) but more to compare the different numerical integration methods.

The scenario I: In this scenario, the sub-division technique is not employed (it means that the value of J mentioned in section 3 is equal to one).

The results of simulations are depicted in Fig. 7 and Table 5. The Fig. 7 shows the results of the three explicit numerical integration methods, which are the explicit Runge-Kutta second-order (RK2), the explicit Runge-Kutta fourth-order (RK4), and the Hybrid Forward Euler with explicit Runge-Kutta second-order (Hybrid). At this simulation condition, the RK2 revealed dramatic oscillations around experimental data. However, the simulation converges despite these high oscillations. The result is less accurate than the ones obtained by the two other explicit methods. The RK4 exposes better stability with less oscillations amplitude, subsequently, the accuracy of results are improved in comparison with the results obtained by RK2. These results confirm the superiority of RK4 in term of enhancing results due to a higher number of evaluation stages. The Hybrid method appears to be the most precise without oscillation and has the best stability condition compared to the two other methods. These results obtained for calculation by finite element code are in agreement with results obtained numerically by [START_REF] Ashour | A new very simple explicit method for the integration of midly stiff ordinary differential equations[END_REF] for ordinary differential equations.

Table 5 presents the results of integration methods. This table illustrates the number of increments to finish the simulation, the automatic incrementation cutbacks, the number of iterations necessary to complete the calculation, CPU time consuming and the CPU time rate. The CPU time rate is obtained by considering the method which consumes the least time as the reference.

The Hybrid is the most efficient method, in this case, it consumes 4.92 (sec) CPU time, thus its CPU time is considered to be the reference, it accomplished the simulation in just 177 increments with zero automatic cutbacks. RK2 and RK4 methods are respectively consuming 2.04 and 1.73 times more than the Hybrid method.

Thus in this scenario, it can be noticed that the Hybrid method has a superior property in terms of accuracy of the result, the stability of calculation, and CPU time efficiency. The scenario II: in this scenario the sub-division technique is used, each time increment is divided by 3 (J = 3).

The rest of simulation conditions mentioned in scenario I are remaining the same. According to the results, all the integration methods are in good agreement with the experimental data : quite similar result in term of stressstrain mechanical behavior are obtained from experimental and numerical results (Fig. 8). Once again, the goal is not to be perfectly predictive but more to compare the different integration processes.

Speaking of time-consuming, Table 6 clearly presents the benefits of sub-division technique: compared to the first scenario (I), the efficiency of different integration methods approaches to the efficient one. The Hybrid method is still the most efficient one. It is worth mentioning that despite sub-division, the total CPU time of all the method decrease in comparison with the first scenario (decrease approximately by half for RK2 or RK4). For the Hybrid method, CPU time decreased from 4.92 to 4.78 (sec), this is a tiny difference for one element but for a complex system with thousands of elements, for example, it could become an important issue. This result led to a third scenario with a finer sub-division to observe the outcome of different methods. The scenario III: the finer sub-division has been designed in this scenario, each increment is divided into 10 sub-increments (J=10). The other conditions are the same as for scenario (I). Simulation results indicate that there is no significant difference between the methods in terms of stress-strain accuracy for relatively high subdivision which means smaller stepsize, as depicted in Fig. 9. The RK4 method became more efficient in CPU time-consuming viewpoint, as is illustrated in Table 7. This result is in agreement with the work of other authors 8 and9. These tables exhibit the number of increments required to complete the simulation, the automatic incrementation cutbacks, the number of iterations, the CPU time cost and CPU time rate to indicate the most efficient schema among employed methods. Besides, a deviation is calculated to illustrate the error between the experimental and the simulation results (the highest level of flowed material corresponds here to the major result).

At 40°C, the Hybrid seems to be the less expensive schema, followed by RK2 then RK4 as depicted in Fig. 11, while at 90°C the CPU time cost for the different schemas changes relatively, it can be noticed that the RK2 becomes the most efficient schema. Concerning accuracy, the RK4 demonstrates a very good precision with the lowest deviation in comparison with other methods as shown in Fig. 12. In conclusion, choosing a method in this example faces a dilemma between accuracy of obtained result and CPU time cost efficiency. On one hand, at sub-increment (J=20), all three explicit methods RK2, RK4 and hybrid method are diverging and running out of calculation without accomplishing the simulation due to unsatisfying plastic consistency at that significant amount of applied load. Thus the results of sub-division (J=20) are not depicted. At the other hand, at the sub-division (J=200), all integration explicit schemas are able to complete the simulations as presented in Table 10. It worth mentioning that, with relatively large integration stepsize the explicit methods are condemned by instability.

Table 10, illustrates the three integration methods at sub-increment divisions (J=200). It can be noticed that the Hybrid method is the most expensive method among employed methods; it achieves the calculation with 494 cutbacks due to instability. While the RK4 has a better stability, it achieves the calculation with only 62 cutbacks, as shown in Fig. 13.

The example of this case illustrates that the amount of applied load on the system changes the stiffness behavior of constitutive model equations. In case I, the Hybrid methods has the best CPU time efficiency but it becomes the most time-consuming method in Case II. This result led to explore other integration schemas than explicit schemas for future works in order to provide an insight into the effect of applied load on the mathematical stiffness behavior of constitutive model equations.

Concluding Remarks

In this work, different explicit numerical integration methods have been investigated to update stress-strain in a mean-stress dependent elastic-viscoplastic constitutive model for SCPs. Three explicit numerical integration schemas (fourth and second order Rung-Kutta, and a hybrid method between Runge-Kutta and Euler method)

have been coded in FORTRAN to be used in a UMAT user subroutine in ABAQUS TM . The numerical integration schemas have been investigated in real engineering structure problem with viscoplastic materials such as thermoplastics. The criteria considered to examine each method are the stability of the method toward convergence, accuracy of the result, satisfying plastic consistency for the elastic-viscoplastic constitutive model, and CPU time cost.

The numerical integration schemas have been implemented into the constitutive model in a flexible modular way, which allows applying these methods to any different constitutive model in future. The implementations are achieved by involving sub-division technique. This technique proves an important impact on the accuracy and convergence of numerical integration methods.

For a tensile test without sub-division, the simulation shows that the Hybrid method demonstrates better performance with very good stability, accuracy and an efficient CPU time-consuming in comparison with two other schemas.

For a tensile test with a sub-division technique, relatively at small integration interval, the integration methods are equivalent in terms of stability, accuracy of result, and CPU time efficiency. It worth mentioning that, the CPU time consuming for both RK2 and RK4 are improved largely with the sub-division technique. In spite of more calculation are performed in the sub-division technique, it consumes less CPU time to complete the simulation.

The experiments and simulation are also driving for creep tests on a structure with different temperature and load conditions at different integration time interval. The RK4 schema exhibit to be very competitive in terms of accuracy, stability and CPU time cost particularly at small integration interval. Results of creep tests on a structure demonstrate that load conditions change the performance of numerical integration schemas; it may be due to the change of mathematical stiffness behavior of constitutive model equations. Therefore, these results initiate the idea of investigating implicit integration schemas for future works in order to provide an insight into the effect of applied load on the mathematical stiffness behavior of constitutive model equations.

J

  variable can describe the irreversible volume change of the solid porous skeleton through the existence of a residual saturated porosity after an unstressed state ( The fluid mass variation into elastic part and viscoplastic part is given by Coussy[START_REF] Coussy | Mechanics of Porous Continua[END_REF], as follows,

C

  denotes the effective Hooke's elasticity acting on the solid network. The associated effective stress states as: free enthalpy potential acting on the free amorphous phase, is homogeneous to a partial pressure

.

.

  Where f is the Von-Mises equivalent stress and s p the mean stress associated with the solid phase. s  is an internal friction parameter which governs the pressure sensitivity of the SCP, s c is a cohesion parameter acting on material resistance to pure hydrostatic tension, s b is a parameter that controls the locking of the yield surface in the (  s eq , s p ) plane.

The 1 L

 1  is the inverse Langevin and   denotes the positive part operator. The authors have considered the evolution of an initially isotropic spherulitic mesostructure into an oriented crystalline structure dependent on the loading direction. For this reason, the parameter 0 c p in Eq. (22) has been introduced to capture the effect of a microstructural texture evolution; it is described in a phenomenological way, through two-phase interaction effects, by means of both parameters ij b and the deviatoric part of E (vp) , as following: denotes the initial value of the residual confining stresses ( 0 c p ).

1 

 1 Elastic modulus (Pa) Poisson's ratio Coupling parameters between the solid network and the free amorphous phase M Biot compressibility modulus of the free amorphous phase contained in the solid network of the SCP (Pa) 0 b Biot's ratio in isotropy case: B = 0 b Plastic Biot's ratio in isotropy case: β =  1 Parameters associated with the yield surface s  Friction ratio of the material s b Control parameter of the yield surface locking, between 0 and 1 s c Cohesion of the material, similar to a triaxial tensile limit (Pa) 00 c p Initial residual confining stresses (Pa) Parameters associated with the viscoplastic flow a Intensity of the isotropic hardening 0 d Initial size of the yield surface s  Hardening parameter associated with viscoplastic volume strain s  Control parameter of the intensity of volume viscoplastic strain  Viscosity parameter equivalent to a long relaxation time (s) s k Viscosity parameter equivalent to a short relaxation time (Pa) Parameters specific to large strains (not activated in this work) ,n Parameters to describe the evolution of the interactions between the two phases at large strains W m Cavitation energy of the material (N.m.kg -1 )

  strain rate. The algorithm allows evaluating intermediate estimation for variables which are depending on stress update at each Gauss point, for the given incremental strains. For pure elastic, the effective stress is directly obtained from the elastic stress rate ) (e s σ  .

  stain rates for first and second estimation,(vp) vol_ E 1 and (vp) vol_ E 2 to the volumetric plastic strain for first and second estimation, the deviatoric viscoplastic strain rate for first and second estimation, (vp) d_ E 1 and (vp) d_ E 2 to the deviatoric viscoplastic strain for first and second estimation.

Fig. 2 .

 2 Fig.2.Flowchart showing only stress-strain integration processes, according to Runge-Kutta second-order

Fig. 6 .

 6 Fig.6. Experimental stress-strain evolution at 23°C and 0.5 mm/min vs. time for the employed semi-crystalline polymer.

Fig. 8 .

 8 Fig.8. Experimental vs. simulation data for different numerical integration methods with (J=3) under tensile test.

Fig. 10 .

 10 Fig.10. Axisymmetric model for creep simulation. (a) Dimension and boundary condition of the model. (b) Typical evolution of a polyethylene-kind material under creep according to the elasto-viscoplastic model presented in Section 2. The red dot corresponds to the maximal creep level recorded during experiments. The case I: The creep tests has been conducted under 21 MPa at the temperatures of 40 and 90°C during 1.5x10 +6 seconds. Each integration time increment is divided into 10 sub-increments. The simulations results are outlined in Table8 and 9. These tables exhibit the number of increments required to complete the simulation, the

Fig. 11 .

 11 Fig.11. Total CPU time cost vs. different integration schemas for creep test. (a) At temperature 40 °C. (b) At temperature 90 °C.

Fig. 12 .

 12 Fig.12. The deviation between experimental and simulation result, for different integration schemas, under creep test. (a) At temperature 40 °C. (b) At temperature 90 °C.

Fig. 14 plotted

 14 Fig.14plotted creep time vs. the increments, it illustrate the capacity of numerical integration schemas to reach the convergence. The RK4 converges rapidly: it needs 6515 increments to complete the calculation vs. 10268 for RK2 and 18708 for the Hybrid method. In conclusion, the RK4 has a better stability property than two other methods.

Fig. 13 .

 13 Fig.13. Total CPU time cost for the different integration schemas, at J=200 under creep test.

Fig. 14 .

 14 Fig.14. Numerical integration schemas rapidity toward convergence at sub-increment J = 200 under creep test.

Table 2 .

 2 Generalized Butcher table for s-stage order explicit Runge-Kutta method.

E

  of the intermediate plastic strain. The number of intermediate effective stresses and its corresponding intermediate strains depend on the number of estimation stages of the employed method. The tensor s σ is the effective stress of the solid matrix, t

			1 σ is an estimation _ s
	of the intermediate effective stress tensor,	(vp _1	)
	 is time increment, is the corresponding plastic strain 1 _ s σ  is the ), the computing of effective stress tensor is the viscoplastic strain, and ) vp ( _1 D rate. After initializing the intermediate effective stress ( intermediate effective stress rate, ) (vp E 1 _ s s σ σ   through (Eq.18) is the next step. The internal variables have been thus updated : ) (vp vol E corresponds to the volumic viscoplastic strain, ) (vp vol E 1 _ is the intermediate volumic viscoplastic strain, and ) (vp vol D is the volumic
	viscoplastic strain rate, while

Table 4 .

 4 Butcher tableau for fourth order Runge-Kutta method.The Runge-Kutta fourth order has four intermediate estimations in order to achieve final estimation at (t i+1 ). The intermediate estimations carried out as mentioned in the Fig.3. The output of each intermediate stage estimation becomes an input for the next intermediate stage evaluation. Besides, at each step, the stepsize has to be adapted

	according to the Table 4. For instance,	and	are evaluated with a time stepsize equal to t  /2J, while	is
	evaluated with a t  /J stepsize.

Table 5 .

 5 The efficiency of numerical integration methods with sub-increment (J=1) for tensile simulation.

	Method	RK2	RK4 Hybrid
	Increments	273	255	177
	Cutbacks	45	39	0
	Iteration	759	514	184
	CPU time (sec)	10.05 8.50	4.92
	CPU time rate	2.04	1.73	1.00

Fig.7. Experimental vs. simulation data for different numerical integration methods at (J=1) under tensile test.

Table 6 .

 6 The efficiency of numerical integration method with sub-increment (J=3) for tensile simulation.

	Method	RK2	RK4	Hybrid
	Increments	177	177	177
	Cutbacks	0	0	0
	Iteration	184	184	184
	CPU time	5.11	4.93	4.78
	CPU time rate	1.069 1.031	1.000

Table 8 .

 8 Performance of numerical integration methods with sub-increment J=10, and temperature 40 °C for creep test.

	Method	RK2 RK4 Hybrid
	Increments	4255 4262	4256
	Cutbacks	0	0	0
	Iteration	4731 4776	4724
	CPU time (sec)	4797 5962	4705
	CPU time rate	1.02 1.267	1.00
	Deviation %	27.63 4.93	27.63

Table 9 .

 9 Performance of numerical integration methods with sub-increment J=10, and temperature 90 °C for creep test.

	Method	RK2 RK4 Hybrid
	Increments	3459 3445	3446
	Cutbacks	2	0	0
	Iteration	3991 3934	3950
	CPU time (sec)	3591 4711	3933
	CPU time rate	1.000 1.312	1.095
	Deviation %	14.55 6.86	14.62

Table 10 .

 10 The efficiency of numerical integration method with sub-increment J=200, and temperature 23 °C under creep test.

	Method	RK2	RK4	Hybrid
	Increments	10268	6515	18708
	Cutbacks	288	62	494
	Iteration	18816 10461	27921
	CPU time (sec)	43230 48891	75636
	CPU time rate	1.00	1.13	1.75

Structural analysis under creep test

Creep experiments have been performed at different loads and temperatures conditions. The highest level of flowed material under the applied creep load has been measured at the end of the experiments. An axisymmetric finite element model has been designed according to the dimension of the corresponding experimental sample as outlined in Fig. 10. The Fig. 10b illustrates the typical flow of the polymer through a gap due to the pressure applied to the polymer. The highest level of flowed material (highlighted in red color on Fig. 10-b) under creep load are analyzed for the different integration methods and are compared to the experimental database.

The stepsize is automatically chosen by ABAQUS TM between 1x10 -5 to 1000s, the maximum number of automatic cutbacks allowed for an increment, and the increasing factor when two consecutive increments converge are set equal to 12 and 1.2 respectively. In this section, the potentials of all numerical explicit integration methods mentioned in Section 2 are investigated for structure calculation under creep test (two conditions of loading are analyzed: case I and case II, depicted in the following). The accuracy of simulation results, the stability of integration schemas, the CPU time efficiency, plastic consistency are still the scope debated of this section.


domain: spi phys sdu

text:
Sashikanth Rapeti

Shakti N Mishra

A Proposal for the Fundamental Structure of South American (GPS Aided) Geo Augmented Navigation System -Space Weather Integrated Monitoring and Surveillance Systems( SAGANS-SWIMS)

Keywords: SAGANS -SWIMS -South American GPS Augmented Navigation System -Space Weather Integrated Monitoring and Surveillance Systems, Space Sandwich, Deep Space Segmentation, GEES segmentation

    A Proposal for the Fundamental Structure of South American (GPS Aided) Geo Augmented Navigation System -Space Weather Integrated Monitoring and Surveillance Systems( SAGANS-SWIMS)

INTRODUCTION

Most of the countries in the Latin American Continent possess their own unique and distinct space satellites on Communication, Earth Observation, Science, Technology and Education. It is worth-noting that all the countries that had their space satellites launched, possess their own unique strategic goals in the space arena in the continent, thereby collaborating with space agencies in either Asia ( China, India etc), ESA [START_REF] Esa | The Global Exploration Strategy: The Framework for Coordination, European Space Agency[END_REF], CNES (FRANCE) , NASA etc. The following table provides the number of satellites for each country which mostly are launched, some of them are yet to be launched at the earliest. Paraguay still needs to launch a satellite of its own. The main areas explored are the Communications, Earth observation, Science, Technology and Education. Every country in the continent of Latin America is equipped with their own unique state of-the-art facilities regarding their space missions. Brazil tops the list in all the three areas [START_REF] Do | O Estado Brasileiro e a Questão das Tecnologias Espaciais: do GOCNAE à AEB[END_REF] having its own space program [START_REF]National Program of Space Activities: PNAE 2012-2021[END_REF], AEB ( BRAZILIAN SPACE AGENCY) [START_REF] De Moraes | Instituições e Agências Brasileiras", Chapter of book "A Conquista do Espaço: do Sputnik à Missão Centenário[END_REF] and INPE [START_REF] Câmara | A pesquisa espacial no Brasil: 50 anos de INPE (1961-2011)[END_REF] (NATIONAL INSTITUTE FOR SPACE RESEARCH) are two premier space research organizations that drive Brazil ahead in the South American continent in the Space sector [START_REF] Santana | O Projeto CBERS de Satélites de Observação da Terra[END_REF], followed by Argentina which are the largest economies in the continent. Brazil ( collaborating with China -CBERS , in terms of Earth Resource satellites) and Argentina have the maximum number of satellites -combined as well as in each specific area individually. Peru launched a satellite which is dedicated to different applications [START_REF] Andina Peru | Nasa launches Peruvian satellite into space[END_REF]. Peru had been very active in developing programs like INCA at UCH, and with nano-satellite and pico-satellite technology in different universities [START_REF] Roman-Gonzalez | INCA PROGRAM FOR DEVELOPING A NANOSATELLITE AT THE UCH. 67th International Astronautical Congress -IAC 2016[END_REF]. Aerospace technology is achieving new frontiers in the country at a good pace [START_REF] Roman-Gonzalez | AEROSPACE TECHNOLOGY IN PERU.66th International Astronautical Congress -IAC 2015[END_REF] [10] [11] [12] [13] [START_REF] Romero | Peruvian University Consortium in the QB50 Project[END_REF]. In this connection, proposals were made regarding small satellite network for the Peruvian region [START_REF] Martín | Design of a Peruvian Small Satellite Network[END_REF] and the nanosatellite 'Chasqui 1' was launched with the help of Russia [START_REF] Info | Peru and Russia allies to launch satellite into space Chasqui I[END_REF]. Colombia entered the space race in the continent with the formation of Colombian Space Commission [17] [START_REF]Decreto Presidencial[END_REF]. The small landlocked country 'Bolivia' had only one satellite to its credit, as aerospace technology is being promoted in the country [START_REF] Vargas-Cuentas | Promotion of Aerospace Technology in Bolivia[END_REF] [START_REF] Vargas-Cuentas | Promotion of Aerospace Technology in Bolivia[END_REF]. Bolivia even had formed its own space agency and launched a satellite mainly for communication purposes. Ecuador, being a small country to the north of Peru is very active in space and has a Civilian Space Agency whose director is Cdr. Ronnie Nader who is the only astronaut of that country [START_REF] Nader | The Ecuadorian Experience in Space: The NEE Satellite Constellation[END_REF]. The Bolivarian Republic of Venezuela also entered the South American Space race [START_REF] Hernández | Current space projects of the Bolivarian Republic of Venezuela[END_REF]. Of course, there are many challenges and opportunities in the space and aerospace arena in the Latin American Continent [START_REF] Gocłowska-Bolek | Latin America Space Research, Challenges and Opportunities[END_REF]. Since, every country in the continent has its own space program, Argentina with the Support of Brazil had put forward a proposal for a common space agency for entire continent [START_REF]Argentina, with Brazilian support proposes a South American Space agency[END_REF]. Many sites provide good details of the various space missions of the countries of the world (of which some examples are provided

)https://eoportal.org/web/eoportal/home , https://space.skyrocket.de/directories/sat_c.htm http://www.earth-observation-directory.com/missions/ and the basic details are outlined in these sites, along with the specific purposes accomplished by those missions.

II. METHODOLOGY

The SAGANS-SWIMS is but a unique project in this regard in terms of integrating the various fields like communication, navigation, earth observation etc., for which every country relies heavily upon technologically, especially in the continent of South America. The continent needs autonomy in the space and aerospace sectors [START_REF] Sánchez | Space Technology Comes to Latin America: Part of the Hemisphere's Road to Autonomy[END_REF]. The detailed review of various space activities in the continent are described in a good paper by Bruno Victorino et al., [START_REF] Victorino | Review of Space Activities in South America[END_REF]. At this crucial juncture, SAGANS -SWIMS will play a much significant role by integrating the Communication, Earth Observation, Science, Technology and Education fields, in the areas such as airborne mapping, ocean floor mapping, agriculture, application of Civil engineering like structural deformations, Vehicle and waypoint navigation, Transit systems in transportation, earthquake and seismic studies, disaster warning systems, natural resources monitoring, marine seismic surveying, precision farming, open pit mining, Oil and natural gas monitoring in the coastal regions of the continent in a more sophisticated manner. Indeed SAGANS-SWIMS can be a new beginning for a South American Space Era [START_REF] Sarli | SOUTH AMERICAN SPACE ERA[END_REF].

The fundamental technical structure of SAGANS-SWIMS is as follows:

A. SPACE SANDWICH ( SPACE SEGMENTATION) :

This segmentation has two parts, completely handled by SWIMS.

1) Deep Space segmentation or Gravitational Equilibrium of Earth Sun (GEES) segmentation:

This deep space segment will be purely handled by SWIMS ( Space weather Integrated Monitoring and Surveillance Systems) , which is a separate module of SAGANS-SWIMS. This segment will collect the data from two space satellite missions -SDEM ( Sun Dynamics Explorer Mission) and UHPE ( Ultra High Energy Particle). These two space satellite missions will station themselves at 1.5 million km from the Earth in a symmetrical way, to the east and to the west of the Earth and will be strictly parallel to the equatorial plane of the Earth, thereby the two spacecrafts will get a good picture of the solar wind interaction with the magnetic field of the Earth at the polar regions. Both the spacecrafts will be at the GEES ( Gravitational Equilibrium of the Earth Sun) east-west points, parallel to the equatorial plane of the Earth as shown in the figure below.

2) Magnetically Turbulent Space (Zone of) Earth Segmentation (MTSE) segmentation: This Space Zone of the

Earth is magnetically active and turbulent all the time and is technically known as Magnetosphere of the Earth which is spread typically up-to an altitude of 60,000 km from the Earth's crust. This Space-zone is active even during quiet geomagnetic disturbance periods. Therefore, a spacecraft known as MEE ( Magnetosphere Earth Explorer) will perform a polar orbit revolution around the Earth perpendicular to the equator passing through the Latin American Continent. This Spacecraft ( satellite) will monitor the magnetosphere region orbiting at an altitude of 56000 km from the Earth' s surface by getting the data from UHPE and SDEM regarding the energetic charged particle influx and its interaction with the Earth's magnetic field lines, the compressing of the field lines by the Copyright @ Dr. Sashikanth Rapeti 2 | P a g e

B. ATMOSPHERE SANDWICH ( ATMOSPHERE SEGMENTATION ):

This segmentation also has two parts:

1) Upper Atmosphere segmentation -UAS: This zone consists of most of the magnetosphere -ionosphere region from the equatorial and polar perspective to monitor the space weather conditions at the Earth's equator and the polar regions in the magnetosphere-ionosphere zone. In addition the polar and equatorial magnetosphere-ionosphere currents are measured which are much important to understand the various anomalies and instabilities like Rayleigh-Taylor, Kelvin Helmholtz etc in the ionosphere and the magnetosphere connection. The Lorentz force along with the electromagnetic energy density is also measured at the polar regions. In this connection, two spacecrafts FALCON EXPLORER ( Polar) and EAGLE EXPLORER ( equatorial) will orbit in the magnetosphere-ionosphere region in the polar and equatorial zones at an altitude of 42000 km from the earth's surface. Equatorial electrojets are given special importance in order to understand the equatorial magnetosphere-ionosphere dynamics.

2) Lower Atmosphere segmentation -LAS: A series of four major space satellite missions need to be launched which will orbit in the lower atmosphere, from thermosphere, ionosphere to the stratosphere-ionosphere boundary region at varying altitudes. The outline is as follows:

a) GEO-SAT Constellation -The GEO-SAT constellation is a series of Earth Observation satellite missions four in number, which will mainly serve the purpose of communication and navigation purposes and airline communication purposes, by integrating the data from the satellites in the GEES segmentation, MTSE segmentation and UAS. These satellites will orbit at a height of 24000 km above the surface of the earth.

b) TERNO-SAT constellation (Terrestrial Network of OBSERVATION Satellites constellation) -The TERNOSAT series of four satellites will orbit the Earth at an altitude of 16000 km above the Earth's surface, will serve the purpose of land mapping, precision farming, forest and natural resources monitoring, agricultural land monitoring, Structural deformations in dams, huge constructions like bridges, various other civil engineering applications etc. c) NEV-SAT Constellation (NETWORK OF ENVIRONMENTAL SATELLITES) -This series of four satellites will orbit at an altitude of 8000 km above the surface of the earth, serving the purposes of Advanced disaster warning systems like earthquake studies, volcanic eruptions, Tsunamis, Open pit mining etc. d) HYDRO-SAT Constellation -This constellation of four satellites will orbit the Earth at an altitude of 4000 km and will monitor the coastal regions of the Latin American continent, pacific and atlantic ocean floor mapping along the coastal length and breadth of the continent. The river water resources, Oil and natural gas reserves etc are also closely monitored for a proper utilization of natural resources.

e) BREEZE Constellation-This constellation of three satellites will monitor the ozone dynamics above the continent of South America

C. TERRESTRIAL SANDWICH ( TERRESTRIAL GROUND SEGMENTATION -TGS) :

. This segmentation also has two parts:

1) Terrestrial Equatorial Ground segmentation : The ground stations will be set up near the equatorial regions in the continent to receive the data from the various constellation of satellites orbiting at various altitudes.

2) Terrestrial Polar Ground segmentation :

The ground stations will be located at the southern most locations of Argentina and Chile to receive the data from the polar satellites of SAGANS.

D. USER SANDWICH ( USER ORIENTED SEGMENTATION -UOS) :

SAGANS-SWIMS is expected to cater to the diversified requirements of users from various fields -academics, scientific and technological , engineering and medical, Disaster studies, Vehicle and way-point navigation, airline communication, governments, defense systems ( military, naval and air-force).

III. THE DIFFERENT PHASES OF SAGANS-SWIMS

A. Initial Phase ( I-Phase):

This Phase consists of the theoretical design of the various space satellite missions, their goals studied, the various parameters which need to be measured are defined, the mathematical analysis will be done, the various technologies to be employed in the sensor and instrument development which will be mounted on the satellite.

B. Testing, Experimentation and Validation Phase ( TEV -

Phase ) : This phase consists of the testing, experimentation and validation of the various instruments to be mounted on the spacecraft in orbit, the software needed to operate the instruments on board.

C. Launch Phase ( L-Phase) :

This phase includes the funding and various options required for the final launch of the spacecrafts and placing them in orbit which includes the strategies being employed.[30]

IV. SATELLITE OR SPACECRAFT PROPULSION SYSTEMS -FUSE -FUEL USAGE STRATEGY AND EXPERIMENTATION

The below are the various options that can be considered for spacecraft propulsion systems and only the outline is given in this paper. Many of the technologies mentioned below are either already existing or some might be in the experimental stages.

A. Electromagnetic and Electrodynamic spacecraft propulsion system (EESP) :

The solar energetic charged particle flux that strikes the spacecraft detectors can be utilized as a source of propulsion systems. In one sense, this is approximately equal to the "Magnetic Sail" propulsion system, but much more can be possible in this type of propulsion. There are many options like Electrodeless Plasma thruster, MPD thruster, Pulsed plasma thruster, Helicon Double layer Thruster, Pulsed inductive thruster which fall under this category.

B. Water (based) propulsion system ( WPS) :

This also can prove to be an effective propulsion system for the satellites in orbit.

C. Hydrazine (based) propulsion system ( HPS) :

This is a well-known spacecraft propulsion system which is currently in use by Ariane Group and others.

D. Green Technologies ( based) propulsion system :

This is also an innovative technology by Launch companies like Ariane Groups.

E. Superfluid helium ( based) propulsion system :

Considering the special properties of Liquid helium-3 and superfluid Helium-4, there can be some possibilities for this kind of propulsion system, and Ariane group is actively working in this area too.

V. SPACECRAFT ORBITAL DYNAMICS This section just outlines the orbital dynamics of spacecrafts, but the details are out of scope of this paper. The orbits will be angular, Quasi-equatorial, Quasi-polar, orbiting only the length and breadth of the South American Continent and need not necessarily be Geo-stationary. (The satellite need not make a complete revolution of the Earth either in equatorial or polar orbit, but will orbit only along the horizontal and vertical size of the continent, without being Geo-stationary at all). Although this can be a radical approach, but not at all impossible, if certain important parameters are taken into account.

VI. RESULTS AND DISCUSSIONS 1. SAGANS-SWIMS is a completely new and

unique project in terms of space weather, communication, navigation, Earth observation etc, which is still in its initial stages. 2. SAGANS-SWIMS can fulfill the need of common space goals and common space agency for entire South American continent as a whole. 3. The minimum and maximum operational timeline for the design, implementation and operation of SAGANS-SWIMS can approximately range from 10 -25 years and is certain to become the most robust communication, navigation and Earth observation system in the continent by integrating the various technologies in the continent. 4. Every country in the continent of Latin America has its own space goals and SAGANS-SWIMS can fulfill these by integrating these space goals. 5. SAGANS-SWIMS will play the major role in shaping the present and future education systems in the continent in the form of SIT -classrooms ( SATELLITE INTEGRATED TECHNOLOGY) from school level to University level. 6. SAGANS-SWIMS will ensure that the safety, security and integrity will never be compromised for any country in the continent under any circumstances. 7. SAGANS-SWIMS will prove to be an effective technological and friendly tool for users in various fields as well as governments and defense systems. 8. SAGANS-SWIMS will facilitate ease of trade and business for various countries inside and outside of the continent. 9. SAGANS-SWIMS will make sure the natural resources of every country are monitored, utilized properly and not fall in the hands of miscreants, by keeping a 'Technological Eagle-eye' on the continent. 10. SAGANS-SWIMS will play the role of an advanced disaster warning system for the entire continent, thereby lots of lives can be saved from the wrath of nature, in the form of earthquakes, Tsunamis etc. 11. SAGANS-SWIMS will be the major boost for agriculture, generating revenue for various countries in the continent in terms of trade with other continents. 12. Special focus will also be given to the coastal area monitoring of the continent in terms of marine seismic surveying, Oil and natural gas monitoring etc. 13. The data received will be real-time and be of a higher accuracy than the traditional GPS systems, which is one of the major goals of SAGANS-SWIMS. 14. The current timeline of this project is to determined yet as the project itself is undergoing a lot of refinements and scrutinizations in its initial stages.

15. The formation of Scientific research teams are currently under progress for this unique big project as of now. 16. Some technical specifications are outlined regarding this unique project. More inputs are necessary from many other scientific researchers in the Latin American continent. 17. The SAGANS-SWIMS is in favor of pico, nano, micro satellites which are cost-effective and technologically robust. 18. Although Cube-Sat technology is considered one of the best for nano satellites, still we can expect some modifications in order to improve efficiency of the satellites. 19. Since SAGANS-SWIMS is a huge project consisting of a total of 24 satellites , and the various sensors, instruments to be mounted, the project will begin with UHPE and SDEM, the spacecrafts belonging to the Deep space segmentation or the GEES Segmentation , thereby going systematically. 20. The technical structure of every space satellite mission will be dealt first and then come the sensors and the instruments to be mounted, then the orbital dynamics and Fuel usage strategy following the Initial, TEV and the Launch phase in a systematic fashion. 21. Some genuine ideas need to be materialized when it comes to important issues like partial funding, step-by step funding etc in this regard.

CONCLUSION AND FUTURE TRENDS

Brazil, Peru, Colombia, Ecuador, Bolivia, Argentina, Uruguay, etc. and some central American countries like Guatemala are some of the foremost countries to accept responsibilities especially regarding space and aerospace technologies as these nations highly recognize and realize the core importance of these advanced technologies for a prosperous society. Since SAGANS -SWIMS offers more features regarding satellites regarding vehicle and way-point navigation, airborne and seafloor mapping, precision farming, monitoring of natural resources, illegal mining, illegal activities, terrorism, smuggling, drug trafficking, forest mapping, utilities and retail industry, monitoring and advanced warning systems of natural calamities like earthquakes, Tsunamis, volcanic eruptions, cadastral surveying, structural deformations and other civil engineering applications and much more etc., the various universities in these countries can actively participate in the preparation of nano-satellites and can begin with any theme listed above like vehicle navigation, weather monitoring, precision farming , seafloor mapping etc.

ACKNOWLEDGMENT

The author is highly indebted to all the dedicated scientists and researchers in Latin America who are working towards these high goals. The author is highly grateful to his co-author Shakti Narayan Mishra for his intense dedication and commitment in the gigantic and long-term SAGANS-project.

The author is intensely grateful to Dr. Avid Roman Gonzalez and his brilliant team of Researchers at INTI Labs, UCH , Peru for their brilliant ideas, dedication, commitment towards the gigantic SAGANS-SWIMS long-term project. Dr. Avid Roman Gonzalez is striving extremely hard to bring all the dedicated scientists, engineers and researchers in Latin America on a common platform for a common space agency and space goals for Latin America as a whole which are clearly evident from the references. Hence, the time is ripe and there is an urgent necessity for all the dedicated scientists, engineers and researchers in Latin America to bring the SA(G)GANS-SWIMS into reality.

Fig. 1 .

 1 Fig. 1. SDEM and UHPE at the equatorial plane on the east and west GEES points.( Figure only for representation purposes only) *courtesy: Google images

TABLE I

 I

	.	NUMBER OF SATELLITES LAUNCHED BY EACH COUNTRY IN LATIN AMERICA
		A. NUMBER OF SATELLITES
	B. COU NTRY	1) COMMUNI CATIONS	2) EARTH OBSERVATIO N	3) SCIENC E, TECHNOLOG
				Y AND
				EDUCATION
	Brazil	21	14 ( series)	13
	Argentina	12	8 ( series)	14
	Peru	-	1	4
	Chile	-	1	3
	Ecuador	-	-	3
	Colombia	-	-	3
	Bolivia	1	-	-
	Venezuela	1	2	-
	Uruguay	-	-	1

* (

Courtesy: https://space.skyrocket.de/directories/sat_c.htm)

(Abbreviation) Journal Name Vol. XXX, No. XXX, 2013