CN109241636B - Finite element-based road surface structure multi-physical field coupling numerical simulation method - Google Patents

Abstract

The invention discloses a finite element-based road surface structure multi-physical field coupling numerical simulation method, which comprises the following steps of: defining mechanical, thermodynamic and hydraulic parameters; defining vehicle load and meteorological functions; adding a solid mechanical module, applying load and setting boundary conditions; adding a porous medium heat transfer module, referring to a meteorological function and setting boundary conditions; adding a Darcy law module, quoting a meteorological function and setting boundary conditions; adding a multi-physical field coupling module to couple the multi-physical field; and calculating and performing post-processing analysis. The invention couples the stress field, the temperature field and the hydraulic field, and comprehensively studies the influence of a plurality of physical fields on the road performance. The method has good guiding significance for the selection and analysis of the road surface structure.

Description

Multi-physics coupling numerical simulation method of pavement structure based on finite element

technical field

The invention relates to the numerical simulation of pavement structure, more specifically, a multi-physics coupling numerical simulation method of pavement structure based on finite element.

Background technique

For a long time, the road performance of asphalt pavement has been one of the research hotspots in the field of road engineering. Researchers can easily analyze the pavement structure with the help of finite element software, thereby saving manpower and material resources.

The change of road performance of pavement structure is often the result of the joint action of multiple physical fields, including stress field, temperature field and hydraulic field. analysis, there must be deviations in the analysis results.

Contents of the invention

In order to avoid the shortcomings of the above-mentioned prior art, the present invention provides a finite element-based multi-physics coupling numerical simulation method for pavement structure, using meteorological data to establish the transient temperature field and hydraulic field of pavement structure, and applying vehicle moving load, Realize the multi-physics coupling numerical simulation of pavement structure, in order to analyze the road performance of pavement structure more accurately, and provide guidance for the selection of pavement materials and structures.

The present invention adopts following technical scheme for solving technical problems:

The present invention is based on the characteristic of multi-physics coupling numerical simulation method of pavement structure in multi-physics coupling finite element software according to the following steps:

Step 1: Define the mechanical parameters, thermodynamic parameters and hydraulic parameters of road materials;

The mechanical parameters include Young's modulus, Poisson's ratio and density, and also include viscoelastic parameters for viscoelastic materials;

Described thermodynamic parameter comprises thermal conductivity, specific heat capacity and thermal expansion coefficient;

The hydraulic parameters include hydraulic conductivity, porosity and Biot elastic parameters;

Step 2: Define the vehicle load function, and define the weather function based on local meteorological data;

Described vehicle load function F (t) is characterized by formula (1), and described vehicle load comprises static load and dynamic load, and described dynamic load is half-sine load:

In the formula (1), A is the load amplitude, t is the time, _t0 is the duration of the load in a single cycle, T is the time interval between two adjacent loads, k=0,1,2...;

The meteorological functions include the daily variation function q(t) of solar radiation represented by formula (2), the daily variation function T _a of atmospheric temperature represented by formula (3), and the daily variation function of rainfall set according to meteorological data ;

in:

为日平均气温，/>

T_m为日气温变化幅度，

为日最高气温，/>

为日最低气温；T₀为初相位；q ₀ is the maximum daily radiation intensity, q ₀ =0.131mQ, m=12/c; Q is the total daily solar radiation; c is the actual sunshine time; ω is the angular frequency;

is the daily average temperature, />

T _m is the daily temperature variation range,

is the daily maximum temperature, />

is the daily minimum temperature; T ₀ is the initial phase;

Step 3: Add a solid mechanics module, specify the material properties of each structural layer for the solid mechanics module, and add boundary conditions;

The governing equation of described solid mechanics module is represented by formula (4):

in:

为位移向量的二阶时间导数；S is the total stress, F is the external force vector, ρ _s is the density of solid material, u is the displacement vector, ▽·S is the divergence of S,

is the second-order time derivative of the displacement vector;

The specified material properties of each structural layer include: specifying the road material of each structural layer as an elastic or viscoelastic material, setting the thermal expansion properties of the road material of each structural layer, and the thermal expansion properties are thermal expansion coefficient and strain reference temperature;

The boundary conditions include: a fixed constraint that the displacement in each direction is 0, a specified displacement that makes the displacement in a certain direction 0, and a boundary load;

Step 4: Add a porous media heat transfer module, specify material properties for the porous media heat transfer module, and define the heat flux of the road surface;

The governing equation of the porous media heat transfer module is characterized by formula (5):

为温度一阶时间导数，u_l为流体速度场，▽T_tem为温度梯度，q为传导热通量，q＝-k_eff▽T_tem，k_eff为有效导热系数，▽·q为q的散度，Q_h为热源；Wherein: ρ _l is fluid density, C _p is heat capacity of fluid at constant pressure, (ρ _l C _p ) _eff is heat capacity of effective volume at constant pressure, T _tem is temperature,

is the first-order time derivative of temperature, u _l is the fluid velocity field, ▽T _tem is the temperature gradient, q is the conduction heat flux, q=-k _eff ▽T _tem , k _eff is the effective thermal conductivity, ▽·q is the Divergence, Q _h is the heat source;

The specified material properties refer to giving the thermodynamic parameters defined in step 1 to the material;

The definition of the heat flux of the path table refers to adding the heat flux interface in the porous media heat transfer module, and inputting the solar radiation diurnal variation function q(t) and the atmospheric temperature diurnal variation function T _a defined in step 2, Among them, the heat flux type of the solar radiation diurnal variation function q _{(t) is generalized inward heat flux, and the heat flux type of the atmospheric temperature diurnal variation function T a is} convective heat flux. The convective heat flux defines the heat transfer coefficient h _c : h _c =3.7v _w +9.4, v _w is the wind speed;

Step 5: Add a Darcy's law module, specify material properties for the Darcy's law module, and add boundary conditions;

The governing equation of described Darcy's law module is represented by formula (6):

为p_l的一阶时间导数，▽·ρ_l为ρ_l的散度，K为水力传导率，g为重力加速度，▽p_l为流体压力梯度，▽D重力方向上的单位向量，Q_m为质量源项；In formula (6), S _p is the storage coefficient, p _l is the pore water pressure,

is the first-order time derivative of p _l , ▽·ρ _l is the divergence of ρ _l , K is the hydraulic conductivity, g is the acceleration of gravity, ▽p _l is the fluid pressure gradient, ▽D is the unit vector in the gravity direction, Q _m is the quality source item;

The specified material properties refer to the hydraulic conductivity and porosity in the hydraulic parameters defined in step 1 are given to the material;

The boundary conditions include pore water pressure and rainwater normal inflow velocity on the boundary;

Step 6: Add the multiphysics coupling module and enter the Biot elastic parameters in it;

The governing equation of the multiphysics coupling module is represented by formula (7):

为ε_vol的一阶时间导数，▽·(S-α_Bp_lI)为S-α_Bp_lI的散度，▽·(ρ_lu_l)为ρ_lu_l的散度；Among them: α _B is the Biot-Willis coefficient, I is the identity matrix, ε _vol is the volume strain,

is the first-order time derivative of ε _vol , ▽·(S-α _B p _l I) is the divergence of S-α _B p _l I, and ▽·(ρ _l u _l ) is the divergence of ρ _l u _l ;

Step 7: Calculate and perform post-processing analysis

The calculation refers to dividing the finite element grid for the constructed pavement structure, and using the finite element method to couple and solve the control equations of each module;

The post-processing analysis includes carrying out stress field analysis, displacement field analysis, temperature field analysis and pore water pressure field analysis on the pavement structure, and obtaining stress cloud maps, displacement cloud maps, temperature cloud maps, and pore water pressure cloud maps, etc., thereby completing multiple pavement structure analysis. Physics coupling analysis.

Compared with the prior art, the beneficial effects of the present invention are reflected in:

1. The present invention couples the pavement structure with the stress field, temperature field and hydraulic field to realize multi-field coupling numerical simulation of the pavement structure, thereby accurately analyzing the pavement performance of the pavement structure and providing good guidance for the selection of pavement materials and structures significance;

2. The present invention considers the influence of the temperature field on the pavement structure and introduces temperature stress, which has good guiding significance for the research on the anti-rutting performance of the pavement structure;

3. The present invention considers the influence of the hydraulic field on the pavement structure, so that the influence of the pore water pressure on the peeling degree of the asphalt can be analyzed;

4. The method for analyzing the pavement structure based on finite elements provided by the present invention enriches the application of finite elements in the field of road engineering.

Description of drawings

Fig. 1 is a flow chart of the method of the present invention;

Fig. 2 is the parametric pavement structure model in the embodiment of the present invention;

Fig. 3 is the change diagram of temperature with time at the middle point of the load region of the calculation results in the embodiment of the present invention;

Fig. 4 is the change diagram of the vertical displacement with time at the middle point of the load region of the calculation results in the embodiment of the present invention;

Fig. 5 is a diagram showing the variation of pore water pressure with time at the middle point of the load region in the calculation results in the embodiment of the present invention.

Labels in the figure: 1 surface layer, 2 base layer, 3 cushion layer, 4 soil foundation;

Detailed ways

Referring to Fig. 1, the finite element-based multi-physics coupling numerical simulation method of pavement structure in this embodiment is carried out in the multi-physics coupling finite element software as follows:

Step 1: Define the mechanical parameters, thermodynamic parameters, and hydraulic parameters of the road material: mechanical parameters include Young's modulus, Poisson's ratio, and density, and viscoelastic parameters for viscoelastic materials; thermodynamic parameters include thermal conductivity, specific heat capacity, and thermal expansion Coefficients; hydraulic parameters include hydraulic conductivity, porosity and Biot elastic parameters. In addition, a variable α _T needs to be defined. α _T is the displacement factor of the time-temperature equivalent principle, which is used to describe the relationship between relaxation time and temperature. It can be defined in the form of WLF equation: log ₁₀ α _T ＝c ₁ T ² +c ₂ T+c ₃ , where c ₁ , c ₂ and c ₃ are constants, and T is the surface layer temperature. The parameters and variables are listed in Table 2.

Step 2: Define the vehicle load function, and define the weather function based on local weather data:

The vehicle load function F(t) is represented by formula (1). The vehicle load includes static load and dynamic load, and the dynamic load can be half-sine load:

In the formula (1), A is the load amplitude, t is the time, _t0 is the duration of the load in a single cycle, T is the time interval between two adjacent loads, k=0,1,2..., k is the non-load integer.

Meteorological functions include the daily variation function q(t) of solar radiation represented by formula (2), the daily variation function T _a of atmospheric temperature represented by formula (3), and the daily variation function of rainfall set according to meteorological data;

in:

为日平均气温，/>

T_m为日气温变化幅度，/>

为日最高气温，/>

为日最低气温；T₀为初相位。q ₀ is the maximum daily radiation intensity, q ₀ =0.131mQ, m=12/c; Q is the total daily solar radiation; c is the actual sunshine time; ω is the angular frequency;

is the daily average temperature, />

T _m is the daily temperature variation range, />

is the daily maximum temperature, />

is the daily minimum temperature; T ₀ is the initial phase.

Due to the limited ability of the road surface to absorb solar radiation, the solar radiation diurnal variation function q(t) needs to be multiplied by the reduction factor, which is generally taken as 0.85; the initial phase T ₀ in formula (3) can be taken as 9, so This makes the formula (3) be calculated from midnight in the morning, which is consistent with the time starting point of a day, so as to be easy to understand.

Step 3: Add a solid mechanics module, specify the material properties of each structural layer for the solid mechanics module, and add boundary conditions:

The governing equation of the solid mechanics module is represented by formula (4):

in:

为位移向量的二阶时间导数。S is the total stress, F is the external force vector, ρ _s is the density of solid material, u is the displacement vector, ▽·S is the divergence of S,

is the second order time derivative of the displacement vector.

Designating the material properties of each structural layer includes: creating a new "Viscoelasticity" interface in the surface layer, selecting "Generalized Maxwell Model", inputting Young's modulus and relaxation time, that is, specifying the surface layer as a viscoelastic material, and the other layers are It is an elastic material; create a new "thermal expansion" interface in each layer, input the temperature, thermal expansion coefficient and strain reference temperature, and generally take the initial temperature of the pavement structure as the strain reference temperature.

Boundary conditions include: set the bottom surface as a fixed constraint that the displacement in all directions is 0, set the right side as a specified displacement in the x direction, set the front and rear sides as a specified displacement in the y direction, and Boundary loads are set in the loading area, and symmetrical boundary conditions can also be added on the left side to simplify the calculation.

Step 4: Add the heat transfer in porous media module, specify the material properties for the heat transfer in porous media module, and define the heat flux of the path table:

The governing equation of the porous media heat transfer module is represented by equation (5):

为温度一阶时间导数，u_l为流体速度场，▽T_tem为温度梯度，q为传导热通量，q＝-k_eff▽T_tem，k_eff为有效导热系数，▽·q为q的散度，Q_h为热源。Wherein: ρ _l is fluid density, C _p is heat capacity of fluid at constant pressure, (ρ _l C _p ) _eff is heat capacity of effective volume at constant pressure, T _tem is temperature,

is the first-order time derivative of temperature, u _l is the fluid velocity field, ▽T _tem is the temperature gradient, q is the conduction heat flux, q=-k _eff ▽T _tem , k _eff is the effective thermal conductivity, ▽·q is the Divergence, Q _h is the heat source.

Specifying material properties means assigning the thermodynamic parameters defined in step 1 to the material.

To define the heat flux of the path table means to add the heat flux interface in the porous media heat transfer module, and input the solar radiation diurnal variation function q(t) and atmospheric temperature diurnal variation function T _a defined in step 2, where the sun The heat flux type of the radiation daily variation function q(t) is generalized inward heat flux, and the heat flux type of the atmospheric temperature daily variation function T _a is convective heat flux. For the convective heat flux to which the atmospheric temperature daily variation function T _a belongs The flux needs to define the heat transfer coefficient h _c : h _c =3.7v _w +9.4, where v _w is the wind speed.

The generalized inward heat flux to which the solar radiation diurnal variation function q(t) belongs and the convective heat flux to which the atmospheric temperature diurnal variation function T _a belongs must be defined on the top surface of the pavement structure. To simplify the calculation, it can also be defined on the left side Add a symmetric boundary condition to the face.

Step 5: Add a Darcy's Law module, specify material properties for said Darcy's Law module, and add boundary conditions:

The governing equation of Darcy's law module is represented by equation (6):

为p_l的一阶时间导数，▽·ρ_l为ρ_l的散度，K为水力传导率，g为重力加速度，▽p_l为流体压力梯度，▽D重力方向上的单位向量，Q_m为质量源项。In formula (6), S _p is the storage coefficient, p _l is the pore water pressure,

is the first-order time derivative of p _l , ▽·ρ _l is the divergence of ρ _l , K is the hydraulic conductivity, g is the acceleration of gravity, ▽p _l is the fluid pressure gradient, ▽D is the unit vector in the gravity direction, Q _m is the quality source item.

Specifying material properties refers to assigning the hydraulic conductivity and porosity in the hydraulic parameters defined in step 1 to the material.

The boundary conditions include the pore water pressure and the normal inflow velocity of rainwater on the boundary. Among them, the pore water pressure is defined on the right side, and the desirable function p _l =-15×(Z+2), Z is the depth of the structure, that is, the pore water pressure changes linearly with the depth; the normal inflow velocity of rainwater is the The defined rainfall diurnal variation function needs to be defined on the top surface.

Step 6: Add the multiphysics coupling module and enter the Biot elasticity parameters in it:

The governing equation of the multiphysics coupling module is represented by equation (7):

为ε_vol的一阶时间导数，▽·(S-α_Bp_lI)为S-α_Bp_lI的散度，▽·(ρ_lu_l)为ρ_lu_l的散度。Among them: α _B is the Biot-Willis coefficient, I is the identity matrix, ε _vol is the volume strain,

is the first-order time derivative of ε _vol , ▽·(S-α _B p _l I) is the divergence of S-α _B p _l I, and ▽·(ρ _l u _l ) is the divergence of ρ _l u _l .

Step 7: Calculate and perform post-processing analysis

Calculation refers to dividing the finite element grid for the constructed pavement structure, and using the finite element method to couple and solve the control equations of each module; post-processing analysis includes stress field analysis, displacement field analysis, temperature field analysis and The analysis of the pore water pressure field obtains the stress cloud map, displacement cloud map, temperature cloud map and pore water pressure cloud map, etc., so as to complete the multi-physics field coupling analysis of the pavement structure.

Simulation process:

Step 1: Set the geometric dimensions of the pavement structure shown in Table 1:

Table 1 Geometric Dimensions (m)

The geometric model of the pavement structure is shown in Figure 2. From top to bottom, they are: surface layer 1, base layer 2, cushion layer 3, and soil foundation 4.

Step 2: Define the mechanical parameters, thermodynamic parameters and hydraulic parameters of each layer as shown in Table 2:

Table 2 Mechanical, thermodynamic and hydraulic parameters

(Continued from Table 2)

According to the method of the present invention, the vehicle load function is defined in the multi-physics coupling finite element software, the meteorological function is defined based on the local meteorological data, the solid mechanics module is added, the load is applied and the boundary conditions are set, the porous medium heat transfer module is added, and the Darcy's law module is added. Add the Darcy's Law module; add the multi-physics field coupling module to couple the multi-physics fields, and through calculation and post-processing analysis, the calculation results in this embodiment shown in Figure 3 change with time at the midpoint of the load region. Fig. 4 shows the variation diagram of the vertical displacement with time at the midpoint of the load region in the calculation results of this embodiment, and the pore water at the midpoint of the load region of the calculation results in the present embodiment as shown in Fig. 5 A graph of pressure over time, from which additional data can be obtained for analysis on demand.

Figure 3 shows that under the dual effects of solar radiation and atmospheric temperature heat exchange, the road surface temperature gradually decreases between 0:00 and 6:00 in the morning, because the solar radiation is almost zero and the temperature is low during this time period; The increase in solar radiation and air temperature causes the road surface temperature to gradually increase between 6 o'clock and 12 o'clock.

Figure 4 shows that under the action of the vehicle load, the vertical displacement of the midpoint of the load area gradually increases, and when the load reaches the peak value, the displacement reaches 0.17mm.

Figure 5 shows that under the action of the vehicle load, when the load reaches its peak value, the pore water pressure rises to 223kPa; after the load is removed, the pore water pressure reaches -50kPa. It is this suction that creates such a strong impact on the asphalt mix that the bitumen is stripped from the aggregate.

Claims (1)

1. A road surface structure multi-physical field coupling numerical simulation method based on finite elements is characterized in that the simulation method is carried out in multi-physical field coupling finite element software according to the following steps:

step 1: defining mechanical parameters, thermodynamic parameters and hydraulic parameters of the road material;

the mechanical parameters comprise Young modulus, poisson's ratio and density, and for the viscoelastic material, the mechanical parameters also comprise viscoelastic parameters;

the thermodynamic parameters include thermal conductivity, specific heat capacity and thermal expansion coefficient;

the hydraulic parameters include hydraulic conductivity, porosity, and Biot elasticity parameters;

step 2: defining a vehicle load function, and defining a meteorological function based on local meteorological data;

the vehicle load function F (t) is characterized by equation (1), the vehicle load includes a static load and a dynamic load, the dynamic load is a half-sinusoidal load:

in the formula (1), A is the load amplitude, t is the time, t ₀ The load acting duration in a single period is shown, T is the time interval between two adjacent loads, k =0,1,2 \ 8230;

the meteorological functions comprise a solar radiation daily variation function q (T) represented by an equation (2), and an atmospheric temperature daily variation function T represented by an equation (3) _a And a rainfall day variation function set according to the meteorological data;

wherein:

q ₀ maximum daily radiation intensity, q ₀ =0.131mq, m =12/c; q is the total daily solar radiation; c is the actual sunshine duration; omega is angular frequency;

is the average daily air temperature of the air,

T _m the change range of the daily air temperature is,

the temperature is the highest temperature of the day,

the daily minimum temperature; t is ₀ Is an initial phase;

and step 3: adding a solid mechanical module, specifying the material properties of each structural layer aiming at the solid mechanical module, and adding boundary conditions;

the governing equation of the solid mechanics module is characterized by equation (4):

wherein:

s is total stress, F is external force vector, rho _s Is the density of the solid material, u is the displacement vector,

is the divergence of the S and is,

is the second time derivative of the displacement vector;

the specifying material properties of the structural layers includes: appointing each structural layer road material to be an elastic or visco-elastic material, and setting the thermal expansion property of each structural layer road material, wherein the thermal expansion property is the thermal expansion coefficient and the strain reference temperature;

the boundary conditions include: fixed constraint that all direction displacements are 0, specified displacement that makes a certain direction displacement be 0, and boundary load;

and 4, step 4: adding a porous medium heat transfer module, specifying material properties for the porous medium heat transfer module, and defining a heat flux of a road surface;

the control equation of the porous medium heat transfer module is characterized by equation (5):

wherein: rho _l Is the density of the fluid, C _p Is constant pressure heat capacity of fluid (rho) _l C _p ) _eff For effective volumetric constant pressure heat capacity, T _tem Is the temperature of the liquid to be treated,

as first time derivative of temperature, u _l In order to be a field of fluid velocity,

for temperature gradients, q is the conducted heat flux,

k _eff in order to be an effective thermal conductivity coefficient,

is the divergence of Q, Q _h Is a heat source;

the specified material properties refer to imparting thermodynamic parameters defined in step 1 into the material;

the step 2 of defining the heat flux of the road surface is to add a heat flux interface in the porous medium heat transfer module and input a solar radiation daily variation function q (T) and an atmospheric temperature daily variation function T defined in the step _a Wherein the heat flux type of the solar radiation daily variation function q (T) is generalized inward heat flux, and the atmospheric temperature daily variation function T _a The type of heat flux of (1) is convective heat flux, which is a function of the daily variation of the atmospheric temperature T _a The associated convective heat flux defines the heat transfer coefficient h _c ：h _c ＝3.7v _w +9.4，v _w Is the wind speed;

and 5: adding a Darcy law module, specifying material properties aiming at the Darcy law module, and adding boundary conditions;

the control equation of the Darcy's law module is characterized by equation (6):

in the formula (6), S _p To store coefficients, p _l In order to obtain the pore water pressure,

is p _l The first time derivative of (a) is,

is rho _l The divergence of (A), K is the hydraulic conductivity, g is the acceleration of gravity,

in order to be a pressure gradient of the fluid,

unit vector in the direction of gravity, Q _m Is a quality source item;

the specified material properties refer to the hydraulic conductivity and porosity in the hydraulic parameters defined in step 1 are imparted into the material;

the boundary conditions include pore water pressure and normal rainwater inflow speed on the boundary;

step 6: adding a multi-physical-field coupling module, and inputting Biot elastic parameters into the multi-physical-field coupling module;

the governing equation of the multi-physical field coupling module is characterized by equation (7):

wherein: alpha is alpha _B Is the Biot-Willis coefficient, I is the identity matrix, ε _vol In order to be a volume strain,

is epsilon _vol The first time derivative of (a) is,

is S-alpha _B p _l The divergence of the beam of I is determined,

is rho _l u _l Divergence of (d);

and 7: calculating and post-processing analysis

The calculation refers to dividing finite element grids aiming at the constructed pavement structure body, and performing coupling solution on control equations of all modules by using a finite element method;

the post-processing analysis comprises stress field analysis, displacement field analysis, temperature field analysis and pore water pressure field analysis of the pavement structure, and a stress cloud picture, a displacement cloud picture, a temperature cloud picture, a pore water pressure cloud picture and the like are obtained, so that the multi-physical field coupling analysis of the pavement structure is completed.

Images