Clay cyclic stability and weakening constitutive model and numerical implementation method

Abstract

The application discloses a kind of clay cycle stability and weakening constitutive model and numerical implementation method, belong to clay characteristic technical field, including the following steps 1, define isotropic hardening ellipse boundary surface, 2, propose mapping center evolution mechanism of servo hardening, calculate image stress by radial mapping rule, 3, construct multi-stage interpolation equation to calculate plastic modulus, 4, based on the stress-strain increment relationship model of incremental constitutive theory under the framework of elastoplasticity, 5, give practical model parameter calibration method, 6, construct a kind of high-precision explicit integration algorithm suitable for the numerical implementation of the proposed boundary surface constitutive model, the application is described by the above model and method Complex mechanical behavior of clay under cyclic loading with fewer parameters, overcome the problem of complex theory and too many parameters of traditional model, realize the dynamic switching of cycle stability and cycle weakening mode, and significantly improve the accuracy and convergence of numerical calculation.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of clay characteristics, and particularly relates to a clay cyclic stability and weakening constitutive model and a numerical implementation method. BACKGROUND

[0002] The rapid development of China's offshore wind power is inseparable from the unique seabed geological conditions - the widespread distribution of clay layers under vast sea areas is not only the key medium supporting the wind turbine foundation, but also the core challenge restricting engineering safety due to its complex mechanical properties. As a typical marine sediment, the marine clay around the pile foundation has low strength, high compressibility, and significant time dependence. The wind turbine needs to continuously withstand the cyclic loads of wind, wave, and flow during long-term operation, and its mechanical response presents complex nonlinear characteristics: cyclic loading can cause soil stiffness degradation, pore water pressure accumulation, and irreversible plastic strain, which in turn leads to foundation settlement and may even cause structural overturning. This characteristic is in sharp contradiction with the trend of China's offshore wind power expanding into the deep sea - the foundation structure needs to withstand various complex environmental loads including seismic loads, and the research demand for the dynamic response of offshore wind turbine clay foundation under cyclic loading is becoming increasingly urgent.

[0003] Marine clay has relatively complex mechanical properties, and the cumulative behavior of clay under different amplitude cyclic loads can be divided into two modes: cyclic stability and cyclic weakening. Most existing constitutive models (such as the Mróz multi-surface model and the CASM-c model) can only describe cyclic stability or cyclic weakening individually and cannot capture the dynamic switching of both modes. This single-mode representation capability does not match the complex soil response under combined loads (such as random wind and wave coupling) in the actual operation of offshore wind turbine foundations, resulting in insufficient accuracy in long-term deformation prediction.

[0004] Secondly, even if some existing boundary surface models can describe the cyclic mechanical properties of clay and couple the cyclic double-mode behavior, their high theoretical complexity and excessive parameter requirements will significantly increase the numerical calculation cost, making it difficult to embed them into commercial finite element software for large-scale engineering simulation. This defect particularly restricts their real-time design and safety assessment applications in deep-sea wind power projects. Based on this, the present application proposes a clay cyclic stability and weakening constitutive model and a numerical implementation method. SUMMARY

[0005] The application aims to provide a clay cyclic stability and weakening constitutive model and a numerical implementation method, which can represent behaviors of soil under cyclic loading, such as hysteresis effect and stiffness nonlinear attenuation, by constructing a dynamic evolution boundary surface with fewer parameters, and can overcome defects of high theoretical complexity, low calculation efficiency and too many required parameters, and can realize dynamic switching of two modes of clay cyclic stability and cyclic weakening, and can solve the stability-efficiency problem of traditional explicit integration in strong nonlinear constitutive calculation by controlling relative error to dynamically adjust the substep size.

[0006] To achieve the above-mentioned purpose, the application provides a clay cyclic stability and weakening constitutive model and a numerical implementation method, comprising the following steps:

[0007] S1, defining an isotropic hardening elliptical boundary surface;

[0008] S2, after defining the boundary surface, combining the evolution mechanism of the mapping center as specified in the stress uniqueness theory, and calculating the image stress based on the radial mapping rule;

[0009] S3, after completing the image stress calculation, constructing a multi-stage interpolation equation to calculate the plastic modulus;

[0010] S4, after determining the plastic modulus, constructing a stress-strain increment relationship model based on the incremental constitutive theory under the elastic-plastic framework;

[0011] S5, after constructing the model, determining the basic parameters through conventional indoor test analysis, and determining the interpolation parameters through fitting test, to complete the model parameter calibration;

[0012] S6, constructing a high-precision explicit integration algorithm suitable for numerical implementation of the proposed boundary surface constitutive model.

[0013] Preferably, the method for defining the isotropic hardening elliptical boundary surface in S1 is as follows:

[0014] S11, calculating the critical state stress ratio M of the elliptical boundary surface, the method being as follows:

[0015]

[0016] wherein represents the internal friction angle;

[0017] S12, constructing the elliptical boundary surface, and the boundary surface equation being as follows:

[0018]

[0019] wherein, wherein p and q respectively represent the average principal stress and shear stress, (p, q) represents the true stress, and p cis the intercept of the boundary surface on the mean principal stress axis, which represents the size of the boundary surface;

[0020] S13, the evolution equation of p c is defined as a hardening parameter, and the evolution equation of p c The evolution equation of the boundary surface is realized by defining p

[0021]

[0022] where e0 is the initial void ratio, λ and κ represent the slopes of the normal consolidation line and the unloading rebound line in the e-lnp space, respectively, represents the plastic strain increment.

[0023] Preferably, the shape and evolution law of the plastic potential surface are also determined by the boundary surface equation in S12 and the isotropic hardening process in S13.

[0024] Preferably, the calculation process of the image stress in S2 is as follows:

[0025] S21, fix the initial mapping center (p pc , q pc ) at the origin, and the line connecting the mapping center and the real stress (p, q) is the intersection point of the boundary surface, which is the image stress point

[0026] S22, according to the uniqueness theory of image stress, that is, when the boundary surface is kinematic hardening, the mapping center and its relative position should remain unchanged, and when stress reversal occurs, the mapping center should be updated to the current real stress point, and the evolution equation of the mapping center is defined as follows:

[0027]

[0028] S23, the image stress is calculated based on the radial mapping rule from the mapping center and the real stress:

[0029]

[0030] where the parameter b ∈ [1, ∞], and the norm quadratic equation about the parameter b is constructed by solving the evolution equation of the mapping center and the boundary surface control equation.

[0031] Preferably, the process of constructing a multi-stage interpolation equation to calculate the plastic modulus in S3 is as follows:

[0032] S31, the plastic modulus of the image stress point is obtained from the boundary surface equation and according to the consistency condition of the boundary surface plastic theory process as follows:

[0033]

[0034] where the loading factor L is given by K is calculated as p the plastic modulus of the true stress point, the operator "<>" is noted as Macaulay brackets, when L is positive <l>= L, and vice versa <l>= 0;

[0035] S32, determining the plastic modulus K of the real stress point S33, determining the plastic modulus K of the real stress point p S34, constructing the interpolation function of K p The interpolation function of K is constructed as follows:

[0036]

[0037] Where (b-1) is used to satisfy the general law of plastic modulus. If b=1, the image stress point coincides with the current stress point, and K=0. If b=∞, stress reversal occurs, and K=∞. The term reflects the influence of the void ratio on the plastic modulus, and h represents a non-negative model parameter, which is determined by the following formula:

[0038]

[0039] Where d represents a damage state variable, and c d represents a model parameter. When the model parameter is positive, it simulates cyclic softening, and when it is negative, it simulates cyclic stability.

[0040] Preferably, the process of constructing the stress-strain increment relationship model in S4 is as follows:

[0041]

[0042] Where dε v represents the volumetric strain increment, dε d represents the shear strain increment, and D ep represents the elastoplastic stiffness matrix, which is calculated as follows:

[0043]

[0044] Where D e represents the elastic stiffness matrix, which is determined by Hooke's law.

[0045] Preferably, the process of calibrating the model parameters in S5 is as follows:

[0046] S51, determining five basic parameters: compression parameter λ, resilience parameter κ, Poisson's ratio v, critical state stress ratio M, and void ratio N corresponding to an isotropic consolidation pressure of 1 kPa through conventional laboratory tests;

[0047] S52, determining two interpolation parameters h0 and c d , the process is as follows:

[0048] S521, the initial cyclic stage accumulates a small amount of plastic strain, and c d has a weak influence on the stress-strain relationship, and h0 and c d the influence effect decoupling processing;

[0049] S522, determining h0 by optimal fitting of the initial stress-strain hysteresis loop;

[0050] S523, determining c d according to whether the soil is in a cyclic stable state or a cyclic weakening state, and then fixing other parameters to determine c d .

[0051] Preferably, the process of establishing the constitutive model numerical implementation algorithm in S6 is as follows:

[0052] S61, setting T = 0 and ΔT = 0, and determining the initial stress state and internal variable of the current increment step through the initial conditions;

[0053] S62, calculating the loading factor L of the current stress state, if L ≤ 0, updating the mapping center to the current stress point, otherwise keeping the relative position of the mapping center and the boundary surface unchanged;

[0054] S63, determining the sub-step strain increment {dε ss} = ΔT {dε}, obtaining the stiffness matrix [D ep {σ}] from the initial stress state, and obtaining the first trial stress increment {dσ1} and the remaining internal variables, the first trial stress calculation formula being as follows:

[0055] {dσ1} = [D ep {σ}] {dε ss};

[0056] S64, constructing the stiffness matrix [D ep {σ+dσ1}] from the stress state {σ+dσ} obtained by the first trial and the updated internal variables, calculating the second trial stress increment {dσ2} and the remaining internal variables, the second trial stress calculation formula being as follows:

[0057] {dσ2} = [D ep {σ+dσ1}] {dε ss};

[0058] S65, calculating the average stress increment {dσ} of the two stress increments:

[0059]

[0060] S66, calculating the relative error, the process being as follows:

[0061]

[0062] Wherein RE represents relative error, SSTOL represents error control value, if RE > SSTOL, the substep length needs to be reduced and returns to S63 to recalculate, and the new substep length is calculated as follows:

[0063] Delta T new = 0.8 [SSTOL / RE] 1 / 2 Delta T

[0064] If RE <= SSTOL, the stress state {sigma} = {sigma} + {d sigma} and the internal variable are updated.

[0065] S67, let T = T + Delta T, the same method as S66 is used to determine the new substep length Delta T new , if T + Delta T new <= 0, let Delta T new = 1-T, return to S63, when T = 1, the increment step calculation is completed, the stiffness matrix is updated with the end stress state and is passed to the analysis software, and the increment step calculation is completed.

[0066] Therefore, the clay cyclic stability and weakening constitutive model and numerical implementation method adopt the above structure, provide a more actual, simple and suitable for engineering construction of clay constitutive model, the boundary surface model can not only accurately describe the mechanical properties of clay, coupled double mode cycle behavior, but also has clear physical meaning and lower parameter calibration cost, at the same time, based on explicit integration method for numerical implementation, by controlling the relative error dynamic adjustment substep length, the balance of numerical simulation precision and efficiency is realized, the calculation performance and result reliability are significantly improved.

[0067] The technical scheme of the present application will be further described in detail below by means of the drawings and examples. BRIEF DESCRIPTION OF DRAWINGS

[0068] Figure 1 It is a boundary surface function and mapping rule schematic diagram of the clay cyclic stability and weakening constitutive model and numerical implementation method of the present application.

[0069] Figure 2 It is a parameter sensitivity analysis schematic diagram of different interpolation parameters h0 of the clay cyclic stability and weakening constitutive model and numerical implementation method of the present application.

[0070] Figure 3 It is a parameter sensitivity analysis schematic diagram of different interpolation parameters c d of the clay cyclic stability and weakening constitutive model and numerical implementation method of the present application.

[0071] Figure 4 It is a parameter sensitivity analysis schematic diagram of different cyclic amplitude q cyc The stress-strain curve simulated by the triaxial undrained shear experiment constitutive model of the lower Georgiakaolin clay is compared with the triaxial undrained shear experiment result and simulation result.

[0072] Figure 5 The different amplitude r of the clay cyclic stability and weakening constitutive model and numerical implementation method of the application is shown in the following table: cyc / S uc The stress-strain curve simulated by the triaxial undrained shear experiment constitutive model of the lower Cloverdale clay is compared with the triaxial undrained shear experiment result and simulation result.

[0073] Figure 6 The numerical implementation method flow chart of the clay cyclic stability and weakening constitutive model and numerical implementation method of the application is shown in the following table:

[0074] Figure 7 The comparison chart of the stress-strain curve of the clay cyclic stability and weakening constitutive model and numerical implementation method of the application in the triaxial drainage unidirectional cycle and bidirectional cycle experiment of Georgiakaolin clay is shown in the following table:

[0075] Figure 8 The cycle number-accumulative vertical strain curve of the clay cyclic stability and weakening constitutive model and numerical implementation method of the application in the process of applying cyclic load on the strip foundation is shown in the following table:

[0076] Figure 9 The local mesh densification of the soil around the pile and the pile-soil contact area and the pore water pressure observation point arrangement schematic diagram of the numerical model for the finite element analysis of the dynamic response of the fan seismic load are shown in the following table:

[0077] Figure 10 a in the formula is the time history evolution curve of the fan single pile settlement, Figure 10 b in the formula is the time history evolution curve of the rotation angle at the sea bed mud surface of the fan single pile;

[0078] Figure 11 a in the formula is the pore water pressure change curve of the soil around the pile under the action of the seismic load, Figure 11 b in the formula is the pore water pressure change curve of the soil around the pile under the action of the seismic load. Figure 11 b in the formula is the pore water pressure change curve of the soil around the pile under the action of the seismic load. DETAILED DESCRIPTION

[0079] To make the purpose, technical scheme and advantages of the embodiments of the application clearer, the technical scheme in the embodiments of the application will be described clearly and completely below with reference to the drawings in the embodiments of the application. Obviously, the described embodiments are part of the embodiments of the application, rather than all the embodiments. The components of the embodiments of the application described and shown in the drawings can be arranged and designed in various different configurations.

[0080] The following detailed description of the embodiments of the application provided in the accompanying drawings is not intended to limit the scope of the application claimed, but merely represents the selected embodiments of the application. Based on the embodiments in the application, all other embodiments obtained by those of ordinary skill in the art without creative labor fall within the scope of protection of the application.

[0081] It should be noted that similar reference numbers and letters represent similar items in the following drawings, and therefore, once an item is defined in one drawing, it need not be further defined and explained in subsequent drawings.

[0082] In the description of the application, it should be noted that the terms "upper", "lower", "inner", "outer", etc. indicate the orientation or positional relationship based on the orientation or positional relationship shown in the drawings, or the orientation or positional relationship in which the product of the application is usually placed, and are only for the convenience of describing the application and simplifying the description, and therefore cannot be understood as indicating or implying that the indicated device or element must have a particular orientation, be constructed and operated in a particular orientation, and therefore cannot be understood as limiting the application.

[0083] In the description of the application, it should also be noted that, unless otherwise explicitly specified and limited, the terms "provided", "mounted", "connected" should be broadly understood, for example, can be fixedly connected, or can be detachably connected, or integrally connected; can be mechanically connected, or can be electrically connected; can be directly connected, or can be indirectly connected through an intermediate medium, or can be the internal communication of two elements. For those of ordinary skill in the art, the specific meaning of the above terms in the application can be understood according to the specific circumstances.

[0084] Some embodiments of the application will be described in detail below with reference to the accompanying drawings. The following embodiments and features in the embodiments can be combined with each other without conflict.

[0085] Embodiment 1

[0086] As shown in the accompanying drawings, the clay cyclic stability and weakening constitutive model and numerical implementation method of the application comprises the following steps: Figures 1-3

[0087] S1, define an isotropic hardening elliptical boundary surface;

[0088] S11, calculate the critical state stress ratio M of the elliptical boundary surface, the method is as follows:

[0089]

[0090] wherein represents the internal friction angle;

[0091] ​S12, constructing the elliptical boundary surface, the boundary surface equation is as follows:

[0092]

[0093] Where p and q represent the average principal stress and shear stress respectively, (p, q) represents the true stress, p c is the intercept of the boundary surface on the average principal stress axis, representing the boundary surface size;

[0094] S13, taking p c as the hardening parameter, the isotropic hardening of the boundary surface is realized by defining the evolution equation of p c , and the evolution equation is:

[0095]

[0096] Where e0 is the initial void ratio, λ and κ represent the slope of the normal consolidation line and the unloading rebound line in the e-lnp space respectively, represents the plastic strain increment.

[0097] Preferably, the shape and evolution law of the plastic potential surface are also determined by the boundary surface equation in S12 and the isotropic hardening process in S13.

[0098] S2, after defining the boundary surface, the evolution mechanism of the mapping center is determined according to the image stress uniqueness theory, and the image stress is calculated based on the radial mapping rule;

[0099] S21, fixing the initial mapping center (p pc , q pc ) at the origin, the intersection of the line connecting the mapping center and the true stress (p, q) and the boundary surface is the image stress point

[0100] S22, according to the image stress uniqueness theory, that is, when the boundary surface is dynamically reinforced, the relative position of the mapping center should be kept unchanged, and when the stress is reversed, the mapping center should be updated to the current true stress point, the evolution equation of the mapping center is as follows:

[0101]

[0102] S23, the image stress is calculated based on the radial mapping rule from the mapping center and the true stress:

[0103]

[0104] Where parameter b ∈ [1, ∞], the norm quadratic equation about parameter b is constructed by solving the mapping center evolution equation and the boundary surface control equation.

[0105] S3, after the stress calculation, construct the multi-stage interpolation equation to calculate the plastic modulus;

[0106] S31, get the plastic modulus of the stress point from the boundary surface equation and according to the consistency condition of the boundary surface plastic theory The process is as follows:

[0107]

[0108]

[0109] Wherein the loading factor L is calculated by K p The plastic modulus of the real stress point, the operator "<>" is recorded as Macaulay bracket, when L is positive <l>= L, and vice versa <l>= 0.

[0110] S32, determining the plastic modulus K of the real stress point according to the plastic modulus of the image stress point S32, determining the plastic modulus K of the real stress point according to the plastic modulus of the image stress point p , constructing an interpolation function for calculating K p is as follows:

[0111]

[0112] where (b-1) is used to satisfy the general law of the plastic modulus, if b=1, the image stress point coincides with the current stress point, at this time if b=∞, stress reversal occurs, at this time The term reflects the influence of the void ratio on the plastic modulus, h represents a non-negative model parameter, which is determined by the following formula:

[0113]

[0114] where d represents a damage state variable, c d represents a model parameter, when the model parameter is positive, the cyclic softening is simulated, and when the model parameter is negative, the cyclic stability is simulated.

[0115] S4, after obtaining the plastic modulus, a stress-strain increment relationship model is constructed based on the incremental constitutive theory under the elastic-plastic framework;

[0116]

[0117] where dε v represents the volume strain increment, dε d represents the shear strain increment, D ep represents the elastic-plastic stiffness matrix, and the calculation formula is as follows:

[0118]

[0119] where D e represents the elastic stiffness matrix, which is determined by Hooke's law.

[0120] S5, after constructing the model, a practical interpolation parameter calibration method is given to determine the model parameters;

[0121] S51, five basic parameters are determined through conventional indoor tests: compression parameter λ, rebound parameter κ, Poisson's ratio v, critical state stress ratio M, and void ratio N corresponding to an isotropic consolidation pressure of 1 kPa;

[0122] S52, two interpolation parameters h0 and c d are determined, and the process is as follows:

[0123] S521, the cumulative plastic strain is minimal in the initial cyclic stage, and c d The influence of h0 and c is decoupled d ;

[0124] S522, determining h0 by optimal fitting of the initial stress-strain hysteresis loop;

[0125] S523, determining the positive and negative of c according to whether the soil is in cyclic stability (hysteresis loop contraction) or cyclic weakening state (hysteresis loop expansion), then fixing other parameters, and determining c by multiple trial calculations. d . d

[0126] S6, constructing a high-precision explicit integration algorithm suitable for numerical implementation of the proposed boundary surface constitutive model:

[0127] S61, let T = 0, ΔT = 0, and determine the initial stress state and internal variable of the current increment step by the initial condition;

[0128] S62, calculate the current stress state loading factor L, if L ≤ 0, update the mapping center to the current stress point, otherwise keep the relative position of the mapping center unchanged;

[0129] S63, determine the sub-step strain increment {dε ss} = ΔT {dε}, obtain the stiffness matrix [D ep {σ}] from the initial stress state, and obtain the first trial stress increment {dσ1} and the remaining internal variables. The first trial stress calculation formula is as follows:

[0130] {dσ1} = [D ep {σ}]{dε ss};

[0131] S64, construct the stiffness matrix [D ep {σ+dσ1}] from the stress state {σ+dσ} obtained by the first trial and the updated internal variables, calculate the next trial stress increment {dσ2} and the remaining internal variables. The second trial stress calculation formula is as follows:

[0132] {dσ2} = [D ep {σ+dσ1}]{dε ss};

[0133] S65, calculate the average stress increment {dσ} of the two stress increments:

[0134]

[0135] S66, calculate the relative error, the process is as follows:

[0136] ​

[0137] where RE represents relative error, SSTOL represents error control value, if RE > SSTOL, the substep needs to be reduced and S63 is returned to recalculate, and the new substep is calculated as follows:

[0138] ΔT new = 0.8 [SSTOL / RE] 1 / 2 ΔT;

[0139] If RE ≤ SSTOL, the stress state {σ} = {σ} + {dσ} and the internal variable are updated.

[0140] S67, T = T + ΔT, the same method as S66 is used to determine the new substep ΔT new , if T + ΔT new ≤ 0, then ΔT new = 1-T, return to S63, when T = 1, the increment step calculation is completed, the stiffness matrix is updated with the end stress state and is passed to the analysis software, and the increment step calculation is completed.

[0141] Example 2

[0142] The process (S1-S4) of establishing the model in the embodiment 1 is extended to a general stress state, and the process is as follows:

[0143] Step 1, the boundary surface equation under the general stress state is established:

[0144]

[0145] where the superscript "-" indicates that the variable is located on the boundary surface, σ and s = σ-pI represent the stress tensor and the partial component of the stress tensor respectively, I represents the unit tensor, p = tr(σ) / 3 represents the average principal stress tensor, r = s / p is recorded as the stress ratio tensor, and the determination of the parameter M is different from that under the triaxial stress state, and the specific determination is as follows:

[0146]

[0147] where c = M c / M e , M c and M e respectively represent the critical state stress ratio under triaxial compression and tension, and θ represents the stress Lode angle;

[0148] The calculation formula of the stress loading direction and the plastic flow direction under the general stress state is as follows:

[0149]

[0150] Step 2, the evolution law of the mapping center under general stress state, like stress calculation and triaxial stress state consistent, the calculation process as follows:

[0151]

[0152]

[0153] The above calculation process and the boundary surface equation are combined and decomposed into only p and s two parts to simplify the solving process of parameter b, the specific formula is as follows:

[0154]

[0155] M 2 (p pc +b(p-p pc )) 2 -M 2 p c (p pc +b(p-p pc ))=A4b 2 +A5b+A6;

[0156]

[0157] A2=3(s-s pc ):s pc ;

[0158]

[0159] A4=M 2 (p-p pc ) 2 ;

[0160] A5=M 2 (2p pc -p c ) 2 (p-p pc );

[0161] A6=M 2 p pc (p pc -p c );

[0162] Where A i (i=1,2,...,6) is the calculation process parameter. The quadratic equation of parameter b can be obtained by combining the above formula:

[0163] (A1+A4)b 2 +(A2+A5)b+(A3+A6)=0;

[0164] Step 3, the calculation process of the plastic modulus of true stress and stress-like stress under general stress state is as follows:

[0165]

[0166]

[0167] Step 4, the stress-strain increment relationship model under general stress state is as follows:

[0168]

[0169] Example 3

[0170] As Figures 4-8 , the calculation program is written according to the model constructed in Example 2 above, and Georgia kaolin clay and Cloverdale clay are selected as experimental sample soils to carry out unidirectional and bidirectional cyclic loading tests under different cyclic amplitudes, so as to verify the rationality of the model construction. Further, strip foundation cyclic loading tests are carried out to further verify the performance of the model in solving boundary value problems.

[0171] Firstly, based on the written calculation program, step S5 is carried out to complete the calibration of basic parameters and interpolation parameters, and each parameter is counted into Table 1.

[0172] Table 1 Parameter statistics table of Example 3

[0173]

[0174] A series of undrained cyclic triaxial tests are carried out on Georgia kaolin clay. Georgia kaolin clay is composed of 62% clay and 38% silt, and its plasticity index is 20%. All experimental samples are first isotropically consolidated under a confining pressure of 345 kPa, and then subjected to cyclic amplitudes q cyc = 136 kPa and 140.7 kPa, respectively, and the results are shown in Figure 4 . Figure 4 The experimental results of a and c in Figure 4 correspond to the simulation results of b and d, respectively, and the degree of agreement is good, so the model can well describe the stress-strain hysteresis loop, cyclic stiffness weakening and stress-strain relationship under different cyclic amplitudes.

[0175] A series of experiments are carried out on Cloverdale clay samples, and the experimental results are used to test the performance of the model. The main components of the clay are 49% clay and 45% silt, and its plasticity index is 24%. The consolidation pressure of the sample is 200 kPa, and the undrained strength S uc = 56 kPa is first determined by monotonic loading. The amplitude r cyc / S uc The model's performance was verified by cyclic triaxial undrained tests at two levels: 0.75 and 0.70. (For example...) Figure 5 As shown, Figure 5 The results of experiment a in the middle and Figure 5 The simulation results of b in the middle Figure 5 The results of experiment c in the middle and Figure 5 The simulation results in the d-axis also show a high degree of agreement, further demonstrating that the model can fully reflect the complex mechanical behavior of clay.

[0176] Based on the sub-incremental step explicit integration method proposed in step S6, a computational subroutine was designed and numerically implemented on a finite element analysis platform. The algorithm flow of the subroutine is detailed below. Figure 6 The developed finite element method was used to simulate a drained triaxial test under cyclic loading to verify the correctness of the subroutine. The error control value SSTOL was set to 1e-6 based on experience. The initial conditions and model parameters were consistent with the previous Georgiakaolin clay experiment. A unidirectional cyclic (q) test was conducted. cyc =50kPa) and bidirectional circulation (q cyc Numerical simulation (e.g., 80 kPa). Figure 7 As shown in Figures a and b, the simulation results from Abaqus and the prediction results from the constitutive model are in good agreement under both loop modes, thus verifying the accuracy of the algorithm design and subroutine writing.

[0177] The constitutive model of this invention was used to test strip foundations under cyclic loading, further verifying the model's performance in solving boundary value problems. The bottom surface of the soil was fixed, the top surface was set as a free displacement surface, and the sides were constrained by horizontal displacement. All boundary surfaces except the top surface were considered undrained interfaces. The soil was discretized using C3D8P elements, resulting in 13266 discretized elements. The strip foundation was modeled as a rigid body connected to the soil surface, with a width of 7.61 cm, a length of 22.9 cm, and a height of 3.81 cm; the model box dimensions were 22.9 cm wide, 91.5 cm long, and 60.7 cm high. Due to the symmetry of the problem, only one-quarter of the model was used for numerical analysis. The soil mechanical parameters were consistent with those of the aforementioned Georgiakaolin clay (Table 1).

[0178] First, a vertical displacement load is applied at the reference point. Because the top boundary can drain water, the pore pressure accumulates below the strip foundation as the number of cycles increases and gradually decreases towards the soil surface. The displacement accumulation rate gradually decreases, and the soil tends to reach a cyclically stable state. Figure 8 The cumulative vertical displacements of the experimental and simulation results were compared. The simulation results were largely consistent with the experimental data, and the numerical simulation results fully demonstrated the practicality of the model in solving boundary value problems.

[0179] Example 4

[0180] As Figures 9-11 shown, the boundary surface constitutive model proposed by the application aims to accurately describe the complex mechanical behavior of marine clay, providing theoretical support for offshore wind power support structure design. Considering the expansion of offshore wind power development to seismic active areas, the safety of wind turbine foundations under cyclic seismic loads faces severe challenges, so this embodiment carries out finite element simulation of the dynamic response of offshore wind turbines under seismic loads to verify the engineering applicability of the model.

[0181] Abaqus is used to establish a three-dimensional finite element model of offshore wind turbines, and the geometric parameters of each part are shown in Table 2. The soil body base boundary is set to Y = Z = 0, and the Y = 0 constraint is applied to the four vertical boundaries. To simulate the seismic response of an infinite horizontal layer, equal displacement constraints are applied to the two vertical boundaries parallel to the YOZ plane. The total model is discretized into 30016 C3D8 elements, and the local mesh is refined in the soil around the pile and the pile-soil contact area (as shown in Figure 9 ), the normal direction of the pile-soil contact is set to hard contact, and the tangential direction uses the Coulomb friction model (friction coefficient 0.3). The tower and single pile are set as steel structures, with an elastic modulus of 210 GPa, a Poisson's ratio of 0.3, and a density of 8500 kg / m 3 . The soil density is 1860 kg / m 3 , and the mechanical properties are considered under undrained conditions, using the aforementioned constitutive model with parameters consistent with the Georgia kaolin clay in Table 1. Two working conditions are set: (a) cyclic stability (c d =-10); (b) cyclic weakening (c d =10). The hydrodynamic pressure is simulated by the added mass method. The seismic input is the Northern California earthquake wave, with a peak acceleration amplitude of 0.5g and a horizontal acceleration applied to the bottom of the foundation in the X direction.

[0182] Table 2 Numerical simulation parameter statistics

[0183]

[0184] Figure 10 The time evolution of the single pile settlement at the seabed mud surface is shown in a. The settlement rate shows a nonlinear characteristic of first increasing and then decreasing, which is consistent with the reported cyclic weakening-stability transition law. Under initial loading, the single pile settlement under the two working conditions develops almost synchronously; with the increase of load duration, the cyclic weakening effect leads to the attenuation of soil stiffness, eventually resulting in the accumulation of larger displacement. Figure 10 The b in the formula reflects the rotation angle of the single pile of the fan at the seabed mud surface, and the time history evolution law when it is in different cyclic behaviors: with the increase of the cycle number, the rotation angle presents a development mode of first rapid rise and then gradual stabilization, and the residual rotation angle is close to zero, which is widely observed in centrifuge tests and numerical simulation; compared with the cyclic steady state, the cyclic weakening weakens the acceleration amplification effect, resulting in a significant decrease in the acceleration response of the seabed mud surface, and then affects the amplitude of the rotation angle of the single pile.

[0185] Considering that the constitutive model is constructed based on the effective stress framework, the dynamic evolution process of the pore water pressure of the soil around the pile can be simulated synchronously, combined with Figure 11 It can be seen that the pore water pressure presents a fluctuating growth characteristic with the loading history, and the closer to the seabed mud surface area, the lower the cumulative amount is; compared with the cyclic steady state, the cyclic weakening will lead to a higher cumulative amount of pore water pressure and a more significant fluctuation amplitude.

[0186] Therefore, the clay cyclic steady and weakening constitutive model and numerical implementation method adopt the above structure, which can represent various complex mechanical behaviors of marine clay under cyclic loading by constructing a dynamically evolving boundary surface with fewer parameters, overcoming the defects of high theoretical complexity and excessive required parameters of traditional constitutive models, and realizing the dynamic switching of the two modes of clay cyclic stability and cyclic weakening; the corresponding numerical implementation method is based on the explicit integration method and adjusts the substep length according to the relative error, realizes the precision-efficiency self-balancing, significantly improves the calculation accuracy and convergence, and provides a numerical implementation scheme with stronger robustness for large-scale geotechnical engineering boundary value problem simulation.

[0187] Finally, it should be noted that: the above examples are only used to illustrate the technical solutions of the present application and not to limit them, although the present application has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that: it can still modify or equivalently replace the technical solutions of the present application, and these modifications or equivalent replacements cannot make the modified technical solutions deviate from the spirit and scope of the technical solutions of the present application.< / l> < / l> < / l> < / l>

Claims

1. A clay cyclic stabilization and weakening constitutive model and its numerical implementation method, characterized in that, Includes the following steps: S1. Define an isotropically hardened elliptical boundary surface; S11. Calculate the stress ratio of the elliptical boundary surface under boundary conditions. The method is as follows: ； in Indicates the angle of internal friction; S12. Construct an elliptical boundary surface. The boundary surface equation is as follows: ； in and They represent the mean principal stress and shear stress, respectively. Represents the actual stress. The intercept of the boundary surface on the mean principal stress axis represents the boundary surface size; S13, with Hardening parameters are defined by... The evolution equation for achieving isotropic hardening of the boundary surface is as follows: ； in The initial void ratio, and These represent the normal consolidation line and the unloading springback line, respectively. Slope in space Indicates the strain increment of a plastic body; S2. After defining the boundary surface, the evolution mechanism of the mapping center is specified by the uniqueness theory of image stress, and then the image stress is calculated based on the radial mapping rule. S3. After completing the image stress calculation, construct a multi-stage interpolation equation to calculate the plastic modulus; S31. Obtain the plastic modulus of the stress point from the boundary surface equation and based on the consistency condition of the boundary surface plasticity theory. The process is as follows: ； ； Among them, loading factor pass Calculations show that The plastic modulus represents the actual stress point, and the operator " "Written as Macaulay brackets, when When it is a positive value ,on the contrary ; S32, Based on the plastic modulus of the stress point Determine the plastic modulus of the actual stress point , construct computing The interpolation function is as follows: ； in To satisfy the general rules of plastic modulus, if If the stress point coincides with the current stress point, then... ,like Then stress reversal occurs, at which point , This item reflects the effect of void ratio on plastic modulus. The nonnegative model parameters are determined by the following formula: ； ； in Represents the damage state variable. This represents the model parameters. When the model parameters are positive, the simulation loop weakens, and when they are negative, the simulation loop is stable. S4. After determining the plastic modulus, construct a stress-strain incremental relationship model based on the incremental constitutive theory under the elastoplastic framework. S5. After constructing the model, the basic parameters are determined through conventional indoor experiments, and the interpolation parameters are determined through fitting and trial calculations to complete the model parameter calibration. S51. Five basic parameters are determined through routine indoor tests: compression parameters springback parameters Poisson's ratio Critical stress ratio And the void ratio corresponding to an isotropic consolidation pressure of 1 kPa ; S52. Determine the two interpolation parameters. and The process is as follows: S521, The cumulative plastic strain is extremely small in the initial cyclic stage and It has a weak effect on the stress-strain relationship. and Decoupling of the impact effects; S522, Determined by optimal fitting of the initial stress-strain hysteresis loop ; S523. Determine whether the soil is in a cyclically stable or cyclically weakened state. The sign of the variable is determined, and then other parameters are fixed, and the result is determined through multiple trial calculations. ; S6. Construct a high-precision explicit integration algorithm suitable for the numerical implementation of the proposed boundary surface constitutive model; S61, Order , The initial stress state and internal variables of this incremental step are determined by the initial conditions. S62. Calculate the current stress state loading factor. ,like Update the mapping center to the current stress point; otherwise, keep the relative position of the mapping center and the boundary surface unchanged. S63. Determine the sub-step strain increment. The stiffness matrix is ​​obtained from the initial stress state. The stress increment of the first trial calculation was obtained. Including other internal variables, the stress calculation formula for the first trial calculation is as follows: ； S64. Stress state obtained from the initial trial calculation Construct the stiffness matrix for the next trial calculation using the updated intrinsic variables. Calculate the stress increment for the next trial calculation. Including other internal variables, the stress calculation formula for the second trial calculation is as follows: ； S65. Calculate the average stress increment between two stress increments. : ； S66. Calculate the relative error, the process is as follows: ； in Indicates relative error. Indicates the error control value, if If so, the substep size needs to be reduced and the calculation returned to S63. The new substep size calculation method is as follows: ； like Then update the stress state. and internal variables; S67, Order The new substep size is determined using the same method as S66. ,like Then let Return to S63, when When the incremental step calculation ends, the stiffness matrix is ​​updated with the final stress state and transmitted to the analysis software, thus completing the incremental step calculation.

2. The clay cyclic stabilization and weakening constitutive model and numerical implementation method according to claim 1, characterized in that: Using the associated flow rule, the shape and evolution of the plastic potential surface are also determined by the boundary surface equation in S12 and the isotropic hardening evolution criterion in S13.

3. The clay cyclic stabilization and weakening constitutive model and numerical implementation method according to claim 2, characterized in that, The calculation process for image stress in S2 is as follows: S21. Set the initial mapping center Fixed at the origin, the mapping center is aligned with the actual stress. The intersection of the line connecting the two sides with the boundary surface is the stress point. ; S22. According to the theory of image stress uniqueness, that is, when the boundary surface undergoes dynamic hardening, the mapping center and its relative position must remain unchanged, and when stress reversal occurs, the mapping center should be updated to the current true stress point. The evolution equation of the mapping center is defined as follows: ； ； S23. Image stress is calculated from the mapping center and the true stress based on the radial mapping rule: ； ； Where parameters By simultaneously solving the evolution equations of the mapping center and the boundary surface governing equations, we can construct a parameter-related equation. The standard quadratic equation is solved.

4. The clay cyclic stabilization and weakening constitutive model and numerical implementation method according to claim 3, characterized in that, The process of constructing the stress-strain increment relationship model in S4 is as follows: ； in Indicates the volumetric strain increment. Represents the increment of shear strain. The elastic-plastic stiffness matrix is ​​represented by the following formula: ； ； in This represents the elastic stiffness matrix, determined by Hooke's Law.

Images