Softening/Hardening Damage Model and Numerical Implementation of Seabed Silt-Steel Interface in Yellow River Underwater Delta

: The interaction between soil and structure is a research hotspot in ocean engineering, and the shear performance of interfaces is an essential factor affecting the bearing capacity of offshore structures. Taking the Yellow River Underwater Delta as the research area, the Softening/Hardening damage model of the silt–steel interface and the determination method of model parameters are proposed based on the statistical damage theory. Through the interface monotonic shear test under the conditions of different normal stress, roughness and water content, the shear mechanical properties and volumetric deformation laws on the silt–steel interface are analyzed, and the damage model parameters are obtained. Finally, a FRIC subroutine for the damage model was developed based on ABAQUS. The research results indicate the following: (1) The interface between silt and steel exhibits two characteristics, softening/hardening and shear shrinkage/expansion, under different conditions. Roughness signiﬁcantly impacts interfacial cohesion, while water content mainly affects the internal friction angle. (2) The softening model based on the classic rock damage model can better simulate the stress–strain relationship of the silt–steel interface under high normal stress and low water content. In contrast, the hardening model based on the classic hyperbola model can better simulate the stress–strain relationship under low normal stress and high water content. The calculated results of the softening/hardening model agree with the experimental results, and the model has 7 parameters. (3) The developed FRIC subroutine can effectively simulate the nonlinear mechanical behavior of the interface between silt and steel. The research results provide a reference for exploring the stability analysis of offshore structures considering interface weakening effects.


Introduction
The ocean contains many mineral, oil and wind energy resources. In today's global energy shortage and deteriorating ecological environment, ocean wind energy, as an emerging strategic resource, is increasingly valued by countries worldwide for its development and utilization [1]. China's coastline is narrow and long, with abundant wind energy reserves that can be developed and utilized at sea. With the support of the "Ocean Power" and "Dual Carbon" policies, the development of marine wind power-related industries is rapid. In 2022, China's offshore wind power increased its lifting capacity by 5.16 million kilowatts, Power" and "Dual Carbon" policies, the development of marine wind power-related industries is rapid. In 2022, China's offshore wind power increased its lifting capacity by 5.16 million kilowatts, accounting for approximately 54% of the global total. The cumulative lifting capacity was 30.51 million kilowatts, an increase of over 20% year-on-year, and it continues to maintain its position as the world's largest [2].
The Yellow River Delta is rich in oil and gas resources and wind energy resources, and the mudflat and shallow sea areas are vast. Under the development strategy of "transformation of new and old kinetic energy", several wind power bases in the south, north and Bohai Peninsula have been planned (Figure 1), with the planned installed capacity of 23.4 million kW, which will significantly promote the large-scale development of offshore wind power in the Yellow River Delta in the next few years [3]. In the shallow sea area of the Yellow River Delta, silty seabed soil is formed in the process of rapid sedimentation, pore water is not effectively discharged and surface sediment is not compacted, so it has the characteristics of high water content and low shear strength [4]. The construction of offshore wind power will change the seabed soil's original stable environment and boundary conditions. Due to the multi-process coupling of complex dynamic loads and the interaction between structures and seabed soil, it is easy to induce the weakening of the soilstructure interface [5], which in turn leads to local/overall failure of the structural foundation, which will seriously affect the safety and stability of wind power foundations. Therefore, studying the mechanical properties of the soil and structural interface in the region is an urgent issue that needs to be addressed. The constitutive interface model is used to describe the relationship between stress and strain during the shear process, analyze the mechanical properties of the interface and is of great significance for studying the interaction between soil and structure [6]. Three standard interface constitutive models exist: elastic, elastic-plastic and damage. Based on the elastic theory of generalized Hooke's law, the hyperbola model has become one of the most popular interface constitutive models because of its few parameters, clear physical meaning, broad applicability to soil types and other advantages. Considering that the rigid plastic model or hyperbola model artificially assumes that the mutual displacement and shear stress in the interface contact zone conform to the inelastic or that nonlinear The constitutive interface model is used to describe the relationship between stress and strain during the shear process, analyze the mechanical properties of the interface and is of great significance for studying the interaction between soil and structure [6]. Three standard interface constitutive models exist: elastic, elastic-plastic and damage. Based on the elastic theory of generalized Hooke's law, the hyperbola model has become one of the most popular interface constitutive models because of its few parameters, clear physical meaning, broad applicability to soil types and other advantages. Considering that the rigid plastic model or hyperbola model artificially assumes that the mutual displacement and shear stress in the interface contact zone conform to the inelastic or that nonlinear cannot fully fit the actual situation, many scholars combine the nonlinear elastic theory with the elastoplastic theory and propose the elastoplastic model of the soil-structure interface. Liu et al. (2014) [7] proposed a constitutive interface model based on critical state soil mechanics and generalized plasticity, which can simulate the monotonic and cyclic 3D behavior of soil-structure interfaces over a wide range of soil density, normal pressure and normal stiffness. Stutz et al. (2016) [8] proposed an enhanced non-plastic interface model that converts existing plastic models with predefined sand limit state surfaces into interface models by reducing stress and tension vectors and redefining tensors. Hosseinali and Toufigh (2018) [9] proposed a new plasticity-based constitutive model that can model strainsoftening behavior in the tangential direction and shear expansion behavior in the normal direction. Saberi et al. (2020) [10] continuously revised and improved a two-sided plastic interface constitutive model by introducing factors such as critical state soil mechanics, particle breakage rate, new failure surfaces and 3D shear coupled loads, enabling the interface model to simulate complex behaviors of various soil-structure interfaces, such as stress hardening and softening, stress path dependence, phase transition, cyclic cumulative shrinkage and stability, stress degradation and particle breakage. Lu et al. (2017) [11] proposed a constitutive model for soil-structure interfaces by combining the tangential stress-strain relationship and normal stress-strain relationship of contact surface soils.
Damage mechanics is the rudiment of the discipline proposed by relevant scholars when studying the deformation and failure of metal materials. Subsequently, Ajcinovic and Silva (1982) [12] first proposed the concept of the "damage factor," which laid the foundation for developing the constitutive model of damage mechanics. The so-called "damage", the gradual creep evolution and ultimate failure of materials under loading is a branch of the solid mechanics of materials. In recent years, relevant scholars have made additional amendments based on metal damage mechanics in geotechnical engineering, considering that the rock and soil still bear part of the shear strength and hydrostatic pressure when damaged. They have developed rock and soil damage mechanics by viewing the mechanical properties of rock and soil and introducing stochastic statistical theory. Referring to the statistical model of damage evolution in rock mechanics, many scholars have also proposed a damage mechanics model suitable for soil-structure contact. Hu and Pu (2004) [13] proposed a damage constitutive model with ten parameters to describe the behavior of rough interfaces. Long et al. (2017) [14] proposed a damage model with only four parameters based on Weibull distribution statistical damage theory, which can effectively describe strain softening behavior. Zhao et al. (2017) [15] adopted the theory of damage mechanics based on the assumption of equal strain distribution, zero elastic region and no coupling of elastic strain and selected the ratio of irreversible volume strain to the maximum irreversible volume strain as the damage factor and established a shear strength and compressive volume strain damage model to describe the shear behavior of the interface between frozen soil and structure. Zhu and Guo (2016) [16] and Chen et al. (2022) [17] derived a constitutive model for frozen soil damage under dynamic loads based on deep learning, which can reflect the relationship between temperature and strain. Sun et al. (2020) [18] proposed an elastic-plastic damage constitutive model based on the overload/overload yield surface by observing the mechanical behavior of frozen sand and frozen saline soil in triaxial compression tests under different confining pressures. Liu et al. (2020) [19] introduced the concept of thermal damage and constructed a thermal damage model for rocks, expanding the research on rock strength theory. Rehman and Zhang (2021) [20] established a three-dimensional elastic-plastic damage model for gravel soil-structure interface considering the shear coupling effect. From this, it can be seen that the evolution of interface strength under shear is analogous to the damage process of metal materials. Assuming that the stress and strain at the interface conform to a specific statistical law based on the existing rock damage constitutive relationship, the damage mathematical model at this interface type is reasonable and feasible by modifying the shear test results.
Based on the research of relevant scholars on statistical damage theory, this paper proposes a softening/hardening damage model for the silt-steel interface in the Yellow River Delta. The model has 7 parameters and was obtained through large-scale interface shear tests. Finally, the rationality of the model is verified. At the same time, its application in the numerical calculation is realized, which is convenient for the subsequent numerical simulation research on the stability of offshore structures, considering the interface weakening effect, and is expected to guide the safety assessment of offshore installations in the Yellow River Delta.

Assumption
Numerous research results have shown that the soil-structure interface exhibits a strain-softening phenomenon during the shear process [21,22]. To derive the softening/hardening damage mathematical model applicable to the silt-steel interface in the Yellow River Delta and describe the two states of hardening and softening, this paper will combine the classic hyperbola model [23] and the classic rock damage model [14] to build.
According to damage mechanics and stochastic statistical theory, the following assumptions are set: (1) Assuming the soil is isotropic and the structure is a fixed rigid body.
(2) Divide the soil-structure interface into several microelements, assuming that there are only two states: undamaged (0) and damaged (1). (3) Assuming that all friction and contact on the soil-structure interface (between soil particles, between soil structure) are unevenly and randomly distributed, the shear strength and shear failure at the interface are also unevenly and randomly distributed. (4) Assuming that the damage process at the soil-structure interface is a continuous process that undergoes an evolution from undamaged to damage, the soil at the interface can be divided into two parts, damaged soil and undamaged soil (Figure 2), and the undamaged soil bears effective stress. River Delta. The model has 7 parameters and was obtained through large-scale interface shear tests. Finally, the rationality of the model is verified. At the same time, its application in the numerical calculation is realized, which is convenient for the subsequent numerical simulation research on the stability of offshore structures, considering the interface weakening effect, and is expected to guide the safety assessment of offshore installations in the Yellow River Delta.

Assumption
Numerous research results have shown that the soil-structure interface exhibits a strain-softening phenomenon during the shear process [21,22]. To derive the softening/hardening damage mathematical model applicable to the silt-steel interface in the Yellow River Delta and describe the two states of hardening and softening, this paper will combine the classic hyperbola model [23] and the classic rock damage model [14] to build.
According to damage mechanics and stochastic statistical theory, the following assumptions are set: (1) Assuming the soil is isotropic and the structure is a fixed rigid body.
(2) Divide the soil-structure interface into several microelements, assuming that there are only two states: undamaged (0) and damaged (1). (3) Assuming that all friction and contact on the soil-structure interface (between soil particles, between soil structure) are unevenly and randomly distributed, the shear strength and shear failure at the interface are also unevenly and randomly distributed. (4) Assuming that the damage process at the soil-structure interface is a continuous process that undergoes an evolution from undamaged to damage, the soil at the interface can be divided into two parts, damaged soil and undamaged soil (Figure 2), and the undamaged soil bears effective stress.

Model Derivation
Based on the above assumptions, the damage process at the soil-structure interface satisfies the following relationship: In the formula, B represents the degree of damage; τ1 and τ2 refer to the shear stress borne by damaged and undamaged soil, respectively.

Model Derivation
Based on the above assumptions, the damage process at the soil-structure interface satisfies the following relationship: In the formula, B represents the degree of damage; τ 1 and τ 2 refer to the shear stress borne by damaged and undamaged soil, respectively. Subsequently, Equation (2) can be transformed into the following form: For undamaged soil N 2 , the stress-strain relationship satisfies linear elasticity, and the shear stress is manifested as: In the formula, γ 2 is the shear strain of undamaged soil; G 2 is the shear modulus of undamaged soil.
For the damaged soil N 1 , when B approaches 1, τ trend τ 1 . The damaged soil describes the mechanical characteristics of the overall interface failure.
Combined with the previous shear test results of silt in the Yellow River Delta [24], the hardening curve's shear stress-strain curve conforms to the traditional hyperbola model; for the softening curve, there is residual friction strength after shear failure. Based on this, it is assumed that the shear stress of the damaged soil mass τ 1 satisfies the following equation: In the formula, γ 1 is the shear strain of the damaged soil; G 1 is the shear modulus of the damaged soil; α is the hyperbola model parameter; τ r is the residual shear stress.
Based on the assumption of equal strain, the following conditions are met: From Equations (3)-(5), a mathematical model for the softening/hardening damage of the silt-steel interface is obtained: Assuming that the interface satisfies the Weibull distribution theory [25]: In the formula, p ƒ is the probability that damage may occur; n and m are parameters of the Weibull distribution.
In damage mechanics, damage degree B can represent the probability of shear failure, that is: The simultaneous Equations (7)-(9) provide a model for the softening/hardening damage of the silt-steel interface based on statistical distribution theory.

Model Parameters
When shear continues, and damage continues to evolve, there are: In the formula, τ max is the maximum shear stress.
Derivation of Equation (7): At the initial shear stage, if the interface is undamaged as a whole, then the initial shear modulus G 0 is the undamaged soil shear modulus G 2 : For hardening curves, when γ→∞, the following equation is satisfied: For the softening type curve, it is first necessary to obtain the peak point of the stressstrain curve, which corresponds to the maximum shear stress before the residual friction stage. The corresponding shear strain at this point is γ m . If the degree of damage is B m , there are the following: The parameters m and n can be obtained by fitting the stress-strain curve of the shear test.

Test Equipment and Plan
In previous studies, the author explored the silt-steel interface's cyclic shear mechanical properties [26], performing monotonic shear tests using the same instruments and materials. The shear tester and steel plate parameters are shown in Figure 3, and the steel plate is set with three different roughness levels (R = 0, 0.025, 0.05 mm). The test soil is a typical silty seabed soil in the Yellow River Delta, with an average particle size of d 50 = 0.03mm and an optimal water content of 16.9%. Three water contents are set (ω = 16%, 20%, 24%), and the soil samples were partially saturated during the experiment. At the same time, three normal stress conditions (σ = 100 kPa, 200 kPa, 300 kPa) were set up. Fifteen sets of undrained fast shear tests were conducted under constant normal load conditions (CNC) to analyze the silt-steel interface's mechanical shear properties and volumetric deformation laws under different conditions (Table 1) and to obtain the softening/hardening damage model parameters. In the shear box, the soil sample is compacted in 5 layers, with each layer measuring 6 cm.

Stress-Strain Rule
(1) Normal stress Figure 4a shows the interface shear stress-strain curves under normal stress of 100 kPa, 200 kPa and 300 kPa, respectively (roughness R = 0.05 and water content ω = 20%). It can be seen that while the roughness of the steel plate and the water content of the soil remain unchanged, the initial shear modulus and maximum shear stress at the interface increase with the increase of normal stress. The shear stress-strain curve under low normal stress (σ = 100 kPa, 200 kPa) exhibits a hardening type, where the maximum shear stress is taken as the shear stress at a shear strain of 50 mm; under high normal stress (σ = 300 kPa), there is an evident softening trend in the shear stress-strain curve, which is a softening type. Here, the peak shear stress is taken as the maximum shear stress. Whether it is a hardening curve or a softening curve, it can always be divided into three stages: the stage of rapid elastic-plastic growth (OA, OD), the stage of slow plastic growth (AB)/plastic softening (DE) and the stage of near failure equilibrium (BC)/residual fraction (EF). Calculated at σ = 300 kPa, the maximum shear stress is 2.44 and 1.35 times that at 100 kPa and 200 kPa, respectively. The increase in shear stress decreases with the rise of normal stress.

Stress-Strain Rule
(1) Normal stress Figure 4a shows the interface shear stress-strain curves under normal stress of 100 kPa, 200 kPa and 300 kPa, respectively (roughness R = 0.05 and water content ω = 20%). It can be seen that while the roughness of the steel plate and the water content of the soil remain unchanged, the initial shear modulus and maximum shear stress at the interface increase with the increase of normal stress. The shear stress-strain curve under low normal stress (σ = 100 kPa, 200 kPa) exhibits a hardening type, where the maximum shear stress is taken as the shear stress at a shear strain of 50 mm; under high normal stress (σ = 300 kPa), there is an evident softening trend in the shear stress-strain curve, which is a softening type. Here, the peak shear stress is taken as the maximum shear stress. Whether it is a hardening curve or a softening curve, it can always be divided into three stages: the stage of rapid elastic-plastic growth (OA, OD), the stage of slow plastic growth (AB)/plastic softening (DE) and the stage of near failure equilibrium (BC)/residual fraction (EF). Calculated at σ = 300 kPa, the maximum shear stress is 2.44 and 1.35 times that at 100 kPa and 200 kPa, respectively. The increase in shear stress decreases with the rise of normal stress.
In addition, the shear strain reaching the inflection point of shear stress rate under different normal stresses also exhibits a certain regularity; the more significant the normal stress, the greater the shear strain required to reach the inflection point of the shear rate. The essence of shear is to perform destructive work, which can be explained by the fact that as the normal stress increases, on the one hand, the lower soil is compressed more   15 300 adopted in this experiment, this situation may be related to the matrix suction of unsaturated soil [30].

Interface Shear Strength
The relationship between peak shear stress/maximum shear stress and normal stress in this experiment was linearly fitted, as shown in Figure 5, with a good fit. The interface failure mode conforms to the Mohr-Coulomb failure criterion, namely: In addition, the shear strain reaching the inflection point of shear stress rate under different normal stresses also exhibits a certain regularity; the more significant the normal stress, the greater the shear strain required to reach the inflection point of the shear rate. The essence of shear is to perform destructive work, which can be explained by the fact that as the normal stress increases, on the one hand, the lower soil is compressed more tightly. On the other hand, the soil particles are in closer contact with the steel surface under pressure, forming a good friction effect, thereby increasing the energy consumption required to reach the critical failure state, leading to a more significant maximum shear stress and rate inflection point shear strain.
(2) Roughness Figure 4b-d shows the shear stress-strain curves on three steel plates with different roughness, R = 0, 0.025 and 0.05mm, under normal stress of 100, 200, and 300 kPa, respectively (water content ω = 20%). The experimental results show that under the same normal stress conditions, different roughness curves overlap in the initial stage; the initial shear stiffness foundation is consistent in the elastic stage. As the roughness increases, shear failure and separation occur successively. The maximum shear stress and shear stress rate strain increase with the roughness increase, indicating that roughness impacts the shear strength of the silt-steel interface under various normal stresses. Analysis of the reasons shows that less soil can be embedded in the groove when the interface roughness is slight. The shear strain is mainly caused by the displacement of the soil inside the groove, and direct sliding failure occurs at the interface, resulting in a smaller displacement of the maximum shear stress and shear stress rate inflection point; when the roughness is high, there is more soil that can be embedded in the groove, and its shear strain is generated by the soil in the interface groove and the soil below the interface. Elastoplastic failure occurs at the interface, and a sliding shear zone appears inside the soil, resulting in a more significant strain of the maximum shear stress and shear stress rate inflection point.
In addition, as the normal stress increases, the shear stress-strain curves under different roughness levels tend to develop towards a hardening type, a folding type and a softening type. Moreover, the curves under different roughness levels become more similar, indicating that the influence of roughness on shear stress decreases accordingly. As the pressure increases, the number of soil particles squeezed into the groove increases, improving the contact surface's non-uniformity. Related scholars have proposed the concepts of "critical roughness" and "ultimate roughness" when studying the effect of roughness on interface shear behavior [27][28][29]. However, due to the simplified sand-filling method used in the roughness setting in this experiment, it is no longer discussed.
(3) Water content The water content directly affects the strength of the soil itself and also affects the strength of its interface. Figure 4e-g shows the shear stress-strain curves on three water content soil samples at 16%, 20% and 24% (with a water content of 15.9%, 19.8% and 23.2% measured after the experiment) under normal stress of 100, 200 and 300 kPa, respectively (R = 0.05 mm). From the experimental results, it can be seen that when the normal stress conditions are the same in the initial stage, the curve trends of the groups (ω = 16%, 20%) closer to the optimal water content of 19.6% are similar. The slope of the curve decreases with the increase in water content. The shear stress-strain curve shows a softening type in the group with low water content. As the water content increases, the shear stress-strain curve transitions from a softening type to a hardening type. The maximum shear stress of the hardening curve and the peak shear stress of the softening curve decrease with the increase in water content. The reason for this is that under low water content conditions, the particle contact of the sample gradually becomes close during the shear process, the soil-water interaction increases, the shear stress increases and the sample exhibits plastic failure; when the water content increases, the thickness of the film water between the silt and the steel structure increases, increasing lubrication, resulting in a decrease in the frictional resistance between the interfaces. An increase in pore water pressure reduces effective stress, ultimately decreasing its shear strength. Due to using a method to control dry density during the soil loading process, the influence of compaction was ruled out.
In addition, under low normal stress conditions (σ = 100 kPa, 200 kPa), the interface with smaller water content exhibits shear failure first. Due to the undrained fast shear adopted in this experiment, this situation may be related to the matrix suction of unsaturated soil [30].

Interface Shear Strength
The relationship between peak shear stress/maximum shear stress and normal stress in this experiment was linearly fitted, as shown in Figure 5, with a good fit. The interface failure mode conforms to the Mohr-Coulomb failure criterion, namely: In the formula, σ is the vertical stress; φ is the apparent angle of internal friction; c is the apparent cohesive force. (1) Roughness Figure 6a,b shows the shear strength and its parameters under different roughness conditions. The cohesion and internal friction angles corresponding to roughness R = 0, 0.025 and 0.05 mm are 2.25 kPa, 4.03 kPa, 5.76 kPa and 11.21°, 12.68° and 13.19°, respectively. As the roughness increases, the cohesion and internal friction angle increase 1.56 times and by 17.67%, respectively, indicating that the roughness greatly affects the interface cohesion. As the roughness increases, the depth of the structural groove deepens and the shear resistance required for soil particles to slip and roll increases; at the same time, the contact area between soil particles and the surface of the structure increases, leading to an increase in the filling of soil in the groove during the shear process, which enhances the bonding effect between the soil and the soil in the groove. At this point, the bonding strength becomes the main part of the shear strength. This pattern is more evident in the experiments of different groove structure steel plates with greater differences in roughness [31]. The friction angle refers to the friction between interfaces. As the roughness increases, the amount of soil filled into the groove increases, resulting in weakened friction between the soil and the smooth part of the structural surface (L2), while the friction between the soil and the soil inside the groove increases, ultimately manifesting as a slight increase in the interface friction angle.
(2) Water content Figure 6c,d shows the shear strength and parameters under different water content conditions. The apparent cohesive force and apparent angle of internal friction corresponding to 16%, 20% and 24% water content are 4.7 kPa, 5.76 kPa, 3.61 kPa and 15.08°, 13.19° and 10.5°, respectively. As the water content increases, the cohesion shows a trend of first increasing and then decreasing, with increases of 22.6% and -59.5%, respectively, and the rate of growth at 16% to 20% is less than the rate of decrease at 20% to 24%. When the water content is low and then increases, the density of the soil sample increases, and at the same time, the adsorption effect of water between soil particles increases the soil cohesion. According to previous research, a peak of cohesion occurs when the water content reaches the plastic limit or the optimal water content [32]. In this study, the peak of cohesion appears near the optimal moisture content (19.6%); as the water content further increases, the thickness of the membrane water between soil particles increases, weak bonding water increases and free water begins to be generated. Soil particles separate un- In the formula, σ is the vertical stress; ϕ is the apparent angle of internal friction; c is the apparent cohesive force.
(1) Roughness Figure 6a,b shows the shear strength and its parameters under different roughness conditions. The cohesion and internal friction angles corresponding to roughness R = 0, 0.025 and 0.05 mm are 2.25 kPa, 4.03 kPa, 5.76 kPa and 11.21 • , 12.68 • and 13.19 • , respectively. As the roughness increases, the cohesion and internal friction angle increase 1.56 times and by 17.67%, respectively, indicating that the roughness greatly affects the interface cohesion. As the roughness increases, the depth of the structural groove deepens and the shear resistance required for soil particles to slip and roll increases; at the same time, the contact area between soil particles and the surface of the structure increases, leading to an increase in the filling of soil in the groove during the shear process, which enhances the bonding effect between the soil and the soil in the groove. At this point, the bonding strength becomes the main part of the shear strength. This pattern is more evident in the experiments of different groove structure steel plates with greater differences in roughness [31]. The friction angle refers to the friction between interfaces. As the roughness increases, the amount of soil filled into the groove increases, resulting in weakened friction between the soil and the smooth part of the structural surface (L 2 ), while the friction between the soil and the soil inside the groove increases, ultimately manifesting as a slight increase in the interface friction angle.
(2) Water content Figure 6c,d shows the shear strength and parameters under different water content conditions. The apparent cohesive force and apparent angle of internal friction corresponding to 16%, 20% and 24% water content are 4.7 kPa, 5.76 kPa, 3.61 kPa and 15.08 • , 13.19 • and 10.5 • , respectively. As the water content increases, the cohesion shows a trend of first increasing and then decreasing, with increases of 22.6% and −59.5%, respectively, and the rate of growth at 16% to 20% is less than the rate of decrease at 20% to 24%. When the water content is low and then increases, the density of the soil sample increases, and at the same time, the adsorption effect of water between soil particles increases the soil cohesion. According to previous research, a peak of cohesion occurs when the water content reaches the plastic limit or the optimal water content [32]. In this study, the peak of cohesion appears near the optimal moisture content (19.6%); as the water content further increases, the thickness of the membrane water between soil particles increases, weak bonding water increases and free water begins to be generated. Soil particles separate under water pressure, and the contribution of structural suction to cohesion decreases, leading to a rapid decrease in cohesion. As the water content increases, the lubrication effect of the water film between the interfaces increases, and the internal friction angle shows a downward trend, with a decrease of 30.37%, indicating that the water content dramatically affects the internal friction angle of the contact surface. leading to a rapid decrease in cohesion. As the water content increases, the lubrication effect of the water film between the interfaces increases, and the internal friction angle shows a downward trend, with a decrease of 30.37%, indicating that the water content dramatically affects the internal friction angle of the contact surface.

Volumetric Deformation Laws
The cross-sectional area of the shear interface in this experiment is a constant value. Therefore, its normal strain is measured to reflect the volume change of the sample during the shear process, which is the shear volumetric deformation laws of the interface.
It can be seen that while the roughness of the steel plate and the water content of the soil remain unchanged, the interfacial deformation patterns presented by different normal stresses are different. Under low-stress (σ = 100 kPa) conditions, the interface exhibits shear dilation. With the increase of normal stress, under medium stress (σ = 200 kPa) conditions, the interfacial deformation shows a pattern of initial shear shrinkage followed by slight shear dilation and overall reduction. The normal stress further increases, and under high-stress (σ = 300 kPa) conditions, the interface exhibits shear shrinkage. However, regardless of its shear shrinkage or shear expansion, the shear strain-normal strain curve can always be divided into three stages: rapid volumetric deformation stage (OA, OD), slow volumetric deformation stage (AB)/volumetric transformation direction stage (DE) and volumetric equilibrium stage (BC, EF). The soil bears a portion of the compaction work during the compaction process. During the shear process, due to the displacement of soil particles, the accumulated cohesive and frictional potential energy during the compaction process is released, resulting in an upward force. When the normal stress is low, this upward force and the normal stress offset the pressure difference, which is positive, thus exhibiting shear dilation; when the normal stress is high, the pressure difference be-

Volumetric Deformation Laws
The cross-sectional area of the shear interface in this experiment is a constant value. Therefore, its normal strain is measured to reflect the volume change of the sample during the shear process, which is the shear volumetric deformation laws of the interface.
It can be seen that while the roughness of the steel plate and the water content of the soil remain unchanged, the interfacial deformation patterns presented by different normal stresses are different. Under low-stress (σ = 100 kPa) conditions, the interface exhibits shear dilation. With the increase of normal stress, under medium stress (σ = 200 kPa) conditions, the interfacial deformation shows a pattern of initial shear shrinkage followed by slight shear dilation and overall reduction. The normal stress further increases, and under high-stress (σ = 300 kPa) conditions, the interface exhibits shear shrinkage. However, regardless of its shear shrinkage or shear expansion, the shear strain-normal strain curve can always be divided into three stages: rapid volumetric deformation stage (OA, OD), slow volumetric deformation stage (AB)/volumetric transformation direction stage (DE) and volumetric equilibrium stage (BC, EF). The soil bears a portion of the compaction work during the compaction process. During the shear process, due to the displacement of soil particles, the accumulated cohesive and frictional potential energy during the compaction process is released, resulting in an upward force. When the normal stress is low, this upward force and the normal stress offset the pressure difference, which is positive, thus exhibiting shear dilation; when the normal stress is high, the pressure difference between the upward force and the normal stress is negative, and the soil particle gap is compressed, thus exhibiting shear shrinkage [33].
of the filled soil, soil samples with lower water content need to withstand more excellent compaction work so that the compacted soil may have over-consolidation properties [35]. At the same time, as the normal stress increases, the degree of bulk deformation changes less, indicating that under low water content conditions, the normal stress increases to a certain extent, and its impact on interface bulk deformation weakens. This is because the more significant the normal stress, the denser the soil compression and the weaker the water sensitivity, making the mechanical properties on the contact surface relatively less affected by the water content. However, high water content (ω = 24%) consistently exhibits shear shrinkage, even under low stress (σ = 100kPa) conditions; the reason for this may be due to the high saturation of the soil sample (Sr = 95.2%) which causes a small amount of soil to be squeezed out during the test process, resulting in shear shrinkage. In addition, under the same normal stress, as the water content increases, the initial slope of the volumetric curve decreases, which has the same regularity as the shear stress-shear displacement curve mentioned above.

Model Validation
According to the shear test results, the silt-steel interface exhibits a hardening curve under low to medium stress (σ = 100 kPa, 200 kPa) and medium to high water content (ω = 20%, 24%) conditions; it is a softening curve under high stress (σ = 300 kPa) and low to medium water content (ω = 16%, 20%) conditions.
For hardening curves, the strength parameter is τmax, G0 and G1, with model parameters m and n; for softening curves, the strength parameters are τmax, γm, τr and G0, and the model parameters are m and n, respectively. The strength parameters were obtained In addition, the reason for the phenomenon of first shrinkage and then shear expansion under medium stress (σ = 200 kPa) may be that the previously compressed voids between soil particles undergo particle dislocation and restore some of the voids with the progress of shear, resulting in a phenomenon of the smaller volume first and then increasing volume.
(2) Roughness Figure 7b-d shows the shear strain-normal strain curves on three steel plates with different roughness, R = 0, 0.025 and 0.05 mm, under normal stress of 100, 200 and 300 kPa, respectively (water content ω = 20%). Under the same normal stress, the interface body becomes more significant as the roughness increases. To be more precise, in the low-stress shear expansion stage, the more significant the roughness, the greater the shear expansion amount; in the high-stress stage, the more significant the roughness, the greater the reduction. The increase in roughness enhances the degree of engagement between soil particles and the interface, leading to an increase in soil particles embedded in rough grooves and an increase in the thickness of the interface shear band [34], increasing volumetric deformation.
(3) Water content Figure 7e-g shows the shear strain-normal strain curves on three water contents soil samples at 16%, 20%, and 24% (with a water content of 15.9%, 19.8%, and 23.2% measured after the experiment) under normal stress of 100, 200, and 300 kPa, respectively (R = 0.05 mm). Under different normal stresses, lower water content (ω = 16%) consistently exhibits a certain degree of shear dilation. Even under high-stress conditions (σ = 300 kPa), shear strain-normal strain curves also exhibit a first shrinkage pattern and then shear dilation. The reason for this phenomenon may be due to the filling method. To control the dry density of the filled soil, soil samples with lower water content need to withstand more excellent compaction work so that the compacted soil may have over-consolidation properties [35]. At the same time, as the normal stress increases, the degree of bulk deformation changes less, indicating that under low water content conditions, the normal stress increases to a certain extent, and its impact on interface bulk deformation weakens. This is because the more significant the normal stress, the denser the soil compression and the weaker the water sensitivity, making the mechanical properties on the contact surface relatively less affected by the water content. However, high water content (ω = 24%) consistently exhibits shear shrinkage, even under low stress (σ = 100 kPa) conditions; the reason for this may be due to the high saturation of the soil sample (S r = 95.2%) which causes a small amount of soil to be squeezed out during the test process, resulting in shear shrinkage. In addition, under the same normal stress, as the water content increases, the initial slope of the volumetric curve decreases, which has the same regularity as the shear stress-shear displacement curve mentioned above.

Model Validation
According to the shear test results, the silt-steel interface exhibits a hardening curve under low to medium stress (σ = 100 kPa, 200 kPa) and medium to high water content (ω = 20%, 24%) conditions; it is a softening curve under high stress (σ = 300 kPa) and low to medium water content (ω = 16%, 20%) conditions. For hardening curves, the strength parameter is τ max , G 0 and G 1 , with model parameters m and n; for softening curves, the strength parameters are τ max , γ m , τ r and G 0 , and the model parameters are m and n, respectively. The strength parameters were obtained through inverse calculation of the shear test results, where G 0 and G 1 were obtained from calculating initial shear modulus in material mechanics and Hooke's law. We selected the softening/hardening test data from the shear test and compared them with the model. Using the custom fitting function of MATLAB software, the m and n values under the above conditions can be obtained through regression analysis, as shown in Table 2. The comparison curves of simulation and experimental results are plotted in Figure 8.   Table 2 shows that the theoretical calculation and shear test results under the four experimental conditions differ slightly for the hardening model. The theoretical values of the maximum shear stress τmax differ by 0.7%, 0.3%, 3.8% and 3.5% compared to the experimental results. From Figure 8a, it can be seen that the theoretical calculation curve has a high degree of agreement with the shear test curve, which can accurately reflect the hardening characteristics. The hardening model in the damage model, based on the classical hyperbola model, can better simulate the stress-strain relationship of the silt-steel interface under conditions of low normal stress and high water content.
From Table 2 Table 2 shows that the theoretical calculation and shear test results under the four experimental conditions differ slightly for the hardening model. The theoretical values of the maximum shear stress τ max differ by 0.7%, 0.3%, 3.8% and 3.5% compared to the experimental results. From Figure 8a, it can be seen that the theoretical calculation curve has a high degree of agreement with the shear test curve, which can accurately reflect the hardening characteristics. The hardening model in the damage model, based on the classical hyperbola model, can better simulate the stress-strain relationship of the silt-steel interface under conditions of low normal stress and high water content.
From Table 2, it can be seen that for the softening model, the theoretical calculation results and shear test results of the four experimental conditions also have slight differences, with the theoretical values of the maximum shear stress τ max varying by 1.5%, 1.6%, 3.4% and 3.6%, respectively, compared to the test results; the theoretical values of residual shear stress τ r differ by 0.8%, 2%, 2.9% and 0.8% compared to the experimental results, and the theoretical values of maximum shear stress displacement γ m differ by 9.7%, 4%, 13.5% and 11.9% compared to the experimental results. From Figure 8b, it can be seen that there is a high degree of overlap between the theoretical calculation curve and the shear test curve. Although there is a slight separation between the curves in the initial shear and residual friction stages, they can show a softening trend. The softening model in the damage model, based on the classic rock damage model, can better simulate the stress-strain relationship at the high normal stress and low moisture content states of the silt-steel interface.
In summary, although the theoretical calculation results do not fully match the shear test results, considering factors such as calculation errors, model assumptions and soil complexity, the difference between the theoretical calculation curve and the shear test curve is within a reasonable range, indicating that the model can reflect the characteristics of softening strain at the silt-steel interface under high normal stress or low water content conditions and hardening strain under low normal stress or high water content conditions. The theoretical model is based on damage mechanics, the model parameters are easy to obtain and the accuracy is appropriate. The theoretical model is suitable for describing the shear failure process of the silt-steel structure interface in the Yellow River Delta.
It should be noted here that on the premise of determining the shear modulus G, since the model parameters m and n in the softening and hardening model are affected by normal stress, roughness and water content simultaneously, they cannot be modified by fitting curves. Considering the complexity of factors affecting the soil-structure interface in natural environments, the model still needs further optimization and improvement.

Numerical Implementation
The derivation and establishment of mathematical interface models aim to match numerical calculation methods and implement them in practical engineering applications. ABAQUS software is commonly used for simulating the stability of offshore structures [36].

Subroutine Programming
To apply the damage model obtained in this paper to the future stability analysis of offshore wind power, FORTRAN language was used, referring to Yan (2020) [37] and Shi (2021) [38], and the subroutine design and code writing were carried out according to the FRIC subroutine format in ABAQUS. Considering the convergence of the calculation results and the saving of calculation time, the default thickness-free element is adopted, which does not consider the stiffness coefficient in the interface direction and only needs to provide the stiffness coefficient in the interface slip direction. Therefore, the interface has no thickness in the calculation, only length [39].
Because ABAQUS determines the slip and detachment of the contact surface, the FRIC subroutine only needs to define the friction characteristics on the interface, that is, select the constitutive relationship. The implementation process of the FRIC subroutine is shown in Figure 9. The subroutine is called once for each load increment step during the primary software calculation process. After the subroutine calculation, it returns to the main program to complete the iterative calculation of this loading step, ultimately achieving the implementation of the soft and hard damage model in ABAQUS.

Subroutine Validation
To verify the reliability of the FRIC subroutine, a simple slider finite element model is set up, as shown in Figure 10a, referring to the method proposed by Lu et al. (2017) [11] and Wang et al. (2012) [40]. The upper part is a 1 m × 1 m ×1 m soil block constructed using CAX4R elements. It should be noted that the Poisson's ratio is set to v = 0 to ignore the influence of lateral deformation on shear. The lower part is a 3 m × 3 m ×1 m rigid block using CAX4I units; the main control surface is a rigid body surface with sparse mesh division, while the subordinate surface is a soil block surface with dense mesh division. The face-to-face discretization method is used to set the normal interaction relationship of its contact surface as the hard contact inherent in ABAQUS, and the tangential interaction is set using the developed FRIC subroutine.
According to the two cases of softening/hardening, the loading is divided into two analysis steps: ① apply normal stresses of 100 kPa (hardening type) and 300 kPa (softening type) to the surface of the soil and establish the contact between the two components; ② set up a soil test with a displacement of 0.05 m in the x-direction, that is, apply a load displacement of U = 0.05m to simulate the shear process. The specific model parameters are shown in Table 3, and the dilation angle (ψ) of soil mass is 0.1°.

Subroutine Validation
To verify the reliability of the FRIC subroutine, a simple slider finite element model is set up, as shown in Figure 10a, referring to the method proposed by Lu et al. (2017) [11] and Wang et al. (2012) [40]. The upper part is a 1 m × 1 m ×1 m soil block constructed using CAX4R elements. It should be noted that the Poisson's ratio is set to v = 0 to ignore the influence of lateral deformation on shear. The lower part is a 3 m × 3 m ×1 m rigid block using CAX4I units; the main control surface is a rigid body surface with sparse mesh division, while the subordinate surface is a soil block surface with dense mesh division. The face-to-face discretization method is used to set the normal interaction relationship of its contact surface as the hard contact inherent in ABAQUS, and the tangential interaction is set using the developed FRIC subroutine.
According to the two cases of softening/hardening, the loading is divided into two analysis steps: 1 apply normal stresses of 100 kPa (hardening type) and 300 kPa (softening type) to the surface of the soil and establish the contact between the two components; 2 set up a soil test with a displacement of 0.05 m in the x-direction, that is, apply a load displacement of U = 0.05 m to simulate the shear process. The specific model parameters are shown in Table 3, and the dilation angle (ψ) of soil mass is 0.1 • .

Conclusions
This article studies the pile-soil interaction of offshore structures. Taking steel and silt in the Yellow River Delta as test materials combined with a monotonic shear test under different normal stress, roughness and water content conditions, the shear mechanical properties of the silt-steel interface were studied. A soil-structure interface softening/hardening damage model based on Weibull distribution was proposed, and the damage model parameters were determined. Finally, the rationality of the model was verified.
(1) Compared to the strain hardening, shear shrinkage deformation and strain softening, the shear expansion deformation behavior of coarse-grained soil and cohesive soil silt exhibits two characteristics: softening/hardening and shear shrinkage/expansion under different conditions. The maximum shear stress increases with normal stress, roughness, and water content increase. Roughness has a significant impact on cohesion. As roughness increases, the cohesion and internal friction angle increase 1.56 times and by 17.67%, respectively; water content mainly affects the internal friction angle. As the water content increases, the internal friction angle decreases, with a decrease of 30.37%. The cohesion shows a trend of increasing and then decreasing, with an increase of 22.6% and −59.5%, respectively, and reaches its peak near optimal water content. (2) The proposed shear damage model has fewer parameters and a clear physical meaning. The strength parameters are τmax, G0, G1, γm, τr, and the model parameters are m and n. The softening model, based on the classic rock damage model, can better simulate the stress-strain relationship of the silt-steel interface under high normal stress and low water content. In contrast, the hardening model based on the classic hyperbola model can better simulate the stress-strain relationship under low normal stress and high water content. After using the model parameters obtained from the test results, the theoretical model curve is in good agreement with the shear test curve, which proves that the proposed contact surface damage model is reasonable. However, the model parameters m and n are simultaneously affected by normal stress, roughness and water content, so they cannot be simply corrected by fitting curves. Considering the complexity of factors affecting the soil-structure interface in real environments, the model still needs further optimization and improvement.  The comparison between the shear stress-strain curve simulated by the subroutine and the theoretical model results is shown in Figure 10b,c. It can be seen from the figure that the simulation results of the subprogram are highly consistent with the calculation results of the constitutive model, which can not only simulate the continuous increase of the strain hardening shear stress but also simulate the strain softening residual friction, which verifies the feasibility of the FRIC subprogram, indicating that the interface shear model established in this paper can be successfully realized in the finite element calculation software ABAQUS by using the subprogram FRIC. This provides a model foundation and technical support for future analysis of the bearing capacity of offshore structures considering interface weakening effects.

Conclusions
This article studies the pile-soil interaction of offshore structures. Taking steel and silt in the Yellow River Delta as test materials combined with a monotonic shear test under different normal stress, roughness and water content conditions, the shear mechanical properties of the silt-steel interface were studied. A soil-structure interface softening/hardening damage model based on Weibull distribution was proposed, and the damage model parameters were determined. Finally, the rationality of the model was verified.
(1) Compared to the strain hardening, shear shrinkage deformation and strain softening, the shear expansion deformation behavior of coarse-grained soil and cohesive soil silt exhibits two characteristics: softening/hardening and shear shrinkage/expansion under different conditions. The maximum shear stress increases with normal stress, roughness, and water content increase. Roughness has a significant impact on cohesion. As roughness increases, the cohesion and internal friction angle increase 1.56 times and by 17.67%, respectively; water content mainly affects the internal friction angle. As the water content increases, the internal friction angle decreases, with a decrease of 30.37%. The cohesion shows a trend of increasing and then decreasing, with an increase of 22.6% and −59.5%, respectively, and reaches its peak near optimal water content.
(2) The proposed shear damage model has fewer parameters and a clear physical meaning. The strength parameters are τ max , G 0 , G 1 , γ m , τ r , and the model parameters are m and n. The softening model, based on the classic rock damage model, can better simulate the stress-strain relationship of the silt-steel interface under high normal stress and low water content. In contrast, the hardening model based on the classic hyperbola model can better simulate the stress-strain relationship under low normal stress and high water content. After using the model parameters obtained from the test results, the theoretical model curve is in good agreement with the shear test curve, which proves that the proposed contact surface damage model is reasonable. However, the model parameters m and n are simultaneously affected by normal stress, roughness and water content, so they cannot be simply corrected by fitting curves. Considering the complexity of factors affecting the soil-structure interface in real environments, the model still needs further optimization and improvement. (3) The numerical simulation application of the damage model was achieved through the FRIC subroutine, providing a mathematical model basis for subsequent research on the stability analysis of offshore structures considering interface weakening effects.