Numerical implementation method of fluid-solid coupling in saturated porous media based on fractional flow model

By constructing and discretizing fractional-order fluid-structure interaction control equations, finite element discrete equations containing degrees of freedom of displacement and pore pressure are generated, solving the stability and scalability problems in the numerical implementation of fractional-order fluid-structure interaction, and realizing stable calculation and multi-field coupling analysis on the ABAQUS platform.

CN121543514BActive Publication Date: 2026-04-10NORTHEASTERN UNIV CHINA
View PDF 2 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Patents(China)
Current Assignee / Owner
Filing Date
2026-01-19
Publication Date
2026-04-10

AI Technical Summary

Technical Problem

Existing numerical methods for fractional-order fluid-structure interaction have shortcomings in terms of stability and multiphysics scalability, especially on the ABAQUS platform, where stable calculations and thermal-fluid-structure multi-field coupling analysis are difficult to achieve.

Method used

Fractional-order fluid-structure interaction control equations are constructed and discretized to generate finite element discrete equations containing degrees of freedom for displacement and pore pressure. User-defined element subroutines are written and numerical solutions are executed on the ABAQUS platform, avoiding the use of the temperature field module.

Benefits of technology

It improves the stability and convergence efficiency of numerical calculation, realizes robust and reliable calculation in complex engineering scenarios, and supports multi-field coupled analysis of thermal-fluid-solid processes.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN121543514B_ABST
    Figure CN121543514B_ABST
Patent Text Reader

Abstract

The application discloses a saturated porous medium fluid-solid coupling numerical implementation method based on a fractional order seepage model and relates to the technical field of numerical simulation. By constructing and discretizing the coupling control equation containing the fractional order derivative, the finite element discrete equation containing the two degrees of freedom of displacement and pore pressure which can be directly coded is generated, and the user-defined unit (UEL) subroutine is compiled accordingly. With the strict element-level discretization and linearization processing, the stability and convergence efficiency of numerical calculation are significantly improved. At the same time, since the pore pressure degree of freedom is directly introduced and solved, the original temperature field module is not occupied, and the ability of the model to be further expanded to real thermal-fluid-solid multi-field coupling analysis is completely retained. The inherent limitation of the traditional UMATHT-based interface scheme is fundamentally solved, the physical expandability and application potential of the model in complex engineering scenarios involving temperature changes are significantly enhanced on the basis of ensuring the robustness and reliability of the calculation process.
Need to check novelty before this filing date? Find Prior Art

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of numerical simulation, and in particular to a saturated porous medium fluid-structure coupling numerical implementation method based on a fractional order seepage model. BACKGROUND

[0002] The saturated porous medium fluid-structure coupling theory is the basis for multi-physical field coupling analysis in geotechnical engineering, petroleum geology and underground engineering, and has key application value in engineering practices such as foundation reinforcement, oil and gas exploitation, nuclear waste disposal and carbon dioxide geological storage. Through the establishment of the constitutive coupling relationship between the pore fluid and the solid skeleton, the theory can co-simulate the seepage process and deformation behavior, thereby realizing the prediction and evaluation of the evolution law of the engineering system.

[0003] In recent years, the theory of fractional calculus has been gradually introduced into the study of saturated porous medium seepage and fluid-structure coupling due to its ability to effectively describe the historical dependence, memory effect and non-local diffusion behavior of materials. By introducing the time fractional derivative, the creep, relaxation and other time-varying processes of the geotechnical medium under load can be described; and the introduction of the spatial fractional derivative can represent the abnormal diffusion phenomenon of the fluid in the heterogeneous porous medium, providing a more accurate mathematical tool for seepage analysis under complex geological conditions.

[0004] However, although the fractional order fluid-structure coupling theory has shown significant advantages in describing complex medium behavior, its efficient and stable numerical implementation at the engineering scale still faces significant challenges. Currently, existing research attempts to develop related computational models on the ABAQUS platform of the general finite element software, mainly by calling its user-defined material subroutine interface (such as UMAT or UMATHT). Although this method has a certain flexibility, it has obvious limitations in practical application: first, the numerical stability is generally insufficient. The implementation scheme based on UMAT or UMATHT usually has sensitivity in the integration algorithm, iterative convergence and fractional operator discretization, which can easily lead to problems such as oscillation, divergence or slow convergence in the calculation process, making it difficult to ensure the reliability of large-scale engineering calculations; second, the physical field expandability is severely restricted. In particular, the scheme using the UMATHT interface essentially uses the mathematical analogy between the seepage equation and the heat conduction equation to indirectly solve the pore pressure degree of freedom by equivalent replacement with the temperature degree of freedom. This shell implementation method can complete the fluid-structure coupling analysis, but it occupies the thermal analysis module of the software, making it impossible to further introduce the real temperature field, thereby blocking the technical path of heat-flow-solid multi-field coupling analysis and greatly limiting the application of the model in engineering scenarios involving temperature changes. SUMMARY

[0005] Therefore, the application provides a saturated porous medium fluid-solid coupling numerical implementation method based on a fractional seepage model to solve the technical problems of insufficient stability and limited multi-physical field expansion in the current fractional fluid-solid coupling numerical implementation.

[0006] In a first aspect, a saturated porous medium fluid-solid coupling numerical implementation method based on a fractional seepage model is provided, and the method comprises the following steps:

[0007] Obtaining geometric parameters, material parameters and boundary conditions of a porous medium to be analyzed, wherein the material parameters include an elastic modulus, a Biot coefficient, a fluid density, a storage coefficient, a fractional seepage coefficient and a fractional order;

[0008] Based on the material parameters, a fractional fluid-solid coupling control equation for describing the interaction between the deformation of the solid skeleton and the fractional seepage of the pore fluid in the porous medium is constructed, wherein the fractional fluid-solid coupling control equation includes a mechanical field control equation and a seepage field control equation based on a fractional derivative;

[0009] Based on a user-defined unit framework, a finite element discrete format derivation is performed on the fractional fluid-solid coupling control equation to generate a finite element discrete equation for compiling a user-defined unit subroutine, wherein the finite element discrete equation includes a first discrete equation derived from the mechanical field control equation and a second discrete equation derived from the seepage field control equation based on the fractional derivative;

[0010] Based on the geometric parameters and the boundary conditions, a numerical calculation model of the porous medium is constructed, and the numerical calculation model is subjected to finite element mesh division;

[0011] Based on the finite element discrete equation and the numerical calculation model, the user-defined unit subroutine is compiled, and a finite element input file containing a user-defined unit configuration is generated;

[0012] The finite element input file is associated with the user-defined unit subroutine, the fluid-solid coupling numerical solution is executed through an ABAQUS platform, and simulation results of the porous medium are obtained, wherein the simulation results include a displacement field, a pore water pressure field and a seepage rate.

[0013] In the scheme realized by the numerical implementation method of saturated porous medium fluid-solid coupling based on the fractional seepage model, the coupled control equation containing the fractional derivative is constructed and discretized to generate the finite element discrete equation containing the two degrees of freedom of displacement and pore pressure which can be directly coded, and the user-defined element subroutine is compiled accordingly. With strict element-level discretization and linearization processing, the stability and convergence efficiency of numerical calculation are significantly improved. At the same time, since the pore pressure degree of freedom is directly introduced and solved, the original temperature field module is not occupied, and the ability of further expanding the model to real thermal-fluid-solid multi-field coupling analysis is completely retained. The inherent limitations of the traditional UMATHT interface scheme are fundamentally solved, the physical expandability and application potential of the model in complex engineering scenarios involving temperature changes are significantly enhanced on the basis of ensuring the robustness and reliability of the calculation process, thereby providing reliable technical support for the fine simulation and safety evaluation of porous media in the fields of geotechnical engineering, groundwater resource assessment, energy geology and the like. BRIEF DESCRIPTION OF DRAWINGS

[0014] Various other advantages and benefits will become apparent to those of ordinary skill in the art upon reading the following detailed description of the preferred embodiments with reference made to the accompanying drawings. The drawings are for purposes of illustration only and are not intended to be limiting in accordance with the present application. In the drawings, like reference numerals have been used to indicate like elements throughout. In the drawings:

[0015] Figure 1 FIG. 1 is a flowchart of the numerical implementation method of saturated porous medium fluid-solid coupling based on the fractional seepage model in an embodiment of the present application;

[0016] Figure 2 FIG. 2 is a comparison chart of numerical calculation result accuracy in an embodiment of the present application;

[0017] Figure 3 FIG. 3 is a comparison chart of the change of numerical calculation convergence efficiency with the analysis step advancing in an embodiment of the present application;

[0018] Figure 4 FIG. 4 is a structural schematic diagram of the numerical implementation device of saturated porous medium fluid-solid coupling based on the fractional seepage model in an embodiment of the present application. DETAILED DESCRIPTION

[0019] In order to make the objectives, technical solutions and advantages of the embodiments of the present application clearer, the technical solutions in the embodiments of the present application will be described below in connection with the drawings of the embodiments of the present application. It should be understood that the drawings in the present application only serve the purpose of illustration and description, and do not serve to limit the protection scope of the present application.

[0020] In addition, it should be understood that the illustrative drawings are not drawn to scale. The flowcharts used in the present disclosure show the operations implemented according to some embodiments of the present disclosure. It should be understood that the operations of the flowcharts can not be implemented in order, and the steps without logical context relationship can be reversed in order or implemented simultaneously. In addition, one or more other operations can be added to the flowcharts or one or more operations can be removed from the flowcharts under the guidance of the present disclosure.

[0021] In addition, the embodiments described in the present disclosure are only some of the embodiments of the present disclosure, not all the embodiments. The components of the embodiments of the present disclosure described and shown in the drawings herein can be arranged and designed in various different configurations. Therefore, the following detailed description of the embodiments of the present disclosure provided in the drawings is not intended to limit the scope of the claimed present disclosure, but only represents selected embodiments of the present disclosure. Based on the embodiments of the present disclosure, all other embodiments obtained by those skilled in the art without creative labor fall within the scope of protection of the present disclosure.

[0022] It should be noted that the term "comprising" will be used in the embodiments of the present disclosure to indicate the presence of the features declared thereafter, but does not exclude the addition of other features. It should also be noted that similar reference numerals and letters represent similar items in the following drawings, so once an item is defined in one drawing, it does not need to be further defined and explained in subsequent drawings.

[0023] The present case will be described in detail below in conjunction with the drawings related to the specification.

[0024] Please refer to Figure 1 The embodiments of the present disclosure provide a saturated porous medium fluid-structure coupling numerical implementation method based on a fractional seepage model, specifically comprising the following steps:

[0025] S10: Obtain the geometric parameters, material parameters and boundary conditions of the porous medium object to be analyzed;

[0026] The material parameters include the elastic modulus, the Biot coefficient, the fluid density, the water storage coefficient, the fractional permeability coefficient and the fractional order.

[0027] It can be understood that the execution subject of the present disclosure can be a saturated porous medium fluid-structure coupling numerical implementation device based on a fractional seepage model, and can also be a terminal or a server, which is not limited here. The embodiments of the present disclosure take the server as the execution subject for example.

[0028] In this step, the porous medium refers to an engineering material and geological body under complex stress, temperature, and time history conditions, and its seepage-deformation behavior exhibits significant non-classical characteristics. The geometric parameters, material parameters, and boundary conditions of the porous medium to be analyzed are collected. Among them, the geometric parameters include the shape, size, and contour of the porous medium object in space; the material parameters include the elastic modulus, Biot's coefficient, fluid density, and water storage coefficient, as well as the fractional order permeability and fractional order number that describe the history dependence and abnormal diffusion characteristics of the material seepage behavior; the boundary conditions include the initial state conditions, mechanical boundary constraints, and seepage boundary constraints.

[0029] S20: Based on the material parameters, a fractional order fluid-solid coupling control equation is constructed to describe the interaction between the solid skeleton deformation and the fractional order seepage of the pore fluid in the porous medium.

[0030] Among them, the fractional order fluid-solid coupling control equation includes a mechanical field control equation and a seepage field control equation based on fractional order derivative.

[0031] In this step, the fractional order fluid-solid coupling control equation is a closed system of differential equations, including two core field equations that are independent of each other but are tightly coupled through the constitutive relationship: the mechanical field control equation and the seepage field control equation based on fractional order derivative. Among them, the mechanical field control equation describes the static equilibrium state of the solid skeleton; the seepage field control equation based on fractional order derivative replaces the traditional Darcy's law and is used to describe the seepage behavior with history dependence and non-local characteristics.

[0032] In an embodiment of the present application, the mechanical field control equation is constructed based on the law of conservation of momentum; the seepage field control equation based on fractional order derivative is constructed based on the law of conservation of mass. Assuming that the porous medium is an isotropic linear elastic body, its stress-strain relationship satisfies the linear elastic constitutive model, and the pore water velocity and pressure satisfy the fractional order seepage model. Then the mechanical field control equation is:

[0033] ;

[0034] Among them, is the stress tensor; is the i-th component of the body force vector;

[0035] ;

[0036] Among them, is the i, j component of the stress tensor; is the i, j component of the effective stress tensor; is the Biot's coefficient; is the pore water pressure; is the Kronecker symbol;

[0037] ;

[0038] ;

[0039] wherein, is the fourth order elastic stiffness tensor of the material; is the k, l component of the strain tensor.

[0040] The seepage field control equation based on the fractional derivative is:

[0041] ;

[0042] wherein, is the fluid density; is the storage coefficient of the porous medium; p is the pore water pressure; t is the time; q is the seepage velocity vector; is the mass source term;

[0043] ;

[0044] wherein, q is the seepage velocity vector; is the fractional order permeability coefficient; is the fluid density; g is the acceleration of gravity; is the fractional derivative operator; is the order of the derivative; a is the start time; t is the current time; is the gradient of the pore water pressure.

[0045] Optionally, the fractional derivative used is the Grünwald-Letnikov (GL) definition, and its specific discrete form is as follows:

[0046] ;

[0047] wherein, is the fractional derivative of order of the function ; is the time step; k is the total number of steps, and i is from 0 to k.

[0048] S30: based on the user-defined unit framework, finite element discrete format derivation is performed on the fractional fluid-solid coupling control equation to generate a finite element discrete equation for writing a user-defined unit subroutine; wherein the finite element discrete equation includes a first discrete equation derived from the mechanical field control equation and a second discrete equation derived from the seepage field control equation based on the fractional derivative.

[0049] In this step, under the user-defined unit framework, starting from the control equation describing the coupling of solid deformation and fractional order seepage of porous media, the finite element method is used to discretize the control equation in space and time, and it is converted into a discrete form that can be calculated numerically. Through systematic discretization derivation, the element-level control equation in the form of matrix and vector is finally obtained. This set of equations clearly gives the complete analytical expression of the core data structure that the user-defined unit subroutine needs to return to the finite element solver, so that the corresponding user-defined unit subroutine can be directly programmed according to the expression. The fractional order fluid-solid coupling theory is embedded in the finite element software platform ABAQUS, realizing stable, efficient and numerical analysis without occupying other physical field degrees of freedom (such as temperature field), and fundamentally overcoming the key technical defects of the traditional UMAT / UMATHT-based method, such as insufficient numerical stability and limited multi-physical field expandability.

[0050] In an embodiment of the present application, a specific discretization scheme is provided, and in S30, the finite element discretization format derivation of the fractional order fluid-solid coupling control equation is performed based on the user-defined unit framework, and the finite element discretization equation for writing the user-defined unit subroutine is generated, which specifically includes the following steps S31-S35:

[0051] S31: Obtain the preset node configuration of the finite element unit;

[0052] The preset node configuration includes the preset number of nodes of the finite element unit, and the preset displacement degree of freedom and the preset pore water pressure degree of freedom of each node.

[0053] In this step, the preset node configuration refers to the unit topology and freedom information that has been clearly specified based on the selected unit type before the numerical simulation is performed. Specifically, it includes: the preset number of nodes, i.e. the total number of nodes contained in this type of unit, for example, the commonly used two-dimensional eight-node quadrilateral unit, and the preset number of nodes is 8; and the preset degree of freedom of each node, i.e. the type and number of independent field variables activated and participating in the calculation at each node. In the fluid-solid coupling problem addressed in the present application, each node is preset with two categories of degrees of freedom: the preset displacement degree of freedom, which is used to describe the spatial motion of the solid skeleton at the node, and in a three-dimensional model, it usually includes the translation degrees of freedom along the x, y and z directions; and the preset pore water pressure degree of freedom, which is used to describe the pressure state of the pore fluid at the node; therefore, a unit with n nodes has a total degree of freedom number of n times the number of displacement degrees of freedom of each node plus the number of pore water pressure degrees of freedom.

[0054] S32: Based on the preset node configuration, construct a unit node variable array.

[0055] In this step, assuming that the number of finite element unit nodes is n, the arrangement order of displacement and pore water pressure in the unit node variable is as follows:

[0056] ;

[0057] Wherein, n is the node variable; e is the unit level; T is the transpose symbol; u, v, w are preset displacement degrees of freedom; is the displacement component of the nth node in the x direction; is the displacement component of the nth node in the y direction; is the displacement component of the nth node in the z direction; p is the preset pore water pressure degree of freedom; is the pore water pressure of the nth node.

[0058] S33: Divide the unit node variable array into a unit node displacement array and a unit node pore water pressure array.

[0059] In this step, in order to facilitate the derivation of the finite element unit stiffness matrix in the subsequent step, the unit node column array is divided into the following two parts:

[0060] Unit node displacement column array:

[0061] ;

[0062] Unit node pore water pressure column array:

[0063] ;

[0064] Wherein, d is displacement.

[0065] S34: Based on the virtual work principle, the finite element discrete format derivation of the mechanical field control equation is carried out under the user-defined unit framework, and the first discrete equation is generated.

[0066] In this step, in order to realize the numerical solution of the mechanical field control equation, it is necessary to convert it from a continuous differential form to a discrete algebraic form. The virtual work principle is adopted as the core mathematical tool in this application. The virtual work principle is a basic principle in mechanics, which states that for a system in equilibrium, the virtual work done by internal forces on any virtual displacement satisfying the kinematic constraints is equal to the virtual work done by external forces on the virtual displacement. Based on this principle, the mechanical field control equation is processed to complete the discretization process, and finally the first discrete equation corresponding to the mechanical field control equation is generated.

[0067] In one embodiment of the present application, a first discrete equation generation scheme is provided, in S34, i.e. under the user-defined unit framework, a first discrete equation is generated based on the virtual work principle for finite element discrete format derivation of the mechanical field control equation, specifically including the following steps S341-S343:

[0068] S341: based on the virtual work principle , the discrete integral equation is obtained:

[0069] ;

[0070] wherein, is the entire volume domain to be solved; is the volume element; S is the boundary; dA is the area element; is the displacement field vector; is the virtual displacement vector; is the strain tensor; is the virtual strain tensor; is the volume force vector acting on the object; is the traction force vector acting on the boundary S.

[0071] In this step, for the continuum in the mechanical field, it is assumed that it is subjected to volume force in the geometric domain , and is subjected to distributed force applied on the outer boundary S. An arbitrary virtual displacement field satisfying the geometric boundary condition is introduced. The mechanical field control equation is weightedly integrated on the calculation domain , and the equivalent integral weak form is obtained by combining the stress boundary condition.

[0072] S342: matrix relationship between virtual displacement, virtual strain and strain and element node variables is established;

[0073] The virtual displacement relationship is:

[0074] ;

[0075] wherein, is the virtual displacement vector at any point in the element e; is the shape function matrix of the element; is the node virtual displacement vector of the element e;

[0076] The virtual strain relationship is:

[0077] ;

[0078] wherein, is the virtual strain vector at any point in element e; B is the geometry matrix of the element;

[0079] The strain relationship is:

[0080] ;

[0081] wherein, is the strain vector at any point in element e.

[0082] S343: Substitute the matrix relationship into the discretized integral equation, and use the element node displacement variation increment arbitrary, eliminate the increment term, and obtain the first discrete equation in the form of finite element discretization:

[0083] ;

[0084] wherein, is the stiffness matrix of element e; is the node displacement vector of element e; is the coupling matrix of element e; is the node pore water pressure vector of element e; is the equivalent node force vector of the mechanical field corresponding to element e

[0085] ;

[0086] wherein, D is the elastic stiffness matrix of the material; B is the geometry matrix of the element;

[0087] ;

[0088] wherein, is the pore pressure influence matrix;

[0089] ;

[0090] wherein, N is the shape function matrix of the unit; F is the traction force vector.

[0091] For steps S342-S343, the standard steps of the finite element method are used to approximate the real displacement field u and the virtual displacement field with the element shape function matrix and the node displacement vector and its variation , that is, and . At the same time, the relationship between strain and node displacement and the virtual strain d. Substitute these relations into the integral weak form above, and introduce the linear elastic constitutive relation and the effective stress principle, and then rearrange and integrate. Since the virtual node displacement is arbitrary, in order to make the virtual work equation always hold, the coefficient in front of it must be zero. Thus, the discretized equation with the element node displacement and the node pore water pressure as the basic unknowns can be directly derived, which is the first discrete equation. This equation reveals the balance relationship between the solid deformation, the pore water pressure and the external load, and is the cornerstone for subsequent assembly of the overall system and numerical solution.

[0092] S35: Under the user-defined element framework, based on the variational principle, the finite element discretization format of the seepage field control equation based on fractional derivative is derived to generate the second discrete equation.

[0093] In this step, in order to realize the numerical solution of the seepage field control equation based on fractional derivative, it needs to be converted into a discrete algebraic system suitable for computer processing. The present application uses the variational principle as the core mathematical tool to complete the discretization process, and finally generates the second discrete equation corresponding to the seepage field control equation based on fractional derivative.

[0094] In an embodiment of the present application, a specific second discrete equation generation scheme is provided, S35, that is, under the user-defined element framework, based on the variational principle, the finite element discretization format of the seepage field control equation based on fractional derivative is derived to generate the second discrete equation, which specifically includes the following steps:

[0095] Construct the functional form of the seepage field control equation based on fractional derivative, and perform spatial interpolation of the pore water pressure based on the shape function;

[0096] The fractional order pressure gradient term in the seepage field control equation based on fractional derivative is time-discretized using the Grünwald-Letnikov finite difference definition, and is expressed as a weighted convolution series of the current time and the historical time;

[0097] Decouple the weighted convolution series into an instantaneous response term and a historical memory term;

[0098] Based on the extremum condition of the variational principle, the tangent stiffness matrix relationship representing the current time permeability characteristics of the element is derived as the second discrete equation.

[0099] In this embodiment, for the continuum problem of the seepage field, the pore water pressure is given on the boundary S to solve the problem of minimizing the following functional I:

[0100] ;

[0101] where I is a functional; is the component of the seepage velocity vector q in the x direction; is the component of the seepage velocity vector q in the y direction; is the component of the seepage velocity vector q in the z direction; is the storage coefficient; p is the pore water pressure; t is time; is the mass source term; is the fluid density;

[0102] The fractional seepage model is expressed as:

[0103]

[0104] where, is the seepage velocity vector within element e; is the component of the fractional permeability in the x direction; is the component of the fractional permeability in the y direction; is the component of the fractional permeability in the z direction; is the gravitational acceleration; is the fractional derivative operator;

[0105] The numerical progressive format of the fractional derivative is:

[0106]

[0107] where, is the value of the fractional derivative of the function f at time t is the discrete time step; i is the time step; k is the total number of time steps; is the weight coefficient; is the value of the function f at the historical time

[0108] The definition of the weight function is:

[0109]

[0110] where, is the weight coefficient; i is the serial index of the weight coefficient, corresponding to the number of steps of time backtracking; is the fractional order;

[0111] The matrix of the fractional seepage model is:

[0112]

[0113] For the finite element unit, the relationship between the spatial derivative of the pore water pressure and the pore water pressure matrix of the unit node is: ​​​​​​

[0114] ;

[0115] where, / 、 / 、 / are the derivatives of shape function matrix with respect to x, y, z direction respectively;

[0116] Substitute the matrix relation into the fractional seepage model matrix, and get the discrete expression of seepage velocity vector:

[0117]

[0118] Take the extreme value of functional I, that is, let , get the second discrete equation corresponding to the seepage field control equation based on fractional derivative:

[0119] ;

[0120] where, is the seepage matrix of unit e; is the equivalent node force vector of seepage field of unit e;

[0121] ;

[0122] .

[0123] In this embodiment, for the seepage constitutive equation containing fractional derivative and the mass conservation equation, first, a functional related to the seepage problem is constructed. The variational principle points out that the distribution of pore water pressure field in the true physical state should make the functional take extreme value. By taking the first variation of the functional and letting it equal to zero, the equivalent integral weak form of the original seepage control equation can be naturally derived. Then, the finite element method is used for spatial discretization. For the fractional derivative term in the equation, the Grünwald-Letnikov definition is used for numerical discretization, which is converted into the weighted sum form of the pressure gradient value at the current and historical time steps. Since the arbitrariness of the variation , to satisfy the variational equation always holds, the sum of all coefficients before must be zero. Through this step and combined with the processing of the time derivative term in the mass conservation equation, the discrete equation with node pore water pressure as the basic unknown quantity is finally derived, that is, the second discrete equation. This equation and the first discrete equation are combined to form a complete discrete equation set describing the fractional fluid-structure coupling problem, which is the direct object of numerical solution.

[0124] S40: Constructing a numerical calculation model of the porous medium based on the geometric parameters and the boundary conditions.

[0125] In this step, after the mathematical construction and discretization derivation of the theoretical model are completed, the engineering implementation phase is entered, and the obtained geometric parameters and boundary conditions are converted into a numerical calculation model that can be recognized and processed by a computer according to the specific engineering problem.

[0126] Specifically, based on the geometric parameters, the physical form of the calculation domain is defined in the finite element pre-processing software; the geometric parameters determine the spatial framework of the model, for example, for foundation settlement analysis, a three-dimensional entity model of soil layer distribution is established according to the survey drawing, and for a simplified verification example, a two-dimensional plane region of specified size can be created. The boundary conditions are applied to the geometric model: the initial conditions are the starting state of the physical process; the mechanical boundary conditions and the seepage boundary conditions are respectively assigned to the corresponding surfaces or edges of the model, which together define the interaction of the system with the external environment and the driving force of the internal evolution. The geometric model is divided into finite elements, and the continuous calculation domain is discretized into a finite number of interconnected subunits using selected element types, forming a grid system composed of nodes and elements. The analysis type, time step solution control parameters, and the above-mentioned geometric, boundary, and grid information are integrated to form a complete numerical calculation model.

[0127] Through the above-mentioned way, the transformation from engineering parameters to standardized digital model is realized, providing complete geometric and boundary data carriers for subsequent generation of finite element input files containing user-defined element information, ensuring the consistency of numerical simulation and engineering practice.

[0128] S50: Compiling user-defined element subroutines based on the finite element discrete equation set and the numerical calculation model, and generating finite element input files containing user-defined element configurations.

[0129] In this step, first, according to the specific calculation formulas of the element stiffness matrix, coupling matrix, and load vector defined by the discrete equation, the subprogram source code conforming to the user-defined element interface specification of ABAQUS is written using Fortran language, and it is converted into an object file that can be called by the solver through the compiler; at the same time, based on the numerical calculation model that has been constructed containing geometric, grid, and boundary conditions, a standard text format input file is generated, and key instructions are embedded in this file to declare the type of user-defined element, node freedom degree, material parameter association, and element set definition, thereby generating a complete input file that can describe specific physical problems and guide the software to call external custom calculation modules.

[0130] Through the above manner, the derived mathematical algorithm and the established digital model are converted into computer executable instructions and structured data files, and deep integration of advanced fractional order theory and industrial standard simulation platform is realized.

[0131] In one embodiment of the present application, a specific UEL subroutine code compiling and finite element input file generation scheme is provided, in S50, that is, based on the finite element discrete equation and the numerical calculation model, the user-defined unit subroutine is compiled, and the finite element input file containing the user-defined unit configuration is generated, specifically including the following steps S51-S53:

[0132] S51: generating a finite element input file based on the numerical calculation model;

[0133] The finite element input file includes the geometric parameters, meshing information, boundary conditions and material parameters of the numerical calculation model.

[0134] In this step, the constructed numerical calculation model is taken as a data source, and the model has integrated complete information such as geometric shape, finite element mesh, boundary constraint and material property. The internal data model is exported to a standardized text format file, that is, a finite element input file, by using the pre-processing function of the finite element software. The file systematically records all elements of the numerical calculation model by using specific keyword syntax, including geometric control point coordinates, element-node connection relationship, definition position and numerical value of various boundary conditions, and preliminary identification of material properties, which constitutes the original data basis for subsequent calculation.

[0135] S52: compiling a user-defined unit subroutine based on the finite element discrete equation;

[0136] The user-defined unit subroutine is used to calculate the stiffness matrix and residual component of the unit.

[0137] In this step, the mathematical algorithm specified by the first discrete equation and the second discrete equation derived is taken as the basis, and a source code is written by using Fortran programming language. The core function of the subroutine is to accurately calculate and return the element stiffness matrix and element residual vector according to the input element node coordinates and material parameters. This subroutine is the core execution module for realizing the calculation function of the fractional order fluid-structure coupling model.

[0138] S53: in the finite element input file, the identification type and degree of freedom of the user-defined unit are declared, the element set using the user-defined unit is defined, and the material parameter list called by the user-defined unit subroutine is written.

[0139] In this step, in order to make the solver correctly call the aforementioned compiled subroutine in the calculation process, the generated original input file must be modified. Specifically, it includes the following steps: by inserting specific control statements, declare a new custom element type using the keyword, assign a unique identifier to it, and specify its node number, geometric dimension, and activated degrees of freedom; using the keyword, redefine the type of part or all of the elements specified in the model as the custom element, thereby forming a calculation set that uses custom elements; using the keyword, associate the material parameters necessary for the calculation with the element set in the form of an ordered list, ensuring that these parameters can be accurately passed to the subroutine. After the above configuration, the original general input file is converted into an instruction file that can not only describe the specific physical problem, but also explicitly instruct the solver to call the external custom algorithm module.

[0140] S60: associate the finite element input file with the user-defined element subroutine, perform fluid-structure coupling numerical solution through the ABAQUS platform, and obtain the simulation results of the porous medium;

[0141] The simulation results include displacement field, pore water pressure field, and seepage rate.

[0142] In this step, the generated input file containing model information and custom element configuration is bound with the compiled user-defined element subroutine, and submitted to the ABAQUS solver for calculation. In this process, the ABAQUS main program calls the external user-defined element subroutine to calculate the stiffness matrix and load contribution of the specified element set when solving the fluid-structure coupling system equation set according to the instructions in the input file, thereby completing the numerical solution of the entire system. The final simulation results include the displacement field of the solid skeleton, the pressure field of the pore fluid, and the seepage rate field reflecting the fractional order seepage characteristics, providing comprehensive quantitative data for engineering analysis and evaluation.

[0143] In an embodiment of the present application, a specific numerical simulation scheme is provided. In S60, the finite element input file is associated with the user-defined element subroutine, and the fluid-structure coupling numerical solution is performed through the ABAQUS platform to obtain the simulation results, which specifically includes the following steps S61-S62:

[0144] S61: create an analysis job in the ABAQUS platform, and associate the finite element input file with the user-defined element subroutine in the analysis job.

[0145] In this step, a new analysis job is created in the job module of ABAQUS. During the setting of this job, the finite element input file that has been configured is associated with the compiled user-defined element subroutine object file by specifying the input file path in the job properties, checking the user subroutine option, and pointing to the corresponding target file path. This ensures that the solver can accurately read the model data containing the custom element instructions and dynamically link to the external computing module that implements the fractional order fluid-structure coupling algorithm when executing the job.

[0146] S62: Submit the analysis job to the ABAQUS solver to perform fluid-structure coupling numerical calculation by calling the user-defined element subroutine, and obtain simulation results including displacement field, pore water pressure field, and seepage rate.

[0147] In this step, after the job creation and association are completed, the analysis job is submitted to the ABAQUS solver, which reads the input file and parses the model information and custom element configuration. During the assembly of the overall system matrix and the balance iteration solution, for the set of elements defined in the file that use the user-defined element, the solver will call the associated UEL subroutine. For each such element, the solver passes node coordinates, current displacement and pressure solution estimates, material parameters, and other data to it, and the UEL subroutine calculates the stiffness matrix and residual vector of the element in real time based on the built-in algorithm derived from the first and second discrete equations, and returns them to the solver. Through this internal and external collaboration, the nonlinear iterative solution of the entire fluid-structure coupling system is completed. After the calculation converges, the solver outputs the result data, and the final simulation results are a complete data field, including: displacement field reflecting the deformation distribution of the solid skeleton in space; pore water pressure field describing the spatial and temporal variation of pore fluid pressure; and seepage rate field, which is directly calculated by the fractional order seepage constitutive model and quantitatively represents the fluid flow velocity with historical dependence and abnormal diffusion characteristics.

[0148] Optionally, the final simulation results are written into a standard result file, which can be visualized, extracted, and analyzed through the post-processing module of ABAQUS, providing comprehensive and quantitative engineering basis for evaluating the mechanical stability, seepage safety, and long-term behavior of porous media structures.

[0149] Through the above method, based on the ABAQUS platform and calling its UEL interface for secondary development, direct numerical solution of the fractional order fluid-structure coupling model is realized, which can independently complete fluid-structure coupling calculation without relying on the temperature module in ABAQUS, has good numerical stability, and provides an effective and feasible implementation path for further expansion of the model to thermal-water-force multi-field coupling analysis.

[0150] In practical application scenarios, the specific implementation of the present application is realized by performing standardized modeling, programming and calculation processes in the ABAQUS software environment. First, in ABAQUS, a numerical model containing geometry, material, boundary conditions and mesh is created according to the conventional modeling steps. Subsequently, a new analysis job is created in the "Job" module, and a standard INP format input file is exported through the "Write Input" command. To ensure that the model can correctly call the user-defined fractional order fluid-structure interaction calculation module, necessary modifications need to be made to the INP file, mainly including: (1) using the corresponding keyword to declare the type of UEL custom element, specifying its node number, geometric dimension and activated degrees of freedom (displacement and pore water pressure); (2) defining the element set that uses the custom element in the file, and specifying which elements will use UEL for calculation; (3) writing all the model parameter lists necessary for UEL subroutine execution in the file, ensuring that the parameter order is consistent with the subroutine reading logic. Through the above configuration, the INP file is converted into an instruction file that can recognize and prepare to call the external UEL subroutine. Subsequently, based on the finite element discretization format of the derived governing equation, the UEL subroutine code is written using Fortran language. The implementation of the subroutine mainly covers the following contents: (1) declaring all necessary variables, arrays and matrices, including data structures for storing historical states; (2) correctly reading the associated material parameters from the INP file; (3) determining the position and weight of the Gauss integration point according to the selected element type, and calculating the corresponding shape function matrix and its derivative; (4) looping at the Gauss integration point, calculating the local stiffness contribution according to the discretization format formula, and integrating to assemble the complete element stiffness matrix (AMATRX); (5) estimating the element residual vector (RHS) based on the current solution. After the code is written, it needs to be compiled into an executable object file using a Fortran compiler compatible with ABAQUS. Finally, return to the ABAQUS environment and create a new analysis job (Job). In the attribute setting of this job, associate the previously modified INP file and the compiled UEL subroutine object file. After submitting this job, the ABAQUS solver will automatically perform the calculation. After the calculation is completed, the output result file (such as.odb file) is post-processed, and the displacement field, pore water pressure field and seepage rate field are extracted and analyzed to verify the correctness and calculation performance of the fractional order fluid-structure interaction model numerical implementation method.

[0151] It can be seen that in the above scheme, by constructing and discretizing the coupled control equation containing fractional derivative, a finite element discrete equation containing two degrees of freedom of displacement and pore pressure is generated, which can be directly coded, and a user-defined unit subroutine is compiled accordingly. With strict element-level discretization and linearization processing, the stability and convergence efficiency of numerical calculation are significantly improved. At the same time, since the pore pressure degree of freedom is directly introduced and solved, the original temperature field module is not occupied, and the ability to further expand the model to real thermal-fluid-solid multi-field coupling analysis is fully retained. The inherent limitations of the traditional UMATHT-based interface scheme are fundamentally solved, the physical expandability and application potential of the model in complex engineering scenarios involving temperature changes are significantly enhanced on the basis of ensuring the robustness and reliability of the calculation process, thereby providing reliable technical support for the fine simulation and safety evaluation of porous media in the fields of geotechnical engineering, groundwater resource assessment, energy geology, etc.

[0152] In actual application scenarios, the application is further described in combination with specific embodiments. Taking a two-dimensional eight-node element as an example, the UEL subroutine code is written, and a numerical model is established in the ABAQUS platform to verify the correctness of the compiled UEL subroutine. The model size is 1.2m*0.2m, the initial condition of pore water pressure is 60Pa, the upper boundary pore pressure is 0Pa, and the remaining boundaries are set as impermeable conditions. In order to verify the results, displacement constraints are applied to the bottom and both sides, so that the model is equivalent to a one-dimensional seepage deformation process in calculation.

[0153] Step 1: Establish the fractional order fluid-solid coupling control equation.

[0154] Step 2: Derive the finite element discrete format of the fractional order fluid-solid coupling control equation, and generate the finite element discrete equation.

[0155] Step 3: Establish the INP input file.

[0156] Specifically, the INP input file is a model information file that ABAQUS reads for user needs to perform finite element simulation, which contains information such as element number of geometric model, node number, node spatial coordinates, model material parameters, model boundary conditions, model initial conditions, model mesh type, model analysis step type, etc. The generation process of the initial INP file is as follows: in ABAQUS, establish a geometric model through the conventional modeling steps, complete the material property definition, assembly and analysis step setting; apply boundary conditions and pre-defined fields to the model, and divide the finite element mesh; create a job in the "Job" module, and export the input file through the "Write Input" command, which will generate the corresponding INP file in the set working directory.

[0157] Step 4: Modify the INP file.

[0158] Specifically, in order to enable the model to correctly call the written UEL subroutine code, the following modifications and supplements need to be made to the generated INP file:

[0159] (1) Declare the UEL custom element type:

[0160] Declare the use of UEL custom elements in the INP file by adding the following keyword statements:

[0161] USER ELEMENT, TYPE=Un, NODES=, COORDINATES=, PROPERTIES=,VARIABLES=

[0162] Data line(s);

[0163] Wherein, TYPE represents the user-defined element type, named Un (n is an integer number); NODES is the number of element nodes; COORDINATES represents the geometric dimension of the element (2 for two-dimensional elements and 3 for three-dimensional elements); PROPERTIES is the number of material property parameters read by the UEL subroutine; VARIABLES is the number of state variables; Dataline(s) is used to define the degree of freedom sequence number of the element. For the fractional order fluid-structure coupling control equation in this application, the displacement and pore water pressure degrees of freedom should be activated, corresponding to ABAQUS degree of freedom numbers 1 (u), 2 (v), and 8 (p) in turn.

[0164] (2) Specify the use of UEL element set:

[0165] When defining the model elements, it is necessary to specify which elements use custom element types. This is achieved through the following key statements:

[0166] ELEMENT, TYPE =Un, ELSET=USER

[0167] Data line(s);

[0168] Wherein, TYPE=Un indicates that the element type used is consistent with the definition described above; ELSET=USER defines an element set containing all custom elements; Data line(s) lists the numbers of each custom element and its corresponding node numbers.

[0169] (3) Input parameters required by UEL subroutine:

[0170] Add the following statements to the file to input the model parameters required by UEL:

[0171] UEL PROPERTY, ELSET=USER

[0172] Data line(s);

[0173] Where, ELSET=USER corresponds to the aforementioned element set name; Data line(s) is the list of material parameters required by UEL, separated by commas between each parameter. The parameter list should be consistent with the order of parameters read by the PROPS array in the UEL subroutine.

[0174] Step 5: Write UEL subroutine code.

[0175] (1) Subroutine structure and variable definition;

[0176] The basic format of the UEL subroutine is as follows:

[0177] SUBROUTINE UEL (…………)

[0178] INCLUDE 'ABA_PARAM.INC'

[0179] DIMENSION……………

[0180] User code lines

[0181] RETURN

[0182] END;

[0183] Where, SUBROUTINE UEL (...) is the standard interface form called by ABAQUS, and the parentheses contain several parameters passed in by the main program; INCLUDE 'ABA_PARAM.INC' is the system predefined parameter file, which contains global constants and precision control definitions; the DIMENSION statement is used to declare the dimensions and storage space of the main variables; User code lines is the UEL code writing area.

[0184] (2) Declare COMMON area;

[0185] COMMON / zone name / matrix name;

[0186] Used to store shared variables such as integration weight functions and historical pore water pressure gradients in the subroutine, so as to be reused in each Gauss integration point calculation.

[0187] (3) Material and model parameter reading;

[0188] UEL reads the material parameters defined in the INP file through the PROPS array, and its calling method is as follows:

[0189] aa = PROPS(1)

[0190] bb = PROPS(2)

[0191] cc = PROPS (3)

[0192] dd = PROPS (4)

[0193] ………………;

[0194] The parameter order must be consistent with the input order under the key statement UEL PROPERTY, ELSET=USER in the INP file to ensure the accuracy of physical quantity transmission.

[0195] (4) Calculate the weight functions. The number of weight functions is the number of analysis steps in the numerical model. Each weight function is stored in a matrix declared in the COMMON area and called later when calculating the fractional derivative.

[0196] (5) Determine the Gaussian integration points and shape function matrix: The cell is a two-dimensional 8-node cell, and its Gaussian integration points and weights are shown in Table 1.

[0197] Table 1

[0198]

[0199] Its shape function is:

[0200]

[0201] in, ; ; ; ; ; ; ; .

[0202] (6) Calculate the shape function and geometric transformation at each Gaussian integration point;

[0203] Specifically, for each Gaussian integration point Perform the following calculations:

[0204] Shape functions and their derivatives: Calculate the shape functions of an 8-node element. and its derivative with respect to the reference coordinates .

[0205] Construct the Jacobi matrix:

[0206] ;

[0207] Construct the spatial coordinate matrix of shape functions:

[0208] ;

[0209] Construct the B matrix:

[0210] ;

[0211] Construct the D matrix, where the two-dimensional elastic stiffness matrix is established based on the plane strain assumption:

[0212] ;

[0213] Construct the matrix :

[0214] ;

[0215] (7) Calculate the local stiffness block at each Gauss integration point and integrate and accumulate:

[0216] Solid mechanics stiffness matrix:

[0217] ;

[0218] ;

[0219] Fluid-structure interaction stiffness matrix:

[0220] ;

[0221] ;

[0222] Porous medium seepage stiffness matrix:

[0223] ;

[0224] ;

[0225] Where tk is the thickness of the two-dimensional model.

[0226] (8) Assemble the element stiffness matrix AMATRX, the size of the stiffness matrix of the two-dimensional eight-node element is 24x24. The assembled matrix is in the form of:

[0227] ;

[0228] (9) Calculate the time and spatial gradient of pore water pressure:

[0229] ;

[0230] ;

[0231] ;

[0232] The spatial gradient of the pore water pressure at the current time is calculated each time, and is recorded in the matrix according to the time step, the Gaussian integration point, the element number and the direction through the matrix declared by the COMMON area.

[0233] (10) Calculate the equivalent spatial gradient of the pore water pressure:

[0234] ;

[0235] ;

[0236] (11) Calculate the local residual term:

[0237] ;

[0238] ;

[0239] (12) Assemble the element residual term:

[0240] .

[0241] wherein, is the weight value of the current Gaussian integration point, is the determinant of the Jacobi matrix.

[0242] Step 6: After the UEL subroutine is programmed and the INP input file is modified, the numerical solution of the fractional order fluid-structure coupling model can be carried out in the ABAQUS software environment. The entire solving process takes the INP file as the core data carrier, and the iterative calculation of the fluid-structure coupling behavior is realized by calling the user-defined unit (UEL). The specific steps are as follows:

[0243] (1) Establish an analysis job: Enter the Job module in the ABAQUS main interface, and create a new analysis job. In the "Model Source" option, select "Input File" and specify the path of the modified INP file. This file contains the geometric information, material parameters, boundary conditions and UEL unit definition of the model.

[0244] (2) Associate the UEL subroutine file: In the new job, select the "User Defined Subroutine" option and specify the path of the programmed UEL subroutine file.

[0245] (3) Submit the calculation job: After confirming that the settings are correct, save the job configuration and click the "Submit" button to start the calculation. ABAQUS will automatically call the UEL subroutine to numerically iterate and solve the user-defined unit part of the model, and complete the coupling calculation.

[0246] (4) Output and Analysis of Calculation Results: After the solution is completed, a result file with the .odb extension will be generated in the working directory. Post-processing can be performed on the solid displacement field, pore water pressure field, fractional-order seepage rate, etc. The reliability of the proposed fractional-order fluid-structure interaction model numerical calculation method is verified by exporting the calculation results and comparing them with existing numerical schemes.

[0247] like Figure 2 The figure shows a comparison of the numerical calculation results when using the UMATHT subroutine and the UEL subroutine of ABAQUS software for numerical solution. The horizontal axis represents time (in seconds), and the vertical axis represents the pore water pressure. By comparing the calculation results based on the finite difference method (FDM) with the results obtained using the UMATHT interface of ABAQUS, it can be seen that the three result curves basically overlap. This verifies the correctness of the numerical simulation method in the embodiments of this application. Figure 3 The figure shows a comparison of the numerical computation convergence efficiency as the analysis step progresses when using the UMATHT and UEL user subroutines in ABAQUS. The horizontal axis represents the analysis step, and the vertical axis represents the number of iterations. The figure shows that when using the UEL interface, the number of iterations required for each analysis step is less than that of the UMATHT-based scheme, indicating that the program has a faster convergence speed. This result further verifies the stability of the numerical simulation method used in this application during the solution process.

[0248] In one embodiment, an apparatus is provided for implementing a numerical method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model. This apparatus corresponds one-to-one with the numerical method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model described in the above embodiments. Figure 4 As shown, the numerical realization device 100 for fluid-structure interaction in saturated porous media based on a fractional-order seepage model includes: an acquisition module 101, a construction module 102, a first generation module 103, a second generation module 104, and a third generation module 105. Detailed descriptions of each functional module are as follows:

[0249] The acquisition module 101 is used to acquire the geometric parameters, material parameters and boundary conditions of the porous medium to be analyzed. The material parameters include elastic modulus, Biot coefficient, fluid density, water storage coefficient, fractional permeability coefficient and fractional order number.

[0250] The construction module 102 is configured to construct a fractional order fluid-solid coupling control equation for describing the interaction between the solid skeleton deformation and the pore fluid fractional order seepage in the porous medium based on material parameters, wherein the fractional order fluid-solid coupling control equation comprises a mechanical field control equation and a seepage field control equation based on a fractional order derivative;

[0251] The first generation module 103 is configured to derive a finite element discrete format of the fractional order fluid-solid coupling control equation based on a user-defined unit framework, and generate a finite element discrete equation for compiling a user-defined unit subroutine, wherein the finite element discrete equation comprises a first discrete equation derived from the mechanical field control equation and a second discrete equation derived from the seepage field control equation based on the fractional order derivative;

[0252] The construction module 102 is further configured to construct a numerical calculation model of the porous medium based on geometric parameters and boundary conditions, and perform finite element mesh division on the numerical calculation model;

[0253] The second generation module 104 is configured to compile the user-defined unit subroutine based on the finite element discrete equation and the numerical calculation model, and generate a finite element input file containing a user-defined unit configuration;

[0254] The third generation module 105 is configured to associate the finite element input file with the user-defined unit subroutine, perform fluid-solid coupling numerical solution through an ABAQUS platform, and obtain simulation results of the porous medium, wherein the simulation results comprise a displacement field, a pore water pressure field and a seepage rate.

[0255] The present application provides a saturated porous medium fluid-solid coupling numerical implementation device 100 based on a fractional order seepage model, which constructs and discretizes the coupling control equation containing the fractional order derivative, generates the finite element discrete equation containing two degrees of freedom of displacement and pore pressure which can be directly coded, and compiles the user-defined unit subroutine accordingly. With strict element-level discretization and linearization processing, the stability and convergence efficiency of numerical calculation are significantly improved. At the same time, since the pore pressure degree of freedom is directly introduced and solved, the original temperature field module is not occupied, and the ability to further expand the model to real thermal-fluid-solid multi-field coupling analysis is fully retained. The inherent limitations of the traditional UMATHT interface scheme are fundamentally solved, the physical expandability and application potential of the model in complex engineering scenarios involving temperature changes are significantly enhanced on the basis of ensuring the robustness and reliability of the calculation process, thereby providing reliable technical support for the fine simulation and safety evaluation of porous media in the fields of geotechnical engineering, groundwater resource assessment, energy geology and the like.

[0256] The specific limitations of the device for numerically implementing fluid-solid coupling of saturated porous media based on a fractional seepage model can refer to the limitations of the method for numerically implementing fluid-solid coupling of saturated porous media based on a fractional seepage model described above, and will not be repeated here. Each module in the device for numerically implementing fluid-solid coupling of saturated porous media based on a fractional seepage model described above can be realized by software, hardware, and combinations thereof, in whole or in part. The above-mentioned modules can be embedded in or independent of the processor in the electronic device in hardware form, or can be stored in the memory in the electronic device in software form, so as to be called and executed by the processor to perform the operations corresponding to each of the above modules.

[0257] In one embodiment, an electronic device is provided, including a memory, a processor, and a computer program stored in the memory and executable on the processor, and the processor implements the method for numerically implementing fluid-solid coupling of saturated porous media based on a fractional seepage model described above when executing the computer program.

[0258] In one embodiment, a computer-readable storage medium is provided, and the computer-readable storage medium stores a computer program, and the computer program is executed by a processor to implement the method for numerically implementing fluid-solid coupling of saturated porous media based on a fractional seepage model described above.

[0259] It should be noted that the functions or steps that the above-mentioned computer-readable storage medium or electronic device can achieve can correspond to the related descriptions of the server side and the client side in the foregoing method embodiments, and to avoid repetition, they will not be described one by one here.

[0260] Those skilled in the art can understand that all or part of the processes in the above-mentioned embodiment methods can be completed by instructing the relevant hardware through a computer program. The computer program can be stored in a non-volatile computer readable storage medium, and when executed, can include the processes of the above-mentioned embodiment methods. Any reference to memory, storage, database or other medium used in the embodiments provided by the present application can include non-volatile and / or volatile memory. Non-volatile memory can include read-only memory (ROM), programmable ROM (PROM), electrically programmable ROM (EPROM), electrically erasable programmable ROM (EEPROM) or flash memory. Volatile memory can include random access memory (RAM) or external cache memory. As an illustration but not limitation, RAM is available in various forms, such as static RAM (SRAM), dynamic RAM (DRAM), synchronous DRAM (SDRAM), double data rate SDRAM (DDR SDRAM), enhanced SDRAM (ESDRAM), synchronous link (Synchlink) DRAM (SLDRAM), memory bus (Rambus) direct RAM (RDRAM), direct memory bus dynamic RAM (DRDRAM), and memory bus dynamic RAM (RDRAM), etc.

[0261] Those skilled in the art can clearly understand that, for the convenience and brevity of description, only the division of the above-mentioned functional units and modules is exemplified, and in actual application, the above-mentioned functions can be completed by different functional units and modules according to needs, that is, the internal structure of the device is divided into different functional units or modules to complete all or part of the functions described above.

[0262] The above-mentioned embodiments are only used to illustrate the technical solutions of the present application, but not limit it. Although the present application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that they can modify the technical solutions recorded in the foregoing embodiments, or make equivalent replacement for part of the technical features. The modification or replacement does not make the essence of the corresponding technical solution deviate from the spirit and scope of the technical solutions of the embodiments of the present application, and should be included in the protection scope of the present application.

Claims

1. A numerical implementation method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model, characterized in that, include: Obtain the geometric parameters, material parameters, and boundary conditions of the porous medium to be analyzed. The material parameters include elastic modulus, Biot coefficient, fluid density, water storage coefficient, fractional permeability coefficient, and fractional order. Based on the material parameters, a fractional-order fluid-structure interaction control equation is constructed to describe the interaction between solid skeleton deformation and fractional-order seepage of pore fluid in porous media. The fractional-order fluid-structure interaction control equation includes a mechanical field control equation and a seepage field control equation based on fractional derivatives. Based on the user-defined element framework, the fractional-order fluid-structure interaction control equations are derived using the finite element discretization format, generating finite element discretization equations for writing user-defined element subroutines. The finite element discretization equations include a first discretization equation derived from the mechanical field control equations and a second discretization equation derived from the seepage field control equations based on fractional derivatives. Based on the geometric parameters and boundary conditions, a numerical calculation model of the porous medium is constructed, and the numerical calculation model is meshed using finite element methods. Based on the finite element discrete equations and the numerical calculation model, compile the user-defined element subroutine and generate a finite element input file containing the user-defined element configuration; The finite element input file is associated with the user-defined unit subroutine, and the fluid-structure interaction numerical solution is executed through the ABAQUS platform to obtain the simulation results of the porous medium, wherein the simulation results include displacement field, pore water pressure field and seepage rate.

2. The numerical implementation method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model according to claim 1, characterized in that, The fractional-order fluid-structure interaction control equations are as follows: The governing equations of the mechanical field: ; in, For stress tensor; Let i be the i-th component of the volume force vector; ; in, Let i be the i-th and j-th components of the stress tensor; Let i be the i-th and j-th components of the effective stress tensor; Biot coefficient; Pore ​​water pressure; The symbol for Kronecker; ; ; in, For the fourth-order elastic stiffness tensor of the material; Let k and l be the kth and lth components of the strain tensor; The seepage field control equations based on fractional derivatives are as follows: ; in, For fluid density; denoted as the water storage coefficient of the porous medium; p is the pore water pressure; t is time; q is the seepage velocity vector. For quality source items; ; Where q is the seepage velocity vector; It is a fractional permeability coefficient; Where is the fluid density; g is the acceleration due to gravity; It is a fractional derivative operator; t is the order of the derivative; a is the initial time; t is the current time. The gradient of pore water pressure; ; in, For function of fractional derivative; is the time step; k is the total number of steps, i iterates from 0 to k.

3. The numerical implementation method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model according to claim 1, characterized in that, The step of deriving the fractional-order fluid-structure interaction control equations using a user-defined element framework and generating finite element discrete equations for writing user-defined element subroutines specifically includes: Obtain the preset node configuration of the finite element element, wherein the preset node configuration includes the preset number of nodes of the finite element element, and the preset displacement degree of freedom and preset pore water pressure degree of freedom of each node; Based on the preset node configuration, construct a unit node variable array: ; Where n is the nodal variable; e is the element level; T is the transpose symbol; u, v, w are the preset displacement degrees of freedom; Let x be the displacement component of the nth node in the x-direction; Let be the displacement component of the nth node in the y-direction; Let be the displacement component of the nth node in the z-direction; p is the preset pore water pressure degree of freedom. Let be the pore water pressure at the nth node; The element node variable array is divided into an element node displacement array and an element node pore water pressure array: ; ; Where d is the displacement; Within the user-defined unit framework, based on the principle of virtual work, the finite element discretization scheme is used to derive the mechanical field control equations, generating the first discrete equation. Within the user-defined unit framework, based on the variational principle, the seepage field control equations based on fractional derivatives are derived using a finite element discretization scheme, generating a second discrete equation.

4. The numerical implementation method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model according to claim 3, characterized in that, The step of deriving the second discrete equation by performing a finite element discretization scheme for the seepage field control equation based on fractional derivatives within a user-defined element framework, based on variational principles, specifically includes: A functional form of the seepage field control equation based on fractional derivatives is constructed, and spatial interpolation of pore water pressure is performed based on shape functions. For the fractional pressure gradient term in the seepage field control equation based on fractional derivatives, the Grünwald-Letnikov finite difference definition is used for time discretization, and it is expressed as a weighted convolution series between the current time and the historical time. The weighted convolution series is decoupled into instantaneous response terms and historical memory terms; Based on the extreme value conditions of the variational principle, the tangent stiffness matrix relationship characterizing the permeability of the unit at the current moment is derived and used as the second discrete equation.

5. The numerical implementation method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model according to claim 3, characterized in that, The finite element discretization equation is: First discrete equation: ; in, Let be the stiffness matrix of element e; Let be the nodal displacement vector of element e; Let be the coupling matrix of element e; Let be the nodal pore water pressure vector of element e; This is the equivalent nodal force vector of the mechanical field corresponding to element e; ; Where D is the elastic stiffness matrix of the material; B is the geometric matrix of the element; ; in, The pore pressure influence matrix; ; Where N is the unit shape function matrix; This is the traction force vector; Second discrete equation: ; in, Let be the seepage matrix of element e; Let be the equivalent nodal force vector of the seepage field in element e; ; 。 6. The numerical implementation method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model according to claim 1, characterized in that, The step of compiling a user-defined element subroutine based on the finite element discrete equations and the numerical calculation model, and generating a finite element input file containing the user-defined element configuration, specifically includes: Based on the numerical calculation model, a finite element input file is generated, wherein the finite element input file includes the geometric parameters, mesh generation information, boundary conditions and material parameters of the numerical calculation model; Based on the finite element discrete equations, a user-defined element subroutine is compiled, wherein the user-defined element subroutine is used to calculate the stiffness matrix and residual components of the element. In the finite element input file, declare the identifier type and degrees of freedom of the user-defined element, define the set of elements using the user-defined element, and write the list of material parameters called by the user-defined element subroutine.

7. The numerical implementation method for fluid-structure interaction in saturated porous media based on a fractional-order seepage model according to claim 1, characterized in that, The step of associating the finite element input file with the user-defined element subroutine and performing fluid-structure interaction numerical solution through the ABAQUS platform to obtain the simulation results of the porous medium specifically includes: Create an analysis job in the ABAQUS platform and associate the finite element input file with the user-defined element subroutine in the analysis job; Submit the analysis job to the ABAQUS solver, and perform fluid-structure interaction numerical calculations by calling the user-defined element subroutine to obtain simulation results including displacement field, pore water pressure field and seepage rate.

Citation Information

Patent Citations

  • Rock-soil body discrete element fluid-solid coupling numerical simulation method based on pore density flow

    CN110263362A

  • Simulation method, device and equipment for predicting flow and deformation of porous medium and medium

    CN120951848A