Asphalt mixture linear and nonlinear viscoelastic model finite element numerical realization method

Abstract

The invention discloses an asphalt mixture linear and nonlinear viscoelastic model finite element numerical value implementation method, and belongs to the technical field of constitutive behaviors of asphalt mixtures. The method aims at solving the problem that no universal and feasible numerical calculation method exists for the linear and nonlinear viscoelastic models of the asphalt mixture at present. The method comprises the following steps: acquiring a complex frequency domain modulus expression of linear and nonlinear viscoelastic models of the asphalt mixture based on Laplace transformation; then, constructing a one-dimensional model stress-strain minimum-order differential expression through equivalent substitution of a time domain and a complex frequency domain; constructing a discrete form model stress updating expression based on differential terms in a GL discrete stress-strain differential form; deducing a three-dimensional model stress updating expression and a Jacobian matrix expression based on a normal or non-normal Poisson's ratio hypothesis and an incremental iteration method; based on a subprogram development interface in the ABAQUS, a user material subprogram is developed by utilizing FORTRAN to realize finite element numerical calculation. According to the method, finite element numerical calculation of the viscoelastic model of the asphalt mixture is realized.

Description

Technical Field

[0001] The invention relates to a finite element numerical realization method for a linear and nonlinear viscoelastic model of an asphalt mixture, and belongs to the technical field of constitutive behavior of asphalt mixtures. Background Art

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

[0003] At present, typical linear and nonlinear viscoelastic models of asphalt mixtures include integer-order viscoelastic models, fractional-order viscoelastic models, modified Burgers models, and five-parameter nonlinear viscoelastic models. Among them, integer-order viscoelastic models are linear viscoelastic models, including Kelvin model, Burgers model, three-parameter solid model, generalized Kelvin model, and generalized Maxwell model. The integer-order model is composed of Newtonian viscosity 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 therefore widely used. However, the integer-order model is difficult to characterize the dynamic mechanical behavior of asphalt mixtures in a wide frequency domain, and it is easily affected by noise during parameter identification, and the identified model parameters may be negative, thus losing physical meaning. The fractional-order viscoelastic model is composed of Abel viscosity pot and spring in series and parallel, which can characterize the mechanical behavior of asphalt mixtures in a wide frequency domain. However, the fractional-order viscoelastic model lacks simple expressions for creep and relaxation moduli, which limits its application in finite element analysis. The modified Burgers model is a type of nonlinear viscoelastic model that replaces the series viscosity pot with an exponential viscosity pot based on the Burgers model. Unlike the Burgers model, in which the creep deformation grows infinitely with time, the modified Burgers model can describe the consolidation effect of asphalt mixtures, and is therefore widely used in pavement rutting analysis. However, the modified Burgers model can only describe stable creep and deceleration creep, but cannot describe accelerated creep. The five-parameter nonlinear viscoelastic model is based on the Burgers model, with a variable cross-section viscosity pot in series. The stress of the viscosity pot can be expressed as the third-order derivative of strain with respect to time, which can describe the three stages of creep deceleration, creep stabilization, and creep acceleration.

[0004] For the above-mentioned different types of linear and nonlinear viscoelastic models, it is necessary to develop finite element numerical algorithms for specific models separately, and the calculation process is relatively complicated and the calculation cost is high. For integer-order viscoelastic models, the exponential algorithm is usually used for calculation. When the model contains nonlinear viscoelastic elements such as Abel fractional-order viscoelasticity, exponential viscoelasticity and variable-section viscoelasticity, it will be difficult to use the exponential algorithm for finite element calculation. Fractional-order viscoelastic models are usually calculated using the Mittag-Leffler function method, which first expresses 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 expressed 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 the approximate calculation formula of the Mittag-Leffler function when calculating, which has a large amount of calculation and the calculation accuracy is affected by the approximate calculation formula. The nonlinear viscoelastic model with exponential viscosity is usually calculated using the implicit Euler method, which uses the constant stress assumption to obtain the expression of the creep strain rate, and solves the creep strain increment through reciprocating iteration in each incremental step, and then solves the stress increment. Therefore, it is only applicable to the case of constant external load and has a large amount of calculation. For the five-parameter nonlinear viscoelastic model used to describe the three-stage creep curve, no research has proposed 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-section viscosity, there is still a lack of effective numerical methods to realize the finite element calculation of such models.

[0005] Therefore, it is urgent to propose a simple and feasible numerical calculation method to apply various linear and nonlinear viscoelastic models of asphalt mixtures to the finite element analysis of asphalt pavements. Summary of the invention

[0006] In view of the problem that there is currently no universal and feasible numerical calculation method for linear and nonlinear viscoelastic models of asphalt mixtures, the present invention provides a finite element numerical implementation method for linear and nonlinear viscoelastic models of asphalt mixtures.

[0007] The present invention provides a finite element numerical realization method for linear and nonlinear viscoelastic models of asphalt mixtures, comprising:

[0008] 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;

[0009] Introducing the equivalent substitution of time domain and complex frequency domain derivatives, constructing the one-dimensional lowest-order differential expression of the total stress and total strain relationship of the model;

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

[0011] Combining the constant Poisson's ratio assumption or the non-constant Poisson's ratio assumption with the incremental iteration method, the stress update expression and Jacobian matrix update expression of the three-dimensional model are constructed;

[0012] Based on the stress update expression and Jacobian matrix update expression of the three-dimensional model, the subroutine development interface in ABAQUS was utilized, and FORTRAN was used to develop the user material subroutine UMAT for linear and nonlinear viscoelastic models. A finite element geometric model was created in ABAQUS, and the material parameters, loads and boundary conditions of the linear and nonlinear viscoelastic models were input after meshing. The user material subroutine UMAT was called to realize the finite element numerical calculation of the linear and nonlinear viscoelastic models of asphalt mixtures.

[0013] The beneficial effects of the present invention are as follows: the method of the present invention first obtains the complex frequency domain modulus expression of the linear and nonlinear viscoelastic models of asphalt mixtures based on Laplace transform; then, through time domain / complex frequency domain equivalent substitution, the lowest-order differential expression of stress and strain of the linear / nonlinear viscoelastic model of asphalt mixture in one-dimensional form is constructed; on this basis, the differential terms in the differential form of stress and strain are discretized based on the Grünwald–Letnikov (GL) method, so as to construct the stress update expression of the linear / nonlinear viscoelastic model of asphalt mixture in discrete form; thereafter, based on the constant / inconstant Poisson's ratio assumption and the incremental iteration method, the stress update expression and Jacobian matrix expression of the linear / nonlinear viscoelastic model of asphalt mixture in three-dimensional form are derived; finally, based on the subroutine development interface in ABAQUS, a user material subroutine (Usermaterial subroutine, UMAT) is developed using FORTRAN to realize the finite element numerical calculation of the linear / nonlinear viscoelastic model of asphalt mixture.

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

[0015] (1) Integer-order viscoelastic models are usually calculated using exponential algorithms, which require that the creep or relaxation modulus of the viscoelastic model has a natural exponential form, and is therefore applicable to linear viscoelastic models containing only Newtonian viscosity. The method of the present invention uses a constitutive equation in differential form, which is not only applicable to linear viscoelastic models containing only Newtonian viscosity, but also to numerical calculations of nonlinear viscoelastic models composed of Newtonian viscosity and other nonlinear viscosity. Applicable nonlinear viscosity includes: Abel fractional viscosity, exponential viscosity and variable cross-section viscosity.

[0016] (2) Fractional viscoelastic models are usually calculated using the Mittag-Leffler function method, which first expresses the time-domain creep compliance or relaxation modulus of the fractional viscoelastic model using the Mittag-Leffler function, and then constructs a stress increment expression. This method is only applicable to simple fractional viscoelastic models whose time-domain creep compliance or relaxation modulus can be expressed using 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 requires the use of the approximate calculation formula of the Mittag-Leffler function during calculation, which results in a large amount of calculation and the calculation accuracy is affected by the approximate calculation formula. The method of the present invention uses a constitutive equation in differential form, which does not require the derivation of the time-domain expression of the model creep compliance or relaxation modulus, and can be applied to finite element numerical calculations of complex fractional viscoelastic models.

[0017] (3) The nonlinear viscoelastic model with exponential viscosity is usually calculated using the implicit Euler method, which uses the constant stress assumption to obtain the expression of the creep strain rate, and solves the creep strain increment through reciprocating iteration in each incremental step, and then solves the stress increment. Therefore, it is only applicable to the case where the external load is constant, and the amount of calculation is large. The method of the present invention does not need to solve the stress increment through reciprocating iteration in each incremental step, and the calculation efficiency is high. In addition, the method of the present invention does not use the assumption of constant external load during 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 with cross-sectional viscosity pots, such as the five-parameter nonlinear viscoelastic model used to describe the three-stage creep curve, the prior art has not yet proposed a finite element numerical calculation method for the model. The method of the present invention can realize the finite element numerical calculation of nonlinear viscoelastic models with cross-sectional viscosity pots.

[0019] (5) When the method of the present invention is used to solve the finite element numerical calculation of linear / nonlinear viscoelastic models containing Newtonian viscosity pot, Abel viscosity pot, exponential viscosity pot or variable cross-section viscosity pot, it has a universal solution step and calculation format, and there is no need to design a special solution algorithm for the viscosity pot type. Therefore, it is more convenient for finite element implementation than other algorithms. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 It is a flow chart of the finite element numerical realization method of the linear and nonlinear viscoelastic model of asphalt mixture according to the present invention;

[0021] Figure 2 It is a flow chart of UMAT development of linear and nonlinear viscoelastic models of asphalt mixtures;

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

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

[0024] Figure 5 It is a schematic diagram of uniaxial loading of the Burgers model;

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

[0026] Figure 7 is the creep curve of the modified fractional Zener model;

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

[0028] Fig. 9 is the modified Burgers model creep curve;

[0029] Fig.10 It is a schematic diagram of the five-parameter nonlinear viscoelastic model;

[0030] Fig.11 is the creep curve of the five-parameter nonlinear viscoelastic model. DETAILED DESCRIPTION

[0031] The following will be combined with the drawings in the embodiments of the present invention to clearly and completely describe the technical solutions in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative work are within the scope of protection of the present invention.

[0032] It should be noted that, in the absence of conflict, the embodiments of the present invention and the features in the embodiments may be combined with each other.

[0033] The present invention will be further described below in conjunction with the accompanying drawings, but is not intended to be a limitation of the present invention.

[0034] Combination Figure 1 and Figure 2 As shown, the present invention provides a finite element numerical realization method of 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] Introducing the equivalent substitution of time domain and complex frequency domain derivatives, constructing the one-dimensional lowest-order differential expression of the total stress and total strain relationship of the model;

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

[0038] Combining the constant Poisson's ratio assumption or the non-constant Poisson's ratio assumption with the incremental iteration method, the stress update expression and Jacobian matrix update expression of the three-dimensional model are constructed;

[0039] Based on the stress update expression and Jacobian matrix update expression of the three-dimensional model, the subroutine development interface in ABAQUS was utilized, and FORTRAN was used to develop the user material subroutine UMAT for linear and nonlinear viscoelastic models. A finite element geometric model was created in ABAQUS, and the material parameters, loads and boundary conditions of the linear and nonlinear viscoelastic models were input after meshing. The user material subroutine UMAT was called to realize the finite element numerical calculation of the linear and nonlinear viscoelastic models of asphalt mixtures.

[0040] This embodiment realizes finite element numerical calculation of an integer-order viscoelastic model with a Newtonian viscoelastic pot (the viscoelastic pot stress is the first-order derivative of the viscoelastic pot strain multiplied by the viscosity), a fractional-order viscoelastic model with an Abel viscoelastic pot (the viscoelastic pot stress is the fractional-order derivative of the viscoelastic pot strain multiplied by the viscosity), a nonlinear viscoelastic model with an exponential viscoelastic pot (the compliance of the viscoelastic pot is an exponential function), and a nonlinear viscoelastic model with a variable-section viscoelastic pot (the viscoelastic pot stress is the second-order or higher integer-order derivative of the viscoelastic pot strain multiplied by the viscosity).

[0041] The method of the present invention 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 1: Combination Figure 3 As shown, the linear viscoelastic model is selected as the Burgers model, and the sub-elements include a spring element 1, a spring element 2, a linear Newtonian viscoelastic element 1, and a linear Newtonian viscoelastic element 2;

[0043] First, based on the Laplace transform, the total stress and total strain relationship 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 update and Jacobian matrix update expressions of the model are derived based on the GL method, and the UMAT subroutine is developed. Finally, a finite element model is constructed based on ABAQUS, and numerical calculations of the model are 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 relationship between the total stress and total strain of the model in the complex frequency domain is established as follows:

[0045]

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

[0047] Simplify formula (1) and introduce the equivalent substitution of time domain and complex frequency domain derivatives with initial value of zero:

[0048]

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

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

[0051]

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

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

[0054]

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

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

[0057]

[0058] Where A and C i , D i , B is the intermediate variable:

[0059]

[0060] Using the constant Poisson's ratio assumption or the non-constant Poisson's ratio assumption, formula (6) is transformed into a three-dimensional form;

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

[0062]

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

[0064] If the extraordinary Poisson's ratio assumption is adopted, A in formula (7) G , A K , B G , B K , and Directly measured by triaxial test;

[0065] Based on the incremental iteration method in ABAQUS, the stress update expression (8) and Jacobian matrix update expression (9) of the three-dimensional model 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 mean stress;

[0068]

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

[0070] Input material parameters in the user material subroutine UMAT, use ABAQUS to create a finite element geometric model, 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 the incremental step and submit the operation to obtain the mechanical response results of the finite element model.

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

[0072] Table 1 Material parameters of linear / nonlinear viscoelastic model

[0073]

[0074] Created using ABAQUS Figure 5 The finite element model shown in the figure adopts the three-dimensional 8-node reduced integration unit C3D8R discretized finite element geometric model, fixes the normal displacement of the back node of the three-dimensional 8-node reduced integration unit in the x direction, and applies the stress load σ in the x direction shown in formula (10): 11 ; Set the incremental step to 0.1s and the analysis step to 100s. Finally, the creep curve of the model in the x direction is as follows: Figure 4 The theoretical solution of the creep curve can be obtained by combining the complex frequency domain expression of formula (1) with formula (10) through Laplace numerical inverse transformation.

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

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

[0077] Embodiment 2: Combination Figure 6 As shown, a modified fractional-order Zener model is selected, and the sub-elements include spring element 1, spring element 2, Abel viscosity pot element 1 and Abel viscosity pot element 2;

[0078] Based on the stress-strain relationship and Laplace transform of the sub-element of the modified fractional-order Zener model, 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 is the modulus in the complex frequency domain, E1 is the modulus value of spring element one, E2 is the modulus value of spring element two, η1 is the viscosity of Abel viscosity element one, η2 is the viscosity of Abel viscosity element two, s is the independent variable in the complex frequency domain, α1 is the fractional order of Abel viscosity element one, α2 is the fractional order of Abel viscosity element two, α2>α1;

[0081] The simplified total stress and total strain relationship of the modified fractional-order Zener model in the complex frequency domain is introduced by introducing the equivalent substitution of the derivatives in the time domain and complex frequency domain with an initial value of zero:

[0082]

[0083] Where L is the Laplace transform, is the α-order derivative operator with respect to time, α is the derivative order, f(t) is the derivative function to be determined, and t is time;

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

[0085]

[0086] Where ε(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 terms in the lowest-order differential expression in one-dimensional form:

[0088]

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

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

[0091]

[0092] Where A and C i , D i , B is the intermediate variable:

[0093]

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

[0095]

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

[0097] If the extraordinary Poisson's ratio assumption is adopted, A in formula (27) G , A K , B G , B K , and Directly measured by triaxial test;

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

[0099]

[0100] 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 mean stress;

[0101]

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

[0103] Input material parameters in the user material subroutine UMAT, use ABAQUS to create a finite element geometric model, 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 the incremental step and submit the calculation to obtain the mechanical response results of the finite element model.

[0104] As an example, based on the developed modified fractional Zener model, the UMAT is created using ABAQUS by inputting material parameters according to Table 1. Figure 5 The finite element model shown in the figure adopts the three-dimensional 8-node reduced integration unit C3D8R discretized finite element geometric model, fixes the normal displacement of the back node of the three-dimensional 8-node reduced integration unit in the x direction, and applies the stress load σ in the x direction shown in formula (10): 11 ; Set the incremental step to 0.1s and the analysis step to 100s. Finally, the creep curve of the model in the x direction is as follows: Figure 7 The theoretical solution of the creep curve can be obtained by combining the complex frequency domain expression of formula (21) with formula (10) through Laplace numerical inverse transformation.

[0105] Example 3: Combination Figure 8 As shown, the modified Burgers model is selected, and the sub-elements include spring element 1, spring element 2, linear viscosity element and exponential viscosity element;

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

[0107]

[0108] In the formula is the modulus in the complex frequency domain, E1 is the modulus value of the spring element 1, a is the parameter 1 of the exponential viscosity element, s is the independent variable in the complex frequency domain, b is the parameter 2 of the exponential viscosity element, η2 is the viscosity of the linear viscosity element, and E2 is the modulus value of the spring element 2;

[0109] The simplified total stress and total strain relationship of the modified Burgers model in the complex frequency domain is introduced by introducing the equivalent substitution of the derivatives in the time domain and complex frequency domain with an initial value of zero:

[0110]

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

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

[0113]

[0114] Where ε(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 terms in the lowest-order differential expression in one-dimensional form:

[0116]

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

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

[0119]

[0120] Where A and C i , D i , B is the intermediate variable:

[0121]

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

[0123]

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

[0125] If the extraordinary Poisson's ratio assumption is adopted, A in formula (37) G , A K , B G , B K , and Directly measured by triaxial test;

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

[0127]

[0128] 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 mean stress;

[0129]

[0130] In the formula is the stress increment in the pq direction at the k+1th time point, and the values ​​of pq are 11, 22, 33, 12, 13, and 23; is the strain increment in the PQ direction at the k+1th time point, and the PQ values ​​are 11, 22, 33, 12, 13, and 23; δ pq , δ PQ , δ pQ and δ qP is the Kronecker symbol;

[0131] Input material parameters in the user material subroutine UMAT, use ABAQUS to create a finite element geometric model, 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 the incremental step and submit the calculation to obtain the mechanical response results of the finite element model.

[0132] As an example, the UMAT based on the developed modified Burgers model is created using ABAQUS by inputting material parameters according to Table 1. Figure 5 The finite element model shown in the figure adopts the three-dimensional 8-node reduced integration unit C3D8R discretized finite element geometric model, fixes the normal displacement of the back node of the three-dimensional 8-node reduced integration unit in the x direction, and applies the stress load σ in the x direction shown in formula (10): 11; Set the incremental step to 0.1s and the analysis step to 100s. Finally, the creep curve of the model in the x direction is as follows: Fig. 9 The theoretical solution of the creep curve can be obtained by combining the complex frequency domain expression of formula (31) with formula (10) through Laplace numerical inverse transformation.

[0133] Embodiment 4: Combination Fig.10 As shown, a five-parameter nonlinear viscoelastic model is selected, and the sub-elements include spring element 1, spring element 2, linear viscoelastic element 1, linear viscoelastic element 2 and variable cross-section viscoelastic element;

[0134] Based on the stress-strain relationship and Laplace transform of the sub-components of the five-parameter nonlinear viscoelastic model, 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 is the modulus in the complex frequency domain, E1 is the modulus value of the spring element 1, E2 is the modulus value of the spring element 2, η1 is the viscosity of the linear viscosity element 1, η2 is the viscosity of the linear viscosity element 2, η3 is the viscosity of the variable cross-section viscosity element, and s is the independent variable in the complex frequency domain;

[0137] The simplified total stress and total strain relationship of the five-parameter nonlinear viscoelastic model in the complex frequency domain is introduced into the time domain and complex frequency domain derivatives with an initial value of zero:

[0138]

[0139] Where L is the Laplace transform, is the α-order derivative operator with respect to time, α is the derivative order, f(t) is the derivative function to be determined, and t is time;

[0140] The lowest-order differential expression for constructing the one-dimensional form of the total stress and total strain relationship of the model is:

[0141]

[0142] Where ε(t) is the total strain and σ(t) is the total stress.

[0143] The Grünwald-Letnikov method shown in formula (44) is used to discretize the differential terms in the lowest-order differential expression in one-dimensional form:

[0144]

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

[0146]

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

[0148]

[0149] Where A and C i , D i , B is the intermediate variable:

[0150]

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

[0152]

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

[0154] If the extraordinary Poisson's ratio assumption is adopted, A in formula (47) G , A K , B G , B K , and Directly measured by triaxial test;

[0155] Based on the incremental iteration method in ABAQUS, the stress update expression (48) and Jacobian matrix update expression (49) of the three-dimensional 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 mean stress;

[0158]

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

[0160] Input material parameters in the user material subroutine UMAT, use ABAQUS to create a finite element geometric model, 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 the incremental step and submit the calculation to obtain the mechanical response results of the finite element model.

[0161] As an example, the UMAT based on the developed modified Burgers model is created using ABAQUS by inputting material parameters according to Table 1. Figure 5 The finite element model shown in the figure adopts the three-dimensional 8-node reduced integration unit C3D8R discretized finite element geometric model, fixes the normal displacement of the back node of the three-dimensional 8-node reduced integration unit in the x direction, and applies the stress load σ in the x direction shown in formula (10): 11 ; Set the incremental step to 0.1s and the analysis step to 100s. Finally, the creep curve of the model in the x direction is as follows: Fig.11 The theoretical solution of the creep curve can be obtained by combining the complex frequency domain expression of formula (41) with formula (10) through Laplace numerical inverse transformation.

[0162] Although the present invention is described herein with reference to specific embodiments, it should be understood that these embodiments are merely examples of the principles and applications of the present invention. It should therefore be understood that many modifications may be made to the exemplary embodiments and that other arrangements may be devised without departing from the spirit and scope of the present invention as defined by the appended claims. It should be understood that the various dependent claims and features described herein may be combined in a manner different from that described in the original claims. It should also be understood that the features described in conjunction with a single embodiment may be used in other described embodiments.

Claims

1. A finite element numerical implementation method for linear and nonlinear viscoelastic models of asphalt mixtures, characterized in that: include, 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; Introducing the equivalent substitution of time domain and complex frequency domain derivatives, constructing the one-dimensional lowest-order differential expression of the total stress and total strain relationship of the model; The Grünwald-Letnikov method is used to discretize the differential terms in the lowest-order differential expression of the one-dimensional form, and the discrete form model stress update expression is constructed; Combining the constant Poisson's ratio assumption or the non-constant Poisson's ratio assumption with the incremental iteration method, the stress update expression and Jacobian matrix update expression of the three-dimensional model are constructed; Based on the stress update expression and Jacobian matrix update expression of the three-dimensional model, the subroutine development interface in ABAQUS was utilized, and FORTRAN was used to develop the user material subroutine UMAT for linear and nonlinear viscoelastic models. A finite element geometric model was created in ABAQUS, and the material parameters, loads and boundary conditions of the linear and nonlinear viscoelastic models were input after meshing. The user material subroutine UMAT was called to realize the finite element numerical calculation of the linear and nonlinear viscoelastic models of asphalt mixtures.

2. The finite element numerical implementation method of linear and nonlinear viscoelastic models of asphalt mixture according to claim 1 is characterized in that: The linear viscoelastic model is selected as the Burgers model, wherein the sub-elements include a spring element 1, a spring element 2, a linear Newtonian viscoelastic element 1, and a linear Newtonian viscoelastic element 2; The relationship between the total stress and total strain of the model in the complex frequency domain is established as: In the formula is the complex frequency domain modulus, E1 is the modulus value of spring element 1, E2 is the modulus value of spring element 2, η1 is the viscosity of linear Newtonian viscosity element 1, η2 is the viscosity of linear Newtonian viscosity element 2, and s is the complex frequency domain independent variable; Introduce the equivalent substitution of time domain and complex frequency domain derivatives with initial value of zero: Where L is the Laplace transform, is the α-order derivative operator with respect to time, α is the derivative order, f(t) is the derivative function to be determined, and t is time; The lowest-order differential expression for constructing the one-dimensional form of the total stress and total strain relationship of the model is: Where ε(t) is the total strain and σ(t) is the total stress.

3. The finite element numerical implementation method of the linear and nonlinear viscoelastic model of asphalt mixture according to claim 2 is characterized in that: The Grünwald-Letnikov method shown in formula (4) is used to discretize the differential terms in the lowest-order differential expression in one-dimensional form: Where Δt is the equidistant time interval, k is the number of equidistant scattered intervals of time t; is the Grünwald coefficient, and The stress update expression of the discrete form model is constructed as follows: Where A and C i , D i , B is the intermediate variable:

4. The finite element numerical implementation method of the linear and nonlinear viscoelastic model of asphalt mixture according to claim 3 is characterized in that: Using the constant Poisson's ratio assumption or the non-constant Poisson's ratio assumption, formula (6) is transformed into a three-dimensional form; If the constant Poisson ratio assumption is adopted, the intermediate variables A and C i , D i , the three-dimensional form conversion relationship of B is: Where I = 1, 2; J = 1, 2; K as a subscript indicates that the variable is related to normal strain or normal stress, G as a subscript indicates that the variable is related to deviatoric strain or deviatoric stress; μ is the constant Poisson's ratio; If the extraordinary Poisson's ratio assumption is adopted, A in formula (7) G , A K , B G , B K , and Directly measured by triaxial test; Based on the incremental iteration method, the stress update expression (8) and Jacobian matrix update expression (9) of the three-dimensional model are obtained: 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 mean stress; In the formula is the stress increment in the pq direction at the k+1th time point, and the values ​​of pq are 11, 22, 33, 12, 13, and 23; is the strain increment in the PQ direction at the k+1th time point, and the PQ values ​​are 11, 22, 33, 12, 13, and 23; δ pq , δ PQ , δ pQ and δ qP is the Kronecker symbol; Input material parameters in the user material subroutine UMAT, use ABAQUS to create a finite element geometric model, 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 the incremental step and submit the calculation to obtain the mechanical response results of the finite element model.

5. The finite element numerical implementation method of linear and nonlinear viscoelastic model of asphalt mixture according to claim 1 is characterized in that: Selecting a modified fractional order Zener model, wherein the sub-elements include a spring element 1, a spring element 2, an Abel viscosity pot element 1 and an Abel viscosity pot element 2; The relationship between the total stress and total strain of the model in the complex frequency domain is established as: In the formula is the modulus in the complex frequency domain, E1 is the modulus value of spring element one, E2 is the modulus value of spring element two, η1 is the viscosity of Abel viscosity element one, η2 is the viscosity of Abel viscosity element two, s is the independent variable in the complex frequency domain, α1 is the fractional order of Abel viscosity element one, α2 is the fractional order of Abel viscosity element two, α2>α1; Introduce the equivalent substitution of time domain and complex frequency domain derivatives with initial value of zero: Where L is the Laplace transform, is the α-order derivative operator with respect to time, α is the derivative order, f(t) is the derivative function to be determined, and t is time; The lowest-order differential expression for constructing the one-dimensional form of the total stress and total strain relationship of the model is: Where ε(t) is the total strain and σ(t) is the total stress.

6. The finite element numerical implementation method of linear and nonlinear viscoelastic models of asphalt mixture according to claim 5 is characterized in that: The Grünwald-Letnikov method shown in formula (24) is used to discretize the differential terms in the lowest-order differential expression in one-dimensional form: Where Δt is the equidistant time interval, k is the number of equidistant scattered intervals of time t; is the Grünwald coefficient, and The stress update expression of the discrete form model is constructed as follows: Where A and C i , D i , B is the intermediate variable: If the constant Poisson ratio assumption is adopted, the intermediate variables A and C i , D i , the three-dimensional form conversion relationship of B is: Where I = 1, 2; J = 1, 2; K as a subscript indicates that the variable is related to normal strain or normal stress, G as a subscript indicates that the variable is related to deviatoric strain or deviatoric stress; μ is the constant Poisson's ratio; If the extraordinary Poisson's ratio assumption is adopted, A in formula (27) G , A K , B G , B K , and Directly measured by triaxial test; Based on the incremental iteration method, the stress update expression (28) and Jacobian matrix update expression (29) of the three-dimensional model are obtained: 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 mean stress; In the formula is the stress increment in the pq direction at the k+1th time point, and the values ​​of pq are 11, 22, 33, 12, 13, and 23; is the strain increment in the PQ direction at the k+1th time point, and the PQ values ​​are 11, 22, 33, 12, 13, and 23; δ pq , δ PQ , δ pQ and δ qP is the Kronecker symbol; Input material parameters in the user material subroutine UMAT, use ABAQUS to create a finite element geometric model, 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 the incremental step and submit the operation to obtain the mechanical response results of the finite element model.

7. The finite element numerical implementation method of linear and nonlinear viscoelastic models of asphalt mixture according to claim 1 is characterized in that: Selecting a modified Burgers model, wherein the sub-elements include a spring element 1, a spring element 2, a linear viscosity element, and an exponential viscosity element; The relationship between the total stress and total strain of the model in the complex frequency domain is: In the formula is the modulus in the complex frequency domain, E1 is the modulus value of the spring element 1, a is the parameter 1 of the exponential viscosity element, s is the independent variable in the complex frequency domain, b is the parameter 2 of the exponential viscosity element, η2 is the viscosity of the linear viscosity element, and E2 is the modulus value of the spring element 2; Introduce the equivalent substitution of time domain and complex frequency domain derivatives with initial value of zero: Where L is the Laplace transform, is the α-order derivative operator with respect to time, α is the derivative order, f(t) is the derivative function to be determined, and t is time; The lowest-order differential expression for constructing the one-dimensional form of the total stress and total strain relationship of the model is: Where ε(t) is the total strain and σ(t) is the total stress.

8. The finite element numerical implementation method of linear and nonlinear viscoelastic models of asphalt mixture according to claim 7 is characterized in that: The Grünwald-Letnikov method shown in formula (34) is used to discretize the differential terms in the lowest-order differential expression in one-dimensional form: Where Δt is the equidistant time interval, k is the number of equidistant scattered intervals of time t; is the Grünwald coefficient, and The stress update expression of the discrete form model is constructed as follows: Where A and C i , D i , B is the intermediate variable: If the constant Poisson ratio assumption is adopted, the intermediate variables A and C i , D i , the three-dimensional form conversion relationship of B is: Where I = 1, 2; K as a subscript indicates that the variable is related to normal strain or normal stress, G as a subscript indicates that the variable is related to deviatoric strain or deviatoric stress; μ is the constant Poisson's ratio; If the extraordinary Poisson's ratio assumption is adopted, A in formula (37) G , A K , B G , B K , and Directly measured by triaxial test; Based on the incremental iteration method, the stress update expression (38) and Jacobian matrix update expression (39) of the three-dimensional model are obtained: 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 mean stress; In the formula is the stress increment in the pq direction at the k+1th time point, and the values ​​of pq are 11, 22, 33, 12, 13, and 23; is the strain increment in the PQ direction at the k+1th time point, and the PQ values ​​are 11, 22, 33, 12, 13, and 23; δ pq , δ PQ , δ pQ and δ qP is the Kronecker symbol; Input material parameters in the user material subroutine UMAT, use ABAQUS to create a finite element geometric model, 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 the incremental step and submit the operation to obtain the mechanical response results of the finite element model.

9. The finite element numerical implementation method of linear and nonlinear viscoelastic models of asphalt mixture according to claim 1 is characterized in that: Select a five-parameter nonlinear viscoelastic model, wherein the sub-elements include a spring element 1, a spring element 2, a linear viscoelastic element 1, a linear viscoelastic element 2, and a variable cross-section viscoelastic element; The relationship between the total stress and total strain of the model in the complex frequency domain is established as: In the formula is the modulus in the complex frequency domain, E1 is the modulus value of the spring element 1, E2 is the modulus value of the spring element 2, η1 is the viscosity of the linear viscosity element 1, η2 is the viscosity of the linear viscosity element 2, η3 is the viscosity of the variable cross-section viscosity element, and s is the independent variable in the complex frequency domain; Introduce the equivalent substitution of time domain and complex frequency domain derivatives with initial value of zero: Where L is the Laplace transform, is the α-order derivative operator with respect to time, α is the derivative order, f(t) is the derivative function to be determined, and t is time; The lowest-order differential expression for constructing the one-dimensional form of the total stress and total strain relationship of the model is: Where ε(t) is the total strain and σ(t) is the total stress.

10. The finite element numerical implementation method of linear and nonlinear viscoelastic models of asphalt mixture according to claim 9, characterized in that: The Grünwald-Letnikov method shown in formula (44) is used to discretize the differential terms in the lowest-order differential expression in one-dimensional form: Where Δt is the equidistant time interval, k is the number of equidistant scattered intervals of time t; is the Grünwald coefficient, and The stress update expression of the discrete form model is constructed as follows: Where A and C i , D i , B is the intermediate variable: If the constant Poisson ratio assumption is adopted, the intermediate variables A and C i , D i , the three-dimensional form conversion relationship of B is: Where I = 1, 2; J = 1, 2, 3; K as a subscript indicates that the variable is related to normal strain or normal stress, G as a subscript indicates that the variable is related to deviatoric strain or deviatoric stress; μ is the constant Poisson's ratio; If the extraordinary Poisson's ratio assumption is adopted, A in formula (47) G , A K , B G , B K , and Directly measured by triaxial test; Based on the incremental iteration method, the stress update expression (48) and Jacobian matrix update expression (49) of the three-dimensional model are obtained: 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 mean stress; In the formula is the stress increment in the pq direction at the k+1th time point, and the values ​​of pq are 11, 22, 33, 12, 13, and 23; is the strain increment in the PQ direction at the k+1th time point, and the PQ values ​​are 11, 22, 33, 12, 13, and 23; δ pq , δ PQ , δ pQ and δ qP is the Kronecker symbol. Input material parameters in the user material subroutine UMAT, use ABAQUS to create a finite element geometric model, 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 the incremental step and submit the operation to obtain the mechanical response results of the finite element model.

Images