Soil-structure interaction finite element analysis method and device

Abstract

The application relates to a soil-structure complex interaction finite element analysis method and device, which comprises the following steps: calling relevant material parameters from a material information database according to material naming; carrying out grid division on a soil body and a structure model, and extracting a material characteristic calculation unit stiffness matrix according to the material information database; determining the actual distance between contact surfaces with plastic displacement according to the plastic displacement between the soil and the structure; calculating the contact node force of each unit through a linear programming calculation method with a contact condition; calculating the unit node displacement and the plastic displacement through an initial stress method, and iteratively obtaining the displacement result, and judging the output strain and stress or re-modifying through a relative difference percentage. Compared with the prior art, the application has higher calculation efficiency for the constitutive of the soil body and the structure, can meet the needs of rapid seismic damage assessment under multiple scenes on one hand, and has strong transplantability on different research problems on the other hand.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of finite element analysis, and relates to a soil-structure complex interaction finite element analysis method and device. BACKGROUND

[0002] The general solution method of nonlinear problems mainly includes variable stiffness method, initial stress method and initial strain method. The variable stiffness method is a commonly used solution method of commercial finite element software at present, and has good universality. However, the unit stiffness matrix needs to be recalculated and combined into the total stiffness matrix at each calculation period, and the calculation efficiency is not high, and the time is long in the application in actual engineering problems. The initial strain method adopts the differential form of the constitutive relation, assumes that the stress is constant in each calculation period, regards the total viscous strain as the initial strain, and calculates the additional force and displacement increment. The stiffness matrix remains unchanged in each period, and the calculation efficiency is higher than that of the variable stiffness method. The initial stress method adopts the differential form of the constitutive relation, and is based on the relaxation modulus of the material. The strain is assumed to be constant in each calculation period, the viscous relaxation stress is regarded as the initial stress, and the additional force and displacement increment are calculated. Similarly, the stiffness matrix remains unchanged in each period. Since the soil is more suitable for relaxation test, the initial stress method has an advantage over the initial strain method in the soil-structure complex interaction finite element analysis.

[0003] The structure and the site are two inseparable factors. The structure is located on or in the site. Regardless of the case, the structure will be affected by the site. This effect is called soil-structure complex interaction effect, that is, the existence of soil will cause the upper part of the site to move differently from the bedrock surface. The response of the upper structure of the site will be very different from that under the assumption of rigid foundation. A large number of studies have shown that soil-structure complex interaction will affect the vibration mode of the structure, and make it produce greater response under dynamic load. Therefore, the rapid finite element simulation of soil-structure complex dynamic interaction is beneficial to rapid and effective earthquake damage assessment work when an earthquake occurs, and is beneficial to the timely issuance of disaster relief decisions, and has great significance for natural disaster emergency response.

[0004] At present, the commonly used constitutive relations of soil in commercial finite element software are viscoelastic constitutive relation and elastoplastic constitutive relation, and the structure is usually simulated by using the elastoplastic constitutive relation of concrete. However, these different constitutive relations are distributed in the tabs under the material properties. If the material constitutive relation needs to be changed, the parameters need to be adjusted and divided again, which leads to the cumbersome operation of modifying the model and the weak reusability of the model.

[0005] Patent CN117313480A discloses a finite element total stiffness matrix integration method for electric power pipe jacking tunnel, which comprises the following steps: firstly, the geometric size and material information data of the target electric power pipe jacking tunnel structure are obtained, the target electric power pipe jacking tunnel structure is discretized, the node coordinate matrix, the unit vector matrix and the node displacement coding matrix are constructed, the target structure is modeled by finite element, the unit stiffness matrix is calculated according to the knowledge of structural mechanics, the total stiffness matrix of the structure is integrated based on the unit stiffness matrix, the total stiffness equation is solved through the total stiffness matrix, and the node displacement array of the target electric power pipe jacking tunnel structure is solved, and the internal force of the target electric power pipe jacking tunnel mechanism is solved. Although the patent discloses a finite element total stiffness matrix integration method, the method cannot consider and solve the complex interaction problem of soil-structure, and the complex interaction of soil-structure has a significant influence on the structural response.

[0006] Patent CN108287945A discloses a deformation calculation method and application technology of foundation soil under large foundation, which comprises the following steps: creating a plane grid system, specifically, dividing the foundation under the test area into a plurality of same subdomains on the base elevation plane, when the number of subdomains is sufficient, the effective additional stress acting on each subdomain can be equivalent to a concentrated force, and the action point of the concentrated force is located at the center of the subdomain; establishing a foundation soil deformation calculation model; calculating the compression modulus Es of the foundation soil and solving the vertical deformation value. Although the patent solves the defect of single-point calculation by dividing the subdomain and applying the equivalent concentrated force, the method cannot consider the contact problem between the foundation and the foundation soil, and cannot meet the needs of seismic damage assessment. SUMMARY

[0007] The purpose of the present application is to overcome at least one of the above-mentioned defects in the prior art, and to provide a soil-structure complex interaction finite element analysis method and device, which has high calculation efficiency for the constitutive of soil and structure, can meet the needs of rapid seismic damage assessment under multiple scenarios on the one hand, and has strong portability on different research problems on the other hand.

[0008] The purpose of the present application can be achieved by the following technical solutions:

[0009] One of the technical solutions of the present application is to provide a soil-structure complex interaction finite element analysis method, which comprises the following steps:

[0010] S1, calling relevant material parameters from a material information database according to material naming;

[0011] S2, meshing the soil body and the structure model, setting node information and element information, and extracting material properties for calculating element stiffness matrix [D] from a material information database and assembling total stiffness matrix [K], and calculating total flexibility matrix [K] -1 ;

[0012] S3, determining the actual distance {δ} between the contact surfaces with plastic displacement according to the plastic displacement {u p} between the soil and the structure;

[0013] S4, calculating the contact node force {P p} of each element by a linear programming calculation method with contact conditions; e ;

[0014] S5, calculating element node displacement {u} and plastic displacement {u p} by the initial stress method;

[0015] S5.1, calculating element node displacement {u} and plastic displacement {u p} of the model by displacement iteration according to the total initial stress node load {P p}, and determining whether the accuracy requirement is met by calculating the relative difference percentage before and after iteration;

[0016] S5.2, if the relative difference percentage meets the accuracy requirement, calculating the strain {ε} and stress {σ} of the model based on the geometric equation (describing the relationship between element displacement and strain) and the physical equation (describing the relationship between element stress and strain), and outputting the strain {ε} and stress {σ};

[0017] S5.3, if the relative difference percentage does not meet the accuracy requirement, extracting the plastic displacement {u p} between the soil and the structure, revising the actual distance {δ} between the contact surfaces in step S3, and performing a new round of iteration calculation, which can improve the accuracy of the actual distance {δ} between the contact surfaces, thereby improving the accuracy of the element node displacement {u} and the plastic displacement {u p}.

[0018] As a preferred technical solution, in step S1, the material information database is established based on one or more of material property tests, literature data sources, and material naming is performed.

[0019] As a preferred technical solution, in step S1, the material property test of the soil body includes one or more of dynamic triaxial test and resonance column test, and the material property test of the structure includes one or more of apparent density test, compressive strength test, and splitting tensile strength test.

[0020] Further, the viscoelastic material in step S2 obtains the element stiffness matrix [D] according to the elasticity mechanics physical equation, and the elasticity mechanics physical equation is as follows:

[0021]

[0022] Wherein, E is the elastic modulus, μ is the poisson's ratio, σ x 、σ y 、σ z is the normal stress of the element, τ xy 、τ yz 、τ zx is the shear stress of the element, ε x 、ε y 、ε z is the normal strain of the element, γ xy 、γ yz 、γ zx is the shear strain of the element.

[0023] Further, the elastoplastic material in step S2 obtains the element stiffness matrix [D] according to the elastoplastic mechanics stress-strain relationship, and the elastoplastic mechanics stress-strain relationship is as follows:

[0024]

[0025] Wherein, E D ' is the equivalent elastic modulus after entering the plasticity, σ ep is the average stress, ε ep is the average strain, σ x 、σ y 、σ z is the normal stress of the element, τ xy 、τ yz 、τ zx is the shear stress of the element, ε x 、ε y 、ε z is the normal strain of the element, γ xy 、γ yz 、γ zx is the shear strain of the element.

[0026] Further, the actual distance {δ} between the contact surfaces in step S3 is calculated according to the following formula:

[0027] {δ}={δ0}+{u p} (3)

[0028] Wherein, {δ} is the actual distance between the contact surfaces, {δ0} is the initial distance between the contact surfaces, {u p} is the plastic displacement between the soil and the structure, and the plastic displacement of the soil body {u p-soil} and the plastic displacement of the structure {up-structure} is 0, and the initial distance {δ0} between the contact surfaces is 0.

[0029] As a preferred technical solution, the initial plastic displacement {u p} between the soil and the structure in step S3 is 0, and the initial distance {δ0} between the contact surfaces is 0.

[0030] Further, in step S4, based on the continuity condition of the contact problem, the balance condition, the contact criterion and the non-negativity constraint condition, the contact node force {P p} e .

[0031] Further, the continuity condition of the contact surface with the plastic displacement is as follows:

[0032] -[F] e {P p} e +a{e}+[I]{δ1}={δ} (4)

[0033] Where [F] e is the flexibility matrix of the contact element, {P p} e is the contact node force of each element, α is the rigid body approach amount, {e} is the unit column vector, [I] is the unit matrix, {δ1} is the actual distance between the contact surfaces after being stressed, and {δ} is the actual distance between the contact surfaces.

[0034] Further, the balance condition of the contact problem is as follows:

[0035] {e} T {P p} e =P (5)

[0036] Where {e} T is the unit row vector, {P p} e is the contact node force of each element, and P is the normal external load.

[0037] The contact criterion is as follows:

[0038]

[0039] Where P p is the single contact node force, and δ k is the distance between the contact nodes.

[0040] The non-negativity constraint condition is as follows:

[0041] P p ≥ 0, δ k ≥ 0, and α ≥ 0 (7)

[0042] where P p is the single contact node force, δ k is the distance between the contact nodes, and α is the rigid body approach.

[0043] Further, the contact node force {P p} of each element is obtained in step S5. e Then, the total initial stress node load {P p} is obtained by superposition.

[0044] Further, the displacement result iteration formula in step S5 is as follows:

[0045] {u} = {u e} = {u p} = [K] -1 {P} + [K] -1 {P p} (8)

[0046] where {u} is the total element node displacement, {u e} is the elastic displacement, {u p} is the plastic displacement, [K] -1 is the total flexibility matrix, P is the normal external load, and {P p} is the initial stress node load.

[0047] As a preferred technical solution, if the relative difference percentage does not meet the accuracy requirement in step S5.3, the strain {ε} and stress {σ} of the model are calculated based on the geometric equation and the physical equation, the element properties of the soil and the structure for the new round of iteration calculation are re-determined based on the elastic strain and the plastic strain in the strain {ε} according to the calculation condition, so as to improve the calculation accuracy, and the corresponding plastic displacement {u p} is obtained according to the plastic strain of the plastic element.

[0048] One of the technical solutions of the present application is to provide a soil-structure complex interaction finite element analysis device, which realizes the method, and comprises a parameter extraction module, a unit division module, an operation module, a solving module, an iteration module and a transmission module, the parameter extraction module generates a material information database for soil and structure models and extracts material parameters, the unit division module divides units for soil and structure models and generates unit stiffness matrix [D], the operation module performs intermediate calculation on data generated by the method itself, the solving module performs optimization calculation on data generated by the method itself under constraint conditions, the iteration module performs convergence comparison on data generated by the method itself under accuracy requirements and returns to the corresponding step for recalculation, and the transmission module inputs and outputs the collected data and data generated by the method itself.

[0049] One of the technical solutions of the present application is to provide a soil-structure complex interaction finite element analysis device, which realizes the method, and comprises a parameter extraction module, a unit division module, an operation module, a solving module, an iteration module and a transmission module, the parameter extraction module generates a material information database for soil and structure models and extracts material parameters, the unit division module divides units for soil and structure models and generates unit stiffness matrix [D], the operation module performs intermediate calculation on data generated by the method itself, the solving module performs optimization calculation on data generated by the method itself under constraint conditions, the iteration module performs convergence comparison on data generated by the method itself under accuracy requirements and returns to the corresponding step for recalculation, and the transmission module inputs and outputs the collected data and data generated by the method itself.

[0050] One of the technical solutions of the present application is to provide a computer readable storage medium, which stores a computer program, and the computer program is executed by a processor to realize the steps of the method.

[0051] One of the technical solutions of the present application is to provide a computer program product, which comprises a computer program, and the computer program is executed by a processor to realize the steps of the method.

[0052] Compared with the prior art, the present application has the following beneficial effects:

[0053] (1) The total stiffness matrix of the soil-structure interaction system based on the initial stress method in the present application does not change in the calculation process, and the purpose of iterative approximation to the final solution is achieved by modifying the force on the right side of the balance equation; the present application considers the characteristics of the balance equation of the soil-structure interaction system in the complex interaction process, and has the advantages of clear theory and controllable calculation accuracy;

[0054] (2) Compared with the variable stiffness method commonly used in existing commercial finite element software, the soil-structure complex interaction finite element analysis method based on the initial stress method developed in the application realizes the coupling iterative calculation of the initial stress method and the contact theory by combining the initial stress method with the contact theory, and improves the calculation efficiency of the model; the application can automatically iterate subsequent calculations according to the required calculation accuracy and material constitutive parameters, effectively improve the calculation accuracy and efficiency of the soil-structure multi-body material interaction, realize the rapid finite element analysis of the large-scale multi-degree-of-freedom soil-structure complex dynamic interaction system, which is beneficial to the rapid damage assessment and decision-making after the earthquake, and is of great significance to the emergency response work of natural disasters;

[0055] (3) The soil-structure complex interaction finite element analysis method based on the initial stress method in the application realizes the rapid adjustment and calculation of multiple interaction problems such as viscoelastic-viscoelastic, viscoelastic-elastic-plastic, and elastic-plastic-elastic-plastic through the encapsulation and visualization of the pre-processing, based on the stiffness matrix of the material, only needs to directly call different material information, through further expansion of the material information library, can give the method strong portability in different research problems, realizes the non-conversion input of the soil material property test and the soil-structure complex interaction finite element analysis calculation parameters, meets the efficient iterative calculation based on the pre-designed calculation accuracy, and significantly improves the calculation accuracy and efficiency. BRIEF DESCRIPTION OF DRAWINGS

[0056] Figure 1 It is a flowchart of the soil-structure complex interaction finite element analysis method in the embodiment of the application;

[0057] Figure 2 It is a boundary condition and load diagram of the finite element model in the embodiment of the application;

[0058] Figure 3 It is a stress distribution diagram after the finite element analysis of the finite element model in the embodiment of the application;

[0059] Figure 4 It is a strain distribution diagram after the finite element analysis of the finite element model in the embodiment of the application. DETAILED DESCRIPTION

[0060] The application will be described in detail below in combination with specific embodiments. The embodiment is implemented on the premise of the technical scheme of the application, and detailed implementation modes and specific operation processes are given, but the protection scope of the application is not limited to the following embodiments.

[0061] Embodiment:

[0062] A soil-structure complex interaction finite element analysis method, which is based on the initial stress method of soil-structure interaction to calculate stress distribution, such as Figure 1The specific steps are as follows:

[0063] S1, based on material property test (dynamic triaxial test, resonance column test, apparent density test, compressive strength test, splitting tensile strength test, etc.), literature data source, establish material information database, and carry out material naming, determine that the soil body adopts Shanghai soft soil material information, the structure adopts C60 concrete material information, and according to the material naming, call the related material parameters from the material information database;

[0064] S1-1, the Shanghai soft soil establishes the material information database in a mixed form based on the dynamic triaxial test and the literature data source. The density of the Shanghai soft soil is 1900kg / m 3 , the elastic modulus is 42MPa, and the Poisson's ratio is 0.33;

[0065] S1-2, the C60 concrete establishes the material information database in a mixed form based on the apparent density test, the compressive strength test, the splitting tensile strength test and the literature data source. The density of the C60 concrete is 2350kg / m 3 , the elastic modulus is 36GPa, and the Poisson's ratio is 0.2;

[0066] S2, as shown in Figure 2 , the model is determined as the left boundary of the structure being fixed, the right boundary of the soil body being applied with a uniform load of 300kN, the soil body and the structure model are meshed, the node information and the element information are set, the material properties are extracted from the material information database, the element stiffness matrix [D] is calculated, the total stiffness matrix [K] is assembled, and the total flexibility matrix [K] -1 is calculated;

[0067] S2-1, the element stiffness matrix calculation method of the viscoelastic material and the elastoplastic material is as follows:

[0068] S2-1-1, the viscoelastic material obtains the element stiffness matrix [D] according to the elasticity mechanics physical equation, and the elasticity mechanics physical equation is as follows:

[0069]

[0070] Wherein, E is the elastic modulus, μ is the Poisson's ratio, σ x , σ y , σ z are the normal stresses of the element, τ xy , τ yz , τ zx are the shear stresses of the element, ε x , ε y , ε z are the normal strains of the element, and γ xy , γ yz , γ zx are the shear strains of the element.

[0071] S2-1-2, the elastic-plastic material obtains the element stiffness matrix [D] according to the stress-strain relationship of elastic-plastic mechanics, and the stress-strain relationship of elastic-plastic mechanics is as follows:

[0072]

[0073] Wherein, E D is the equivalent elastic modulus after entering the plastic stage, σ ep is the average stress, ε ep is the average strain, σ x , σ y , σ z is the normal stress of the element, τ xy , τ yz , τ zx is the shear stress of the element, ε x , ε y , ε z is the normal strain of the element, γ xy , γ yz , γ zx is the shear strain of the element.

[0074] S2-2, the soil body is set as a viscoelastic material, and a total of 20000 elements are divided, the soil body element stiffness matrix [D soil ] is constructed, the soil body element stiffness matrix [D soil ] is assembled into the total stiffness matrix [K soil ] based on the element number, and the total flexibility matrix [K soil ] -1 is calculated.

[0075] S2-3, the structure is set as an elastic-plastic material, and a total of 20000 elements are divided, the structure element stiffness matrix [D structure ] is constructed, the structure element stiffness matrix [D structure ] is assembled into the total stiffness matrix [K structure ] based on the element number, and the total flexibility matrix [K structure ] -1 is calculated.

[0076] S3, according to the plastic displacement {u p} between the soil and the structure, the actual distance {δ} between the contact surfaces with the plastic displacement is determined.

[0077] S3-1, the actual distance {δ} between the contact surfaces is calculated according to the following formula:

[0078] {δ}={δ0}+{u p} (3)

[0079] Among them, {δ} is the actual distance between the contact surfaces, {δ0} is the initial distance between the contact surfaces, {u p} is the plastic displacement between soil and structure, and {u p-soil} and the plastic displacement of the structure {u p-structure} sum;

[0080] S3-2. Soil is a viscoelastic material and does not have plastic displacement. Therefore, the plastic displacement between soil and structure {u p} is the plastic displacement of the structure {u p-structure};

[0081] S3-3. Initial plastic displacement between soil and structure {u p} is 0, and the initial distance between the contact surfaces {δ0} is 0;

[0082] S4. Calculate the contact node force {P p} e ;

[0083] S4-1. Based on the continuity conditions, equilibrium conditions, contact criteria and non-negativity constraints of the contact problem, the contact node force {P p} e ;

[0084] S4-1-1. The continuity conditions for the contact surface with plastic displacement are as follows:

[0085] -[F] e {P p} e +α{e}+[I]{δ1}={δ} (4)

[0086] Among them, [F] e is the flexibility matrix of the contact element, {P p} e is the contact node force of each element, α is the rigid body approach, {e} is the unit column vector, [I] is the unit matrix, {δ1} is the actual distance between the contact surfaces after the force is applied, and {δ} is the actual distance between the contact surfaces;

[0087] S4-1-2-1. The equilibrium conditions based on the contact problem are as follows:

[0088] {e} T {P p} e =P (5)

[0089] Among them, {e} T is a unit row vector, {P p} eP is the contact node force of each element, P is the normal external load;

[0090] S4-1-2-2, the contact criterion is as follows:

[0091]

[0092] P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount; p P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount; k P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount;

[0093] S4-1-2-3, the non-negativity constraint condition is as follows:

[0094] P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount; p P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount; k P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount;

[0095] P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount; p P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount; k P is the contact node force of each element, δ is the distance between the contact nodes, and α is the rigid body approach amount;

[0096] S5, the initial stress method is used to calculate the element node displacement {u} and plastic displacement {u p};

[0097] S5.1, the contact node force {P p} of each element is obtained e , and then the total initial stress node load {P p} is obtained by superposition;

[0098] S5.2, according to the total initial stress node load {P p}, the element node displacement {u} and plastic displacement {u p} of the model are obtained by displacement result iteration, and whether the accuracy requirement is met is determined by calculating the relative difference percentage before and after iteration;

[0099] S5.2-1, the displacement result iteration formula is as follows:

[0100] {u}={u e}+{u p}=[K] -1 {P}+[K] -1 {P p} (8)

[0101] {u} is the total element node displacement, {u e} is the elastic displacement, {u p} is the plastic displacement, [K] -1 is the total flexibility matrix, P is the normal external load, and {P p} is the initial stress node load;

[0102] S5.2-2, the soil is a viscoelastic material, and there is no plastic displacement, and the displacement result iteration calculation unit node displacement {u soil} is obtained.

[0103] S5.2-3, the structure is a viscoelastic material, and the displacement result iteration calculation unit node displacement {u structure} and plastic displacement {u p-structure} are obtained.

[0104] S5.2-4, the relative difference percentage of the unit node displacement {u soil} of the soil before and after iteration is calculated, and the relative difference percentage of the unit node displacement {u structure} and plastic displacement {u p-structure} of the structure before and after iteration is calculated.

[0105] S5.3, if the relative difference percentage meets the accuracy requirement, the strain {ε} and stress {σ} of the model are calculated based on the geometric equation (describing the relationship between the displacement and strain of the unit) and the physical equation (describing the relationship between the stress and strain of the unit), and the strain {ε} and stress {σ} are output.

[0106] S5.3-1, if the relative difference percentage meets the accuracy requirement, the strain {ε} and stress {σ} of the soil and the structure are calculated based on the geometric equation and the physical equation, respectively.

[0107] S5.4, if the relative difference percentage does not meet the accuracy requirement, the plastic displacement {u p} between the soil and the structure is extracted, the actual distance {δ} between the contact surfaces is corrected in step S3, and a new round of iteration calculation is performed. The accuracy of the actual distance {δ} between the contact surfaces can be improved through iteration, so as to improve the accuracy of the unit node displacement {u} and the plastic displacement {u p}.

[0108] S5.4-1, if the relative difference percentage does not meet the accuracy requirement, the strain {ε} and stress {σ} of the soil and the structure are calculated based on the geometric equation and the physical equation, respectively, and the unit properties of the soil and the structure for the new round of iteration calculation are determined again according to the elastic strain and plastic strain in the strain {ε}, so as to improve the accuracy of the calculation, and the corresponding plastic displacement {u p} is obtained according to the plastic strain of the plastic unit.

[0109] As shown in Figure 3 and Figure 4 , in the calculation result of the embodiment, the maximum stress of the soil is 0.64 MPa, the maximum strain is 4.05x10 -3 , the maximum stress of the structure is 0.90 MPa, and the maximum strain is 1.05x10-5 The calculation time of the embodiment is 49s. The average calculation error of stress is-3.8%, the average calculation error of strain is-3.3%, the calculation time is reduced by 26.5%, the stress and strain of the soil and structure calculated by the embodiment are close to the actual results, and the calculation time is also shortened to a certain extent compared with the conventional analysis method.

[0110] The above description of the embodiments is to facilitate those skilled in the art to understand and use the application. Those skilled in the art can easily make various modifications to the embodiments, and apply the general principles described herein to other embodiments without creative labor. Therefore, the application is not limited to the above embodiments, and the improvements and modifications made by those skilled in the art without departing from the scope of the application should be within the protection scope of the application.

Claims

1. A soil-structure interaction finite element analysis method, characterized by, The method comprises the following steps: S1, calling relevant material parameters from a material information database according to material naming; S2, meshing the soil and structure model, setting node information and element information, and extracting material properties for calculating element stiffness matrix [D] from the material information database and assembling total stiffness matrix [K], and calculating total flexibility matrix [K] -1 ; S3. determining the actual distance {δ} between the contact surfaces with plastic displacement {u p} between the contact surfaces with plastic displacement {u S4. Calculate the contact node force {P of each element by linear programming method with contact condition p} e ; S5, the initial stress method calculation unit node displacement {u} and plastic displacement {u p}; S5.1、According to the total initial stress node load {P p} and the plastic displacement {u p} of the model, the element node displacement {u} of the model is obtained by iteration of displacement results. The relative difference percentage before and after calculation is compared to determine whether the accuracy requirement is met. S5.2, if the relative difference percentage meets the accuracy requirement, calculating the strain {ε} and stress {σ} of the model, and outputting the strain {ε} and stress {σ}; S5.3, if the relative difference percentage does not satisfy the accuracy requirement, extract the plastic displacement {u} between the soil and the structure p} in step S3, the actual distance {δ} between the contact surfaces is revised again, and a new round of iterative calculation is performed.

2. The soil-structure interaction finite element analysis method according to claim 1, wherein In step S2, the viscoelastic material obtains the element stiffness matrix [D] according to the elasticity mechanics physical equation, and the elasticity mechanics physical equation is as follows: where E is the modulus of elasticity, μ is the Poisson's ratio, σ x , σ y , σ z is the normal stress of the unit, τ xy , τ yz , τ zx is the shear stress of the unit, ε x , ε y , ε z is the normal strain of the unit, γ xy , γ yz , γ zx is the shear strain of the unit.

3. The method of claim 1, wherein the method is characterized by: In step S2, the elastoplastic material obtains the element stiffness matrix [D] according to the elastoplastic mechanics stress-strain relationship, and the elastoplastic mechanics stress-strain relationship is as follows: where E' is the equivalent elastic modulus after entering plasticity, σ is the average stress, ε is the average strain, σ is the normal stress of the unit, τ is the shear stress of the unit, ε is the normal strain of the unit, and γ is the shear strain of the unit. D ep ep x y z xy yz zx x y z xy yz zx where E' is the equivalent elastic modulus after entering plasticity, σ is the average stress, ε is the average strain, σ is the normal stress of the unit, τ is the shear stress of the unit, ε is the normal strain of the unit, and γ is the shear strain of the unit.​​​​​​​​​​​​​​ 4. The soil-structure interaction finite element analysis method according to claim 1, wherein In step S3, the actual distance {δ} between the contact surfaces is calculated according to the following formula: {δ} = {δ0} + {u p} (3) where {δ} is the actual distance between the contact surfaces, {δ0} is the initial distance between the contact surfaces, {u p} is the plastic displacement between the soil and the structure, and {u p-soil} is the plastic displacement of the soil and {u p-structure} is the plastic displacement of the structure.

5. The soil-structure interaction finite element analysis method according to claim 1, wherein In step S4, the contact node force {P of each element is solved by the simplex method based on the continuity condition of the contact problem, the equilibrium condition, the contact criterion and the non-negativity constraint condition. p} e .

6. The soil-structure interaction finite element analysis method according to claim 5, wherein The continuous condition of the contact surface with plastic displacement is as follows: -[F] e {P p} e + α{e} + [I] {δ1} = {δ} (4) where [F] = [K] {δ} + {F} (1) e is the flexibility matrix of the contact elements, {P p} e is the contact node force of each element, a is the rigid body approach, {e} is the unit column vector, [I] is the identity matrix, {δ1} is the actual distance between the contact surfaces after the force is applied, and {δ} is the actual distance between the contact surfaces.

7. The method of claim 5, wherein the method is characterized by: The balance condition based on the contact problem is as follows: {e} T {P p} e = P (5) where {e} T is the unit row vector, {P p} e is the contact node force of each element, and P is the normal external load. The contact criterion is as follows: where P p is the force at a single contact node, δ k is the distance between contact nodes; The non-negativity constraint condition is as follows: P p ≥0,δ k ≥0,α≥0 (7) where P p is the force at a single contact node, δ k is the distance between contact nodes, and a is the rigid body approach.

8. The soil-structure interaction finite element analysis method according to claim 1, wherein The contact node forces {P of each element in step S5 p} e The total initial stress node loads {P are then obtained by superposition p} 9. The soil-structure interaction finite element analysis method according to claim 1, wherein In step S5, the displacement result iteration formula is as follows: {u} = {u e} + {u p} = [K] -1 {P} + [K] -1 {P p} (8) where {u} is the total nodal displacement, {u e} is the elastic displacement, {u p} is the plastic displacement, [K] -1 is the total flexibility matrix, P is the normal external load, {P p} is the initial stress nodal load.

10. A finite element analysis device for soil-structure complex interaction, characterized in that: The device realizes the method in any one of claims 1 to 9, and comprises a parameter extraction module, an element division module, an operation module, a solving module, an iteration module and a transmission module. The parameter extraction module generates a material information database and extracts material parameters for a soil body and a structure model. The element division module divides elements for the soil body and the structure model and generates an element stiffness matrix [D]. The operation module performs intermediate calculation on data generated by the method itself. The solving module performs optimization calculation on data generated by the method itself under constraint conditions. The iteration module performs convergence comparison on data generated by the method itself under accuracy requirements and returns to a corresponding step for recalculation. The transmission module inputs and outputs collected data and data generated by the method itself.

Images