Finite element and artificial boundary coupled wave field computation method and system

Abstract

The application provides a wave field calculation method and system coupling finite elements and artificial boundaries, and relates to the technical field of underground structure anti-seismic technology. In view of the problems that false reflection is prone to be generated in the artificial boundary of the finite element, the calculation efficiency of the boundary element method is low, the material nonlinearity is difficult to handle, and coupling is difficult, the application divides a calculation domain into an external infinite domain and an internal finite domain, solves the external infinite domain by using a fast multipole boundary element method, obtains artificial boundary wave field data, converts the wave field data into equivalent node forces on boundary nodes of the internal finite domain by using an artificial boundary substructure method, and performs non-linear dynamic response analysis on the internal finite domain subjected to the equivalent node forces by using the finite element method, so that the defects of low calculation efficiency and difficult nonlinearity handling are overcome, the ability of the finite element method in handling complex nonlinearity is retained, the unification of high precision and high calculation efficiency is realized, and the stability and reliability of seismic dynamic response analysis are improved.

Description

Technical Field

[0001] This invention relates to the field of seismic resistance technology for underground structures, and in particular to a wave field calculation method and system that couples finite element method and artificial boundary. Background Technology

[0002] In the seismic design and safety evaluation of underground structures, seismic dynamic response analysis is a crucial step. To accurately assess the mechanical behavior of underground structures under seismic loading, numerical methods are typically used to simulate the propagation process of seismic waves in infinite or semi-infinite media and their complex dynamic interactions with underground structures.

[0003] Currently, numerical simulation methods for infinite domain wave problems mainly include the finite element method and the boundary element method. However, these methods have the following shortcomings in practical applications: The finite element method requires setting artificial boundaries at the boundary of the computational domain to simulate the radiation damping effect in an infinite domain. However, traditional artificial boundaries are mostly based on simplified wave theory, which has limited accuracy in simulating complex incident wave fields (such as oblique incidence and topographic effects). They are difficult to accurately reflect the back propagation effect of external scattered waves, and are prone to producing false boundary reflections, leading to distorted calculation results.

[0004] Although the boundary element method is applicable to infinite domain fluctuation problems, its coefficient matrix is ​​an asymmetric full matrix, and the storage requirements and computational cost increase sharply with the problem size, resulting in low computational efficiency. At the same time, the traditional boundary element method is difficult to handle material nonlinearity problems, and when coupled with the finite element method, there are difficulties such as variable interpolation and matrix property mismatch, resulting in poor solution stability.

[0005] Therefore, there is an urgent need for a wave field simulation method that can both efficiently simulate infinite-domain fluctuations and ensure computational accuracy under complex wave fields, in order to overcome the aforementioned technical deficiencies. Summary of the Invention

[0006] To address the problems of spurious reflections easily generated by artificial boundaries in the finite element method (FEM), low computational efficiency, difficulty in handling material nonlinearity, and coupling difficulties in existing technologies, this invention proposes a wave field calculation method that couples finite element analysis with artificial boundaries, specifically including the following steps: A wave field calculation method coupled with finite element method and artificial boundary, characterized by comprising: Obtain the computational domain and the seismic incident wave field; The computational domain is divided into an outer infinite domain and an inner finite domain based on an artificial boundary, wherein the artificial boundary is set outside the nonlinear influence range of the structure-surrounding rock system. The boundary wave field data on the artificial boundary is obtained by solving the external infinite domain based on the fast multipole boundary element method and the seismic incident wave field. The boundary wave field data is transformed into equivalent nodal forces acting on the boundary nodes of the internal finite domain based on the artificial boundary substructure method. The equivalent nodal force is applied as a seismic input load to the boundary nodes of the internal finite domain. The finite element method is used to perform nonlinear dynamic response analysis on the internal finite domain to which the equivalent nodal force is applied, and the seismic dynamic response results of the internal finite domain are obtained.

[0007] Furthermore, the determination of the artificial boundary includes: Calculate the extent of the plastic zone of the surrounding rock under seismic loading based on the geometric dimensions, burial depth, and mechanical properties of the underground structure. Obtain the safety distance outside the plastic zone of the surrounding rock; The artificial boundary is determined based on the safety distance so that the medium on the artificial boundary satisfies the linear elasticity assumption.

[0008] Furthermore, after determining the artificial boundary based on the safety distance, the method further includes: The region within the artificial boundary is discretized using finite element methods to obtain a finite element mesh, wherein the finite element mesh is used for nonlinear dynamic response analysis. The artificial boundary is discretized based on the boundary element model to obtain the boundary element mesh; Based on the earthquake incident wave field and the finite element method, the finite element dynamic equations corresponding to the finite element mesh are solved to obtain the first nodal displacement and the first nodal surface force of each node on the artificial boundary. Based on the earthquake incident wave field and the fast multipole boundary element method, the boundary integral equation corresponding to the boundary element mesh is solved to obtain the second nodal displacement and the second nodal surface force at the corresponding position on the artificial boundary. Calculate the displacement difference between the displacement of the first node and the displacement of the second node; Calculate the vector sum of the forces on the first and second nodal surfaces; If both the displacement difference and the vector sum satisfy the corresponding preset thresholds, it is determined that the propagation of the seismic incident wave field at the artificial boundary conforms to physical laws.

[0009] Furthermore, if at least one of the displacement difference and the vector sum does not satisfy a corresponding preset threshold, the method further includes: The position offset and offset direction of the artificial boundary are determined based on the distribution characteristics of the displacement difference on the artificial boundary, and the position of the artificial boundary is adjusted based on the position offset and offset direction. The mesh refinement region is determined based on the vector and the coordinates of the extreme points on the artificial boundary, and the finite element mesh and boundary element mesh of the mesh refinement region are refined and densified. Based on the adjusted artificial boundary and the refined and densified finite element mesh and boundary element mesh, the displacement difference and vector sum are recalculated until the recalculated displacement difference and the vector sum both meet the corresponding preset threshold.

[0010] Furthermore, the solution to the external infinite domain based on the fast multipole boundary element method and the seismic incident wave field includes: Boundary integral equations for the external infinite domain are established based on the boundary element grid of artificial boundaries and the seismic incident wave field. The kernel function in the boundary integral equation is expanded and recursively calculated using a fast multipole expansion algorithm to obtain the displacement vector and surface force vector of each node on the artificial boundary. The displacement vector and surface force vector are used as the boundary wave field data.

[0011] Furthermore, the method of converting the boundary wavefield data into equivalent nodal forces acting on the boundary nodes of the internal finite domain based on the artificial boundary substructure method includes: Obtain the displacement vector and surface force vector of the artificial boundary nodes in the boundary wave field data; Construct a transformation matrix from boundary wave field data to equivalent nodal forces based on the stiffness matrix and mass matrix in the artificial boundary substructure method; Substituting the displacement vector and the surface force vector into the transformation matrix, the equivalent nodal force vector on each boundary node is obtained. The equivalent node force vectors are sorted based on the node numbers to obtain the equivalent node force load sequence.

[0012] Furthermore, the equivalent nodal force is applied as a seismic input load to the boundary nodes of the internal finite domain, and the finite element method is used to perform nonlinear dynamic response analysis on the internal finite domain to obtain the seismic dynamic response results of the internal finite domain, including: Based on the time points of the ground motion time history, the equivalent nodal force load sequence is applied one by one to the boundary nodes corresponding to the internal finite domain. At each time point, the finite element dynamic equations are solved based on the time integration algorithm to obtain the displacement, velocity, and acceleration responses of all nodes in the internal finite domain; The displacement, velocity, and acceleration responses of all nodes at all time points are sorted in chronological order to form the seismic dynamic response results.

[0013] A wave field calculation system coupled with finite element method and artificial boundary, the system employing a wave field calculation method coupled with finite element method and artificial boundary as described in any of the preceding claims, specifically including the following modules: The acquisition module is used to acquire the computational domain and the seismic incident wave field; A partitioning module, connected to the acquisition module, is used to partition the computational domain into an outer infinite domain and an inner finite domain based on an artificial boundary, wherein the artificial boundary is set outside the nonlinear influence range of the structure-surrounding rock system; The solution module, connected to the partitioning module, is used to solve the external infinite domain based on the fast multipole boundary element method and the seismic incident wave field to obtain the boundary wave field data on the artificial boundary. A conversion module, connected to the solution module, is used to convert the boundary wave field data into equivalent nodal forces acting on the boundary nodes of the internal finite domain based on the artificial boundary substructure method. The analysis module, connected to the partitioning module and the conversion module, is used to apply the equivalent nodal force as a seismic input load to the boundary nodes of the internal finite domain, and to perform nonlinear dynamic response analysis on the internal finite domain to which the equivalent nodal force is applied using the finite element method, thereby obtaining the seismic dynamic response results of the internal finite domain.

[0014] Compared with the prior art, the beneficial effects of the present invention are as follows: Firstly, the computational domain is divided into an outer infinite domain and an inner finite domain. The outer infinite domain is solved using the fast multipole boundary element method to obtain artificial boundary wavefield data. The wavefield data is then transformed into equivalent nodal forces on the boundary nodes of the inner finite domain using the artificial boundary substructure method. The finite element method is then used to perform nonlinear dynamic response analysis on the inner finite domain with the applied equivalent nodal forces. This overcomes the shortcomings of low computational efficiency and difficulty in nonlinear processing, while retaining the ability of the finite element method to handle complex nonlinearities. This achieves a balance between high accuracy and high computational efficiency, and improves the stability and reliability of seismic dynamic response analysis. Secondly, by setting the artificial boundary outside the nonlinear influence range of the structure-surrounding rock system, the medium on the artificial boundary always satisfies the linear elastic assumption, effectively avoiding the distortion of calculation results caused by the nonlinearity of the medium at the artificial boundary. By calculating the range of the plastic zone of the surrounding rock and reserving a safe distance, the artificial boundary setting is not too conservative, which would lead to an unnecessary expansion of the internal finite domain calculation scale. Thus, the overall calculation efficiency is improved while ensuring the calculation accuracy. Attached Figure Description

[0015] To more clearly illustrate the specific embodiments of the present invention or the technical solutions in the prior art, the drawings used in the description of the specific embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained from these drawings without creative effort.

[0016] Figure 1This is a flowchart illustrating a wave field calculation method that couples finite element analysis with artificial boundaries. Figure 2 This is a structural block diagram of a wave field calculation system that couples finite element analysis and artificial boundaries. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below. Obviously, the described embodiments are only a part of the embodiments of this invention, and not all of them. Based on the embodiments of this invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this invention.

[0018] The specific embodiments of the present invention will be described below.

[0019] To address the issues of spurious reflections easily generated by artificial boundaries in the finite element method (FEM), low computational efficiency, difficulty in handling material nonlinearities, and coupling difficulties in the FEM, this invention divides the computational domain into an outer infinite domain and an inner finite domain. The outer infinite domain is solved using the fast multipole FEM to obtain artificial boundary wavefield data. The wavefield data is then transformed into equivalent nodal forces on the boundary nodes of the inner finite domain using the artificial boundary substructure method. The FEM is then applied to perform nonlinear dynamic response analysis on the inner finite domain under these equivalent nodal forces. This approach overcomes the shortcomings of low computational efficiency and difficulty in handling nonlinearities while retaining the FEM's ability to handle complex nonlinearities. It achieves a balance between high accuracy and high computational efficiency, improving the stability and reliability of seismic dynamic response analysis.

[0020] Example 1 like Figure 1 As shown, this invention proposes a wave field calculation method coupled with finite element method and artificial boundary, which specifically includes the following steps: Step S1: Obtain the computational domain and the seismic incident wave field; In this embodiment, the geometric range of the computational domain is determined according to the analysis requirements of the actual engineering project. The computational domain includes the underground structure and the surrounding rock area extending outward from the outer edge of the underground structure. The surrounding rock area is a continuous medium, and its extension range should be large enough to reduce the error introduced by subsequent artificial boundaries. The geometric dimensions of the computational domain are determined based on the structural characteristic dimensions, the wavelength corresponding to the dominant seismic wave frequency, and the wave velocity of the surrounding rock. For example, the distance from the boundary of each direction of the computational domain to the outer edge of the structure is not less than 1 to 2 times the wavelength of the dominant seismic wave.

[0021] In this embodiment, the incident wave field of the earthquake can be obtained from the ground motion data recorded by the seismic station. For example, the site conditions, such as the soil type, overburden thickness, and shear wave velocity, can be selected to match the strong earthquake records of the actual engineering site. If no suitable records are available, simulations can be performed based on the seismic wave propagation theory. For example, the synthetic ground motion method can be used to synthesize artificial ground motion time histories that conform to the seismic hazard level of the actual engineering site based on the target response spectrum and power spectral density function. The incident wave field includes acceleration, velocity, or displacement time histories.

[0022] Step S2: Divide the computational domain into an outer infinite domain and an inner finite domain based on an artificial boundary, wherein the artificial boundary is set outside the nonlinear influence range of the structure-surrounding rock system; In this embodiment, the structure-surrounding rock system is an interactive system formed by the underground structure and its surrounding rock and soil.

[0023] The determination of the artificial boundary includes: calculating the range of the plastic zone of the surrounding rock under seismic action based on the geometric dimensions, burial depth and mechanical properties of the underground structure; obtaining the safety distance outside the range of the plastic zone of the surrounding rock; and determining the artificial boundary based on the safety distance so that the medium on the artificial boundary satisfies the linear elastic assumption.

[0024] In this embodiment, the extent of the plastic zone of the surrounding rock can be determined by nonlinear finite element analysis. For example, a two-dimensional or three-dimensional finite element model including the underground structure and the surrounding rock can be established. The constitutive model of the surrounding rock, such as the Mohr-Coulomb model or the Drucker-Prag model, can be input, and an incident seismic wave field can be applied. The boundary of the plastic strain or yield region can be determined by elastoplastic time history analysis. This boundary is the extent of the plastic zone of the surrounding rock. Alternatively, the extent of the plastic zone of the surrounding rock can be calculated based on empirical formulas or semi-analytical methods. For example, based on the size and burial depth of the underground structure, the strength parameters of the surrounding rock, and the seismic intensity parameters, the formulas provided in the seismic design code for geotechnical engineering or related research results can be used to calculate the approximate extent of the plastic zone of the surrounding rock. For example, for a circular tunnel, the radius of the plastic zone can be calculated using the modified Fenner formula or the Kastner formula.

[0025] In this embodiment, the safety distance can be determined based on engineering experience, specification requirements, or sensitivity analysis results. For example, the safety distance can be a certain percentage of the maximum size of the plastic zone, such as 10%–30%, or a fixed distance value, such as 5 meters or 10 meters. It can also be determined through parametric studies. For example, after initially determining the plastic zone, multiple simulations can be performed by changing the safety distance to observe whether the stress-strain response of the medium at the artificial boundary satisfies the linear elastic assumption, and the minimum safety distance that satisfies the linear elastic assumption can be selected as the optimal distance to avoid unnecessary expansion of the computational domain.

[0026] In this embodiment, the range of the plastic zone of the surrounding rock can be extended outward by a safe distance, and the extended boundary can be used as an artificial boundary. When the plastic zone of the surrounding rock is an irregular shape, its outer envelope can be uniformly extended outward by a safe distance, or the radius of the artificial boundary can be the distance from the farthest point of the plastic zone of the surrounding rock to the center plus the safe distance, with the center of the structure as the origin. The artificial boundary can also be determined by iterative or optimization methods. After the artificial boundary position is initially determined, a preliminary dynamic response analysis is performed to check whether the stress-strain state at the artificial boundary meets the linear elastic assumption. When nonlinearity occurs in a local area, the artificial boundary of the area is moved outward by a preset step size along the outer normal direction. For example, the current boundary point is moved by 5% to 10% of the distance between the current boundary point and the center of the structure, or a fixed distance value, such as 1 meter to 2 meters, to further move it away from the nonlinear region. The above outward movement, analysis, and checking process is repeated until all boundary points meet the linear elastic assumption. The artificial boundary position at this time is the final position.

[0027] In this embodiment, by setting the artificial boundary outside the nonlinear influence range of the structure-surrounding rock system, the medium on the artificial boundary always satisfies the linear elastic assumption, avoiding the distortion of calculation results caused by the nonlinearity of the medium at the artificial boundary. By calculating the range of the plastic zone of the surrounding rock and reserving a safe distance, the artificial boundary setting is not set too conservatively, which would lead to an unnecessary expansion of the internal finite domain calculation scale. Thus, the overall calculation efficiency is improved while ensuring the calculation accuracy.

[0028] After determining the artificial boundary based on the safety distance, the method further includes: discretizing the region within the artificial boundary using finite element methods to obtain a finite element mesh, wherein the finite element mesh is used for nonlinear dynamic response analysis; discretizing the artificial boundary based on the boundary element model to obtain a boundary element mesh; solving the finite element dynamic equations corresponding to the finite element mesh based on the seismic incident wave field and the finite element method to obtain the first nodal displacement and the first nodal surface force of each node on the artificial boundary; solving the boundary integral equations corresponding to the boundary element mesh based on the seismic incident wave field and the fast multipole boundary element method to obtain the second nodal displacement and the second nodal surface force at the corresponding position on the artificial boundary; calculating the displacement difference between the first nodal displacement and the second nodal displacement; calculating the vector sum of the first nodal surface force and the second nodal surface force; and determining that the transmission of the seismic incident wave field at the artificial boundary conforms to physical laws when both the displacement difference and the vector sum satisfy the corresponding preset thresholds.

[0029] Finite element discretization involves dividing a continuous finite domain into a finite number of non-overlapping elements. These elements and their nodes approximate the physical behavior of the continuum. In this embodiment, for geometrically regular regions, structured meshing methods such as quadrilateral or hexahedral meshes can be used for finite element discretization. The mesh density is determined based on the shortest wavelength of the seismic wave; for example, the element size is 1 / 8 to 1 / 5 of the shortest wavelength. The element type is determined based on the analysis dimension and physical assumptions, such as plane stress elements, plane strain elements, or solid elements. For geometrically complex regions, unstructured meshing methods such as triangular or tetrahedral meshes can be used for finite element discretization. The mesh is automatically generated using algorithms such as Delaunay triangulation or Advancing Front. After meshing, material properties, boundary conditions, and load cases are assigned to the elements in the mesh to form the finite element model of the internal finite domain.

[0030] In this embodiment, linear or quadratic isoparametric elements can be used to discretize the artificial boundary, thereby dividing the artificial boundary into a series of boundary elements. Each element contains several nodes, which are used to describe the displacement and surface force on the boundary. Higher-order boundary element elements, such as cubic or quartic elements, can also be used to improve the accuracy of boundary discretization. The boundary element model includes a boundary element mesh, boundary integral equations, basic solutions, and solution algorithms. In this embodiment, the boundary element model is used to simulate the mechanical response of an external infinite domain at the artificial boundary.

[0031] In this embodiment, the seismic incident wave field is applied as an excitation load to the finite element mesh. For example, it can be applied to the artificial boundary nodes through equivalent nodal forces or to the bottom of the finite element model through a base input method. Explicit or implicit time integration algorithms, such as the Newmark-β method or the central difference method, are used to solve the finite element dynamic equations, thereby obtaining the first nodal displacement and the first nodal surface force of the artificial boundary nodes. Alternatively, the frequency domain finite element method can be used to solve the finite element dynamic equations. The seismic incident wave field is subjected to a Fourier transform, and the finite element equations are solved in the frequency domain. Then, the displacement and surface force responses in the time domain, i.e., the first nodal displacement and the first nodal surface force, are obtained through an inverse Fourier transform.

[0032] It should be noted that the finite element equation is the structural dynamic equilibrium equation obtained by discretizing the internal finite domain within the artificial boundary using finite elements. The specific form of the finite element dynamic equation and its standard solution method are well known to those skilled in the art, and will not be elaborated here.

[0033] In this embodiment, the seismic incident wave field can be applied to the boundary element model as a free field response or equivalent load, and the solution process of the boundary integral equation can be accelerated by using a fast multipole expansion algorithm, thereby obtaining the second nodal displacement and the second nodal surface force of the artificial boundary node; in addition, other acceleration algorithms such as H-matrix or adaptive cross approximation can be used to solve the boundary integral equation, thereby obtaining the second nodal displacement and the second nodal surface force on the artificial boundary.

[0034] The displacement difference is an indicator for evaluating whether the displacement continuity at the artificial boundary meets the physical laws. In this embodiment, the displacement vectors of the corresponding nodes on the artificial boundary can be subtracted point by point to obtain the displacement difference vector of each node. The magnitude or a norm of the displacement difference vector, such as the L2 norm, can be calculated as the overall displacement difference. Alternatively, the relative error between the displacement of the first node and the displacement of the second node can be calculated. For example, the displacement difference can be divided by the magnitude of a reference displacement (such as the maximum displacement or the average displacement) to obtain the displacement difference index.

[0035] According to physical principles, at an artificial boundary, the surface forces applied to the boundary by the internal finite domain and the surface forces applied to the boundary by the external infinite domain should be in balance, that is, their vector sum should approach zero. In this embodiment, the surface force vectors of the corresponding nodes on the artificial boundary can be directly added point by point to obtain the surface force vector sum of each node, and the magnitude or norm of the force vector sum can be calculated as the overall surface force vector sum. Alternatively, the relative ratio of the surface force vector sum to a certain reference surface force (such as the maximum surface force or the average surface force) can be calculated to assess the degree of force balance.

[0036] In this embodiment, the preset threshold can be determined based on engineering experience, the accuracy requirements of numerical simulation, or the results of previous parameter studies. For example, the preset threshold for displacement difference can be a percentage of the maximum displacement, and the preset threshold for the sum of surface force vectors can be a percentage of the maximum surface force. An adaptive threshold setting method can also be used to dynamically adjust the preset threshold based on factors such as the characteristics of the computational domain, the frequency components of seismic waves, or the grid density.

[0037] In this embodiment, the region within the artificial boundary is discretized using the finite element method, and the first nodal displacement and first nodal surface force on the artificial boundary are obtained based on the seismic incident wave field and the finite element method. This provides data on the response of the artificial boundary from the perspective of the internal finite domain. The artificial boundary is discretized based on the boundary element model, and the second nodal displacement and second nodal surface force at the corresponding position on the artificial boundary are obtained using the fast multipole boundary element method. This provides data on the response of the artificial boundary from the perspective of the external infinite domain. By calculating the displacement difference and vector sum, the displacement continuity and force balance characteristics at the artificial boundary are evaluated. When both the displacement difference and vector sum meet the preset threshold, it can be determined that the transmission of the seismic incident wave field at the artificial boundary conforms to the physical laws. This ensures the rationality of the artificial boundary position and the accuracy of the discrete mesh, avoiding distortion of subsequent wave field calculations and nonlinear dynamic response analysis results caused by improper artificial boundary settings or insufficient mesh discretization.

[0038] If at least one of the displacement difference and the vector sum does not meet the corresponding preset threshold, the method further includes: determining the position offset and offset direction of the artificial boundary based on the distribution characteristics of the displacement difference on the artificial boundary, and adjusting the position of the artificial boundary based on the position offset and offset direction; determining the mesh refinement region based on the extreme point coordinates of the vector sum on the artificial boundary, and refining the finite element mesh and boundary element mesh of the mesh refinement region; recalculating the displacement difference and the vector sum based on the adjusted artificial boundary and the refined finite element mesh and boundary element mesh, until the recalculated displacement difference and the vector sum both meet the corresponding preset threshold.

[0039] In this embodiment, the distribution characteristics reflect whether the position of the artificial boundary is reasonable and the degree and direction of its deviation from the reasonable position. For example, when the displacement difference shows a systematic positive or negative deviation in a certain area, it indicates that the artificial boundary is too close to or too far from the actual infinite domain boundary in that area. The distribution characteristics of the displacement difference are determined by numerical analysis or empirical judgment to determine the distance and direction of the artificial boundary to be moved. For example, the least squares method can be used to fit the distribution curve of the displacement difference to find its centroid or trend line, thereby determining the overall offset and direction. Alternatively, the location of the maximum or minimum displacement difference can be analyzed. When the maximum value appears in a certain area, it indicates that the artificial boundary in that area is too close and should be moved outward. When the minimum value appears in a certain area, it indicates that the artificial boundary in that area is too far and should be moved inward. The moving distance can be an estimated value of the ratio of the displacement difference to the displacement gradient at the extreme point, or it can be gradually adjusted according to a preset step size, such as 5% to 10% of the distance between the boundary point and the center of the structure.

[0040] To determine the position offset and offset direction, the geometric position of the artificial boundary can be modified. For example, the artificial boundary can be moved to a new position by means of translation, scaling or local deformation to ensure that the medium on the artificial boundary satisfies the linear elastic assumption and is outside the nonlinear influence range of the structure-surrounding rock system.

[0041] In this embodiment, the spatial location of the point where the value of the vector sum reaches its maximum or minimum is found, that is, the coordinates of the extreme point of the vector sum on the artificial boundary. The extreme point indicates that the gradient change of the seismic incident wave field in these areas is drastic, or there is stress concentration, indicating that the current grid discretization accuracy is insufficient to capture these details.

[0042] In this embodiment, one or more local regions where the mesh density needs to be increased are delineated based on the coordinates of the extreme points of the vector sum, i.e., the mesh refinement region is determined. For example, a preset range can be extended outward from the extreme point as the refinement region. The preset range is determined based on the shortest wavelength of the seismic wave, the typical size of the current mesh, the gradient attenuation distance of the vector sum, or a fixed empirical value. For example, it can be 1 / 4 to 1 / 2 of the shortest wavelength of the seismic wave, or 3 to 5 times the typical size of the current mesh. The boundary of the refinement region is also determined based on the gradient change of the vector sum around the extreme point. Within the determined mesh refinement region, the number of elements in the finite element mesh and the boundary element mesh is increased, and the element size is reduced, thereby improving the discretization accuracy of the region. For example, methods such as quadtree / octree subdivision, local mesh re-division, or adaptive meshing techniques are used to locally increase the mesh resolution without affecting the mesh density of other regions, so as to accurately capture the details of the seismic incident wave field.

[0043] In this embodiment, after adjusting the position of the artificial boundary and refining the mesh, the steps of calculating the first node displacement, the second node displacement, the first node surface force, the second node surface force, the displacement difference, and the vector sum on the artificial boundary are performed again until the new displacement difference and vector sum are all within the corresponding preset threshold. The preset threshold is set according to engineering experience or accuracy requirements.

[0044] In this embodiment, by analyzing the distribution characteristics of displacement difference on the artificial boundary, the degree and direction of deviation of the current artificial boundary position are determined, thereby adjusting the position of the artificial boundary in a targeted manner to ensure that it is always outside the nonlinear influence range of the structure-surrounding rock system and meets the linear elasticity assumption of the medium. This avoids false reflections and calculation distortions caused by improper artificial boundary settings. By identifying the extreme point coordinates of vectors on the artificial boundary, areas with insufficient mesh discretization accuracy are located, and only the finite element mesh and boundary element mesh of these local areas are refined and densified, avoiding unnecessary mesh densification of the entire computational domain and reducing the amount of computation to a certain extent.

[0045] Step S3: Solve the external infinite domain based on the fast multipole boundary element method and the seismic incident wave field to obtain the boundary wave field data on the artificial boundary; Specifically, a boundary integral equation for an external infinite domain is established based on the boundary element mesh of the artificial boundary and the seismic incident wave field; a fast multipole expansion algorithm is used to expand and recursively calculate the kernel function in the boundary integral equation to obtain the displacement vector and surface force vector of each node on the artificial boundary; the displacement vector and surface force vector are used as the boundary wave field data.

[0046] In this embodiment, a direct boundary integral equation can be used to decompose the seismic incident wave field into a free field and a scattered field. Based on the Helmholtz integral equation and the Sommerfeld radiation condition, an integral equation involving only displacement and surface force on the artificial boundary can be derived. The derivation process and specific form of this boundary integral equation are well known to those skilled in the art and will not be elaborated here. Alternatively, the weighted residual method can be used to perform weighted integration of the wave equation on the artificial boundary and obtain the boundary integral equation using integration by parts.

[0047] Among them, the fast multipole expansion algorithm is a technique to solve the problem of full-rank coefficient matrix and large computational cost in the boundary element method. It significantly reduces the computational complexity by performing multipole expansion and local expansion of far-field interactions and using recursive relations, thereby improving computational efficiency. In this embodiment, the fast multipole boundary element method is extended from solving conventional seismic wave scattering problems to being applied to the finite element-boundary element coupled framework to solve the wavefield response of the external infinite domain, so that it can directly provide the required boundary wavefield input for the artificial boundary substructure method.

[0048] In this embodiment, a classic fast multipole expansion algorithm can be used to divide the nodes on the artificial boundary into boxes of different levels. The interaction between far-field boxes is expanded in multiple stages, and the interaction between near-field boxes is calculated directly. Efficient calculation is achieved through upward aggregation and downward distribution. Alternatively, a generalized fast multipole algorithm can be used to further optimize the computational performance by expanding and approximating the kernel function. In this algorithm, the displacement vector and surface force vector are direct results obtained by the boundary element method, representing the motion state and force conditions of each point on the artificial boundary.

[0049] In this embodiment, the displacement vector and surface force vector contain the necessary information for the propagation of the seismic incident wave field on the artificial boundary. For example, the calculated displacement vector and surface force vector are stored as a two-dimensional array or matrix, where each row or column corresponds to the displacement and surface force components of a boundary node at different time steps.

[0050] It should be noted that in this embodiment, a one-way coupling strategy is adopted, that is, it is assumed that because the artificial boundary is far enough away, the secondary scattered waves generated by the internal finite field have been greatly attenuated when they propagate to the boundary, and their influence on the boundary wave field can be ignored.

[0051] In this embodiment, a boundary integral equation for an external infinite domain is established based on an artificial boundary element mesh and the seismic incident wave field. Using the already discretized artificial boundary mesh and the known input seismic wave field, a control equation conforming to the wave propagation law of the infinite domain is constructed, thus ensuring that the solution model does not deviate from the setting of the actual problem. A fast multipole expansion algorithm is used to expand and recursively calculate the kernel function in the boundary integral equation, thereby quickly classifying and processing the kernel function interactions between different nodes in a hierarchical manner. The interaction terms that originally required full matrix storage are transformed into an approximate sparse storage form, reducing the storage requirements and computational load of the solution process. This solves the problem of the sharp decline in computational efficiency after the scale of the traditional boundary element method increases. The calculated displacement vector and surface force vector are used as boundary wave field data without additional interpolation and transformation processing, which simplifies the solution process and ensures the original accuracy of the boundary wave field data.

[0052] Step S4: Based on the artificial boundary substructure method, the boundary wave field data is converted into equivalent nodal forces acting on the boundary nodes of the internal finite domain; Specifically, the displacement vectors and surface force vectors of the artificial boundary nodes in the boundary wave field data are obtained; a transformation matrix from the boundary wave field data to the equivalent nodal forces is constructed based on the stiffness matrix and mass matrix in the artificial boundary substructure method; the displacement vectors and surface force vectors are substituted into the transformation matrix to obtain the equivalent nodal force vectors on each boundary node; the equivalent nodal force vectors are sorted based on the node number to obtain the equivalent nodal force load sequence.

[0053] In this embodiment, the boundary wave field data in step S103 is obtained to obtain the displacement vector and surface force vector. When the boundary element mesh does not completely coincide with the boundary node mesh of the internal finite domain, the displacement and surface force data obtained by the boundary element solution can be mapped to the boundary nodes of the internal finite domain by interpolation methods, such as shape function interpolation, radial basis function interpolation or kriging interpolation, so as to obtain the corresponding displacement vector and surface force vector.

[0054] The transformation matrix can be constructed based on the dynamic stiffness matrix of the artificial boundary substructure. In this embodiment, the artificial boundary substructure is discretized using the finite element method to obtain its stiffness and mass matrices. Based on the dynamic equilibrium relationship in the substructure method, the transformation matrix is ​​constructed using the stiffness and mass matrices. For example, the displacement vector and surface force vector in the boundary wave field data are used as inputs. By solving the dynamic equilibrium equations of the artificial boundary substructure, the displacements and surface forces at the boundary nodes are transformed into equivalent nodal forces. The construction of the transformation matrix can be completed in the frequency domain or the time domain. In the frequency domain... The transformation matrix involves the combined operation of the stiffness matrix, mass matrix, and excitation frequency. For each preset frequency point, the dynamic stiffness matrix at that frequency is formed according to the dynamic equilibrium relationship, which is the frequency domain transformation matrix. In the time domain, the transformation matrix is ​​obtained by substituting the stiffness matrix and mass matrix into the state space equation or recursive scheme. For example, the time step is discretized using the central difference method or the Newmark-β method to derive the recursive coefficients between the equivalent nodal force at the current time and the displacement and surface force at the current and historical times. Arranging these coefficients into a matrix form yields the recursive transformation matrix in the time domain.

[0055] It should be noted that the specific values ​​of the transformation matrix can be pre-calculated using the cohesive method of the artificial boundary substructure. The cohesive method refers to eliminating the degrees of freedom of the nodes inside the substructure and retaining only the degrees of freedom of the artificial boundary nodes, thereby obtaining the equivalent stiffness matrix, equivalent mass matrix and equivalent damping matrix of the boundary nodes.

[0056] In this embodiment, the obtained displacement vector and surface force vector are multiplied with the constructed transformation matrix to calculate the equivalent nodal force vector acting on each boundary node. For some complex nonlinear or time-varying systems, the transformation process needs to be carried out through iterative solution or step-by-step integration to ensure mechanical balance and energy conservation in each time step or iteration step, thereby obtaining the accurate equivalent nodal force.

[0057] In this embodiment, a standard sorting algorithm, such as quicksort, mergesort, or heapsort, can be used to rearrange the equivalent node force vectors according to the node number, so that they form an ordered load sequence according to the node order of the internal finite field.

[0058] In this embodiment, the finite element method no longer directly relies on traditional artificial boundary conditions. Instead, it uses the equivalent nodal forces generated by the combined processing of the fast multipole boundary element method and the artificial boundary substructure method as the seismic input load. Since the fast multipole boundary element method directly solves the total wave field (the sum of incident and scattered waves) of the external infinite domain and the transformation process is a faithful mapping, the equivalent nodal forces include all the wave information of the external infinite domain. This allows the nonlinear analysis of the internal finite domain to be based on the accurate external wave field input, thereby improving the overall calculation accuracy to a certain extent.

[0059] Step S5: Apply the equivalent nodal force as a seismic input load to the boundary nodes of the internal finite domain, and perform nonlinear dynamic response analysis on the internal finite domain to which the equivalent nodal force is applied, thereby obtaining the seismic dynamic response results of the internal finite domain.

[0060] Specifically, based on the time points of the seismic motion time history, the equivalent nodal force load sequence is applied one by one to the boundary nodes corresponding to the internal finite domain; at each time point, the finite element dynamic equation is solved based on the time integration algorithm to obtain the displacement, velocity, and acceleration responses of all nodes in the internal finite domain; the displacement, velocity, and acceleration responses of all nodes at all time points are sorted in chronological order to form the seismic dynamic response results.

[0061] In this embodiment, a ground motion time history is preset, which includes multiple discrete time points within the duration of the earthquake action. The equivalent nodal force load sequence is associated with these discrete time points. During the numerical simulation, the equivalent nodal force loads at the corresponding time points are extracted and applied to the boundary nodes of the internal finite domain one by one according to the order of the discrete time points. An event-driven mechanism can also be used, which triggers the corresponding equivalent nodal force load application operation when the simulation time step reaches a certain preset time point, thereby ensuring the accuracy and timing of the load application.

[0062] Among them, the finite element dynamic equation is the structural dynamic equilibrium equation obtained by discretizing the internal finite domain using finite element methods. Its specific form is well known to those skilled in the art and will not be elaborated here. The role of the time integration algorithm is to discretize the time domain, transforming the continuous dynamic equation into a series of algebraic equations for solution, thereby obtaining the displacement, velocity, and acceleration of the structure at different time points. For example, the time integration algorithm can use the Newmark-β method, selecting algorithm parameters β and γ according to preset accuracy and stability requirements, such as taking β=0.25 and γ=0.5, and calculating the displacement, velocity, and acceleration step by step at each time step. The time integration algorithm can also use the central difference method; the time integration algorithm can also use the Wilson-θ method.

[0063] In this embodiment, during the calculation process, the displacement, velocity, and acceleration response data of all nodes in the internal finite domain obtained at each time point can be stored in a pre-allocated data structure, such as a multidimensional array, list, or matrix. These data are sorted and concatenated according to their corresponding time indices to obtain the response results. Alternatively, a database management system can be used to store the node response data at each time point as records and index and sort them using timestamps to facilitate subsequent data queries.

[0064] In this embodiment, the equivalent nodal force load sequence is applied one by one to the boundary nodes corresponding to the internal finite domain based on the time points of the seismic motion time history. This sequential application of loads according to the time nodes not only matches the ordered sequence of equivalent nodal force loads but also conforms to the actual process of seismic motion input, avoiding problems such as load mismatch and chaotic application. At each time point, the finite element dynamic equation is solved based on the time integration algorithm to obtain the displacement, velocity, and acceleration responses of all nodes in the internal finite domain. The displacement, velocity, and acceleration responses of all nodes at all time points are sorted in chronological order to form the seismic dynamic response results, which facilitates the subsequent seismic design and safety evaluation of underground structures and solves the problem that scattered data is difficult to apply directly.

[0065] In this embodiment, the finite element method does not need to be iteratively coupled or synchronously solved with the boundary element method. Instead, it completes the internal nonlinear analysis independently. In this process, the wave field calculation of the outer infinite domain and the nonlinear dynamic analysis of the inner finite domain are decoupled in time. That is, the boundary wave field data on the artificial boundary is first calculated by the fast multipole boundary element method, and then the boundary wave field data is transformed into equivalent nodal forces by the artificial boundary substructure method. Finally, the finite element method independently completes the time history response solution under the action of the equivalent nodal forces. The whole process avoids the problem of repeated interactive iterations required in the traditional finite element-boundary element coupled method, and greatly improves the computational efficiency.

[0066] Example 2 like Figure 2 As shown, this invention also proposes a wave field calculation system coupled with finite element method and artificial boundary, using a wave field calculation method coupled with finite element method and artificial boundary as described in any of Embodiment 1, including the following modules: The acquisition module is used to acquire the computational domain and the seismic incident wave field; A partitioning module, connected to the acquisition module, is used to partition the computational domain into an outer infinite domain and an inner finite domain based on an artificial boundary, wherein the artificial boundary is set outside the nonlinear influence range of the structure-surrounding rock system; The solution module, connected to the partitioning module, is used to solve the external infinite domain based on the fast multipole boundary element method and the seismic incident wave field to obtain the boundary wave field data on the artificial boundary. A conversion module, connected to the solution module, is used to convert the boundary wave field data into equivalent nodal forces acting on the boundary nodes of the internal finite domain based on the artificial boundary substructure method. The analysis module, connected to the partitioning module and the conversion module, is used to apply the equivalent nodal force as a seismic input load to the boundary nodes of the internal finite domain, and to perform nonlinear dynamic response analysis on the internal finite domain to which the equivalent nodal force is applied using the finite element method, thereby obtaining the seismic dynamic response results of the internal finite domain.

[0067] To address the problems of spurious reflections easily generated by artificial boundaries in the finite element method (FEM), low computational efficiency, difficulty in handling material nonlinearity, and coupling difficulties in existing technologies, this invention divides the computational domain into an outer infinite domain and an inner finite domain. The outer infinite domain is solved using the fast multipole boundary element method to obtain artificial boundary wavefield data. The wavefield data is then transformed into equivalent nodal forces on the boundary nodes of the inner finite domain using the artificial boundary substructure method. The FEM is then applied to perform nonlinear dynamic response analysis on the inner finite domain under the applied equivalent nodal forces. This overcomes the shortcomings of low computational efficiency and difficulty in handling nonlinearity, while retaining the FEM's ability to handle complex nonlinearities. It achieves a balance between high accuracy and high computational efficiency, improving the stability and reliability of seismic dynamic response analysis.

[0068] 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 the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the technical solutions of the embodiments of the present invention.

Claims

1. A wave field calculation method coupled with finite element method and artificial boundary, characterized in that, include: Obtain the computational domain and the seismic incident wave field; The computational domain is divided into an outer infinite domain and an inner finite domain based on an artificial boundary, wherein the artificial boundary is set outside the nonlinear influence range of the structure-surrounding rock system. The boundary wave field data on the artificial boundary is obtained by solving the external infinite domain based on the fast multipole boundary element method and the seismic incident wave field. The boundary wave field data is transformed into equivalent nodal forces acting on the boundary nodes of the internal finite domain based on the artificial boundary substructure method. The equivalent nodal force is applied as a seismic input load to the boundary nodes of the internal finite domain. The finite element method is used to perform nonlinear dynamic response analysis on the internal finite domain to which the equivalent nodal force is applied, and the seismic dynamic response results of the internal finite domain are obtained. The determination of the artificial boundary includes: Calculate the extent of the plastic zone of the surrounding rock under seismic loading based on the geometric dimensions, burial depth, and mechanical properties of the underground structure. Obtain the safety distance outside the plastic zone of the surrounding rock; The artificial boundary is determined based on the safety distance, so that the medium on the artificial boundary satisfies the linear elasticity assumption. After determining the artificial boundary based on the safety distance, the method further includes: The region within the artificial boundary is discretized using finite element methods to obtain a finite element mesh, wherein the finite element mesh is used for nonlinear dynamic response analysis. The artificial boundary is discretized based on the boundary element model to obtain a boundary element mesh. Based on the earthquake incident wave field and the finite element method, the finite element dynamic equations corresponding to the finite element mesh are solved to obtain the first nodal displacement and the first nodal surface force of each node on the artificial boundary. Based on the earthquake incident wave field and the fast multipole boundary element method, the boundary integral equation corresponding to the boundary element mesh is solved to obtain the second nodal displacement and the second nodal surface force at the corresponding position on the artificial boundary. Calculate the displacement difference between the displacement of the first node and the displacement of the second node; Calculate the vector sum of the forces on the first and second nodal surfaces; If both the displacement difference and the vector sum satisfy the corresponding preset thresholds, it is determined that the propagation of the seismic incident wave field at the artificial boundary conforms to physical laws.

2. The wave field calculation method coupled with finite element method and artificial boundary as described in claim 1, characterized in that, If at least one of the displacement difference and the vector sum does not satisfy a corresponding preset threshold, the method further includes: The position offset and offset direction of the artificial boundary are determined based on the distribution characteristics of the displacement difference on the artificial boundary, and the position of the artificial boundary is adjusted based on the position offset and offset direction. The mesh refinement region is determined based on the vector and the coordinates of the extreme points on the artificial boundary, and the finite element mesh and boundary element mesh of the mesh refinement region are refined and densified. Based on the adjusted artificial boundary and the refined and densified finite element mesh and boundary element mesh, the displacement difference and vector sum are recalculated until the recalculated displacement difference and the vector sum both meet the corresponding preset threshold.

3. The wave field calculation method coupled with finite element method and artificial boundary as described in claim 1, characterized in that, The solution to the external infinite domain based on the fast multipole boundary element method and the seismic incident wave field includes: Boundary integral equations for the external infinite domain are established based on the boundary element grid of artificial boundaries and the seismic incident wave field. The kernel function in the boundary integral equation is expanded and recursively calculated using a fast multipole expansion algorithm to obtain the displacement vector and surface force vector of each node on the artificial boundary. The displacement vector and surface force vector are used as the boundary wave field data.

4. A wave field calculation method coupled with finite element method and artificial boundary as described in claim 1 or 3, characterized in that, The method based on artificial boundary substructure transforms the boundary wavefield data into equivalent nodal forces acting on the boundary nodes of the internal finite domain, including: Obtain the displacement vector and surface force vector of the artificial boundary nodes in the boundary wave field data; Construct a transformation matrix from boundary wave field data to equivalent nodal forces based on the stiffness matrix and mass matrix in the artificial boundary substructure method; Substituting the displacement vector and surface force vector into the transformation matrix, the equivalent nodal force vector on each boundary node is obtained. The equivalent node force vectors are sorted based on the node numbers to obtain the equivalent node force load sequence.

5. The wave field calculation method coupled with finite element method and artificial boundary as described in claim 4, characterized in that, The equivalent nodal force is applied as a seismic input load to the boundary nodes of the internal finite domain. The finite element method is then used to perform nonlinear dynamic response analysis on the internal finite domain to obtain the seismic dynamic response results of the internal finite domain, including: Based on the time points of the ground motion time history, the equivalent nodal force load sequence is applied one by one to the boundary nodes corresponding to the internal finite domain. At each time point, the finite element dynamic equations are solved based on the time integration algorithm to obtain the displacement, velocity, and acceleration responses of all nodes in the internal finite domain; The displacement, velocity, and acceleration responses of all nodes at all time points are sorted in chronological order to form the earthquake dynamic response results.

6. A wave field calculation system coupled with finite element method and artificial boundary, characterized in that, The system employs a wave field calculation method coupled with finite element and artificial boundary as described in any one of claims 1 to 5, specifically including the following modules: The acquisition module is used to acquire the computational domain and the seismic incident wave field; A partitioning module, connected to the acquisition module, is used to partition the computational domain into an outer infinite domain and an inner finite domain based on an artificial boundary, wherein the artificial boundary is set outside the nonlinear influence range of the structure-surrounding rock system; The solution module, connected to the partitioning module, is used to solve the external infinite domain based on the fast multipole boundary element method and the seismic incident wave field to obtain the boundary wave field data on the artificial boundary. A conversion module, connected to the solution module, is used to convert the boundary wave field data into equivalent nodal forces acting on the boundary nodes of the internal finite domain based on the artificial boundary substructure method. The analysis module, connected to the partitioning module and the conversion module, is used to apply the equivalent nodal force as a seismic input load to the boundary nodes of the internal finite domain, and to perform nonlinear dynamic response analysis on the internal finite domain to which the equivalent nodal force is applied using the finite element method, thereby obtaining the seismic dynamic response results of the internal finite domain.

Images