NUMERICAL SIMULATION OF VISCOELASTIC MATERIALS

Purpose. The main goal of this paper is to develop the numerical model of viscoelastic material using finite element method (FEM). The model was applied to asphalt-aggregate mixtures. Additionally the obtained numerical results with use of different FEM software were compared with experimental data. Methodology. In order to perform the investigation, the numerical specimen was built within FEM software. Material of the specimen was assumed to be viscoelastic. Viscoelastic materials are characterized by a combination of elastic behavior, which stores energy during deformation, and viscous behavior, which dissipates energy during deformation. It was assumed that the behavior of the material corresponds to generalized Maxwell model. The model consists of a spring element in parallel with a number of spring and dashpot Maxwell elements. Generalized Maxwell model consisting of 5 ele-


Introduction
The constitutive behaviour of many engineering materials combines with its elastic, viscous, plastic and fracturing response.Determining the parameters of such materials within a wide range of loadings is extremely complicated.Even in such a case when plasticity and fracture are not considered, material viscosity causes temperature and strain rate dependence of the stiffness, what is of great importance analysing asphaltic materials [1,2] or polymers.For that reason, temperature and strain rate are sometimes integrated in the elastic parameters of the models.
Another problem arising from application of non-elastic materials in the process of structural design is to perform appropriate experimental test in order to calibrate material parameters.While such tests are performed, we need to use special optimization algorithms which allow obtaining the best curve fitting results.This problem is especially important in case of cyclic tests used for viscoelastic materials [5, 7, 9-11, 14, 15].

Purpose
The main goal of this paper is to develop the numerical model of viscoelastic material using finite element method.The model will be applied to asphalt-aggregate mixtures.Additionally the obtained numerical results with use of different FEM software will be compared with experimental data.

Methodology
In order to perform the investigation, the numerical specimen was built within FEM software.Material of the specimen was assumed to be viscoelastic.Viscoelastic materials are characterized by a combination of elastic behavior, which stores energy during deformation, and viscous behavior, which dissipates energy during deformation.The elastic behavior is rate-independent and represents the recoverable deformation due to mechanical loading.The viscous behavior is rate-dependent and represents dissipative mechanisms within the material [3,4,8].A wide range of materials such as polymers, glassy materials, soils, biologic tissue, and textiles exhibit viscoelastic behavior.Material is viscoelastic if its stress response consists of an elastic part and viscous part.Upon application of a load, the elastic response is instantaneous while the viscous part changes over time.Generally, the stress function of viscoelastic materials can be defined in an integral form.Within the context of small strain theory, the constitutive equation for isotropic viscoelastic material can be written as follows: , where e -denotes deviatoric part of the strain, ∆represents volumetric part of the strain, G(t)represents shear relaxation kernel function, K(t)bulk relaxation kernel function, t -current time, τpast time and I -is unit tensor.The kernel functions are represented in terms of Prony series, which assumes that the kernel functions can be equivalently expressed as The integral function can recover the elastic behaviour at the limits of very slow and very fast loadings.Here, 0 G and 0 K denote shear and bulk moduli at the fast load limit (i.e. the instantaneous moduli), while G ∞ and K ∞ are the moduli at the slow limit.The elasticity parameters input correspond to those of the fast load limit.Moreover, by admitting, the deviatoric and volumetric parts of the stress are assumed to follow different relaxation behaviour.The number of Prony terms for shear G n and for volumetric behavior K n need not be the same, nor do the relaxation times G i τ and K i τ .
The Prony representation has a prevailing physical meaning in that it corresponds to the solution of the classical differential model (the parallel Maxwell model) of viscoelasticity.This physical rooting is the key to understand the extension of the above constitutive equations to largedeformation cases as well as the appearance of the time-scaling law (for example, pseudo time) at the presence of time-dependent viscous parameters.
At the same time, generalized Maxwell material model was offered.The following Fig. 1 shows one dimensional representation of a generalized Maxwell solid.It consists of a spring element in parallel with a number of spring and dashpot Maxwell elements.The development of material model was accomplished with software based on FEM.In order to solve the problem the numerical specimen of asphalt pavement was built.The specimen and its coordinate system are shown in Fig. 2, a, while the loading scheme is shown in Fig. 2, b.
Mechanical properties of the material correspond to the properties of real asphaltic material.It In order to develop the numerical model of the material the data based on [13,14] was used.Corresponding parameters of Maxwell model are shown in Table 1.In Table 1, volumetric (bulk elastic modulus) and deviatoric (shear elastic modulus) parts of the Maxwell elements are given.Moreover, corresponding viscosities are shown for each modulus.In order to implement appropriate loading and boundary conditions and to eliminate rigid motion, the 1/8 part of the specimen was investigated (see Fig. 3, a).Corresponding boundary and symmetry conditions were applied to reduced model (see Fig. 3, b).

Optimized data of Maxwell model
The specimen was loaded kinematically applying a displacement in z-direction on the upper part of the cube (Fig. 3, b).The magnitude of the displacement was determined in such a way the strain in this direction to be equal 5 • 10 -5 .Thus, the small deformations theory can be assumed.For our specimen, the displacement of the upper part of the cube was determined to be equal 5 • 10 -5 m.Moreover, the loading was taken with different frequencies corresponding to formula: sin( ) u A t = ω .The following frequencies were proposed: f = 0.1, 0.5, 1.0, 5.0, 10.0, 25.0 Hz (see Table 2).Here, A denotes amplitude of the load (A = 5 • 10 -5 m), ω is circular frequency, 2 f ω = π .

a b
Fig. 3. Investigated scheme Table 2 shows frequencies and corresponding values of dynamic moduli (E dyn ) and phase angles (φ) obtained in experiments.In the paper, the comparison of results is shown obtained with ANSYS and LS-DYNA software, based on finite element method.

Findings
As the result of calculations the stress and strain state of the asphalt specimen were obtained.It is necessary to note that the aim of this part was to develop numerical model of asphaltic material and to compare mechanical characteristics of the model with experimental data.To fulfill the task we had to calculate the dynamic elastic modulus and phase angle.
We can note that stress and strain at such type of loading can be written as: ( ) ( ) where σ 0 , ε 0 (ω) -amplitude values of stress and strain.
To determine dynamic elastic modulus and phase angle we have to write some relationships.So, to define the modulus we used such formulas: Thus, to calculate the dynamic modulus, we have to know the amplitude values of stress and strain in the model.
In order to visualize the relations between stress and strain during numerical experiment the history graphs were built.The diagram of stress and strain changes for load case with 0.1 Hz frequency is shown in Fig. 4.

Fig. 4. Diagrams of stress and strain
To obtain the needed values we can use the diagrams of stress and strain as functions of time (Fig. 4).The amplitude values of stresses and strains are shown in the table 3. Also, obtained dynamic modulus are written in the table 3.
It is necessary to note that the characteristics achieve their maximal values at different time.So, maximal stress in the specimen appears early than maximal strain.To calculate phase angle we will use this difference.
To calculate phase angle, we also can use the diagrams (Fig. 4).On the diagrams we can see that there is the difference in time between maximal values of stresses and strains.It is equaled to the phase angle.To obtain the angle value it is necessary to determine the difference between maximal values of stresses and strains on the time-axis and take into account the scale of the time.Also, we need to use the formula: The difference as time between values of stresses and strains in the model and obtained values of the angle are shown in the table 4.   So, we can write main results (elastic dynamic modulus and phase angle) in the table 5.The obtained results of investigation can be given as amplitude-phase characteristics.The storage modulus E′ is situated on the x-axis, and loss one E′′ is situated on the y-axis, according to relations:

Numerical results (ANSYS)
The modulus were calculated based on the formulas.Obtained results are shown in the table 6.Also, in the table 6 experimental results are written.
Using data showed in the table 6, we can obtain the diagram of the modulus to compare experimental and numerical results (Fig. 5).On the Fig. 5 we can see good agreement between experimental and numerical (ANSYS) data.
To estimate the different software, we will give the data obtained in LS-DYNA [14].The values are showed in the table 7.
The comparison of experimental and numerical (LS-DYNA) data is shown in Fig. 6.Comparing the diagrams (Figs. 5, 6) for numerical results obtained with ANSYS and LS-DYNA software we can see that the results of ANSYS calculations are better satisfied to experimental data.Let us to estimate the values of the corresponding modulus and phase angles obtained with different software.To do it, the values were written to common table 8.
As we can see from the table 8, the differences between the values are not essential.To confirm it the ratios of results obtained in ANSYS and LS-DYNA were determined and written to the table In the table 9, the ratios of elastic dynamic modulus and phase angles are shown.The relations of ANSYS results to LS-DYNA ones were calculated.It is necessary to note that every value is clothed to 1.It points on similarity of the results.
We can note that differences in dynamic modulus are small, and in phase angle are bigger.So, for the elastic dynamic modulus ratios the averaged difference in the results does not exceed 1%.As an exception, maximal difference is about 2% for load case with 10 Hz frequency.Minimal difference is 0.15% for load case with 0.1 Hz frequency.
If we speak about differences between phase angle results we can note such ratios.The maximal difference is 10.9% for load case with 25 Hz frequency.And minimal difference is 1.6% for load case with 0.1 Hz frequency.
So, the results of calculations for different software at load case with 0.1 Hz frequency coincided most exactly.

Originality and practical value
In the article, viscoelastic material model based on generalized Maxwell scheme was developed using ANSYS software.In order to compare the numerical model of the material with real asphalt-aggregate mix data, the values of dynamic modulus and phase angle were selected.The procedure of the calculation was explained in [14].Moreover, in this work the experimental data are given.To visualize the results the diagram of storage modulus vs. loss modulus was built (so called Cole-Cole graph) (Fig. 5).Analyzing the diagram one can see good correspondence of numerical material model to experimental data.
Moreover, the numerical data obtained with ANSYS and LS-DYNA software were compared.The compared data are written in Table 8.In Figs. 5, 6 the diagrams of ANSYS and LS-DYNA results are shown.
It should be emphasized that FEM gives the possibility of determining stresses and strains for asphalt pavement non-elastic models what is of great importance using mechanistic design procedures.

Conclusions
Comparing numerical results obtained with ANSYS and LS-DYNA software we can see that the results of ANSYS calculations are better satisfied to experimental data.
Analysis of data are written to the table 8 showed.Increasing of load frequency causes in increasing of elastic dynamic modulus and decreasing of phase angle.It can be noted that dynamic modulus increase from 1 to 6 load cases.But, values of phase angles decrease from 1 to 5 load cases and increase at 6 load case.Such behavior of the values is fulfilled for every investigated software.
Also, we can see that for load cases 1-5 dynamic modulus obtained with ANSYS software is bigger than the modulus obtained with LS-DYNA one.For load case 6, the LS-DYNA modulus is bigger than ANSYS one.
For values of phase angle, we note that phase angle is bigger for ANSYS software at load cases 1-3.And it is smaller for load cases 4-6.
It can be noted that relationships between load frequency and obtained values (elastic dynamic modulus and phase angle) are nonlinear.By the way, at little load frequency the velocity of value changing is bigger than at big one.So, most significant changes occur at load cases 1-3 (frequency 0.1; 0.5; 1.0 Hz).