Method for constructing multi-scale coupled freeze-thaw-dp concrete phase field fracture model
By using a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model, the problem of the influence of freeze-thaw cycles on concrete fracture behavior was solved, and accurate prediction and dynamic simulation of concrete fracture behavior under freeze-thaw conditions were achieved, overcoming the problems of scale disconnect and lack of physical basis for parameters in traditional models.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- Filing Date
- 2026-06-08
- Publication Date
- 2026-07-03
AI Technical Summary
Existing technologies struggle to accurately characterize the effects of freeze-thaw cycles on concrete fracture behavior, particularly the dynamic impact of tensile-compressive asymmetric fracture behavior and freeze-thaw damage on the fracture process. Furthermore, the lack of quantitative correlation between microstructural features and macroscopic fracture behavior leads to insufficient prediction accuracy.
A multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model is adopted. The micro-macro damage mapping is established through asymptotic homogenization theory. The Drucker-Prager criterion is introduced to describe tensile-compressive asymmetric fracture. The series coupling relationship between freeze-thaw damage and mechanical damage is established through strain equivalence assumption. A model of crack initiation, propagation and failure of concrete under freeze-thaw conditions is constructed.
It achieves accurate prediction of concrete fracture behavior under freeze-thaw conditions, overcomes the problem of empirical value taking of characteristic length in traditional models, can naturally simulate the initiation, propagation and merging process of cracks, reduces numerical complexity and improves prediction accuracy.
Smart Images
Figure CN122333928A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of durability assessment and fracture mechanics simulation technology for concrete structures in cold regions, and in particular to a method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model. Background Technology
[0002] Freeze-thaw cycles are a major cause of degradation in concrete structures in cold regions. Concrete structures such as bridges, dams, tunnels, and ports operating in cold regions not only bear complex mechanical loads but also experience diurnal and seasonal freeze-thaw cycles. These cycles cause repeated freezing and thawing of pore water within the concrete, resulting in volume expansion and contraction. This triggers the initiation, propagation, and even penetration of microcracks, ultimately leading to a decrease in structural strength, stiffness, and durability. Therefore, accurately predicting the fracture behavior of concrete under freeze-thaw conditions is of significant theoretical and engineering value for the durability design and life assessment of engineering structures in cold regions.
[0003] Extensive research has been conducted by scholars both domestically and internationally on freeze-thaw damage to concrete, primarily focusing on the description of performance degradation at the material level. Descriptions of material performance degradation under a single environmental influence mainly focus on the empirical degradation patterns of macroscopic performance indicators. When concrete that has undergone freeze-thaw degradation continues to bear loads, its fracture behavior exhibits more complex characteristics, not only displaying typical quasi-brittle fracture properties but also accompanied by strength degradation, stiffness reduction, and increased ductility due to freeze-thaw damage. Existing freeze-thaw damage models struggle to reflect the dynamic impact of freeze-thaw damage on the fracture process, and the coupling mechanism between freeze-thaw degradation and mechanical fracture has not yet been systematically studied.
[0004] In the field of fracture simulation, existing research also has significant limitations and is independent of freeze-thaw degradation research. Phase-field fracture models have received widespread attention in recent years because they can naturally simulate the entire process of crack initiation, propagation, bifurcation, and merging without pre-setting crack paths. Existing research models mainly focus on simulating fracture behavior under ambient temperature conditions and do not consider the deterioration effect of freeze-thaw cycles on material properties. Traditional fracture models cannot accurately characterize the significant tensile-compressive asymmetric fracture behavior of quasi-brittle materials such as concrete: brittle fracture under tension and ductile-shear failure under compression. Especially after freeze-thaw degradation, the dynamic evolution of key mechanical parameters of concrete, such as tensile-compressive strength ratio, internal friction angle, and cohesion, with the number of freeze-thaw cycles is difficult to describe, resulting in severely insufficient prediction accuracy. In addition, most existing fracture models are macroscopic phenomenological models, and their constitutive parameters lack clear microscopic physical meaning, making it difficult to establish a quantitative correlation between the microstructural characteristics of concrete and macroscopic fracture behavior.
[0005] Existing technologies for simulating concrete freeze-thaw-fracture coupling have the following core shortcomings:
[0006] (1) Mechanism separation: Freeze-thaw degradation studies mainly focus on macroscopic empirical attenuation laws, while fracture mechanics models are mostly designed for ambient temperature environments. The two are independent of each other in terms of theoretical framework and numerical implementation, and lack an effective coupling mechanism. The interaction mechanism between freeze-thaw damage and stress-induced damage is still unclear, and it is difficult to reflect the dynamic amplification effect of freeze-thaw damage on the fracture process.
[0007] (2) Constitutive limitations: Existing fracture models are difficult to accurately characterize the tensile-compressive asymmetric fracture behavior of concrete. Although the Drucker-Prager (DP) criterion can describe pressure sensitivity, it has not yet been coupled with freeze-thaw degradation. At the same time, freeze-thaw cycles cause significant changes in key parameters such as the tensile-compressive strength ratio, internal friction angle, and cohesion of concrete. Existing models cannot reflect the dynamic evolution of these parameters with the number of freeze-thaw cycles.
[0008] (3) Scale disconnect: Most existing models are macroscopic phenomenological models, and their constitutive parameters lack clear microscopic physical meaning, making it difficult to establish a physical relationship between microscopic structural features and macroscopic fracture behavior. In addition, the grid dependency problem in phase-field fracture simulation has not been fundamentally solved, and there is a lack of regularization schemes based on physical mechanisms. Summary of the Invention
[0009] This invention provides a method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model, which integrates freeze-thaw degradation, tensile-compressive asymmetric fracture, and multi-scale damage mapping. It can accurately predict the fracture behavior of concrete under freeze-thaw conditions, providing a theoretical basis and numerical tool for durability assessment in cold regions.
[0010] To solve the above-mentioned technical problems, the technical solution proposed by this invention is as follows:
[0011] A method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model includes the following steps:
[0012] Step S1: Based on the asymptotic homogenization theory, establish the micro-macro damage mapping of concrete and derive the expression for the phase field characteristic length related to the microstructure.
[0013] Step S2: Introduce the Drucker-Prager criterion into the phase field fracture framework to decompose the total strain energy density into a damage-driven part and a protected part to describe the tensile-compressive asymmetric fracture behavior of concrete.
[0014] Step S3: Based on the strain equivalence assumption, establish the series coupling relationship between freeze-thaw damage and mechanical damage. By introducing a freeze-thaw damage factor into the phase field evolution equation to correct the strain energy release rate, the dynamic influence of freeze-thaw degradation on the mechanical fracture process is simulated.
[0015] Step S4: Based on steps S1-S3, construct the control equations of the coupled damage model, and perform numerical discretization and solution using the finite element method to predict the entire process of crack initiation, propagation and failure of concrete structures under freeze-thaw conditions.
[0016] A further improvement to the above technical solution is as follows:
[0017] Preferably, in step S1, the phase field characteristic length is represented as a strain rate-related function, with the following expression:
[0018]
[0019] In the formula, The phase field characteristic length, For quasi-static characteristic length, The current strain rate; Characteristic strain rate is an intrinsic parameter of the material; For scaling parameters, and ; It is the triaxial stress coefficient, and .
[0020] Preferably, in step S2, the total strain energy density is decomposed into a damage-driving part and a protected part:
[0021]
[0022] In the formula, The total strain energy density; Damage-driven strain energy density; The protected strain energy density; For strain tensor; As a damage variable, , Indicates no loss. Indicates complete destruction; This is the degradation function.
[0023] Preferably, in step S2, the damage-driven strain energy density is constructed based on the Drucker-Prager criterion. The expression is:
[0024]
[0025] In the formula, It is the elastic modulus; The main part of the brackets in Macaulay; For the stress tensor in the undamaged state, , The elastic tensor is in a lossless state; The second deviatoric stress invariant is based on the stress tensor of the lossless state. This represents the positive part of the first stress invariant based on the stress tensor of the undamaged state; These are parameters related to the internal friction angle.
[0026] Preferably, in step S2, the protected strain energy density is determined based on the Drucker-Prager criterion. The expression is:
[0027]
[0028] In the formula, For the protected strain energy density, Bulk modulus These are parameters related to the internal friction angle. Shear modulus It is the elastic modulus; The negative part of the brackets in Macaulay. The main part of the brackets in Macaulay; As the first strain invariant, This is the second partial strain invariant.
[0029] Preferably, in step S3, stiffness reduction is used to couple freeze-thaw damage and mechanical damage based on the strain equivalence assumption, and the effective stiffness matrix is expressed as:
[0030]
[0031] In the formula, This represents the effective stiffness matrix after freeze-thaw and mechanical coupling damage. The elastic tensor is in a lossless state; Let be the mechanical damage degradation function. It is a freeze-thaw damage factor. For freeze-thaw damage variables, For damage variables.
[0032] Preferably, the coupled damage model in step S4 consists of the following governing equations:
[0033] Momentum balance equation:
[0034]
[0035] In the formula, For the total spatial gradient operator, For stress tensor, The density of force per unit volume;
[0036] Stress-strain-damage relationship:
[0037]
[0038] In the formula, For stress tensor, For the lossless elastic tensor, For freeze-thaw damage variables, As a damage variable, For strain tensor;
[0039] Geometric equations:
[0040]
[0041] In the formula, For strain tensor, Let be the displacement gradient tensor. For displacement field, This is a transpose operation;
[0042] Phase field evolution equation:
[0043]
[0044] In the formula, To be related to the number of freeze-thaw cycles The relevant fracture energy, As a damage variable, The characteristic length of the phase field; For the Laplace operator, For historical variables.
[0045] Preferably, in step S4, the displacement field and phase field are solved alternately using a separate iteration method. Specifically, in the current increment step, the phase field is fixed first, and the displacement field is solved based on the momentum balance equation; then the displacement field is fixed, and the phase field is solved based on the phase field evolution equation; the above alternating process is repeated until both the displacement residual and the phase field residual satisfy the preset convergence tolerance; in each increment step, the linearized discrete equation system is solved using Newton-Raphson iteration.
[0046] Preferably, in step S4, the numerical solution of the coupled damage model is achieved by writing the UMAT user material subroutine of the Abaqus software. The subroutine performs the following operations in each increment step: reads the current strain increment and freeze-thaw damage state variables, updates the stress and Jacobian matrices according to the coupled constitutive model, calls the phase field solution submodule to update the historical variables and phase field variables, and stores the state variables for use in the next increment step.
[0047] The multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model construction method provided by this invention has the following advantages compared with the prior art:
[0048] (1) The multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model construction method of the present invention is the first to unify the asymptotic homogenization multi-scale theory, the Drucker-Prager quasi-brittle fracture criterion, and the freeze-thaw damage evolution theory into the same phase field framework. By establishing a quantitative mapping relationship between microstructural parameters and macroscopic fracture parameters, each parameter in the model has a clear physical meaning, fundamentally overcoming the shortcomings of traditional phase field models where characteristic length values are empirically determined and lack physical basis. At the same time, the strain energy decomposition scheme based on the Drucker-Prager criterion can naturally distinguish between tensile brittle fracture and compressive ductile shear failure of concrete, without the need to manually set fracture mode judgment criteria, so that the model can maintain physical consistency under different stress states.
[0049] (2) The multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model construction method of the present invention dynamically reduces the effective stiffness of the material by the freeze-thaw damage factor and updates the fracture energy, strength parameters and Drucker-Prager parameters simultaneously, thereby realizing the accurate simulation of the amplification effect of freeze-thaw degradation on the mechanical fracture process.
[0050] (3) The multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model construction method of the present invention adopts a series coupling scheme based on the strain equivalence assumption, and directly applies the freeze-thaw damage factor to the lossless elastic tensor, which significantly reduces the numerical complexity of the multi-field coupled model. By writing the UMAT user material subroutine of Abaqus software, the proposed model is seamlessly embedded in the commercial finite element platform. In the solution process, the displacement field and phase field are updated alternately and iteratively, which results in stable convergence and high computational efficiency. The model can naturally simulate the entire process of crack initiation, propagation, bifurcation, merging and diffuse damage without pre-setting crack paths, thus avoiding mesh sensitivity. Attached Figure Description
[0051] Figure 1 This is a schematic diagram comparing the length of the microstructure with the length of the macroscopic feature in this invention.
[0052] Figure 2 The diagram shows the micro-cracks and intact areas of the damaged material, where (a) represents the undeformed state and (b) represents the deformed state.
[0053] Figure 3 This is a flowchart for the freeze-thaw cycle test in the experimental verification.
[0054] Figure 4 The shear damage cloud maps were used to verify the different freeze-thaw cycles in the experiment. (a) represents 0 freeze-thaw cycles, (b) represents 25 freeze-thaw cycles, (c) represents 50 freeze-thaw cycles, and (d) represents 75 freeze-thaw cycles.
[0055] Figure 5 To experimentally verify the load-displacement curves of the model under freeze-thaw damage.
[0056] Figure 6 To experimentally verify the stress-strain curve of concrete under freeze-thaw cycles. Detailed Implementation
[0057] The following provides a detailed description of specific embodiments of the present invention. It should be understood that the specific embodiments described herein are for illustrative and explanatory purposes only and are not intended to limit the scope of the invention.
[0058] The multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model construction method of the present invention includes the following steps:
[0059] Step S1: Based on the asymptotic homogenization theory, establish the micro-macro damage mapping of concrete and derive the expression for the phase field characteristic length related to the microstructure.
[0060] A physical relationship between the microstructural characteristics of concrete and its macroscopic fracture parameters is established. By introducing the scale separation assumption and the asymptotic homogenization method, the elastic tensor of the microscopic non-homogeneous material is averaged into a macroscopic effective constitutive relation, and an expression for the characteristic length of the phase field with clear microscopic physical meaning is derived, solving the problem of empirically determined characteristic length values and lack of physical basis in traditional phase field models.
[0061] S1-1, the basic assumptions of multiscale analysis and scale separation.
[0062] The microstructure of concrete (aggregates, pores, microcracks, etc.) exhibits significant scale differences compared to its macrostructure. To accurately describe this multi-scale characteristic, the scale separation assumption is introduced, meaning that the phase field feature length is much smaller than the macroscopic feature length, such as... Figure 1 As shown. The scale separation parameter is defined as:
[0063] (1)
[0064] In the formula, It is the characteristic length of the phase field, which is related to microstructural features such as microcrack size and aggregate particle size; The macroscopic characteristic length represents the size of the structural component or the loading wavelength scale; is the scale separation parameter. This ensures that the micro and macro scales are mathematically separable.
[0065] A dual-scale coordinate system is established based on the scale separation assumption: macroscopic coordinates. Describes macroscopic location, microscopic coordinates Describe the relative position within the RVE. Represents scale separation parameters, such as Figure 1 As shown.
[0066] According to the asymptotic homogenization theory, the displacement field can be expressed as a function of macroscopic and microscopic coordinates, and the parameters can be separated by scale. Expand into an asymptotic series:
[0067] (2)
[0068] In the formula, Representation depends on macroscopic coordinates and micro coordinates The total displacement field, The macroscopic uniform displacement component describes the overall deformation of the material and mainly relies on macroscopic coordinates. ; These are first-order microscopic fluctuation displacement components, caused by microscopic structural inhomogeneities, and dependent on macroscopic coordinates. and micro coordinates ; For second-order and higher-order small quantities. This represents the scale separation parameter.
[0069] Under the assumption of small deformation, the strain tensor is defined as the symmetric part of the displacement gradient:
[0070] (3)
[0071] In the formula, The components representing the strain tensor are macroscopic coordinates. and micro coordinates The function; and Let be the components of the displacement vector, where the subscripts indicate the direction; and For macroscopic spatial coordinate components; subscript and It is a free index, with a value range of [value range missing]. , respectively corresponding to the Cartesian coordinate system , , Three directions.
[0072] Considering the relationship between microscopic and macroscopic coordinates, the spatial gradient operator can be decomposed into two parts: macroscopic gradient and microscopic gradient. ,in Represents the total spatial gradient operator. Represents the macroscopic gradient operator. Let represent the micro-gradient operator. Using the chain rule, substituting equation (2) into equation (3) and rearranging, we can obtain the multi-scale decomposition form of the strain tensor:
[0073] (4)
[0074] In the formula, Represents the total strain tensor; It represents the zero-order strain tensor and is the dominant term in the asymptotic expansion; This represents the first-order strain tensor, corresponding to the strain contributed by microscopic fluctuation displacements. Multiplying it by the scale separation parameter ε results in a first-order correction term for the strain. It is a second-order or higher-order small quantity.
[0075] Zero-order strain tensor It can be further broken down into:
[0076] (5)
[0077] In the formula, For macroscopic average strain, it depends only on macroscopic coordinates. The zero-order macroscopic displacement field represents the average deformation of the material as a whole; Microscopic fluctuation strain reflects local strain disturbances caused by microstructural inhomogeneities. Let be a first-order microscopic fluctuation displacement field, representing local displacement fluctuations relative to the macroscopic average motion, which satisfies periodic boundary conditions on the RVE boundary; where , , , Represents the zeroth-order macroscopic uniform displacement component, subscript and For free indicators; , This represents the first-order micro-fluctuation displacement component. , Represents the components of microscopic spatial coordinates. , For macroscopic spatial coordinate components, subscript and For free indicators.
[0078] To ensure compatibility between Representative Volume Elements (RVEs) and simplify the solution of microscopic problems, periodic boundary conditions are applied to the microscopic fluctuation displacements. On the relative boundaries of the RVEs, the microscopic fluctuation displacements should satisfy periodicity, while the surface forces should satisfy anti-periodicity. The mathematical expression of the periodic boundary conditions is as follows:
[0079] (6)
[0080] In the formula, This represents the value of the micro-fluctuation displacement after translating a periodic vector; This represents the first-order micro-fluctuation displacement component, that is, the value of the micro-fluctuation displacement at the original position; Let be the periodic vector of the RVE, which describes the periodicity of the microstructure. Equation (6) not only ensures the displacement continuity between adjacent RVEs, but also provides the necessary mathematical constraints for the subsequent derivation of the effective elastic tensor.
[0081] S1-2, an effective constitutive relation based on asymptotic homogenization.
[0082] At the microscale, assuming the material within the RVE is linearly elastic, the micro-stress and micro-strain satisfy a linear constitutive relationship, and the micro-equilibrium equation is:
[0083] (7)
[0084] In the formula, It is the micro-stress tensor, which reflects the stress state at the micro-scale. Let be the microscopic elastic tensor, which is about the microscopic coordinates. The function reflects the non-uniformity of the microstructure; It is the microscopic strain tensor.
[0085] For isotropic materials:
[0086] (8)
[0087] In the formula, For the microelastic tensor, This is Lamé's first constant, which varies with microscopic position; This refers to the shear modulus, which varies with microscopic position. Let Kronecker function be used when hour ,otherwise ; Indicates when hour, ,otherwise Indicates when hour, ,otherwise ; Indicates when hour, ,otherwise ; Indicates when hour, ,otherwise ; Indicates when hour, ,otherwise .
[0088] By using asymptotic homogenization methods, microscale inhomogeneities can be averaged down to the macroscale, yielding macroscopically effective constitutive relations. The core of this process is solving for the characteristic displacement field. It represents the deformation response of the microstructure under unit macroscopic strain. The characteristic displacement field satisfies the following governing equations:
[0089] (9)
[0090] In the formula, For macroscopic effective elasticity tensor; For RVE volume, For RVE fields; Let be the characteristic displacement field, be the characteristic displacement function, and satisfy the periodic boundary conditions. This represents the value of the macroscopic strain component. Indicates the direction index of the characteristic displacement function; , Represents the microelastic tensor under different boundary conditions. The microscopic coordinate components are represented. Equation (9) condenses the geometric and material information of the microstructure into macroscopic constitutive parameters, realizing cross-scale information transmission.
[0091] Macroscopic stress is defined as the volume average of microscopic stress on the RVE, and its expression is:
[0092] (10)
[0093] In the formula, Represents the macroscopic stress tensor. It is the micro-stress tensor, which describes the stress state at the micro-scale.
[0094] Substituting equation (7) and the solution of the characteristic displacement field, we obtain the macroscopic stress-strain relationship:
[0095] (11)
[0096] In the formula, For macroscopic strain tensor, and , Indicators of Freedom The macroscopic average strain It is the macroscopic effective elastic tensor. That is, the non-uniform material is equivalent to a uniform anisotropic material, and its stiffness is determined by equation (9).
[0097] S1-3, strain rate-dependent phase field characteristic length.
[0098] Under dynamic loading conditions, the fracture behavior of the material exhibits a significant rate-dependent characteristic. Based on Griffith's fracture theory and fragment size law, a strain rate-dependent phase field characteristic length is introduced, expressed as follows:
[0099] (12)
[0100] In the formula, The phase field characteristic length, The quasi-static characteristic length; Given the current strain rate, , Represents the components of the strain rate tensor; Characteristic strain rate is an intrinsic parameter of the material; The scaling parameter controls the rate sensitivity intensity. ; It is the triaxial stress coefficient, and Equation (12) reflects the influence of complex stress states on the characteristic length. The physical meaning of Equation (12) is that as the strain rate increases, the characteristic length decreases and the damage becomes more localized, which is consistent with experimental observations that the fracture surface of the material is smoother and the fragment size is smaller under high strain rates.
[0101] Quasi-static characteristic length It is determined by both the material's energy parameters and microstructural characteristics. Based on the classical crack band theory and considering the microstructural characteristics of concrete, this invention proposes the following expression:
[0102] (13)
[0103] In the formula, For quasi-static characteristic length, It is the elastic modulus; Fracture energy represents the energy required to generate a crack per unit area. The critical strength is the material's breaking strength. The average aggregate particle size is the average size of the aggregate in the concrete. Porosity , is the volume fraction of pores in concrete.
[0104] Stress state has a significant impact on damage evolution, especially under complex loading conditions. To reflect this effect, a triaxial stress coefficient is introduced. , defined as the ratio of mean stress to equivalent stress:
[0105] (14)
[0106] In the formula, For average stress, The equivalent stress is expressed as:
[0107] (15)
[0108] (16)
[0109] In the formula, It is the trace of the stress tensor, that is, the sum of the normal stress components; , , These are the normal stress components in the x, y, and z directions, respectively; For the deviatoric stress tensor, , Let be the components of the stress tensor, where the subscripts are... and It is a free index, with a value range of 1, 2, 3, when When represents the normal stress component, when Time represents the shear stress component; For average stress, Let Kronecker function represent the expression when... hour ,otherwise The range of values for the triaxial stress coefficient is as follows: Positive values represent triaxial tensile stress, negative values represent triaxial compressive stress, and zero values represent pure shear stress. In damage evolution, By influencing the characteristic length to adjust the degree of damage localization, stress state-related fracture behavior description can be achieved.
[0110] Step S2 introduces the Drucker-Prager criterion to decompose the total strain energy density into a damage-driven part and a protected part to describe the tensile-compressive asymmetric fracture behavior of concrete.
[0111] By introducing the Drucker-Prager criterion into the phase field framework, the total strain energy density is decomposed into a damage-driven part and a protected part. This allows tensile and shear strain energies exceeding the threshold to drive damage evolution, while pure compressive strain energy is protected, thus reproducing the typical characteristics of tensile brittle fracture and compressive ductile shear failure of concrete.
[0112] S2-1, Total Strain Energy Density Decomposition.
[0113] like Figure 2 As shown in (a) and (b), quasi-brittle materials such as concrete exhibit significant tension-compression asymmetry, exhibiting brittle fracture under tension and ductile shear failure under compression. To accurately describe this characteristic, this invention decomposes the total strain energy density into a damage-driven portion and a protected portion:
[0114] (17)
[0115] In the formula, The total strain energy density represents the elastic strain energy stored per unit volume. The damage-driven strain energy density is the part that drives damage evolution, abbreviated as... ; The protected strain energy density, the portion unaffected by damage, is abbreviated as... ; As a damage variable, , Indicates no loss. Indicates complete destruction; This is a degradation function that describes the decrease in material stiffness as damage occurs; For strain tensor.
[0116] A quadratic degradation function is employed to ensure the well-posedness of the damage evolution problem, and a numerical stability parameter is introduced:
[0117] (18)
[0118] In the formula, As a damage variable, It is a numerically stable parameter, and .
[0119] S2-2 introduces the Drucker-Prager criterion.
[0120] The Drucker-Prager criterion (DP criterion) is expressed as follows:
[0121] (19)
[0122] In the formula, This is the Drucker-Prager destructive function; Indicates the elastic state. Indicates destruction; Let be the first stress invariant, and , The trace of the stress tensor is the sum of the three normal stress components. Indicates compressive stress. Indicates tensile stress; It is the second deviatoric stress invariant, and , It is the deviatoric stress tensor; These are parameters related to cohesion. These are parameters related to the internal friction angle. The Drucker-Prager criterion describes the effect of pressure on material strength by introducing a first stress invariant: compressive stress increases strength, while tensile stress decreases strength.
[0123] parameter and Uniaxial tensile strength and uniaxial compressive strength Confirmed, we can obtain:
[0124] (20)
[0125] (twenty one)
[0126] For concrete materials, there are usually ,therefore This demonstrates the stress-enhancing effect.
[0127] Based on the Drucker-Prager criterion, the damage-driven strain energy density is constructed. The expression is:
[0128] (twenty two)
[0129] In the formula, It is the elastic modulus; The positive part of the brackets in Macaulay ensures that damage is only driven by the portion exceeding the destruction threshold. For the stress tensor in the undamaged state, , The elastic tensor is in a lossless state; The second deviatoric stress invariant is based on the stress tensor of the lossless state. This represents the positive part of the first stress invariant based on the stress tensor of the undamaged state; These are parameters related to the internal friction angle.
[0130] Protected strain energy density This corresponds to the elastic strain energy that can be stored even under complete damage. By solving the partial differential equations derived from thermodynamic consistency and considering the Drucker-Prager criterion as a constraint, we can obtain:
[0131] (twenty three)
[0132] In the formula, The protected strain energy density, abbreviated as ; Bulk modulus , Shear modulus ; Poisson's ratio, It is the elastic modulus; The negative part of the brackets in Macaulay. The main part of the brackets in Macaulay; As the first strain invariant, , , , Let be the normal strain components of the strain tensor in the x, y, and z directions, respectively; For the second partial strain invariant, ,in For the partial strain tensor, , For the total strain tensor, Let Kronecker function be used.
[0133] In equation (23), the first term Corresponding to the pure compressive strain energy density, where the negative part of the Macaulay brackets Ensure that only when the volume is compressed (i.e. When this term is not zero, the pure compressive strain energy can still be stored after the material is completely damaged, and therefore is not affected by the damage; the second term Corresponding to the shear strain energy density corrected by the Drucker-Prager criterion, where The positive value portion of the first strain invariant, outer layer This ensures that the term contributes energy only when the correction term is positive, reflecting the strengthening effect of pressure on shear strength, i.e., compressive stress increases shear resistance, while tensile stress decreases shear resistance.
[0134] S2-3, stress-strain relationship after damage.
[0135] Taking the partial derivative of equation (17), we can obtain the stress-strain relationship under the damage state:
[0136] (twenty four)
[0137] In the formula, For stress tensor, abbreviated as ; For damage-driven strain energy density, abbreviation; For the protected strain energy density, abbreviation; For strain tensor.
[0138] Substituting equations (22) and (23) into the equations and rearranging them, we can obtain the stress tensor. Explicit expression:
[0139] (25)
[0140] In the formula, For damage-driven strain energy density, For the protected strain energy density, For strain tensor.
[0141] This expression indicates that the stress after damage consists of two parts: one part is the undamaged stress after reduction by the degradation function, and the other part is the compensation term of the protected stress.
[0142] Based on the variational principle, the total potential energy functional is varied, and historical variables are introduced to ensure the irreversibility of damage evolution, resulting in the phase field evolution equation:
[0143] (26)
[0144] In the formula, This is the fracture energy; The characteristic length of the phase field; For the Laplace operator, , The sign of the second-order partial derivative; For macroscopic spatial coordinate components, ; For historical variables, it is defined as the maximum value of the damage-driven strain energy density during the loading history: , for Damage-driven strain energy density at time step For historical time variables, The current time. Equation (26) This represents the surface energy of the crack's resistance to damage. This represents the driving force of strain energy release on damage.
[0145] Step S3: Based on the strain equivalence assumption, establish the series coupling relationship between freeze-thaw damage and mechanical damage. By introducing a freeze-thaw damage factor into the phase field evolution equation to correct the strain energy release rate, the dynamic influence of freeze-thaw degradation on the mechanical fracture process is simulated.
[0146] The attenuation law of concrete material parameters under freeze-thaw cycles was established, and the series coupling of freeze-thaw damage and mechanical damage was realized through the strain equivalence assumption. Freeze-thaw damage variables were defined, and a linear evolution model of elastic modulus, fracture energy, strength parameters, and Drucker-Prager parameters with the number of freeze-thaw cycles was established. The freeze-thaw degradation effect was introduced into the phase field fracture framework by reducing the effective stiffness matrix, realizing a deep integration of environmental effects and mechanical fracture.
[0147] S3-1, Evolution of material parameters under freeze-thaw cycles.
[0148] Freeze-thaw cycle test procedure as follows Figure 3As shown, within the linear accumulation framework of the UMAT user material subroutine (i.e., the Fortran program interface for user-defined material constitutive models) in Abaqus software, this invention employs a linear decay model to describe the elastic modulus with respect to the number of freeze-thaw cycles. Changes:
[0149] (27)
[0150] In the formula, For experience Elastic modulus after one freeze-thaw cycle The initial elastic modulus, The maximum attenuation coefficient of the elastic modulus. ; The limit is the number of freeze-thaw cycles. The computational efficiency of this linear model is better than that of the exponential model. The damage to concrete under freeze-thaw cycles is shown in Table 1.
[0151] Table 1. Damage to concrete under freeze-thaw cycles
[0152]
[0153] The fracture energy is also subject to freeze-thaw cycle degradation. Drawing on the decay law of elastic modulus, the decay model of fracture energy can be expressed as:
[0154] (28)
[0155] In the formula, For experience Fracture energy after one freeze-thaw cycle The initial fracture energy, The attenuation coefficient of the fracture energy. For freeze-thaw damage variables.
[0156] Compressive strength and tensile strength are the most basic mechanical properties of concrete. The decay law of compressive strength with freeze-thaw cycles can be expressed as follows:
[0157] (29)
[0158] In the formula, For experience Compressive strength after one freeze-thaw cycle Initial compressive strength, The attenuation coefficient of compressive strength, , For freeze-thaw damage variables.
[0159] Tensile strength reduction:
[0160] (30)
[0161] In the formula, For experience Tensile strength after one freeze-thaw cycle This represents the initial tensile strength. The attenuation coefficient of tensile strength, The linear decay model used in formulas (27) to (30) and the freeze-thaw damage variable This maintains consistency with the linear accumulation model, making it easy to implement in the UMAT user material subroutine.
[0162] Drucker-Prager parameters and The changes with freeze-thaw cycles can be indirectly derived from the decrease in intensity. For parameters related to cohesion, For parameters related to the internal friction angle. Substituting equations (29) and (30) into equations (20) and (21), we get:
[0163] (31)
[0164] (32)
[0165] In the formula, For history Drucker-Prager parameters after one freeze-thaw cycle , For experience Drucker-Prager parameters after one freeze-thaw cycle , For experience Compressive strength after one freeze-thaw cycle For experience Tensile strength after one freeze-thaw cycle.
[0166] This evolutionary pattern reflects the impact of freeze-thaw degradation on the material's pressure sensitivity. As the number of freeze-thaw cycles increases, A gradual decrease indicates a decline in cohesion; while The change reflects the effect of the change in the tensile-compressive strength ratio on the friction angle.
[0167] To quantify the degree of material degradation caused by freeze-thaw cycles, a freeze-thaw damage variable is defined. as follows:
[0168] (33)
[0169] In the formula, For experience Elastic modulus after one freeze-thaw cycle This represents the initial elastic modulus. Damage factor. and As shown in Table 2, the number of freeze-thaw cycles That is, the number of cycles.
[0170] Table 2 Damage factors of concrete under freeze-thaw cycles
[0171]
[0172] Freeze-thaw damage variables satisfy , This indicates no freeze-thaw damage. This indicates complete freeze-thaw damage.
[0173] Combining the linear decay model (27), the evolution rate of freeze-thaw damage with the number of cycles is:
[0174] (34)
[0175] In the formula, The maximum attenuation coefficient of the elastic modulus. For the number of freeze-thaw cycles, This represents the maximum number of freeze-thaw cycles.
[0176] Freeze-thaw damage accumulates at a constant rate, and this simplified form allows for direct calculation in the UMAT user material subroutine using the ratio of the current analysis time to the freeze-thaw cycle.
[0177] S3-2, Coupled Damage Evolution Theory.
[0178] Based on the strain equivalence assumption in damage mechanics, this invention employs a simplified stiffness reduction method to couple freeze-thaw damage with mechanical damage. Freeze-thaw damage variables. The mechanical response is indirectly influenced by reducing the initial stiffness matrix of the material, without explicitly modifying the driving force terms in the phase-field evolution equations. This simplification significantly reduces the complexity of numerical implementation and effectively reflects the weakening effect of freeze-thaw degradation on the macroscopic mechanical properties of the material under small deformation quasi-static conditions. Specifically, the effective stiffness matrix is expressed as:
[0179] (35)
[0180] In the formula, To account for the effective stiffness matrix after considering freeze-thaw and mechanical coupling damage; The elastic tensor is in a lossless state; This represents the mechanical damage degradation function. It is a freeze-thaw damage factor.
[0181] Step S4: Based on steps S1 to S3, construct the control equations of the coupled damage model, and perform numerical discretization and solution using the finite element method to predict the entire process of crack initiation, propagation and failure of concrete structures under freeze-thaw conditions.
[0182] A complete set of governing equations coupling freeze-thaw damage and mechanics was constructed, including momentum balance equations, geometric equations, constitutive equations, and phase field evolution equations. The governing equations were spatially discretized using the finite element method, and the displacement field and phase field were solved alternately using a separation iterative method. Numerical solutions were implemented using the UMAT user material subroutine in Abaqus software, ultimately outputting predicted results such as damage contour maps, stress-strain curves, load-displacement curves, and crack propagation paths.
[0183] S4-1, Construct a coupling damage model.
[0184] Taking into account momentum balance, geometric relationships, constitutive relations, and phase field evolution, this invention proposes a coupled damage model consisting of the following governing equations:
[0185] (1) Momentum balance equation (quasi-static):
[0186] (36)
[0187] In the formula, For the total spatial gradient operator, For stress tensor, This is the density of force per unit volume.
[0188] (2) Stress-strain-damage relationship:
[0189] (37)
[0190] In the formula, For the lossless elastic tensor, This represents the mechanical damage degradation function. It is a freeze-thaw damage factor. For freeze-thaw damage variables, As a damage variable, For strain tensor.
[0191] In the incremental implementation of the UMAT user material subroutine, the stress-strain-damage relationship is discretized in an incremental form:
[0192] (38)
[0193] In the formula, Let be the stress increment tensor, representing the change in stress within the current increment step; For the strain increment tensor, it represents the change in strain within the current increment step.
[0194] The stress at the current increment step is obtained by summing the stress at the previous increment step and the current stress increment, i.e. , This is the stress tensor at the end of the current increment step; This is the stress tensor at the end of the previous increment step, i.e., the stress at the beginning of the current increment step; The time step size of the current increment step. This is the current time.
[0195] (3) Geometric equations (small deformation assumption):
[0196] (39)
[0197] In the formula, For strain tensor, Let be the displacement gradient tensor. For displacement field, This is for the transpose operation.
[0198] (4) Phase field evolution equation:
[0199] (40)
[0200] In the formula, To be related to the number of freeze-thaw cycles Related fracture energy; These are damage variables, i.e., phase field variables; The characteristic length of the phase field; For the Laplace operator, , The sign of the second-order partial derivative. For macroscopic spatial coordinate components ( ); These are historical variables used to ensure the irreversibility of damage evolution.
[0201] This set of equations constitutes a strongly coupled set of nonlinear partial differential equations, which needs to be solved numerically.
[0202] To ensure the boundary value problem is well-posed, appropriate boundary conditions must be applied to the boundaries. The displacement boundary conditions and the traction force boundary conditions are as follows:
[0203] (41)
[0204] (42)
[0205] In the formula, For displacement field, Given the boundary displacement values, To specify the boundary of the displacement, To specify the boundary of the traction force, Indicated on the boundary, used to specify the boundary area where boundary conditions are applied; For stress tensor; It is the unit vector of the outward normal.
[0206] For phase field variables, homogeneous Neumann boundary conditions are typically used, indicating that the normal derivative of the damage at the boundary is zero.
[0207] (43)
[0208] In the formula, Phase field variables The gradient represents the rate of change of the damage variable in space; Let be the outward normal unit vector. To find the solution domain.
[0209] This condition ensures a natural transition of damage evolution at the boundary.
[0210] To facilitate finite element discretization, the governing equations are transformed into a weak form. The weak form of the momentum equation is:
[0211] (44)
[0212] In the formula, To solve for the domain, This is the symmetrical part of the virtual displacement gradient. This is a virtual displacement. For volume, For area, The density of force per unit volume. To specify the boundary of the traction force, For stress tensor, This represents the displacement field.
[0213] Weak form of the phase field equation:
[0214] (45)
[0215] In the formula, To solve for the domain, To be related to the number of freeze-thaw cycles The relevant fracture energy, As a damage variable, The phase field characteristic length, For the total spatial gradient operator, This is a virtual displacement. For historical variables, For volume.
[0216] Equations (44) and (45) form the weak formal basis for solving coupled problems.
[0217] S4-2, numerical discretization and solution are performed using the finite element method.
[0218] The finite element method is used to discretize the weak form spatially. The solution domain is then defined. Divided into a finite number of elements, the displacement field and phase field are approximated by nodal value interpolation:
[0219] (46)
[0220] (47)
[0221] In the formula, Let be the displacement field, representing the displacement field within a range. In this formula, it represents the displacement field at any point within the element. Number the cell nodes with indexes. ; This represents the number of unit nodes. For the displacement field The shape function of each node. For the phase field Shape functions of nodes; For the first The degrees of freedom of displacement of each node, i.e., the nodal displacement values; For the first Phase field degrees of freedom of each node (node damage variable values); These are damage variables, i.e., phase field variables.
[0222] The strain tensor can be represented by the nodal displacements through a strain-displacement matrix:
[0223] (48)
[0224] In the formula, For the first The strain-displacement matrix of each node is composed of the partial derivatives of the shape functions; For strain tensor, As a damage variable, The number of unit nodes, Number the cell nodes with indexes. For the first The displacement degrees of freedom of each node.
[0225] Substituting the discrete approximation into the weak form, we obtain the displacement residual vector:
[0226] (49)
[0227] In the formula, For the corresponding to the first The displacement residual vector with each degree of freedom. For the first The strain-displacement matrix of each node. For the displacement field The shape function of each node. Represents the stress tensor. For volume, To solve for the domain, The density of force per unit volume. To specify the boundary of the traction force, For transpose operation, For the current time, For area.
[0228] Substituting the discrete approximation into the weak form, we obtain the phase field residual vector:
[0229] (50)
[0230] In the formula, For the corresponding to the first The phase field residual vector with each phase field degree of freedom. For the first The phase field gradient-shape function matrix at each node is composed of the partial derivatives of the phase field shape functions, i.e. , For macroscopic spatial coordinate components, subscript For free indicators, ; For transpose operation, To solve for the domain, To be related to the number of freeze-thaw cycles The relevant fracture energy, The phase field characteristic length, For the total spatial gradient operator, For historical variables, For volume, For the phase field The shape function of each node. For phase field variables, Phase field variables The gradient.
[0231] The goal of solving the discrete system of equations is to make the residual vector zero.
[0232] The Newton-Raphson iterative method is used to solve the nonlinear equations:
[0233] (51)
[0234] In the formula, The algorithm-consistent tangent stiffness tensor is obtained by linearizing the constitutive relation; The displacement-displacement tangent stiffness matrix is... Number the cell nodes with indexes; For the first The displacement degrees of freedom of each node, Number the cell nodes with indexes. For the corresponding to the first The displacement residual vector with each degree of freedom. For the first The strain-displacement matrix of each node. For the first The strain-displacement matrix of each node. For volume, This is for the transpose operation.
[0235] Similarly, the displacement-phase field stiffness matrix can be derived. Phase field-phase field stiffness matrix The formula is as follows:
[0236] (52)
[0237] (53)
[0238] In the formula, For the corresponding to the first The displacement residual vector with each degree of freedom. For the first The strain-displacement matrix of each node. For the first The phase field degrees of freedom of each node, For phase field variables; For freeze-thaw damage variables; The elastic tensor is in a lossless state; For strain tensor; For the phase field The shape function of each node. For the phase field The shape function of each node. To solve for the domain, This is a transpose operation; For the corresponding to the first The phase field residual vector with each phase field degree of freedom. To be related to the number of freeze-thaw cycles The relevant fracture energy, The phase field characteristic length, For the first The phase field gradient-shape function matrix of each node. For the first The phase field gradient-shape function matrix of each node. For historical variables;
[0239] Assemble the displacement and phase field degrees of freedom into a global vector, and the nonlinear equations can be written as follows:
[0240] (54)
[0241] In the formula, The displacement-displacement tangent stiffness matrix has the following component form. As defined in formula (51), represents the derivative of the displacement degree of freedom with respect to the displacement residual, and represents the derivative of the displacement degree of freedom with respect to the displacement residual; Let be the displacement-phase-field coupling stiffness moment, and let be the derivative of the phase field degrees of freedom with respect to the displacement residual; Let be the phase-displacement coupling stiffness matrix, representing the derivative of the displacement degrees of freedom with respect to the phase residual; Let be the phase field-phase field tangent stiffness matrix, and let represent the derivative of the phase field degrees of freedom with respect to the phase field residuals. This is the vector of displacement degree of freedom increments; This is the phase field degree of freedom increment vector; The displacement residual vector; Let be the phase field residual vector.
[0242] In each iteration step, solve for the increment. and Then update the displacement and phase field:
[0243] (55)
[0244] (56)
[0245] In the formula, For the first The displacement degree of freedom vector after the next iteration (updated value); For the first The displacement degree of freedom vector at the next iteration (current value); For the first The phase field degree of freedom vector after the next iteration (updated value); For the first Phase field degree of freedom vector at the next iteration (current value); Number the Newton-Raphson iteration steps. ; This is the vector of displacement degree of freedom increments; This is the phase field degree of freedom increment vector.
[0246] Iterate until the norm of the residual vector is less than the preset tolerance.
[0247] Experimental verification:
[0248] (a) Simulation of direct shear crack propagation under freeze-thaw cycles.
[0249] A direct shear specimen model was established, with pre-defined shear band regions on both sides of the shear plane to induce crack initiation and propagation at predetermined locations. Boundary conditions were set as follows: Bottom constraint: All degrees of freedom were fixed in the lower region of the specimen; Top loading: A vertical displacement control load was applied to the upper region at a loading rate of [value missing]. Shear: Apply a horizontal displacement to a 24 mm long region on the left edge to simulate shear action.
[0250] The model was spatially discretized using the finite element method. The model employed approximately 80,000 four-node plane strain quadrilateral elements, with mesh refinement in the expected crack propagation region to ensure that the element feature size was less than half the phase field feature length, guaranteeing accurate crack path capture. Material parameters were calibrated based on experimental data, as shown in Table 3.
[0251] Table 3 Material parameters for the direct shear model
[0252]
[0253] Figure 4 (a) to (d) show the damage contour plots of specimens with different freeze-thaw cycles. In the figures, (a) n=0, (b) n=25, (c) n=50, and (d) n=75. The color bars on the left represent the damage variables. The range of values, blue No damage, red Complete destruction; NT11 is the node number identifier; the red line on the right is the crack propagation path, damage variable. The contour lines. It can be seen that with the number of freeze-thaw cycles... As the damage increases, the damaged area gradually expands from the tip of the pre-set shear band into the interior of the specimen, and the high-damage area ( The area of ) increases significantly. When At this time, the damage is localized within the main shear zone, exhibiting a typical brittle shear failure morphology. At that time, the damage was almost spread throughout the entire shear region, and the specimen changed from a single macroscopic crack propagation to a failure mode dominated by diffuse damage.
[0254] This evolutionary pattern can be explained from the perspective of the coupled damage model of this invention as follows: freeze-thaw damage variables The accumulation of these factors leads to a decrease in the effective stiffness of the material and Attenuation. Reduced stiffness decreases the elastic strain energy generated under the same external load, but the attenuation of fracture energy also weakens the material's ability to resist crack propagation. Under the combined effect of these two factors, damage is no longer concentrated in a single shear band, but is distributed in a more diffuse manner. Macroscopically, this manifests as a decrease in shear strength and a shift in the failure mode from brittle to ductile.
[0255] Figure 5 The load-displacement curves under different freeze-thaw cycles are shown in Table 4. The main results are as follows:
[0256] (1) Unfrozen specimens ( The peak load is 74.56 N; The temperature dropped to 69.05 N, a decrease of 7.39%. The temperature dropped to 61.74 N, a decrease of 19.19%. The temperature dropped to 58.54 N, a decrease of 21.49%. This trend is consistent with the freeze-thaw damage variable. The cumulative effect is consistent with the observation that freeze-thaw degradation continuously weakens the shear strength of the material.
[0257] (2) The peak displacement of the unfrozen specimen was 0.089 mm; The thickness decreased to 0.087 mm. The thickness was reduced to 0.080 mm. The thickness decreased to 0.077 mm. This is due to freeze-thaw damage causing a degradation in material stiffness and a decrease in elastic modulus. The strength decreases, meaning the material reaches its peak strength and fails with minimal deformation.
[0258] (3) The post-peak curve of the unfrozen specimen drops sharply, indicating brittle fracture characteristics; as the number of freeze-thaw cycles increases, the post-peak curve declines more gradually, and the ductility characteristics are enhanced. This change is due to the microcrack network induced by freeze-thaw cycles, which makes the damage distribution inside the material more diffuse and the crack propagation path more tortuous, thereby consuming more fracture energy.
[0259] The above evolutionary pattern indicates that, in the concrete phase field fracture model of this invention, the freeze-thaw damage factor... Directly reducing the effective stiffness of the material leads to a decrease in the stress level under the same external load; simultaneously, the fracture energy caused by freeze-thaw cycles... Attenuation and strength parameter degradation together weaken the material's ability to resist crack propagation. Therefore, although the driving force term in the phase field evolution equation (Equation (40)) Even without explicitly including freeze-thaw damage factors, freeze-thaw degradation still affects fracture behavior through two pathways: stiffness reduction and energy parameter decay. Macroscopically, this manifests as a comprehensive degradation effect of reduced peak load, decreased peak displacement, and enhanced post-peak ductility.
[0260] Table 4. Numerical values of the direct shear model under freeze-thaw damage.
[0261]
[0262] (II) Verification of concrete crack propagation under freeze-thaw cycles.
[0263] To verify the effectiveness of the concrete phase-field fracture model proposed in this invention, which couples freeze-thaw degradation with the Drucker-Prager criterion, a three-dimensional finite element model of a concrete specimen consistent with experimental results was established. Boundary conditions were set as follows: a fixed constraint was applied to the bottom, and a displacement of 0.05 mm was applied to the top. The model was spatially discretized using eight-node hexahedral linear reduced integral elements. Mesh refinement was performed in the expected crack propagation region, ensuring that the element feature size was less than half the phase-field feature length, thus guaranteeing accurate crack path capture.
[0264] As the number of freeze-thaw cycles (n) increases, the high-stress zone gradually shrinks, the peak stress continuously decreases, and the stress gradient tends to flatten. Specifically, this manifests as follows:
[0265] (1) When At that time, the load-bearing area in the middle of the specimen showed a significant high stress concentration zone, with the peak equivalent stress being [value missing]. The stress in the x-direction is The stress in the y direction is The stress in the z-direction is At this point, the freeze-thaw damage is relatively minor, and the material still retains good stress transfer capabilities.
[0266] (2) When At that time, the peak equivalent stress dropped to The stress in the x-direction, y-direction, and z-direction decreased to [values to be inserted here]. , and The significant shrinkage of the high-stress zone indicates that the freeze-thaw induced microcrack network has begun to weaken the stress transfer efficiency of the material, leading to stress redistribution.
[0267] (3) When At that time, the peak equivalent stress further decreased to The stress in the x-direction, y-direction, and z-direction decreased to [values to be inserted here]. , and At this point, the high-stress zone is limited to a local area, the stress distribution is more uniform, and the gradient change tends to be gentle, reflecting that freeze-thaw damage has led to a significant degradation in the overall stiffness of the material and the disruption of the stress transmission path.
[0268] The above evolutionary pattern can be explained by the coupled damage model proposed in this invention: freeze-thaw damage factor The introduction of freeze-thaw cycles reduces the effective stiffness of the material, resulting in a decrease in the stress level under the same external load. Simultaneously, the dispersed microcrack network induced by freeze-thaw cycles weakens the stress transfer capacity to specific regions, ultimately manifesting as a combined deterioration effect of reduced ultimate stress peak, shrinkage in high-stress areas, and homogenization of stress distribution. The numerical simulation results are qualitatively consistent with experimental observations, verifying the effectiveness of the model in describing the stress response of concrete after freeze-thaw degradation.
[0269] Figure 6 The stress-strain curves of concrete specimens under different freeze-thaw cycles are shown in Table 5, which presents the corresponding peak stress and peak strain values. The stress-strain curves of concrete specimens under different freeze-thaw cycles are presented in Table 5. With the increase of , the curve shape exhibits the following evolutionary characteristics:
[0270] (1) Peak stress drops to (decline Peak strain reduced to ; Peak stress drops to (decline Peak strain is ; Peak stress further decreased (decline Peak strain reduced to Meanwhile, the slope of the initial straight segment of the stress-strain curve gradually decreases, indicating that freeze-thaw cycles cause a significant degradation in the material's elastic modulus, and the material reaches its peak strength with smaller deformations.
[0271] (2) The peak curve of the unfrozen specimen drops sharply, showing typical brittle failure characteristics; the peak curve of the frozen specimen drops more gradually, and the ductility is enhanced. This is because the microcrack network induced by freeze-thaw makes the damage distribution more diffuse and the crack propagation path more tortuous, thus consuming more energy during the failure process.
[0272] (3) The above evolutionary pattern shows that the freeze-thaw damage factor The introduction of this effect reduces the effective stiffness of the material, causing it to enter the damage evolution stage under lower stress and strain. Macroscopically, this manifests as a comprehensive deterioration effect, including reduced peak stress, decreased peak strain, stiffness degradation, and enhanced post-peak ductility. The numerical simulation results agree well with the experimental data, verifying the effectiveness of the model in describing the mechanical response of concrete after freeze-thaw degradation.
[0273] Table 5 Numerical values of the model under freeze-thaw damage
[0274]
[0275] This invention addresses the fracture behavior of concrete structures in cold regions under the coupled effects of freeze-thaw cycles and mechanical loads. It proposes a physically consistent multi-scale phase-field fracture model and systematically establishes a unified description method for freeze-thaw degradation, tensile-compressive asymmetric fracture mechanisms, and multi-scale damage mapping. Experimental verification results are as follows:
[0276] (1) This invention utilizes an effective stiffness matrix By coupling freeze-thaw damage with mechanical damage and combining a linear freeze-thaw damage accumulation model, a deep integrated description of environmental degradation and mechanical fracture is achieved. This scheme achieves an effective balance between computational efficiency and physical rationality, and is easy for UMAT user material subroutines to implement.
[0277] (2) Numerical simulation was performed on the specimens by direct shearing. The results showed that the number of freeze-thaw cycles increased with the number of freeze-thaw cycles. Increase, when At that time, the damage was almost widespread throughout the shear region, and the failure mode changed from single macroscopic crack propagation to diffuse damage dominance. The force-displacement curves show that compared to the unthawed specimen (peak load 74.56 N), The peak load decreased to 58.54 N, a reduction of 21.5%; the peak displacement decreased from 0.089 mm to 0.077 mm, a decrease of 13.5%; the post-peak curve changed from a steep drop to a gentler curve, indicating enhanced ductility. This evolutionary pattern verifies the continuous weakening effect of freeze-thaw damage accumulation on the shear resistance of the material.
[0278] (3) The model was independently verified through concrete freeze-thaw cycle tests and ultrasonic testing. The stress-strain curves showed that the peak stress of the unfrozen specimen was... The peak strain is After 75 freeze-thaw cycles, the peak stress decreased to (Decreased by 23.3%), peak strain decreased to (Decrease of 14.7%). Ultrasonic testing showed that the signal amplitude was 1.95mV before freeze-thaw cycles, decreased by 80.87% after 50 freeze-thaw cycles, and by 84.10% after 75 freeze-thaw cycles, with a significant difference in the time-domain waveform compared to the unfrozen state. These verifications collectively confirm the effectiveness and reliability of the model in describing the mechanical response of concrete after freeze-thaw degradation.
[0279] In summary, the phase-field fracture model proposed in this invention, which couples freeze-thaw degradation with the Drucker-Prager criterion, provides a physically based theoretical description and numerical tool for the fracture behavior of concrete structures in cold regions under the coupled effects of freeze-thaw cycles and loads. Future development can further extend this model to predict the performance of complex three-dimensional structures, multiaxial stress states, and long-term service capabilities, and provide theoretical support for freeze-thaw durability design and life assessment.
[0280] The above embodiments are merely preferred examples of the present invention and are not intended to limit the present invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the present invention. Therefore, any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention should fall within the protection scope of the present invention.
Claims
1. A method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase-field fracture model, characterized in that, Includes the following steps: Step S1: Based on the asymptotic homogenization theory, establish the micro-macro damage mapping of concrete and derive the expression for the phase field characteristic length related to the microstructure. Step S2: Introduce the Drucker-Prager criterion into the phase field fracture framework to decompose the total strain energy density into a damage-driven part and a protected part to describe the tensile-compressive asymmetric fracture behavior of concrete. Step S3: Based on the strain equivalence assumption, establish the series coupling relationship between freeze-thaw damage and mechanical damage. By introducing a freeze-thaw damage factor into the phase field evolution equation to correct the strain energy release rate, the dynamic influence of freeze-thaw degradation on the mechanical fracture process is simulated. Step S4: Based on steps S1-S3, construct the control equations of the coupled damage model, and perform numerical discretization and solution using the finite element method to predict the entire process of crack initiation, propagation and failure of concrete structures under freeze-thaw conditions.
2. The method according to claim 1, wherein, In step S1, the phase field characteristic length is represented as a strain rate-related function, as shown in the following expression: ； wherein is a phase-field characteristic length, is a quasi-static characteristic length, is a current strain rate; is a characteristic strain rate, which is a material intrinsic parameter; is a scaling parameter, and ; is a triaxial stress coefficient, and .
3. The method according to claim 2, wherein, In step S2, the total strain energy density is decomposed into the damage-driving part and the protected part: ； wherein is the total strain energy density; is the damage driven strain energy density; is the protected strain energy density; is the strain tensor; is the damage variable, , denotes no damage, denotes complete damage; is the degradation function.
4. The method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model according to claim 3, characterized in that, In step S2, the damage-driven strain energy density is constructed based on the Drucker-Prager criterion. The expression is: ； In the formula, It is the elastic modulus; The main part of the brackets in Macaulay; For the stress tensor in the undamaged state, , The elastic tensor is in a lossless state; The second deviatoric stress invariant is based on the stress tensor of the lossless state. This represents the positive part of the first stress invariant based on the stress tensor of the undamaged state; These are parameters related to the internal friction angle.
5. The method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model according to claim 4, characterized in that, In step S2, the protected strain energy density is determined based on the Drucker-Prager criterion. The expression is: ； In the formula, For the protected strain energy density, Bulk modulus These are parameters related to the internal friction angle. Shear modulus It is the elastic modulus; The negative part of the brackets in Macaulay. The main part of the brackets in Macaulay; As the first strain invariant, This is the second partial strain invariant.
6. The method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model according to claim 5, characterized in that, In step S3, stiffness reduction is used based on the strain equivalence assumption to achieve coupling between freeze-thaw damage and mechanical damage. The effective stiffness matrix is expressed as: ； In the formula, This represents the effective stiffness matrix after freeze-thaw and mechanical coupling damage. The elastic tensor is in a lossless state; Let be the mechanical damage degradation function. It is a freeze-thaw damage factor. For freeze-thaw damage variables, For damage variables.
7. The method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model according to claim 6, characterized in that, The coupled damage model in step S4 consists of the following governing equations: Momentum balance equation: ； In the formula, For the total spatial gradient operator, For stress tensor, The density of force per unit volume; Stress-strain-damage relationship: ； In the formula, For stress tensor, For the lossless elastic tensor, For freeze-thaw damage variables, As a damage variable, For strain tensor; Geometric equations: ； In the formula, For strain tensor, Let be the displacement gradient tensor. For displacement field, This is a transpose operation; Phase field evolution equation: ； In the formula, To be related to the number of freeze-thaw cycles The relevant fracture energy, As a damage variable, The characteristic length of the phase field; For the Laplace operator, For historical variables.
8. The method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model according to claim 7, characterized in that, In step S4, the displacement field and phase field are solved alternately using the separation iteration method. Specifically, in the current increment step, the phase field is fixed first, and the displacement field is solved based on the momentum balance equation; then the displacement field is fixed, and the phase field is solved based on the phase field evolution equation; the above alternating process is repeated until the displacement residual and the phase field residual both meet the preset convergence tolerance; in each increment step, the linearized discrete equation system is solved using Newton-Raphson iteration.
9. The method for constructing a multi-scale coupled freeze-thaw-DP criterion concrete phase field fracture model according to claim 8, characterized in that, In step S4, the numerical solution of the coupled damage model is realized by writing the UMAT user material subroutine of Abaqus software. The subroutine performs the following operations in each increment step: read the current strain increment and freeze-thaw damage state variables, update the stress and Jacobian matrices according to the coupled constitutive model, call the phase field solution submodule to update the historical variables and phase field variables, and store the state variables for use in the next increment step.