Finite element numerical implementation method of linear and nonlinear viscoelastic model of asphalt mixture

Abstract

The application relates to a finite element numerical implementation method of a linear and nonlinear viscoelastic model of asphalt mixture, and belongs to the technical field of constitutive behavior of asphalt mixture. The application aims at the problem that there is no general and feasible numerical calculation method for the linear and nonlinear viscoelastic model of asphalt mixture. The application obtains the complex frequency domain modulus expression of the linear and nonlinear viscoelastic model of asphalt mixture based on Laplace transform; then, the lowest-order differential expression of stress and strain of a one-dimensional model is constructed through equivalent substitution between the time domain and the complex frequency domain; the stress updating expression of a discrete form model is constructed based on the differential term in the differential form of the GL discrete stress and strain; the stress updating expression and the Jacobian matrix expression of a three-dimensional model are derived based on the constant or non-constant Poisson ratio assumption and the incremental iteration method; and the finite element numerical calculation is realized by developing a user material subroutine by using FORTRAN based on the subroutine development interface in ABAQUS. The application realizes the finite element numerical calculation of the viscoelastic model of asphalt mixture.

Description

TECHNICAL FIELD

[0001] The present application relates to a finite element numerical implementation method of linear and nonlinear viscoelastic models of asphalt mixture, and belongs to the technical field of constitutive behavior of asphalt mixture. BACKGROUND

[0002] Asphalt mixture is a typical viscoelastic material, and its mechanical behavior is closely related to temperature and load frequency, and conforms to the principle of time-temperature equivalence. Researchers obtain the viscoelastic parameter values of asphalt mixture through creep tests or dynamic modulus tests, and construct linear and nonlinear viscoelastic models of asphalt mixture to characterize the mechanical behavior of asphalt mixture under different loading modes.

[0003] At present, typical linear and nonlinear viscoelastic models of asphalt mixture include integer order viscoelastic model, fractional order viscoelastic model, modified Burgers model and five-parameter nonlinear viscoelastic model. Among them, the integer order viscoelastic model, i.e. linear viscoelastic model, includes Kelvin model, Burgers model, three-parameter solid model, generalized Kelvin model and generalized Maxwell model. The integer order model is composed of Newtonian viscous pot and spring in series and parallel, which can be directly converted into Prony series, which is convenient for finite element numerical calculation, and is widely used. However, the integer order model is difficult to characterize the dynamic mechanical behavior of asphalt mixture in a wide frequency range, and it is easy to be affected by noise in the parameter identification process, and the model parameters identified may be negative, thus losing physical meaning. The fractional order viscoelastic model is composed of Abel viscous pot and spring in series and parallel, which can characterize the mechanical behavior of asphalt mixture in a wide frequency range. However, the fractional order viscoelastic model lacks simple creep and relaxation modulus expressions, which limits its application in finite element analysis. The modified Burgers model is a class of nonlinear viscoelastic models based on the Burgers model, which replaces the series viscous pot with an exponential viscous pot. Unlike the Burgers model, the creep deformation of the modified Burgers model increases indefinitely with time, so it is widely used in pavement rut analysis. However, the modified Burgers model can only describe stable creep and deceleration creep, and cannot describe acceleration creep. The five-parameter nonlinear viscoelastic model is a series of variable cross-section viscous pot based on the Burgers model. The stress of the viscous pot can be expressed as the third derivative of strain with respect to time, so it can describe the three stages of creep deceleration, stable creep and accelerated creep.

[0004] For the above different types of linear and nonlinear viscoelastic model, it is necessary to develop a finite element numerical algorithm for a specific model alone, and the calculation process is more complex and the calculation cost is higher. For integer order viscoelastic model, the exponential algorithm is usually used for calculation. When the model contains nonlinear viscoelastic elements such as Abel fractional order viscoelastic, exponential viscoelastic and variable cross-section viscoelastic, it will be difficult to use the exponential algorithm for finite element calculation. The fractional order viscoelastic model usually uses the Mittag-Leffler function method for calculation, which first represents the time domain creep compliance or relaxation modulus of the fractional order viscoelastic model with the Mittag-Leffler function, and then constructs the stress increment expression. This method is only applicable to simple fractional order viscoelastic models whose time domain creep compliance or relaxation modulus can be represented by the Mittag-Leffler function. For complex fractional order viscoelastic models, it is difficult to derive the time domain expression of their creep compliance or relaxation modulus. In addition, this method requires the use of an approximate calculation formula of the Mittag-Leffler function during calculation, which has a large amount of calculation and the calculation accuracy is affected by the approximate calculation formula. The nonlinear viscoelastic model with exponential viscoelastic element usually uses the implicit Euler method for calculation, which uses the constant stress assumption to obtain the expression of the creep strain rate, and solves the creep strain increment in each increment step by reciprocal iteration, and then solves the stress increment, so it is only applicable to the case where the external load is constant, and the calculation amount is large. For the five-parameter nonlinear viscoelastic model used to describe the three-stage creep curve, there is no research to propose a finite element numerical calculation method for this model. In addition, for other types of nonlinear viscoelastic models composed of integer order, fractional order and variable cross-section viscoelastic elements, there is still a lack of effective numerical method to realize the finite element calculation of this type of model.

[0005] Therefore, it is urgent to propose a simple and feasible numerical calculation method to use various types of asphalt mixture linear and nonlinear viscoelastic models for asphalt pavement finite element analysis. SUMMARY

[0006] In order to solve the problem that there is no general and feasible numerical calculation method for asphalt mixture linear and nonlinear viscoelastic model, the present application provides a finite element numerical implementation method for asphalt mixture linear and nonlinear viscoelastic model.

[0007] The finite element numerical implementation method for asphalt mixture linear and nonlinear viscoelastic model provided by the present application comprises,

[0008] Based on the stress-strain relationship of each sub-element in the linear and nonlinear viscoelastic model and the Laplace transform, the total stress and total strain relationship formula in the complex frequency domain is established;

[0009] The time domain and complex frequency domain derivative equivalent substitution is introduced, and the one-dimensional form of the lowest order differential expression of the total stress and total strain relationship formula of the model is constructed;

[0010] The Grunwald-Letnikov method is used to discretize the differential term in the lowest order differential expression of the one-dimensional form, and a stress update expression of the discrete form model is constructed;

[0011] Combined with the normal Poisson ratio assumption or the very Poisson ratio assumption and the incremental iteration method, a stress update expression and a Jacobian matrix update expression of the three-dimensional form model are constructed;

[0012] Based on the stress update expression and the Jacobian matrix update expression of the three-dimensional form model, a user material subroutine (UMAT) of the linear and nonlinear viscoelastic model is developed by using the subroutine development interface in ABAQUS and FORTRAN; after a finite element geometric model is created in ABAQUS, the linear and nonlinear viscoelastic model material parameters, loads and boundary conditions are input after meshing, and the user material subroutine UMAT is called, so that the linear and nonlinear viscoelastic model finite element numerical calculation of the asphalt mixture is realized.

[0013] The method of the present application first obtains the complex frequency domain modulus expression of the linear and nonlinear viscoelastic model of the asphalt mixture based on the Laplace transform; then, the lowest order differential expression of the stress and strain of the linear / nonlinear viscoelastic model of the asphalt mixture in one-dimensional form is constructed through time domain / complex frequency domain equivalent substitution; on this basis, the differential term in the differential form of the stress and strain is discretized based on the Grunwald-Letnikov (GL) method, so that the stress update expression of the discrete form of the linear / nonlinear viscoelastic model of the asphalt mixture is constructed; thereafter, based on the normal / very Poisson ratio assumption and the incremental iteration method, the stress update expression and the Jacobian matrix expression of the three-dimensional form of the linear / nonlinear viscoelastic model of the asphalt mixture are derived; finally, based on the subroutine development interface in ABAQUS, the user material subroutine (UMAT) is developed by using FORTRAN, so that the finite element numerical calculation of the linear / nonlinear viscoelastic model of the asphalt mixture is realized.

[0014] Compared with the prior art, the method of the present application has the following advantages:

[0015] (1) The integer order viscoelastic model usually uses an exponential algorithm for calculation, which requires that the creep or relaxation modulus of the viscoelastic model has a natural exponential form, and is therefore suitable for the linear viscoelastic model containing only a Newtonian dashpot. The method of the present application uses a differential form constitutive equation, which can not only be applied to the linear viscoelastic model containing only a Newtonian dashpot, but also be used for the numerical calculation of the nonlinear viscoelastic model composed of a Newtonian dashpot and other nonlinear dashpots. The applicable nonlinear dashpots include: Abel fractional order dashpot, exponential dashpot and variable cross-section dashpot.

[0016] (2) Fractional viscoelastic model usually uses the Mittag-Leffler function method to calculate, which firstly uses the Mittag-Leffler function to express the time-domain creep compliance or relaxation modulus of the fractional viscoelastic model, and then constructs the stress increment expression. This method is only applicable to simple fractional viscoelastic models whose time-domain creep compliance or relaxation modulus can be expressed by the Mittag-Leffler function. For complex fractional viscoelastic models, it is difficult to derive the time-domain expression of their creep compliance or relaxation modulus. In addition, this method needs to use the approximate calculation formula of the Mittag-Leffler function in calculation, and the calculation amount is large, and the calculation accuracy is affected by the approximate calculation formula. The method of the present application uses the constitutive equation in differential form, which does not need to derive the time-domain expression of the model creep compliance or relaxation modulus, and can be applied to the finite element numerical calculation of complex fractional viscoelastic models.

[0017] (3) The nonlinear viscoelastic model with exponential viscopot has usually used the implicit Euler method for calculation, which uses the constant stress assumption to obtain the expression of the creep strain rate, and solves the creep strain increment in each increment step by reciprocal iteration, and then solves the stress increment, so it is only applicable to the case where the external load is constant, and the calculation amount is large. The method of the present application does not need to solve the stress increment by reciprocal iteration in each increment step, and has high calculation efficiency. In addition, the method of the present application does not use the assumption of constant external load in calculation, so it can be applied to the finite element numerical calculation of nonlinear viscoelastic models under various loads.

[0018] (4) For nonlinear viscoelastic models containing cross-sectional viscopot, such as five-parameter nonlinear viscoelastic models used to describe three-stage creep curves, the prior art has not yet proposed a finite element numerical calculation method for the model. The method of the present application can realize the finite element numerical calculation of nonlinear viscoelastic models containing cross-sectional viscopot.

[0019] (5) The method of the present application is used for solving the finite element numerical calculation of linear / nonlinear viscoelastic models containing Newton viscopot, Abel viscopot, exponential viscopot or variable cross-sectional viscopot, which has a general solving step and calculation format, and does not need to design a special solving algorithm for the type of viscopot, so it is more convenient for finite element implementation than other algorithms. BRIEF DESCRIPTION OF DRAWINGS

[0020] Figure 1 is a flow chart of the finite element numerical implementation method of the linear and nonlinear viscoelastic model of the asphalt mixture described in the present application;

[0021] Figure 2 is a UMAT development flow chart of the linear and nonlinear viscoelastic model of the asphalt mixture;

[0022] Figure 3 is a schematic diagram of the Burgers model;

[0023] Figure 4 It is the creep curve of the Burgers model;

[0024] Figure 5 This is a schematic diagram of a single-axis loading of the Burgers model;

[0025] Figure 6 This is a schematic diagram of the modified fractional Zener model;

[0026] Figure 7 It is a correction of the creep curve of the fractional Zener model;

[0027] Figure 8 This is a schematic diagram of the modified Burgers model;

[0028] Figure 9 It is a modified Burgers model creep curve;

[0029] Figure 10 This is a schematic diagram of a five-parameter nonlinear viscoelastic model;

[0030] Figure 11 It is a creep curve of a five-parameter nonlinear viscoelastic model. Detailed Implementation

[0031] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0032] It should be noted that, unless otherwise specified, the embodiments and features described in the present invention can be combined with each other.

[0033] The present invention will be further described below with reference to the accompanying drawings, but this should not be construed as limiting the invention.

[0034] Combination Figure 1 and Figure 2 As shown, this invention provides a finite element numerical implementation method for linear and nonlinear viscoelastic models of asphalt mixtures, including,

[0035] Based on the stress-strain relationship and Laplace transform of each sub-element in the linear and nonlinear viscoelastic models, the relationship between the total stress and total strain of the model in the complex frequency domain is established.

[0036] By introducing equivalent substitutions for time-domain and complex frequency-domain differentiation, a one-dimensional lowest-order differential expression for the relationship between total stress and total strain in the model is constructed.

[0037] The Grünwald-Letnikov method is used to discretize the differential term in the lowest order differential expression of one-dimensional form, and a stress updating expression of the discrete form model is constructed.

[0038] The three-dimensional form model stress updating expression and Jacobian matrix updating expression are constructed by combining the constant Poisson ratio assumption or the very high Poisson ratio assumption with the incremental iteration method.

[0039] Based on the three-dimensional form model stress updating expression and Jacobian matrix updating expression, a user material subroutine UMAT of the linear and nonlinear viscoelastic model is developed by using the FORTRAN language and the subroutine development interface in ABAQUS; a finite element geometric model is created in ABAQUS, and after meshing, the material parameters, load and boundary conditions of the linear and nonlinear viscoelastic model are input, and the user material subroutine UMAT is called, so that the linear and nonlinear viscoelastic model finite element numerical calculation of the asphalt mixture is realized.

[0040] The embodiment realizes the finite element numerical calculation of the linear and nonlinear viscoelastic model including the integer order viscoelastic model with Newtonian dashpot (the dashpot stress is the first order derivative of the dashpot strain multiplied by the viscosity), the fractional order viscoelastic model with Abel dashpot (the dashpot stress is the fractional order derivative of the dashpot strain multiplied by the viscosity), the nonlinear viscoelastic model with exponential dashpot (the compliance of the dashpot is an exponential function), and the nonlinear viscoelastic model with variable cross-section dashpot (the dashpot stress is the integer order derivative higher than two of the dashpot strain multiplied by the viscosity).

[0041] The method of the application is described below in combination with four types of linear / nonlinear viscoelastic models. The four types of viscoelastic models include: Burgers model (linear viscoelastic model), fractional order modified Zener model (nonlinear viscoelastic model), modified Burgers model (nonlinear viscoelastic model), and five-parameter nonlinear viscoelastic model (nonlinear viscoelastic model).

[0042] Embodiment one: in combination with Figure 3 As shown, the linear viscoelastic model is selected as the Burgers model, and the sub-elements include spring element one, spring element two, linear Newtonian dashpot element one and linear Newtonian dashpot element two.

[0043] Firstly, based on the Laplace transform, the total stress and total strain relationship formula in the complex frequency domain of the model is obtained, and the differential constitutive form of the model is further obtained. On this basis, the stress updating and Jacobian matrix updating expressions of the model are derived based on the GL method, and the UMAT subroutine is developed. Finally, the finite element model is constructed based on ABAQUS, and the numerical calculation of the model is carried out. The specific steps are as follows:

[0044] Based on the stress-strain relationship and Laplace transform of the Burgers model sub-element, the total stress and total strain relationship formula in the complex frequency domain of the model is established as:

[0045]

[0046] wherein is the complex frequency domain modulus, E1 is the modulus value of the spring element one, E2 is the modulus value of the spring element two, η1 is the viscosity of the linear Newtonian viscous element one, η2 is the viscosity of the linear Newtonian viscous element two, and s is the complex frequency domain independent variable;

[0047] The formula (1) is simplified, and the time domain and the complex frequency domain derivative equivalent substitution with the initial value of zero are introduced:

[0048]

[0049] wherein L is the Laplace transform, is the derivative operator with respect to time, α is the derivative order, f(t) is the function to be derived, and t is the time;

[0050] The one-dimensional form of the lowest order differential expression of the total stress and total strain relationship of the model is constructed as:

[0051]

[0052] wherein ε(t) is the total strain, and σ(t) is the total stress.

[0053] The Grünwald-Letnikov method shown in the formula (4) is used to discretize the differential term in the one-dimensional form of the lowest order differential expression:

[0054]

[0055] wherein Δt is the equidistant time interval, and k is the number of equidistant scattering intervals of the time t; is the Grünwald coefficient, and

[0056] The stress update expression of the discrete form model is constructed as:

[0057]

[0058] wherein A, C G , D K , and B are intermediate variables:

[0059]

[0060] The formula (6) is converted into a three-dimensional form by using the constant Poisson's ratio assumption or the non-constant Poisson's ratio assumption.

[0061] If the constant Poisson's ratio assumption is used, the three-dimensional form conversion relationship of the intermediate variables A, C G , D K , and B is:

[0062]

[0063] where I = 1, 2; J = 1, 2; K as a superscript, indicates that the variable is related to the normal strain or normal stress, G as a superscript, indicates that the variable is related to the deviatoric strain or deviatoric stress; μ is a constant Poisson's ratio;

[0064] If the assumption of a constant Poisson's ratio is adopted, A G , A K , B G , B K , and are directly measured by triaxial tests;

[0065] Based on the incremental iteration method in ABAQUS, the three-dimensional form of the stress updating expression (8) and the Jacobian matrix updating expression (9) are obtained:

[0066]

[0067] where oo represents the normal stress or normal strain component in different directions, and wu represents the shear stress or shear strain component in different directions; is the average strain, is the average stress;

[0068]

[0069] where is the stress increment in the pq direction at the k+1 time point, and pq takes the values of 11, 22, 33, 12, 13, and 23; is the strain increment in the PQ direction at the k+1 time point, and PQ takes the values of 11, 22, 33, 12, 13, and 23; δ pq , δ PQ , δ pQ , and δ qP are the Kronecker symbols;

[0070] The material parameters are input in the user material subroutine UMAT, a finite element geometric model is created using ABAQUS, three-dimensional solid elements are selected to discretize the finite element geometric model, and loads and boundary conditions are applied to the discretized finite element geometric model; after setting the incremental step, the operation is submitted to obtain the mechanical response results of the finite element model.

[0071] As an example, based on the developed UMAT of the Burgers model, the material parameters of the corresponding variables in Table 1 are input;

[0072] Table 1 Material parameter table of linear / nonlinear viscoelastic models

[0073]

[0074] Create using ABAQUS Figure 5 The finite element model shown adopts a three-dimensional 8-node reduced integral element C3D8R discrete finite element geometric model. The normal displacement of the back nodes of the three-dimensional 8-node reduced integral element in the x-direction is fixed, and the stress load σ in the x-direction shown in Equation (10) is applied. 11 The increment step was set to 0.1s, and the analysis step to 100s. The resulting creep curve in the x-direction of the model is shown below. Figure 4 As shown. The theoretical solution of its creep curve can be obtained by combining the complex frequency domain expression of equation (1) with that of equation (10) and then performing an inverse Laplace numerical transform.

[0075] σ 11 =σ0(t(U(t)-U(t-t0))+U(t-t0)) (10).

[0076] In the formula, σ0 is the stress load amplitude, σ0 = 1 MPa, t0 is the time for linear load growth, t0 = 1 s; and U is the unit step function.

[0077] Example 2: Combination Figure 6 As shown, a modified fractional Zener model is selected, and the sub-elements include spring element one, spring element two, Abel sticky pot element one, and Abel sticky pot element two;

[0078] Based on the stress-strain relationship of the sub-elements of the modified fractional Zener model and the Laplace transform, the relationship between the total stress and total strain of the model in the complex frequency domain is established as follows:

[0079]

[0080] In the formula Let be the modulus in the complex frequency domain, E1 be the modulus of spring element one, E2 be the modulus of spring element two, η1 be the viscosity of Abel sticky pot element one, η2 be the viscosity of Abel sticky pot element two, s be the independent variable in the complex frequency domain, α1 be the fractional order of Abel sticky pot element one, α2 be the fractional order of Abel sticky pot element two, and α2 > α1.

[0081] The simplified stress and strain relationships of the modified fractional Zener model in the complex frequency domain are obtained by introducing equivalent substitutions in the time and complex frequency domains with initial values ​​of zero:

[0082]

[0083] In the formula, L is the Laplace transform. is the α-order derivative operator with respect to time, where α is the order of the derivative, f(t) is the function to be differentiated, and t is time;

[0084] The one-dimensional form of the lowest order differential expression of the model of the total stress and total strain relationship is as follows:

[0085]

[0086] wherein ε(t) is the total strain, and σ(t) is the total stress.

[0087] The Grünwald-Letnikov method shown in formula (24) is used to discretize the differential term in the one-dimensional form of the lowest order differential expression:

[0088]

[0089] wherein Δt is the equidistant time interval, and k is the number of equidistant dispersion intervals of time t; is the Grünwald coefficient, and

[0090] The stress updating expression of the discrete form model is constructed as follows:

[0091]

[0092] wherein A, C i , D i , and B are intermediate variables:

[0093]

[0094] If a constant Poisson's ratio is assumed, the three-dimensional form conversion relationship of the intermediate variables A, C i , D i , and B is as follows:

[0095]

[0096] wherein I=1, 2; J=1, 2; K is used as an index, indicating that the variable is related to the normal strain or normal stress, and G is used as an index, indicating that the variable is related to the deviatoric strain or deviatoric stress; and μ is the constant Poisson's ratio.

[0097] If a non-constant Poisson's ratio is assumed, A G , A K , B G , B K , and are directly measured through a triaxial test.

[0098] Based on the incremental iteration method in ABAQUS, the three-dimensional form of the stress updating expression (28) and the Jacobian matrix updating expression (29) of the model are obtained as follows:

[0099]

[0100] where oo represents different directional normal stress or normal strain components, and wu represents different directional shear stress or shear strain components; is the average strain, is the average stress;

[0101]

[0102] where is the stress increment in the pq direction at the k+1 time point, and pq takes the values 11, 22, 33, 12, 13, 23; is the strain increment in the PQ direction at the k+1 time point, and PQ takes the values 11, 22, 33, 12, 13, 23; pq , δ PQ , δ pQ , and δ qP is the Kronecker symbol;

[0103] The material parameters are input in the user material subroutine UMAT, a finite element geometric model is created using ABAQUS, three-dimensional solid elements are selected to discretize the finite element geometric model, and the discretized finite element geometric model is subjected to a load and boundary conditions; an increment step is set and the operation is submitted to obtain the mechanical response results of the finite element model.

[0104] As an example, based on the developed UMAT of the modified fractional Zener model, the material parameters in Table 1 are input, a finite element model as shown in Figure 5 is created using ABAQUS, a three-dimensional 8-node reduced integration element C3D8R is used to discretize the finite element geometric model, the normal displacement of the back nodes of the three-dimensional 8-node reduced integration element in the x direction is fixed, and the stress load σ 11 in the x direction shown in equation (10) is applied; the increment step is set to 0.1 s, the analysis step is 100 s, and finally the x direction creep curve of the model is obtained as shown in Figure 7 The theoretical solution of the creep curve can be obtained by Laplace numerical inverse transformation from equation (21) combined with the complex frequency domain expression of equation (10).

[0105] Example Three: in combination with Figure 8 as shown, the modified Burgers model is selected, and the sub-element includes a spring element one, a spring element two, a linear dashpot element, and an exponential dashpot element;

[0106] Based on the stress-strain relationship of the modified Burgers model sub-element and the Laplace transformation, the total stress and total strain relationship formula in the complex frequency domain is:

[0107]

[0108] where E1 is the modulus value of the spring element one, a is the parameter one of the exponential dashpot element, s is the complex frequency domain argument, b is the parameter two of the exponential dashpot element, η2 is the viscosity of the linear dashpot element, E2 is the modulus value of the spring element two;

[0109] The total stress and total strain relationship of the simplified complex frequency domain modified Burgers model is obtained, and the initial value zero time domain and complex frequency domain derivative equivalent substitution is introduced:

[0110]

[0111] In the formula, L is the Laplace transform, is the α-order derivative operator with respect to time, α is the derivative order, f(t) is the function to be derived, and t is time;

[0112] The one-dimensional form of the lowest order differential expression of the total stress and total strain relationship of the model is constructed as:

[0113]

[0114] In the formula, ε(t) is the total strain, and σ(t) is the total stress.

[0115] The Grünwald-Letnikov method shown in formula (34) is used to discretize the differential term in the one-dimensional form of the lowest order differential expression:

[0116]

[0117] In the formula, Δt is the equidistant time interval, and k is the number of equidistant scattering intervals of time t; is the Grünwald coefficient, and

[0118] The discrete form of the model stress update expression is constructed as:

[0119]

[0120] In the formula, A, C i , D i , and B are intermediate variables:

[0121]

[0122] If the constant Poisson ratio assumption is adopted, the three-dimensional form conversion relationship of the intermediate variables A, C i , D i , and B is:

[0123]

[0124] In the formula, I = 1, 2; K is a subscript indicating that the variable is related to normal strain or normal stress, G is a subscript indicating that the variable is related to deviatoric strain or deviatoric stress; μ is the constant Poisson's ratio;

[0125] If the assumption of a very Poisson's ratio is adopted, A in formula (37) G A K B G B K , and Measured directly through triaxial testing;

[0126] Based on the incremental iteration method in ABAQUS, the stress update expression (38) and Jacobian matrix update expression (39) of the three-dimensional formal model are obtained:

[0127]

[0128] In the formula, oo represents the normal stress or normal strain components in different directions, and wu represents the shear stress or shear strain components in different directions; For average strain, The average stress;

[0129]

[0130] In the formula Let pq be the stress increment in the pq direction at the (k+1)th time point, where pq takes values ​​of 11, 22, 33, 12, 13, and 23. Let δ be the strain increment in the PQ direction at the (k+1)th time point, where PQ takes values ​​of 11, 22, 33, 12, 13, and 23; pq δ PQ δ pQ and δ qP The symbol for Kronecker;

[0131] Input material parameters in the user material subroutine UMAT, create a finite element geometric model using ABAQUS, select three-dimensional solid elements to discretize the finite element geometric model, apply loads and boundary conditions to the discretized finite element geometric model, set incremental steps and submit the calculation to obtain the mechanical response results of the finite element model.

[0132] As an example, based on the developed modified Burgers model, UMAT is created using ABAQUS by inputting material parameters as shown in Table 1. Figure 5 The finite element model shown adopts a three-dimensional 8-node reduced integral element C3D8R discrete finite element geometric model. The normal displacement of the back nodes of the three-dimensional 8-node reduced integral element in the x-direction is fixed, and the stress load σ in the x-direction shown in Equation (10) is applied. 11The increment step was set to 0.1s, and the analysis step to 100s. The resulting creep curve in the x-direction of the model is shown below. Figure 9 As shown. The theoretical solution of its creep curve can be obtained by combining the complex frequency domain expression of equation (31) with that of equation (10) through the inverse Laplace numerical transform.

[0133] Example 4: Combination Figure 10 As shown, a five-parameter nonlinear viscoelastic model is selected, and the sub-components include spring element one, spring element two, linear sticky pot element one, linear sticky pot element two, and variable cross-section sticky pot element.

[0134] Based on the stress-strain relationship of the sub-elements of the five-parameter nonlinear viscoelastic model and the Laplace transform, the relationship between the total stress and total strain of the model in the complex frequency domain is established as follows:

[0135]

[0136] In the formula E1 is the modulus in the complex frequency domain, E2 is the modulus of spring element one, η1 is the viscosity of linear sticky pot element one, η2 is the viscosity of linear sticky pot element two, η3 is the viscosity of the variable cross-section sticky pot element, and s is the independent variable in the complex frequency domain.

[0137] The simplified stress and strain relationships of the five-parameter nonlinear viscoelastic model in the complex frequency domain are equivalently replaced by time-domain and complex frequency-domain differentiation with initial values ​​of zero:

[0138]

[0139] In the formula, L is the Laplace transform. is the α-order derivative operator with respect to time, where α is the order of the derivative, f(t) is the function to be differentiated, and t is time;

[0140] The lowest-order one-dimensional differential expression for the relationship between total stress and total strain in the model is as follows:

[0141]

[0142] In the formula, ε(t) is the total strain and σ(t) is the total stress.

[0143] Discretize the differential terms in the lowest-order differential expression in one-dimensional form using the Grünwald-Letnikov method shown in formula (44):

[0144]

[0145] In the formula, Δt is the equidistant time interval, and k is the number of equidistant intervals of time t; The Grünwald coefficient, and

[0146]

[0147] The stress update expression of the discrete form model is constructed as:

[0148]

[0149] where A, C i , D i , and B are intermediate variables:

[0150]

[0151] If the constant Poisson's ratio assumption is adopted, the three-dimensional form conversion relationship of the intermediate variables A, C i , D i , and B is:

[0152]

[0153] where I = 1, 2; J = 1, 2, 3; K as a subscript, indicates that the variable is related to the normal strain or normal stress, and G as a subscript, indicates that the variable is related to the deviatoric strain or deviatoric stress; μ is the constant Poisson's ratio;

[0154] If the non-constant Poisson's ratio assumption is adopted, A G , A K , B G , B K , and are directly measured through triaxial tests;

[0155] Based on the incremental iteration method in ABAQUS, the stress update expression (48) and the Jacobian matrix update expression (49) of the three-dimensional form model are obtained:

[0156]

[0157] where oo represents the normal stress or normal strain component in different directions, and wu represents the shear stress or shear strain component in different directions; is the average strain, is the average stress;

[0158]

[0159] where is the stress increment in the pq direction at the k+1 time point, and pq takes values of 11, 22, 33, 12, 13, and 23; is the strain increment in the PQ direction at the k+1 time point, and PQ takes values of 11, 22, 33, 12, 13, and 23; pq , δ PQ , δpQ and δ qP is the Kronecker symbol.

[0160] The material parameters are input in the user material subroutine UMAT, the finite element geometric model is created by ABAQUS, the three-dimensional solid element is selected to discretize the finite element geometric model, the load and boundary conditions are applied to the discretized finite element geometric model, the operation is submitted after setting the increment step, and the mechanical response result of the finite element model is obtained.

[0161] As an example, the UMAT based on the developed modified Burgers model is used, the material parameters in Table 1 are input, the finite element model shown in FIG. 1 is created by ABAQUS, the three-dimensional 8-node reduced integration element C3D8R is used to discretize the finite element geometric model, the normal displacement of the back node of the three-dimensional 8-node reduced integration element in the x direction is fixed, and the stress load σ 11 in the x direction shown in formula (10) is applied. Figure 5 The increment step is set to 0.1 s, the analysis step is set to 100 s, and finally the x direction creep curve of the model is obtained as shown in FIG. 2. The theoretical solution of the creep curve can be obtained by the Laplace numerical inverse transformation of formula (41) combined with the complex frequency domain expression of formula (10). Figure 11

[0162] While the application has been described with reference to particular embodiments, it will be understood that the examples are merely illustrative of the principles and applications of the present application. It will be understood that various modifications can be made to the illustrative embodiments, and other arrangements can be devised without departing from the spirit and scope of the present application as defined by the appended claims. It will be understood that the features of the various embodiments can be combined with each other, in different ways than those described. It will be understood that features described with respect to one embodiment can be used in other embodiments.​

Claims

1. A method for numerical implementation of linear and nonlinear viscoelastic models of asphalt mixture by finite elements, characterized in that, The application relates to a linear and nonlinear viscoelastic model finite element numerical implementation method for asphalt mixture. Based on stress-strain relations of each sub-element in linear and nonlinear viscoelastic models and Laplace transform, a total stress and total strain relation formula in complex frequency domain is established; A one-dimensional form lowest-order differential expression of the total stress and total strain relation formula of the model is constructed by introducing time domain and complex frequency domain derivation equivalent substitution; A discrete form model stress updating expression is constructed by discretizing differential terms in the one-dimensional form lowest-order differential expression by using a Grunwald-Letnikov method; A three-dimensional form model stress updating expression and Jacobian matrix updating expression are constructed by combining a constant Poisson ratio assumption or a non-constant Poisson ratio assumption and an incremental iteration method; Based on the three-dimensional form model stress updating expression and Jacobian matrix updating expression, a user material subroutine UMAT of the linear and nonlinear viscoelastic model is developed by using a FORTRAN language through a subroutine development interface in ABAQUS; a finite element geometric model is created in ABAQUS, and after meshing, material parameters, loads and boundary conditions of the linear and nonlinear viscoelastic model are input, and the user material subroutine UMAT is called, so that linear and nonlinear viscoelastic model finite element numerical calculation of asphalt mixture is realized.

2. The method of claim 1, wherein, The linear viscoelastic model is a Burgers model, and the sub-elements include spring element one, spring element two, linear Newtonian dashpot element one and linear Newtonian dashpot element two. The total stress and total strain relation formula in complex frequency domain is established as follows: wherein E1is the modulus value of spring element one, E2is the modulus value of spring element two, η1is the viscosity of linear Newtonian viscous element one, η2is the viscosity of linear Newtonian viscous element two, and s is the complex frequency domain argument; The time domain and complex frequency domain derivation equivalent substitution with an initial value of zero is introduced as follows: where L is a Laplace transform, is an αth derivative operator with respect to time, α is the order of the derivative, f(t) is the function to be differentiated, and t is time. The one-dimensional form lowest-order differential expression of the total stress and total strain relation formula of the model is constructed as follows: In the formula, epsilon (t) is total strain, and sigma (t) is total stress.

3. The linear and nonlinear viscoelastic model finite element numerical implementation method for asphalt mixture according to claim 2, wherein the differential terms in the one-dimensional form lowest-order differential expression are discretized by using a Grunwald-Letnikov method shown in formula (4): The discrete form model stress updating expression is constructed as follows: where Δt is the equidistance time interval, k is the number of equidistance dispersion intervals at time t; is the Grunwald coefficient, and The three-dimensional form conversion of formula (6) is performed by using a constant Poisson ratio assumption or a non-constant Poisson ratio assumption. wherein A, C i , D i , B are intermediate variables:

4. The method of claim 3, wherein, In the formula, I=1, 2; J=1, 2; K is an index indicating that a variable is related to normal strain or normal stress, and G is an index indicating that a variable is related to deviatoric strain or deviatoric stress; and mu is a constant Poisson ratio. If the constant Poisson's ratio is assumed, the three-dimensional form conversion relationship of the intermediate variables A, C i , D i , B is as follows: Based on the incremental iteration method, the three-dimensional form model stress updating expression (8) and the Jacobian matrix updating expression (9) are obtained: If the very Poisson's ratio assumption is adopted, A in equation (7) is given by G , K , G , K , and are directly measured through triaxial tests; Material parameters are input in the user material subroutine UMAT, a finite element geometric model is created by using ABAQUS, a three-dimensional solid element is selected to discretize the finite element geometric model, loads and boundary conditions are applied to the discretized finite element geometric model, an incremental step is set, and then operation is submitted to obtain mechanical response results of the finite element model. where oo represents the different directional normal stress or normal strain components and wu represents the different directional shear stress or shear strain components; is the average strain, is the average stress; wherein is the stress increment in the pq direction at the k+1 time point, pq taking the values 11, 22, 33, 12, 13, 23; is the strain increment in the PQ direction at the k+1 time point, PQ taking the values 11, 22, 33, 12, 13, 23; δ pq , δ PQ , δ pQ and δ qP are Kronecker symbols; The modified fractional Zener model is selected, and the sub-elements include spring element one, spring element two, Abel dashpot element one and Abel dashpot element two.

5. The method of claim 1, wherein, The total stress and total strain relation formula in complex frequency domain is established as follows: The time domain and complex frequency domain derivation equivalent substitution with an initial value of zero is introduced as follows: wherein E1is the modulus value of the spring element one, E2is the modulus value of the spring element two, η1is the viscosity of the Abel dashpot element one, η2is the viscosity of the Abel dashpot element two, s is the complex frequency domain argument, α1is the fractional order of the Abel dashpot element one, and α2is the fractional order of the Abel dashpot element two, with α2> α1. The one-dimensional form lowest-order differential expression of the total stress and total strain relation formula of the model is constructed as follows: where L is the Laplace transform, is the αth derivative operator with respect to time, α is the order of the derivative, f(t) is the function to be differentiated, and t is time. ​ where ε(t) is the total strain and σ(t) is the total stress.

6. The method of claim 5, wherein, The Grünwald-Letnikov method shown in formula (24) is used to discretize the differential term in the one-dimensional form of the lowest order differential expression of the model: where Δt is the equidistance time interval, k is the number of equidistance dispersion intervals at time t; is the Grunwald coefficient, and The stress updating expression of the discrete form of the model is constructed as: wherein A, C i , D i , B are intermediate variables: If the constant Poisson's ratio is assumed, the three-dimensional form conversion relationship of the intermediate variables A, C i , D i , B is as follows: where I=1, 2; J=1, 2; K is a subscript, indicating that the variable is related to the normal strain or the normal stress, G is a subscript, indicating that the variable is related to the deviatoric strain or the deviatoric stress; μ is a constant Poisson's ratio; If the very Poisson's ratio assumption is adopted, A in equation (27) G , K , G , K , and are directly measured through triaxial tests; Based on the incremental iteration method, the stress updating expression (28) and the Jacobian matrix updating expression (29) of the three-dimensional form of the model are obtained: where oo represents the different directional normal stress or normal strain components and wu represents the different directional shear stress or shear strain components; is the average strain, is the average stress; wherein is the stress increment in the pq direction at the k+1 time point, pq taking the values 11, 22, 33, 12, 13, 23; is the strain increment in the PQ direction at the k+1 time point, PQ taking the values 11, 22, 33, 12, 13, 23; δ pq , δ PQ , δ pQ , and δ qP are Kronecker symbols; The material parameters are input in the user material subroutine UMAT, a finite element geometric model is created by using ABAQUS, three-dimensional solid elements are selected to discretize the finite element geometric model, and the discrete finite element geometric model is subjected to loads and boundary conditions; after setting the incremental step, the operation is submitted to obtain the mechanical response results of the finite element model.

7. The method of claim 1, wherein, The modified Burgers model is selected, and the sub-elements include spring element one, spring element two, linear dashpot element and exponential dashpot element; The relationship between the total stress and the total strain of the model in the complex frequency domain is: wherein E1 is the modulus value of the spring element one, a is the parameter one of the exponential dashpot element, s is the complex frequency domain argument, b is the parameter two of the exponential dashpot element, η2 is the viscosity of the linear dashpot element, and E2 is the modulus value of the spring element two. The initial value of zero is introduced into the time domain and the complex frequency domain to perform derivative equivalent substitution: where L is a Laplace transform, is an αth derivative operator with respect to time, α is the order of the derivative, f(t) is the function to be differentiated, and t is time. The one-dimensional form of the lowest order differential expression of the relationship between the total stress and the total strain of the model is constructed as: where ε(t) is the total strain and σ(t) is the total stress.

8. The method of claim 7, wherein, The Grünwald-Letnikov method shown in formula (34) is used to discretize the differential term in the one-dimensional form of the lowest order differential expression of the model: where Δt is the equidistance time interval, k is the number of equidistance dispersion intervals at time t; is the Grunwald coefficient, and The stress updating expression of the discrete form of the model is constructed as: wherein A, C i , D i , B are intermediate variables: If the constant Poisson's ratio is assumed, the three-dimensional form conversion relationship of the intermediate variables A, C i , D i , B is as follows: where I=1, 2; K is a subscript, indicating that the variable is related to the normal strain or the normal stress, G is a subscript, indicating that the variable is related to the deviatoric strain or the deviatoric stress; μ is a constant Poisson's ratio; If the very Poisson's ratio assumption is adopted, A in equation (37) G , A K , B G , B K , and are directly measured through triaxial tests; Based on the incremental iteration method, the stress updating expression (38) and the Jacobian matrix updating expression (39) of the three-dimensional form of the model are obtained: where oo represents the different directional normal stress or normal strain components and wu represents the different directional shear stress or shear strain components; is the average strain, is the average stress; wherein is the stress increment in the pq direction at the k+1 time point, pq taking the values 11, 22, 33, 12, 13, 23; is the strain increment in the PQ direction at the k+1 time point, PQ taking the values 11, 22, 33, 12, 13, 23; δ pq , δ PQ , δ pQ , and δ qP are Kronecker symbols; The material parameters are input in the user material subroutine UMAT, a finite element geometric model is created by using ABAQUS, three-dimensional solid elements are selected to discretize the finite element geometric model, and the discrete finite element geometric model is subjected to loads and boundary conditions; after setting the incremental step, the operation is submitted to obtain the mechanical response results of the finite element model.

9. The method of claim 1, wherein, The five-parameter nonlinear viscoelastic model is selected, and the sub-elements include spring element one, spring element two, linear dashpot element one, linear dashpot element two and variable cross-section dashpot element; The relationship between the total stress and the total strain of the model in the complex frequency domain is established as: wherein E1is the modulus value of spring element one, E2is the modulus value of spring element two, η1is the viscosity of linear dashpot element one, η2is the viscosity of linear dashpot element two, η3is the viscosity of the variable cross-section dashpot element, and s is the complex frequency domain argument. The initial value of zero is introduced into the time domain and the complex frequency domain to perform derivative equivalent substitution: where L is a Laplace transform, is an αth derivative operator with respect to time, α is the order of the derivative, f(t) is the function to be differentiated, and t is time. The one-dimensional form of the lowest order differential expression of the relationship between the total stress and the total strain of the model is constructed as: where ε(t) is the total strain and σ(t) is the total stress.

10. The method of claim 9, wherein, The Grünwald-Letnikov method shown in formula (44) is used to discretize the differential term in the one-dimensional form of the lowest order differential expression of the model: where Δt is the equidistance time interval, k is the number of equidistance dispersion intervals at time t; is the Grunwald coefficient, and The stress updating expression of the discrete form of the model is constructed as: wherein A, C i , D i , B are intermediate variables: If the constant Poisson's ratio is assumed, the three-dimensional form conversion relationship of the intermediate variables A, C i , D i , B is as follows: where I=1, 2; J=1, 2, 3; K is a subscript, indicating that the variable is related to the normal strain or the normal stress, G is a subscript, indicating that the variable is related to the deviatoric strain or the deviatoric stress; μ is a constant Poisson's ratio; If the very Poisson's ratio assumption is adopted, A in equation (47) is given by G , K , G , K , and are directly measured through triaxial tests; Based on the incremental iterative method, the stress updating expression(48) and the Jacobian matrix updating expression(49) of the three-dimensional form model are obtained: where oo represents the different directional normal stress or normal strain components and wu represents the different directional shear stress or shear strain components; is the average strain, is the average stress; wherein is the stress increment in the pq direction at the k+1 time point, pq taking the values 11, 22, 33, 12, 13, 23; is the strain increment in the PQ direction at the k+1 time point, PQ taking the values 11, 22, 33, 12, 13, 23; δ pq , δ PQ , δ pQ , and δ qP are Kronecker symbols; The material parameters are input in the user material subroutine UMAT, and the finite element geometric model is created by ABAQUS. The three-dimensional solid element is selected to discretize the finite element geometric model. The load and boundary conditions are applied to the discretized finite element geometric model. The operation is submitted after setting the incremental step, and the mechanical response results of the finite element model are obtained.

Images