A numerical implementation method of thermo-hydro-elastic coupling based on fractional viscoelastic model

Abstract

The application discloses a kind of based on fractional order viscoelastic model's thermal fluid-solid coupling numerical implementation method, belong to finite element model numerical implementation technical field, including the establishment of containing equilibrium equation, fractional order Merchant constitutive equation, fluid continuity equation based on Darcy's law and energy equation based on Fourier heat transfer law thermal fluid-solid coupling control equation;Fractional order derivative in constitutive equation is dispersed and approximated using Grünwald-Letnikov definition fractional order derivative, the finite element discrete format of control equation is deduced, and unit stiffness matrix and residual matrix are obtained;In finite element software platform, based on unit stiffness matrix and residual matrix, numerical solution is carried out by user-defined unit interface construction unit calculation module.The application can independently determine the coupling term in multi-field coupling control equation, and is suitable for numerical simulation of actual engineering model with complex geometry and boundary conditions.

Description

Technical Field

[0001] This invention relates to the field of numerical implementation technology of finite element models, and in particular to a numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model. Background Technology

[0002] The fractional Merchant model is a constitutive model describing the deformation of viscoelastic media. Existing numerical developments of the fractional Merchant model involving multi-field coupling are mostly based on the following two approaches: (1) The governing equations of the model are simplified to a one-dimensional case and solved by Laplace transform or finite difference programming. Although this approach can complete the model calculation and parameter analysis, it lacks application capability in actual engineering with complex geometry and boundary conditions.

[0003] (2) The calculation code for the fractional-order Merchant model was written in the finite element software ABAQUS using UMAT (User Material Subroutine) and embedded into the software's default multi-field coupled control equation framework. Although this approach can solve the limitations of the above approach in practical engineering applications, since UMAT can only define material constitutive relations, the overall framework and coupling terms of the multi-field coupled control equations are still determined by the ABAQUS built-in solver. Users cannot adjust the coupling terms in the control equations according to the physical mechanism of the actual problem.

[0004] The fundamental reason for the aforementioned limitations is that the UMAT interface is designed solely to define material constitutive relations, and its input / output parameters do not involve the overall architecture of the multi-field coupled governing equations; while directly modifying the ABAQUS core solver code is neither feasible nor permitted. Therefore, a numerical solution scheme is needed that can autonomously determine the coupling terms of the multi-field governing equations and can be applied to the simulation of actual engineering models. Summary of the Invention

[0005] The purpose of this invention is to provide a numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model, so as to solve the problems mentioned in the background art.

[0006] To achieve the above objectives, this invention provides a numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model, comprising the following steps: Step 1: Establish fractional-order thermal-fluid-structure interaction (TFS) control equations. These equations include equilibrium equations, fractional-order Merchant constitutive equations, fluid continuity equations based on Darcy's law, and energy equations based on Fourier's heat transfer law. The fractional derivative terms in the fractional-order Merchant constitutive equations are discretized and approximated. The finite element discretization scheme for the fractional-order TFS control equations is then derived, resulting in an element stiffness matrix assembled from multiple stiffness sub-matrices and a residual matrix assembled from multiple residual sub-matrices. Step 2: In the finite element software platform, based on the element stiffness matrix and residual matrix obtained in Step 1, an element calculation module is constructed through a user-defined element interface to numerically solve the displacement field, pore water pressure field, and temperature field of the target porous media engineering structure.

[0007] Preferably, in step one, the discrete approximation uses the fractional derivative as defined by Grünwald-Letnikov. The definition of the Grünwald-Letnikov fractional derivative is: ; in, From arrive of Fractional derivative operators, Let be the function to be differentiated. For time step, For fractional order, For historical time step indexing, For the total number of historical items, The coefficients are binomial coefficients; The discrete form of the Grünwald-Letnikov fractional derivative is approximated as: ; in, For the current moment, This is the current time step number. For the first The weighting coefficients for each historical time step For function exist The value of each time step; In step one, the derivation of the finite element discretization scheme uses the integration by parts method.

[0008] Preferably, in step one, the fractional-order thermo-fluid-structure interaction governing equations include the following four sets of equations: Equilibrium equations: ; in, These are components of the stress tensor. , For tensor subscripts, Biot coefficient, Pore ​​water pressure, For Kroneck symbol, when hour, ,when hour, , These are the components of the volume force vector; Fractional Merchant constitutive equations: ; in, and For the elastic modulus parameter, The viscosity-time coefficient, From 0 to of Fractional derivative operators, Poisson's ratio, These are the components of the strain tensor. The coefficient of thermal expansion of a solid is... This refers to the change in temperature. Fluid continuity equation based on Darcy's law: ; in, The water storage coefficient of porous media, For time, Permeability coefficient, For fluid density, It is the acceleration due to gravity. For the Laplace operator, For volumetric strain, The coefficient of thermal expansion of the fluid. For temperature; Energy equation based on Fourier's law of heat transfer: ; in, For the effective volumetric heat capacity of porous media, For fluid density, Specific heat capacity of the fluid This is the Darcy velocity vector. For gradient operators, The thermal conductivity coefficient, The thermal expansion factor of solids This is a reference temperature.

[0009] Preferably, in step one, when deriving the finite element discretization scheme using the integration by parts method, trial functions of the displacement field, pore water pressure field, and temperature field are introduced to obtain the following weak form equations: Weak form of the equilibrium equation: ; in, To solve for the volume of the region, For virtual strain tensor, and For the elasticity matrix, For strain tensor, For the vector form of the Kronecker symbol, For equivalent stress, To account for the equivalent strain due to temperature effects, , For elasticity model parameters, For virtual displacement vectors, The volume force vector, superscript Indicates matrix transpose; Weak form of the fluid continuity equation: ; in, For virtual pore water pressure, For fluid specific weight; The weak form of the energy equation: ; in, This is a virtual temperature.

[0010] Preferably, the element stiffness matrix obtained in step one is assembled from the following stiffness sub-matrices: The stiffness submatrix corresponding to the equilibrium equations is: ; in, The stiffness submatrix corresponds to the displacement degrees of freedom. This is the stiffness submatrix coupled with displacement and pore water pressure. This is the stiffness submatrix coupled with displacement and temperature. For geometric matrices, For the elasticity matrix, For the shape function matrix, To solve for the region; The stiffness submatrix corresponding to the fluid continuity equation is: ; in, The stiffness submatrix is ​​the coupling between pore water pressure and displacement. This is the stiffness submatrix corresponding to the degree of freedom of pore water pressure. This is the stiffness submatrix coupled with pore water pressure and temperature; The stiffness submatrix corresponding to the energy equation is: ; in, This is the stiffness submatrix coupled with temperature and displacement. The stiffness submatrix is ​​a coupling matrix between temperature and pore water pressure. This is the stiffness submatrix corresponding to the temperature degree of freedom.

[0011] Preferably, the residual matrix obtained in step one is assembled from the residual submatrix of the equilibrium equation, the residual submatrix of the fluid continuity equation, and the residual submatrix of the energy equation.

[0012] Preferably, in step two, based on the actual coupling requirements of the target multi-field coupling problem, one or more stiffness sub-matrices are selectively activated or disabled when assembling the element stiffness matrix; when assembling the element stiffness matrix, the matrix elements corresponding to the disabled stiffness sub-matrices are set to zero; when assembling the residual matrix, one or more residual sub-matrices are selectively activated or disabled; when assembling the residual matrix, the matrix elements corresponding to the disabled residual sub-matrices are set to zero.

[0013] Preferably, in step two, when constructing the unit calculation module, a historical variable storage area is set up to store the historical weight coefficients of each Gaussian integral point at each time step. Historical stress tensor and historical strain tensor Historical weighting coefficients, historical stress tensors, and historical strain tensors are all accessed and indexed by element number, Gaussian integration point number, and time step number.

[0014] Preferably, in step two, before performing numerical solutions through the user-defined element interface in the finite element software platform, the step also includes an input file generation and modification step: A numerical model is established in the finite element software platform, displacement boundary conditions, pore water pressure boundary conditions and temperature boundary conditions are applied, finite element mesh generation is completed, and the input file is exported. Add a user element declaration statement to the input file to activate the displacement degree of freedom, pore water pressure degree of freedom, and temperature degree of freedom. Create a set of user cells in the input file, specifying the cell region to be calculated by the cell calculation module built using the user-defined cell interface; Add a parameter input statement to the input file, and input the material parameters in the order they are read from the parameter array in the user-defined unit interface.

[0015] Preferably, the displacement degree of freedom corresponds to the translational degree of freedom along three orthogonal spatial directions in the finite element software platform, the pore water pressure degree of freedom is the fluid pressure degree of freedom, and the temperature degree of freedom is the thermodynamic temperature degree of freedom.

[0016] Therefore, the present invention employs the above-mentioned numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model, which has the following beneficial effects: (1) This invention uses the UEL (User Element) interface of ABAQUS for secondary development. Users can selectively activate or disable one or more stiffness sub-matrices and residual sub-matrices when assembling the stiffness matrix and residual matrix of the element according to the physical mechanism of the actual problem, thereby flexibly adjusting the coupling terms in the multi-field coupling equation and overcoming the defect of the existing UMAT scheme that cannot independently determine the coupling terms. (2) This invention establishes a finite element discretization scheme for the three-dimensional fractional-order thermo-fluid-structure interaction control equations based on the finite element method, and embeds it into the ABAQUS platform through a UEL subroutine, enabling numerical simulation of thermo-fluid-structure interaction for actual engineering structures with complex geometries and boundary conditions. This overcomes the limitation that existing Laplace transform or finite difference schemes are only applicable to one-dimensional simplified cases; (3) This invention establishes a complete process of fractional-order thermal-fluid-structure interaction control equation → finite element discretization → UEL program implementation, providing a general numerical implementation framework for the application of fractional-order Merchant constitutive model in thermal-fluid-structure interaction analysis of saturated porous media. (4) The ABAQUS platform of this invention implements a complete numerical solution process (INP file generation → modification → UEL subroutine compilation → calculation submission → result export and verification). By comparing and verifying the results with existing numerical implementation schemes, the reliability of the proposed method is ensured.

[0017] The technical solution of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0018] Figure 1 This is a flowchart of a numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to the present invention. Figure 2 This is a diagram showing the verification results of the fractional-order viscoelastic constitutive model in an embodiment of the present invention. Figure 3 The figure shows the verification results of the thermal-fluid-structure interaction solution framework of the present invention. (a) compares the numerical solution and analytical solution of displacement response at different heights, (b) compares the numerical solution and analytical solution of pore water pressure response at different heights, and (c) compares the numerical solution and analytical solution of temperature response at different heights. Detailed Implementation

[0019] The following detailed description of embodiments of the invention provided in the accompanying drawings is not intended to limit the scope of the claimed invention, but merely illustrates selected embodiments of the invention. All other embodiments obtained by those skilled in the art based on the embodiments of the invention without inventive effort are within the scope of protection of the invention.

[0020] Example like Figure 1 As shown, this invention provides a numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model, comprising the following steps: Step 1: Establish fractional-order thermal-fluid-structure interaction (TFS) control equations. These equations include equilibrium equations, fractional-order Merchant constitutive equations, fluid continuity equations based on Darcy's law, and energy equations based on Fourier's heat transfer law. The fractional derivative terms in the fractional-order Merchant constitutive equations are discretized and approximated. The finite element discretization scheme for the fractional-order TFS control equations is then derived, resulting in an element stiffness matrix assembled from multiple stiffness sub-matrices and a residual matrix assembled from multiple residual sub-matrices.

[0021] Assuming the porous medium is isotropic, the stress-strain relationship satisfies the fractional-order Merchant constitutive model, the pore water velocity and pressure obey Darcy's law, and the heat flow velocity and temperature obey Fourier's law of heat conduction. The fractional-order thermo-fluid-structure interaction governing equations consist of the following four sets of equations.

[0022] Equilibrium equations: ; in, These are components of the stress tensor. , For tensor subscripts, Biot coefficient, Pore ​​water pressure, For Kroneck symbol, when hour, ,when hour, , The components of the volume force vector are denoted by ; this equation describes the force balance relationship of porous media under stress.

[0023] Fractional Merchant constitutive equations: ; in, and For the elastic modulus parameter, The viscosity-time coefficient, From 0 to of Fractional derivative operators, Poisson's ratio, These are the components of the strain tensor. The coefficient of thermal expansion of a solid is... The equation represents the temperature change and describes the stress-strain relationship of a fractional-order viscoelastic material, where the fractional derivative term reflects the material's time memory characteristics.

[0024] Fluid continuity equation based on Darcy's law: ; in, The water storage coefficient of porous media, For time, Permeability coefficient, For fluid density, It is the acceleration due to gravity. For the Laplace operator, For volumetric strain, The coefficient of thermal expansion of the fluid. For temperature; this equation describes the mass conservation relationship of porous fluid in porous media.

[0025] Energy equation based on Fourier's law of heat transfer: ; in, For the effective volumetric heat capacity of porous media, For fluid density, Specific heat capacity of the fluid This is the Darcy velocity vector. For gradient operators, The thermal conductivity coefficient, The thermal expansion factor of solids The reference temperature is used; this equation describes the heat conservation relationship in porous media.

[0026] The fractional derivative has multiple definitions; this embodiment uses the Grünwald-Letnikov (GL) definition for the discrete approximation. The Grünwald-Letnikov fractional derivative is defined as follows: ; in, From arrive of Fractional derivative operators, Let be the function to be differentiated. For time step, For fractional order, For historical time step indexing, For the total number of historical items, The coefficients are binomial coefficients; The discrete form of the Grünwald-Letnikov fractional derivative is approximated as: ; in, For the current moment, This is the current time step number. For the first The weighting coefficients for each historical time step can be determined according to the fractional order. The pre-computed storage is equal in number to the total number of analysis steps in the numerical model. For function exist The value of each time step.

[0027] Substituting the definition of the fractional derivative into the fractional Merchant constitutive equation, we can derive the equation. First, rearranging the fractional Merchant constitutive equation, neglecting thermal expansion, we get: ; Substituting the definition of the fractional derivative into the above equation and rearranging, we get: ; in, For average strain, , For volumetric strain, , , , , These are coefficients obtained after algebraic operations, used to simplify the expression of equations, and are intermediate results in mathematical derivation.

[0028] .

[0029] After rearranging, the stress-strain expression using the fractional-order Mercat constitutive equation can be obtained again: ; After considering thermal strain, the above equation becomes: ; in, , , , , .

[0030] For equivalent change, For equivalent stress, To account for the equivalent strain of temperature effects, since the calculation of fractional derivatives requires the accumulation of historical variables, which have no specific physical meaning, they are just abbreviations of long formulas to make the formulas look simpler.

[0031] Substituting the above stress-strain expressions into the equilibrium equations, and combining them with the fluid continuity equation and the energy equation, the governing equations for the thermal-fluid-structure interaction of saturated porous media considering the fractional-order Merchant constitutive model can be rewritten as follows: .

[0032] In deriving the finite element discretization scheme of the governing equations, the integration by parts method is employed. It is assumed that the number of nodes in the finite element is... The order of displacement, pore water pressure, and temperature among the element node variables is as follows: ; To facilitate the subsequent derivation of the finite element stiffness matrix, the element node array is divided into the following three parts: Element nodal displacement array: ; Unit node pore water pressure array: ; Unit node temperature array: ; in, , , These are the finite element elements. Each node Directional displacement, For the finite element element Pore ​​water pressure at each node, For the finite element element Temperature of each node.

[0033] Introducing trial functions (virtual displacement vectors) of displacement field, pore pressure field, and temperature field. Virtual pore water pressure and virtual temperature The governing equations are discretized using the Galerkin method. After integration by parts, the following weak-form equations are obtained.

[0034] Weak form of the equilibrium equation: ; in, To solve for the volume of the region, For virtual strain tensor, and For the elasticity matrix, For strain tensor, For the vector form of Kronecker notation, in the three-dimensional case, it can be represented as: , For equivalent stress, To account for the equivalent strain due to temperature effects, , For elasticity model parameters, For virtual displacement vectors, The volume force vector, superscript Indicates matrix transpose; Weak form of the fluid continuity equation: ; in, For virtual pore water pressure, For fluid specific weight; The weak form of the energy equation: ; in, This is a virtual temperature.

[0035] In finite element discretization, the element displacement field Using shape function matrix and nodal displacement matrix Approximation; pore water pressure and temperature Using shape function matrices respectively and the corresponding node vectors and Approximation: , , ; in, It is a geometric matrix.

[0036] Substituting the above discrete unit equations into the weak form equations, we get: ; ; ; Due to dummy variables , , Having arbitrary characteristics, after eliminating the dummy variable, the above equation becomes: ; ; ; Based on the above equation, the relationship between the stiffness matrix and the residual matrix can be obtained: ; The stiffness submatrix corresponding to the equilibrium equations is: ; in, The stiffness submatrix corresponds to the displacement degrees of freedom. This is the stiffness submatrix coupled with displacement and pore water pressure. This is the stiffness submatrix coupled with displacement and temperature. For the elasticity matrix, For the shape function matrix, To solve for the region; The stiffness submatrix corresponding to the fluid continuity equation is: ; in, The stiffness submatrix is ​​the coupling between pore water pressure and displacement. This is the stiffness submatrix corresponding to the degree of freedom of pore water pressure. This is the stiffness submatrix coupled with pore water pressure and temperature; The stiffness submatrix corresponding to the energy equation is: ; Where is the stiffness submatrix coupled with temperature and displacement, is the stiffness submatrix coupled with temperature and pore water pressure, and is the stiffness submatrix corresponding to the temperature degree of freedom.

[0037] The residual matrix is ​​derived from the residual submatrix of the equilibrium equation. residual submatrix of fluid continuity equation and the residual submatrix of the energy equation It is assembled from the weak forms of the equilibrium equation, the fluid continuity equation, and the energy equation, excluding the stiffness term (including the contribution of the history term and the boundary condition term).

[0038] The residual submatrix of the equilibrium equation The integral expression is: ; in, This is the node displacement vector at the current time step. This is the node pore water pressure vector at the current time step. This represents the temperature change at the current time step. Let be the stress tensor at the i-th historical time step.

[0039] The integral expression for the residual submatrix of the fluid continuity equation is: ; in, This is the nodal pore water pressure vector from the previous time step. Let be the nodal displacement vector of the previous time step. This is the node temperature vector at the current time step. This is the node temperature vector from the previous time step.

[0040] Energy equation residual submatrix The integral expression is: ; In the above derivation, corresponding terms can be added or deleted according to the actual coupling terms of the governing equations to adjust the coupling situation of the multi-field coupling equations. For example, when the effect of temperature change on the mechanical response can be ignored, the intermediate... and The relevant terms are set to zero; when the coupling between pore pressure change and mechanical response can be ignored, it is possible to... It will be set to zero. This flexible ability to adjust coupling terms is one of the important improvements of this embodiment compared to existing UMAT schemes.

[0041] Step two: In the finite element software platform, based on the element stiffness matrix and residual matrix obtained in step one, an element calculation module is constructed through a user-defined element interface to numerically solve the displacement field, pore water pressure field, and temperature field of the target porous media engineering structure. This embodiment is based on the ABAQUS platform and achieves the numerical solution of the above-mentioned fractional viscoelastic model of thermo-fluid-structure interaction through secondary development by calling the UEL interface.

[0042] Before performing numerical solutions through the user-defined element interface in the finite element software platform, the process includes generating and modifying the input file. The input file (i.e., the ABAQUS Input File, INP file) defines the basic information of the numerical model and serves as the input carrier for ABAQUS to call user subroutines for solution. The INP file includes element and node numbers, node spatial coordinates, material parameters, analysis step settings, boundary conditions, initial conditions, and mesh generation information. The specific process for generating and modifying the INP file is as follows.

[0043] The process of generating the initial INP file is as follows: (1) Establish the geometric model to be analyzed in the ABAQUS graphical modeling environment; (2) Apply displacement boundary conditions, pore water pressure boundary conditions and temperature boundary conditions or other predefined field variables to the model, and complete the finite element mesh generation; (3) Create an analysis job in the Job module of ABAQUS, and generate the corresponding INP file through the input file export function ("Write Input" command) for subsequent UEL subroutines to call.

[0044] In order to realize custom multi-field coupling analysis based on UEL subroutines, the generated INP file needs to be extended and modified so that it can correctly identify user-defined elements and related parameters. The specific modification operations are as follows: (1) Add user element declaration statements (*USER ELEMENT statement) to the INP file, and activate the displacement, pore water pressure and temperature degrees of freedom through the corresponding statements, corresponding to ABAQUS degree of freedom number 1 ( u ), 2 ( v ), 3 ( w ), 8 ( p (2) In order to specify the element region for UEL calculation in the model, create a user element set (UEL element set) in the INP file; (3) Add UEL parameter input statements to the INP file and write the required model parameters. Separate the parameters with commas and input them in the order of reading the material parameter array (PROPS array) in the UEL subroutine.

[0045] Regarding the writing of UEL subroutine code, firstly, following the ABAQUS help documentation, the standard interface form of the UEL subroutine, SUBROUTINE UEL(...), is established using Fortran, with parentheses containing several parameters passed from the main program. Since fractional-order models have historical memory characteristics, it is necessary to save the stress and strain data at each Gaussian integration point at each time step. In this embodiment, a common area is used in the UEL framework to store the weighting coefficients, stresses, and strains required for fractional-order thermo-fluid-structure interaction model calculations.

[0046] The matrix storing the weighting coefficients is a one-dimensional matrix, while the stress and strain matrices are four-dimensional matrices. UEL reads the material parameters defined in the INP file through the PROPS array; the parameter order must match the order of the corresponding parameters in the parameter input statements added to the INP file. The number of weighting functions is equal to the number of analysis steps in the numerical model. Each weighting function is stored in a matrix declared in the COMMON area and is called later when calculating fractional derivatives.

[0047] When constructing the unit calculation module, due to the historical memory characteristic of the fractional-order model, it is necessary to save the historical data at each Gaussian integration point at each time step. Specifically, a historical variable storage area is set up, using the COMMON area to store the historical weight coefficients, historical stress tensors, and historical strain tensors at each Gaussian integration point at each time step. The matrix storing the historical weight coefficients is a one-dimensional matrix, and the matrix storing the historical stress and strain tensors is a four-dimensional matrix. The historical weight coefficients, historical stress tensors, and historical strain tensors are all indexed by unit number, Gaussian integration point number, and time step number; that is, the position of each set of data in the storage matrix is ​​determined by the unit number, Gaussian integration point number, and time step number.

[0048] In this embodiment, the unit is a three-dimensional eight-node element, and its Gaussian integration points and weights are as follows: Table 1 Gaussian integration points and weight vectors

[0049] At the corresponding Gaussian integration points, the shape function formulas are as follows: ; ; ; ; ; ; ; ; in, These represent the coordinates of any point within the element in the local coordinate system, along with the coordinates of the three directions. Substituting the coordinates of the corresponding Gaussian integration point into the above shape function formula yields the shape function matrix at the corresponding Gaussian integration point. Construct the displacement shape function matrix from the shape functions. The pore water pressure shape function matrix and the temperature shape function matrix are both .

[0050] for example: ; .

[0051] Geometric matrix The construction method is as follows: For each Gaussian integration point, first construct the Jacobian matrix. : ; The Jacobian matrix is ​​composed of the partial derivatives of the physical coordinates of the unit nodes with respect to their local coordinates. Then, the derivatives of the shape functions with respect to the physical coordinates are calculated using the inverse of the Jacobian matrix. ; Finally, the geometric matrix is ​​constructed using the derivatives of the shape functions with respect to the physical coordinates. , .

[0052] Multiplying the shape function matrix by the corresponding nodal variable matrix yields the element pore water pressure and temperature; multiplying the geometric matrix... Multiplying the strain tensor by the nodal displacement matrix yields the element strain tensor. The data from this strain tensor is then recorded in the strain storage matrix declared in the COMMON area, with the position of this data set within the storage matrix determined by the element number, node number, and time step.

[0053] Elasticity matrix , The construction method is as follows: Calculate based on the parameters read from the INP file. ; Constructing the three-dimensional elastic stiffness matrix: ; Will Compared to Poisson Substituting into the above matrix yields the following result. , .

[0054] Based on the historical data in the strain storage matrix and the data in the weighting function matrix, the equivalent strain is calculated using the following formula: ; in, Calculated based on the difference between the temperature at the current time step and the temperature at the previous time step.

[0055] At the initial time step, i.e., when the time step count is 1, the stress is calculated using the following formula: ; The data is recorded in the stress storage matrix declared in the COMMON area.

[0056] When the time step count is greater than 1, the equivalent stress is first calculated based on the stress storage matrix: ; Then according to the matrix And the stress is calculated using the weighting function: ; The data is recorded in the stress storage matrix declared in the COMMON area.

[0057] Calculate the local stiffness matrix at each Gaussian integration point and integrate and sum them to assemble the element stiffness matrix: The local stiffness submatrix corresponding to the equilibrium equations is: ; The local stiffness submatrix corresponding to the mass equation (fluid continuity equation) is: ; The local stiffness submatrix corresponding to the energy equation is: ; By performing Gaussian integrals on each local stiffness submatrix and summing them, the global stiffness submatrix is ​​obtained: ; ; ; in, Number the Gaussian integration points. For the first The weights of each Gaussian integration point (each integration point in a 3D eight-node element has a weight of 1) is the determinant value of the Jacobian matrix.

[0058] The stiffness matrix of the three-dimensional eight-node element is 40×40, and the stiffness matrix of the assembled element is as follows: .

[0059] In the discrete form of the Gaussian integral points, each residual submatrix is ​​calculated as follows: ; in, The first The local residual submatrix calculated at each Gaussian integration point.

[0060] At each Gaussian integration point, the local residual submatrix The specific expression is: ; ; ; The unit residual matrix is ​​assembled from the above three residual submatrices: .

[0061] The numerical solution of the fractional-order thermo-fluid-structure interaction model in the ABAQUS software environment is performed as follows: (1) Create an analysis job: Create a new calculation task in the ABAQUS analysis environment and select to create the analysis model by input file, and then import the modified INP file; (2) Load user subroutines: Configure the user subroutine file in the calculation task settings interface, specify the path of the written UEL subroutine file, and associate the written UEL program with the current analysis task; (3) Submit the calculation job: After confirming that the above settings are correct, save the job configuration and submit the calculation; (4) Output and analyze the calculation results: After the solution is completed, an output database (ODB) format result file will be generated in the working directory. By exporting the calculation results, the reliability of the proposed numerical calculation method is verified by comparing them with the calculation results of the existing numerical implementation scheme.

[0062] The verification of the numerical method proposed in this embodiment can be accomplished by comparing it with analytical solutions. However, for thermo-fluid-structure interaction models considering fractional-order viscoelasticity, there is currently a lack of complete analytical solutions that can be directly used for comparison. Therefore, this embodiment adopts a degenerate verification method to verify the accuracy of the proposed numerical method. Specifically, the verification process is divided into two parts: on the one hand, the coupling model is degenerated into a pure mechanical model that does not consider pore pressure and temperature changes, which is used to verify the correctness of the fractional-order viscoelastic constitutive model; on the other hand, the fractional-order viscoelastic constitutive model is degenerated into an elastic constitutive model, transforming the model into a classical thermo-fluid-structure interaction model, which is used to verify the accuracy of each coupling term and the overall solution framework.

[0063] In the first part of the verification, the model degenerates into a purely mechanical problem with no pore pressure and no temperature change. To verify the correctness of the fractional-order viscoelastic model, a cube model with a side length of 1m is built in ABAQUS and divided into only one element. Since the model size is 1m, under uniaxial loading, its axial displacement and axial strain are numerically equal. Therefore, the axial strain response of the model can be directly characterized by the axial displacement of the top node. During the calculation, a rectangular pulse tensile stress in the Z direction is applied to the top of the element, and its mathematical expression is as follows: ; Under the above loading conditions, the strain of the model The response conforms to the following analytical solution expression: ; in, , It is the Mittag-Leffler function.

[0064] The strain response obtained from the numerical model was compared with the analytical solution, and the two showed good agreement, indicating that the fractional-order viscoelastic constitutive integral scheme used in this embodiment and its implementation in the ABAQUS user element are correct. The results are as follows: Figure 2 As shown.

[0065] In the second part of the verification, the fractional-order viscoelastic constitutive model is degenerated into a linear elastic constitutive model. Simultaneously, the thermal convection coupling terms in the energy equation are ignored, transforming the established model into a classical thermo-fluid-solid coupled model. This further verifies the correctness of the realization of mechanical deformation, pore pressure diffusion, temperature conduction, and their mutual coupling effects in the coupled control equations. Specifically, a high-resolution model is established in ABAQUS. A 1.2m cylindrical model was constructed, with displacement constraints applied to the bottom and sides of the model, as well as constraints on the fluid and heat fluxes at the corresponding boundaries. Pore water pressure and temperature boundaries were applied to the top of the model to create pore pressure and temperature gradients within the model, thereby inducing a thermal-fluid-structure interaction response.

[0066] Under these degradation conditions, the displacement, pore water pressure, and temperature response of the model can all be obtained analytically from classical thermo-fluid-structure interaction theory, and their expressions are as follows: ; ; ; in, , and Let Lamé be the constant. ， ， ， ， ， ， ， ， To simplify the parameters of the formula, see below: ; ; ; ; ; ; ; ; ; ; ; in, This is the initial pore pressure.

[0067] The displacement, pore water pressure, and temperature variations over time at model heights of 0.4m, 0.6m, and 0.8m were extracted and compared with the corresponding analytical solutions. The results show that the numerical solutions for displacement, pore water pressure, and temperature at different heights are in good agreement with the analytical solutions. The results are as follows: Figure 3 As shown in (a), (b), and (c), the verification results from the two parts above together demonstrate the accuracy and reliability of the numerical scheme proposed in this embodiment.

[0068] Therefore, this invention employs the aforementioned numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model. This method can autonomously determine the coupling terms in the multi-field coupling control equations and can perform numerical simulations of actual engineering models with complex geometry and boundary conditions in finite element software. This method solves the technical problems of existing UMAT schemes being unable to autonomously determine coupling terms and one-dimensional simplified schemes being unsuitable for complex engineering applications, providing a general numerical implementation framework for multi-field coupling analysis of fractional-order viscoelastic materials.

[0069] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit them. Although the present invention has been described in detail with reference to preferred embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the technical solutions of the present invention, and these modifications or equivalent substitutions cannot cause the modified technical solutions to deviate from the spirit and scope of the technical solutions of the present invention.

Claims

1. A numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model, characterized in that, Includes the following steps: Step 1: Establish fractional-order thermal-fluid-structure interaction (TFS) control equations. These equations include equilibrium equations, fractional-order Merchant constitutive equations, fluid continuity equations based on Darcy's law, and energy equations based on Fourier's heat transfer law. The fractional derivative terms in the fractional-order Merchant constitutive equations are discretized and approximated. The finite element discretization scheme for the fractional-order TFS control equations is then derived, resulting in an element stiffness matrix assembled from multiple stiffness sub-matrices and a residual matrix assembled from multiple residual sub-matrices. Step 2: In the finite element software platform, based on the element stiffness matrix and residual matrix obtained in Step 1, an element calculation module is constructed through a user-defined element interface to numerically solve the displacement field, pore water pressure field and temperature field of the target porous medium engineering structure. Based on the actual coupling requirements of the target multi-field coupling problem, one or more stiffness submatrices are selectively activated or disabled when assembling the element stiffness matrix. When assembling the element stiffness matrix, the matrix elements corresponding to the disabled stiffness submatrix are set to zero; when assembling the residual matrix, one or more residual submatrices are selectively activated or disabled. When assembling the residual matrix, the matrix elements corresponding to the disabled residual submatrix are set to zero.

2. The numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to claim 1, characterized in that, In step one, the discrete approximation uses the fractional derivative as defined by Grünwald-Letnikov. The definition of the Grünwald-Letnikov fractional derivative is: ； in, From arrive of Fractional derivative operators, Let be the function to be differentiated. For time step, For fractional order, For historical time step indexing, The total number of historical items, The coefficients are binomial coefficients; The discrete form of the Grünwald-Letnikov fractional derivative is approximated as: ； in, For the current moment, This is the current time step number. For the first The weighting coefficients for each historical time step For function exist The value of each time step; In step one, the finite element discretization scheme is derived using the integration by parts method.

3. The numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to claim 2, characterized in that, In step one, the fractional-order thermo-fluid-structure interaction governing equations include the following four sets of equations: Equilibrium equations: ； in, These are components of the stress tensor. , For tensor subscripts, Biot coefficient, Pore ​​water pressure, For Kroneck symbol, when hour, ,when hour, , These are the components of the volume force vector; Fractional Merchant constitutive equations: ； in, and For the elastic modulus parameter, The viscosity-time coefficient, From 0 to of Fractional derivative operators, Poisson's ratio, These are the components of the strain tensor. The coefficient of thermal expansion of a solid is _____. This refers to the change in temperature. Fluid continuity equation based on Darcy's law: ； in, The water storage coefficient of porous media, For time, Permeability coefficient, For fluid density, It is the acceleration due to gravity. For the Laplace operator, For volumetric strain, The coefficient of thermal expansion of the fluid. For temperature; Energy equation based on Fourier's law of heat transfer: ； in, For the effective volumetric heat capacity of porous media, For fluid density, Specific heat capacity of the fluid This is the Darcy velocity vector. For gradient operators, The thermal conductivity coefficient, The thermal expansion factor of solids, This is a reference temperature.

4. The numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to claim 3, characterized in that, In step one, when deriving the finite element discretization scheme using the integration by parts method, trial functions for the displacement field, pore water pressure field, and temperature field are introduced, resulting in the following weak form equations: Weak form of the equilibrium equation: ； in, To solve for the volume of the region, For virtual strain tensor, and For the elasticity matrix, For strain tensor, For the vector form of the Kronecker symbol, For equivalent stress, To account for the equivalent strain due to temperature effects, , For elasticity model parameters, For virtual displacement vectors, The volume force vector, superscript Indicates matrix transpose; Weak form of the fluid continuity equation: ； in, For virtual pore water pressure, For fluid specific weight; The weak form of the energy equation: ； in, This is a virtual temperature.

5. The numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to claim 4, characterized in that, The element stiffness matrix obtained in step one is assembled from the following stiffness sub-matrices: The stiffness submatrix corresponding to the equilibrium equations is: ； in, The stiffness submatrix corresponds to the displacement degrees of freedom. This is the stiffness submatrix coupled with displacement and pore water pressure. This is the stiffness submatrix coupled with displacement and temperature. For geometric matrices, For the elasticity matrix, For the shape function matrix, To solve for the region; The stiffness submatrix corresponding to the fluid continuity equation is: ； in, The stiffness submatrix is ​​the coupling between pore water pressure and displacement. This is the stiffness submatrix corresponding to the degree of freedom of pore water pressure. This is the stiffness submatrix coupled with pore water pressure and temperature; The stiffness submatrix corresponding to the energy equation is: ； in, This is the stiffness submatrix coupled with temperature and displacement. The stiffness submatrix is ​​a coupling matrix between temperature and pore water pressure. This is the stiffness submatrix corresponding to the temperature degree of freedom.

6. The numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to claim 5, characterized in that: The residual matrix obtained in step one is assembled from the residual submatrix of the equilibrium equation, the residual submatrix of the fluid continuity equation, and the residual submatrix of the energy equation.

7. The numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to claim 1, characterized in that: In step two, when constructing the unit calculation module, a historical variable storage area is set up to store the historical weight coefficients of each Gaussian integral point at each time step. Historical stress tensor and historical strain tensor Historical weighting coefficients, historical stress tensors, and historical strain tensors are all accessed and indexed by element number, Gaussian integration point number, and time step number.

8. The numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to claim 1, characterized in that, Step two, before performing numerical solutions through the user-defined element interface in the finite element software platform, also includes steps for generating and modifying input files: A numerical model is established in the finite element software platform, displacement boundary conditions, pore water pressure boundary conditions and temperature boundary conditions are applied, finite element mesh generation is completed, and the input file is exported. Add a user element declaration statement to the input file to activate the displacement degree of freedom, pore water pressure degree of freedom, and temperature degree of freedom. Create a set of user cells in the input file, specifying the cell region to be calculated by the cell calculation module built using the user-defined cell interface; Add a parameter input statement to the input file, and input the material parameters in the order they are read from the parameter array in the user-defined unit interface.

9. The numerical implementation method for thermo-fluid-structure interaction based on a fractional-order viscoelastic model according to claim 1, characterized in that: The displacement degree of freedom corresponds to the translational degree of freedom along three orthogonal spatial directions in the finite element software platform, the pore water pressure degree of freedom is the fluid pressure degree of freedom, and the temperature degree of freedom is the thermodynamic temperature degree of freedom.

Images