Soil-rock aggregate deformation-seepage coupling simulation method

Abstract

The invention discloses an earth-rock aggregate deformation-seepage coupling simulation method, and relates to the technical field of coupling simulation, and the method comprises the steps: building a polygon matched with the shape of an earth-rock aggregate according to the actual geometric characteristics of the earth-rock aggregate, taking the polygon information as a model boundary, and inputting corresponding material parameters and control parameters; establishing a numerical manifold porous medium model; applying boundary conditions to the numerical manifold porous medium model by adopting a modified Nitsche method, and establishing a Nitsche type symmetric weak form of a hydraulic coupling problem; searching the cutting unit on the boundary and the interface, and performing ghost penalty stabilization processing on the cutting unit; and calculating the displacement and pore water pressure of each manifold unit in the numerical manifold porous medium model subjected to ghost penalty stabilization treatment under the hydraulic coupling effect. According to the method, the hydraulic coupling solution of the complex area can be carried out on the unmatched grids.

Description

Technical Field

[0001] The present invention relates to the technical field of coupled simulation, and more particularly to a method for simulating the deformation-seepage coupling of soil-rock mixtures. Background Art

[0002] So far, many numerical simulation methods have been proposed to solve the hydro-mechanical coupling problem, such as the finite element method, the boundary element method (BEM), the meshless methods (MMs), etc. To better handle the hydro-mechanical coupling problem of porous media in complex regions, several mesh-incompatible methods have been proposed by different scholars based on the finite element method, among which the representative ones are the numerical manifold method (NMM) and the generalized / extended finite element method (GFEM / XFEM). In addition, there are some other mesh-incompatible methods such as the cut finite element method (Cut-FEM), the finite cell element method (FCM), and the isogeometric analysis.

[0003] Currently, the numerical calculation methods for the hydro-mechanical coupling problem are diverse and complex. When using the traditional finite element method for calculation, it is required that the mesh adapts to various boundaries and interfaces, facing huge challenges in terms of preprocessing efficiency and stability. Therefore, some scholars have proposed some mesh-independent finite element methods. As a representative of the mesh-incompatible methods, the numerical manifold method has shown great advantages in solving the hydro-mechanical coupling problem involving complex geometric models and has achieved some results. However, nowadays, the numerical manifold method still faces some common problems of mesh-based methods, such as the imposition of boundary / interface conditions, the ill-conditioning of the coefficient matrix caused by mesh cutting, etc. The existing research pays relatively little attention to these problems, restricting the application of mesh-incompatible methods in solving the hydro-mechanical coupling problem in complex regions.

[0004] Therefore, how to solve the hydro-mechanical coupling in complex regions on incompatible meshes is an urgent problem to be solved by those skilled in the art. Summary of the Invention

[0005] In view of this, the present invention provides a method for simulating the deformation-seepage coupling of soil-rock mixtures, which can realize the hydro-mechanical coupling solution in complex regions on incompatible meshes.

[0006] To achieve the above object, the present invention adopts the following technical solutions:

[0007] A method for simulating the deformation-seepage coupling of soil-rock mixtures, comprising:

[0008] Establish a polygon that matches the actual geometric characteristics of the soil-rock mixture, use the polygon information as the model boundary, and input the corresponding material parameters and control parameters to establish a numerical manifold porous medium model;

[0009] Apply boundary conditions to the numerical manifold porous medium model using the modified Nitsche method to establish the Nitsche-type symmetric weak form of the hydro-mechanical coupling problem;

[0010] Retrieve the cut elements on the boundary and at the interface, and perform ghost penalty stabilization on the cut elements;

[0011] Calculate the displacements and pore water pressures of each manifold element in the numerical manifold porous medium model after ghost penalty stabilization under the action of hydro-mechanical coupling.

[0012] Preferably, applying boundary conditions to the numerical manifold porous medium model using the modified Nitsche method to establish the Nitsche-type symmetric weak form of the hydro-mechanical coupling problem includes:

[0013] The Nitsche-type weak form of the equilibrium equation at the material interface, find the solution (u, p) ∈ U u ×U p , such that for any v ∈ V u , the following equation holds:

[0014] a u (v, u) + b(v, p) = l u (v);

[0015]

[0016] The Nitsche-type weak form of the mass conservation equation, find the solution (u, p) ∈ U u ×U p , such that for any q ∈ V p , the following equation holds:

[0017] -b 1 (u, q) + a M (q, p) + a p (p, q) = l p (q);

[0018]

[0019] Where V u , U u are the test function and trial function spaces of displacement, v p , U pis the test function and trial function space for pore water pressure, α is the Biot coefficient, κ is the hydraulic conductivity tensor, v is the trial function related to displacement, is the divergence of the displacement test function, represents the first-order derivative of the displacement vector u with respect to time t, {q} represents the average value of the pore water pressure test function, M represents the Biot modulus, represents the first-order derivative of the pore water pressure p with respect to time t, represents the gradient of the pore water pressure, represents the gradient of the pore water pressure test function, Γ f represents the Neumann boundary that prescribes the normal flux, Γ p represents the Dirichlet boundary that prescribes the pore pressure, is defined on the Dirichlet boundary Γ p is the pore pressure on it, u is the solid skeleton displacement vector, p is the pore water pressure, q is the trial function related to the pore water pressure, ε(u) is the strain tensor with respect to u, σ′(v) and σ′(u) are the effective stress tensors with respect to v and u respectively, Ω represents the entire problem domain, Γ u represents the displacement boundary, Γ t is the force boundary condition, Γ I is the material interface, n is the unit outer normal vector on the boundary, b is the body force vector, is defined on the boundary Γ t is the known surface force vector on it, is defined on the boundary Γ u is the known displacement vector on it, and are the jump terms with respect to v and u on the interface respectively, {σ′(u)} and {σ′(v)} are the effective stresses on both sides of the interface, {p} is the average term of the water pressure, h is the characteristic length of the element, and are the penalty parameters.

[0020] Preferably, the retrieval of the cut elements on the boundary and the interface and the ghost penalty stabilization treatment of the cut elements specifically include:

[0021] At the boundary, three ghost penalty stabilization terms are added to the inner edges of each cut element; at the interface, three ghost penalty stabilization terms are added to the inner edges of each virtual element respectively as follows:

[0022]

[0023] Among them, represents the j-th normal partial derivative of the displacement test function v i and is defined as Denotes the j-th normal partial derivative of the pore water pressure function p, defined as Denotes the j-th normal partial derivative of the pore water pressure test function q, defined as h F Is the size of the edge F, Denotes the (2j - 1)-th power of the size of the edge F; Denotes the (2j - 1)-th power of the size of the edge F; Denotes the (2j + 1)-th power of the size of the edge F; Is the user-defined penalty parameter. The virtual element is formed by cutting an element through the interface, and after decomposition of the cut element, two elements located on both sides of the interface are formed; F Γ Is the set of internal edges of the cut element, defined as: The virtual element is formed by cutting an element through the interface, and after decomposition of the cut element, two elements located on both sides of the interface are formed; F Γ Is the set of internal edges of the cut element, defined as:

[0024] F Γ ={F: F = E 1 ∩E 2 , E 1 ≠E 2 , E 1 ∈T Γ or E 2 ∈T Γ};

[0025] Where d is the spatial dimension, p u Is the order of displacement interpolation, Is the given penalty parameter, h F Is the size of the edge F, E1, E2 represent two different elements, T Γ Represents the boundary cut element;

[0026] Denotes the j-th normal partial derivative of the displacement component u i , defined as:

[0027]

[0028] The square brackets [·] denote the difference across the edge F:

[0029] [u i = u i |F l - u i |F r ;

[0030] After introducing the ghost penalty stabilization term, the weak form of the system becomes: Find the solution For any The following holds:

[0031]

[0032] Wherein:

[0033]

[0034] v h is the approximate function of the displacement test function, is the j-th power of the unit normal vector, F l represents the edge of the left virtual element, F r represents the edge of the right virtual element, u h is the displacement function, p h is the pore water pressure function, q h is the pore water pressure test function, and are the displacement test function and the trial function space, and are the pore water pressure test function and the trial function space; i s (v h , u h )、 and i p (q h , p h ) are the ghost penalty stabilization terms, and the superscript h is the approximate function.

[0035] Preferably, after the ghost penalty stabilization treatment of the cutting element, it further includes time discretization of the weak form introducing the ghost penalty stabilization term, and its formula is:

[0036]

[0037] Wherein, X represents u and p, Δt is the time step, u is the solid skeleton displacement vector, and p is the pore water pressure;

[0038] After time discretization, the fully discrete weak form is obtained: for n≥1, given Find the solution such that for any the following equation holds:

[0039]

[0040] Wherein, n and n - 1 respectively represent the current time step and the previous time step, and the t in n represents the corresponding value at the n-th time step, p h,n-1 、p h,n respectively represent the finite element approximate solutions of the pore water pressure at the previous time step and the current time step, Denote the prescribed displacement on the Dirichlet boundary Γ u at the n-th time step.

[0041] Preferably, the calculation specifically includes the displacements and pore water pressures of each manifold element in the numerical manifold porous media model under the action of hydro-mechanical coupling after the ghost penalty stabilization treatment:

[0042] Global approximation of displacement:

[0043] [u h (x,t)] = [u x u y T = N u (x)d u (t);

[0044] Global approximation of water pressure:

[0045] p h (x,t) = N p (x)d p (t);

[0046] The matrix form expression of strain is:

[0047] [ε h (x,t)] = [ε xx (x,t) ε yy (x,t) γ xy (x,t)] T = B u (x)d u (t);

[0048] The matrix form expression of pore water pressure gradient is:

[0049]

[0050] The expression of volumetric strain is:

[0051]

[0052] The calculation expressions of displacement and pore water pressure are:

[0053] Kd = F;

[0054] where N u (x) and N p (x) are the shape function matrices of displacement and pore water pressure respectively, u h (x,t) is the finite element approximation function of displacement, p h (x,t) is the finite element approximation function of pore water pressure,​ is the pore water pressure gradient, are its components in the x - direction and y - direction respectively, d u (t) is the nodal vector of displacement, d p (t) is the nodal vector of pore water pressure, ε h (x, t) is the finite - element approximation function of strain, and its components are expressed as ε xx (x, t), ε yy (x, t), γ xy (x, t), B u (x) is the strain matrix, B p (x) is the pore water pressure gradient matrix, B vol is the volumetric strain matrix, K is the system stiffness matrix, d is the nodal displacement vector, and F is the nodal force vector.

[0055] Preferably, the material parameters include the elastic modulus, Poisson's ratio, permeability, seepage flow rate, permeability coefficient, cross - sectional area, water pressure difference between upstream and downstream, seepage length, hydraulic gradient, fluid viscosity, pore pressure gradient, Biot modulus, and Biot - Willis constant of each medium.

[0056] Through the above - mentioned technical solutions, compared with the prior art, the present invention discloses a method for simulating the deformation - seepage coupling of soil - rock mixtures. The Nitsche method is applied to and modified for the equilibrium equation and the mass conservation equation respectively. By introducing additional terms at the boundaries and interfaces, the accurate imposition of the essential boundary conditions of displacement and water pressure at the boundaries and interfaces is achieved. The introduction of three ghost - penalty stabilization terms reduces the ill - conditioning problem of the coefficient matrix caused by small cuts. BRIEF DESCRIPTION OF THE DRAWINGS

[0057] In order to more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the following will briefly introduce the drawings required for the description of the embodiments or the prior art. Obviously, the drawings in the following description are only the embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on the provided drawings.

[0058] Figure 1 is the problem domain diagram with an interface provided by the present invention;

[0059] Figure 2(a) is the physical region and mesh diagram with a material interface provided by the present invention;

[0060] Figure 2(b) is the cut - element and internal - edge diagram provided by the present invention;

[0061] Figure 3 is the flow - chart of the method steps provided by the present invention;

[0062] Figure 4 The finite element interpolation diagram provided by the present invention. Specific implementation manners

[0063] Next, the technical solutions in the embodiments of the present invention will be clearly and completely described in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without making creative efforts shall fall within the protection scope of the present invention.

[0064] The embodiments of the present invention disclose a method for simulating the deformation-seepage coupling of soil-rock mixtures, as Figure 3 shown, including:

[0065] Establish a polygon that matches the actual geometric characteristics of the soil-rock mixture, use the polygon information as the model boundary, and input the corresponding material parameters and control parameters to establish a numerical manifold porous medium model;

[0066] Apply boundary conditions to the numerical manifold porous medium model by using the modified Nitsche method to establish the Nitsche-type symmetric weak form of the hydro-mechanical coupling problem;

[0067] Retrieve the cut elements on the boundary and the interface, and perform ghost penalty stabilization processing on the cut elements;

[0068] Calculate the displacements and pore water pressures of each manifold element in the numerical manifold porous medium model after the ghost penalty stabilization processing under the hydro-mechanical coupling action.

[0069] When the ratio of the actual area of the manifold element to the area of the mathematical element (referred to as the cut ratio, denoted by η) is very small, this manifold element is called a cut element. For different problems, η has different values and needs to be determined according to the actual situation.

[0070] Among them, a polygon that matches the actual geometric characteristics of the soil-rock mixture is established, the polygon information is imported into the MATLAB system as the model boundary, and it is embedded in a polygon or polyhedron mesh, or other partially overlapping closed regions of any shape. Thus, a weakly discontinuous porous medium model composed of a finite number of manifold elements that do not overlap with each other and whose union constitutes the entire problem domain is established and the corresponding material parameters and control parameters are input.

[0071] In a specific embodiment, the material parameters include the elastic modulus, Poisson's ratio, permeability, seepage flow rate, permeability coefficient, cross-sectional area, upstream and downstream water pressure difference, seepage length, hydraulic gradient, fluid viscosity, pore pressure gradient, Biot modulus, and Biot-Willis constant of each medium.

[0072] In a specific embodiment, for the problem of the material interface, to obtain the weak form of its governing equation, for each subdomain Ω i (i = 1, 2), the following function spaces are defined:

[0073]

[0074] where d is the spatial dimension, and respectively represent the test function and trial function spaces of the displacement, and respectively represent the test function and trial function spaces of the pore pressure, Ω i represents the sub-region, Ω i ×[0, T] represents the Cartesian product of the sub-region and the time interval [0, T], H 1 (Ω i ) is the first-order Sobolev space, u is the solid skeleton displacement vector, and p is the pore water pressure.

[0075] For the problem of the material interface, to obtain the weak form of its governing equation, for each subdomain Ω i (i = 1, 2), the following function spaces are redefined:

[0076]

[0077] where d is the spatial dimension, and respectively represent the test function and trial function spaces of the displacement, and respectively represent the test function and trial function spaces of the pore pressure. The global function space is defined as:

[0078]

[0079] Define the weighted average of the function v at the interface:

[0080] {v} = c 1 v 1 + c 2 v 2 (1.3)

[0081] where c 1 and c 2 are the weights, and c 1 + c 2 = 1. Selecting appropriate weight values is beneficial to the stability and convergence of the algorithm, and can avoid the ill-conditioning of the unit stiffness matrix and the global stiffness matrix caused by large parameter differences.

[0082] Regarding the elasticity problem, the weighting coefficient values are reasonably designed, and the expression is as follows:

[0083]

[0084] Among them,

[0085] η i =λ i +2G i (1.5)

[0086] λ i is the Poisson's ratio of the material, and G i represents the shear modulus.

[0087] In a specific embodiment, retrieving the cutting elements at the boundaries and interfaces and performing ghost penalty stabilization on the cutting elements specifically includes:

[0088] At the boundaries, three ghost penalty stabilization terms are added to the inner edges of each cutting element; at the interfaces, three ghost penalty stabilization terms are respectively added to the inner edges of each virtual element as follows:

[0089]

[0090] Among them, the virtual element is formed by cutting through the interface, and after the cut element is decomposed, two elements located on both sides of the interface are formed, as shown in Figure 2(b); F Γ is the set of inner edges of the cut element, defined as:

[0091] F Γ ={F:F = E 1 ∩E 2 , E 1 ≠E 2 , E 1 ∈T Γ or E 2 ∈T Γ};

[0092] d is the spatial dimension, p u is the order of displacement interpolation, is the given penalty parameter, h F is the size of the edge F;

[0093] represents the j-th normal derivative of the displacement component u i , defined as:

[0094]

[0095] The square brackets [·] represent the difference between both sides of the edge F:

[0096] [u i =ui |F l -u i |F r ;

[0097] After introducing the ghost penalty stabilization term, the weak form of the system becomes: Find the solution For any the following holds:

[0098]

[0099] where:

[0100]

[0101] v h is the approximate function of the displacement test function, u h is the displacement function, p h is the pore water pressure function, q h is the pore water pressure test function, and are the displacement test function and trial function spaces, and are the pore water pressure test function and trial function spaces.

[0102] To facilitate the presentation of the process of the ghost penalty stabilization scheme, consider a two-dimensional porous medium problem domain Ω with a material interface as shown Figure 1 in which it is divided into two different sub-domains Ω I and Ω 1 and Ω 2 .

[0103] The boundary conditions for the coupled solid and liquid phases are given independently. For the solid phase, the boundary conditions include the prescribed displacement defined on the Dirichlet boundary Γ u and the external force load defined on the Neumann boundary Γ t , and Γ u , Γ t satisfy Γ u ∪Γ t = Γ,

[0104]

[0105] For the fluid, the boundary conditions include the pore pressure defined on the Dirichlet boundary Γ p and the normal flux defined on the Neumann boundary Γ f , and similarly Γ p ∪Γ f = Γ, These boundary conditions are expressed as:

[0106]

[0107] The initial conditions are the displacements and pore pressures specified at the initial time \(t = 0\), expressed as:

[0108]

[0109] where \(u(x, 0)\) is the displacement function, \(u 0 (x)\) is its value at \(t = 0\), \(p(x, 0)\) is the pore water pressure function, and \(p 0 (x)\) is its value at \(t = 0\).

[0110] For any material interface \(\Gamma I \), the displacements and water pressures on both sides are continuous, and the stresses and fluxes are in balance. The interface conditions are expressed as:

[0111]

[0112] where \(n\) represents the unit normal vector pointing from \(\Omega 1 \) to \(\Omega 2 \), and \( \) represents the jump at the interface, defined by the following equation:

[0113]

[0114] where \(V u \) and \(U u \) represent the trial function space and test function space of the displacement respectively, \(V p \) and \(U p \) represent the trial function space and test function space of the pore pressure respectively, \(Q l \) and \(Q m \) represent the spaces of the \(l\)-th and \(m\)-th order polynomials respectively, \(E\) represents the element after the discretization of the problem domain space, and \(T h \) represents the set of elements.

[0115] The equilibrium equation of the fully coupled Biot consolidation model is:

[0116]

[0117] The mass conservation equation is:

[0118]

[0119] Take the trial function \(v\in V u \), multiply both sides of Equation (1.12) by \(v\) and integrate over the problem domain. The weak form equations established include:

[0120] Establish the Nitsche-type symmetric weak form of the equilibrium equations for the material interface and find the solution \((u, p)\in U\) u \(\times U\) p such that for any \(v\in V\) u the following holds:

[0121] \(a\) u (v, u)+b(v, p)=l u (v); (1.13)

[0122]

[0123] Establish the Nitsche-type symmetric weak form of the mass conservation equations and find the solution \((u, p)\in U\) u \(\times U\) p such that for any \(q\in V\) p the following holds:

[0124] -b 1 (u, q)+a M (q, p)+a p (p, q)=l p (q); (1.15)

[0125]

[0126] where \(V\) u and \(U\) u are the test function and trial function spaces for displacements, \(v\) p and \(U\) p are the test function and trial function spaces for pore water pressures, \(\alpha\) is the Biot coefficient, \(\kappa\) is the hydraulic conductivity tensor, \(v\) is the trial function related to the displacement, is the divergence of the displacement test function, denotes the first-order derivative of the displacement vector \(u\) with respect to time \(t\), \(\{q\}\) denotes the average value of the pore water pressure test function, \(M\) denotes the Biot modulus, denotes the first-order derivative of the pore water pressure \(p\) with respect to time \(t\), denotes the gradient of the pore water pressure, denotes the gradient of the pore water pressure test function, \(\Gamma\) f denotes the Neumann boundary where the normal flux is prescribed, \(\Gamma\) p denotes the Dirichlet boundary where the pore pressure is prescribed, is the pore pressure defined on the Dirichlet boundary \(\Gamma\) p , \(u\) is the solid skeleton displacement vector, \(p\) is the pore water pressure, \(q\) is the trial function related to the pore water pressure, \(\varepsilon(u)\) is the strain tensor with respect to \(u\), \(\sigma'(v)\) and \(\sigma'(u)\) are the effective stress tensors with respect to \(v\) and \(u\) respectively, \(\Omega\) denotes the entire problem domain, \(\Gamma\)u Denotes the displacement boundary, Γ t Is the force boundary condition, Γ I Is the material interface, n is the unit outer normal vector on the boundary, b is the body force vector, Is defined on the boundary Γ t Of the known surface force vector, Is defined on the boundary Γ u Of the known displacement vector, And Are the jump terms of v and u on the interface respectively, {σ′(u)} and {σ′(v)} are the effective stresses on both sides of the interface, {p} is the average term of the water pressure, h is the characteristic length of the element, And Are the penalty parameters.

[0127] Perform spatial discretization on the above weak form equation:

[0128] In the case where the mesh completely covers the physical domain, standard finite element interpolation is used for displacement and water pressure. Specifically, the 9-node quadrilateral element (Q9Q4 shown in the following figure, 9-node element interpolation for displacement and 4-node element interpolation for water pressure) is used. Among them, the black dots represent the nodes for displacement interpolation, and the circles represent the nodes for water pressure interpolation, as Figure 4 Shown.

[0129] Construct a biquadratic interpolation for displacement using all nodes. For any node P i (x i ,y i ), the corresponding shape function is:

[0130]

[0131] When constructing the water pressure interpolation, only the corner nodes 1, 2, 3, and 4 are used, and the corresponding shape functions are:

[0132]

[0133] Thus, the finite element weak form equation of spatial discretization can be obtained:

[0134] The discrete equation of the Nitsche type weak form of the equilibrium equation containing the material interface is: Find the solution (u h ,p h )∈U u ×U p , such that for any v∈V u , the following equation holds:

[0135] a u (v h ,u h )+b(v h,p h )=l u (v h )

[0136] in:

[0137]

[0138] The Nitsche-type weak form of the mass conservation equation is discrete: Find the solution (u, p)∈U u ×U p , so that for any q∈V p , the following formula holds:

[0139] -b 1 (u h ,q h )+a M (q h ,p h )+a p (p h ,q h )=l p (q h )

[0140] in:

[0141]

[0142] Among them, v h is the approximate function of the displacement test function, u h is the displacement function, p h is the pore water pressure function, q h is the pore water pressure test function, and is the displacement test function and the test function space and It is the pore water pressure test function and the test function space.

[0143] In a specific embodiment, in order to solve the problem of cutting units, it is also necessary to apply the ghost penalty method to the weak form. In order to better describe the ghost penalty method, the following notation is first introduced:

[0144] Consider the problem domain Ω surrounded by a quasi-uniform grid T h Full coverage, T h = {E}, E is a regular quadrilateral unit, and The boundary cutting unit is denoted as T Γ , The inner edge set of the cutting unit is F Γ :

[0145] F Γ ={F:F=E1 ∩E 2 ,E 1 ≠E 2 ,E 1 ∈T Γ or E 2 ∈T Γ} (1.17)

[0146] First, for the solid stiffness, introduce the first stabilization term i s (u, v), and add it to the weak form of the equilibrium equation, and the specific expression is as follows:

[0147]

[0148] where d is the spatial dimension, corresponding to 2 and 3 for two-dimensional and three-dimensional problems respectively, p u is the order of displacement interpolation, is the given penalty parameter, h F is the size of the edge F, represents the j-th normal derivative of the displacement component u i and is defined as:

[0149]

[0150] The square brackets [·] represent the difference on both sides of the edge F:

[0151]

[0152] When the solution domain is cut by material interfaces (such as Γ I ) in Fig. 2(a) or discontinuous surfaces such as internal cracks, the region is separated by the material interface. At this time, each element cut by the interface is decomposed into two virtual elements on both sides of the interface. For example, the green-filled elements in Fig. 2(b), and their inner edges are marked in red. The specific implementation scheme is to add three ghost penalty stabilization terms i s (u, v), i p (q, p) and

[0153] respectively on the inner edges of the two-sided elements. By introducing stabilization terms at these positions, the numerical instability caused by interface cutting can be effectively addressed. In addition, no additional stabilization terms are required on the outer boundaries of the grids on each side after cutting. After introducing the ghost penalty terms, the weak form of the system becomes: find the solution For any the following holds:

[0154]

[0155] where v h is the approximate function of the displacement test function, uh is the displacement function, p h is the pore water pressure function, q h is the pore water pressure test function, and are the displacement test function and the trial function space and are the pore water pressure test function and the trial function space.

[0156] Among them, represents the j-th normal partial derivative of the displacement test function v i and is defined as represents the j-th normal partial derivative of the pore water pressure function p and is defined as represents the j-th normal partial derivative of the pore water pressure test function q and is defined as h F is the size of the edge F, represents the (2j - 1)-th power of the size of the edge F; represents the (2j - 1)-th power of the size of the edge F; represents the (2j + 1)-th power of the size of the edge F; is the user-defined penalty parameter. The virtual element is formed by cutting the interface. After the cut element is decomposed, two elements located on both sides of the interface are formed, as shown in Fig. 2(b); F Γ is the set of inner edges of the cut element and is defined as:

[0157] F Γ ={F: F = E 1 ∩E 2 ,E 1 ≠E 2 ,E 1 ∈T Γ or E 2 ∈T Γ}

[0158] F represents the inner edge of the cut element, E i represents the element (i = 1, 2), d is the space dimension, p u is the order of displacement interpolation, is the given penalty parameter, h F is the size of the edge F, T Γ represents the boundary cut element.

[0159] represents the j-th normal partial derivative of the displacement component u i and is defined as:

[0160]

[0161] The square brackets [·] denote the difference on both sides of edge F:

[0162]

[0163] It is obtained that

[0164] A u (v, u)+b(v, p) = l u (v)

[0165] A u (v, u) = a u (v, u)+i s (v, u)

[0166] Introduce a stabilization term into the left - hand side of the weak - form equation of the mass - conservation equation:

[0167]

[0168] where p p is the order of the interpolation of the water pressure, and are the given penalty parameters, and thus the weak form with the introduced stabilization term can be obtained

[0169]

[0170] v h is the approximate function of the displacement test function, is the j - th power of the unit normal vector, F l denotes the edge of the left - hand - side virtual element, F r denotes the edge of the right - hand - side virtual element, u h is the displacement function, p h is the pore - water pressure function, q h is the pore - water pressure test function, and are the displacement test function and trial - function spaces, and are the pore - water pressure test function and trial - function spaces; i s (v h , u h ) and with i p (q h , p h ) are the ghost - penalty stabilization terms, and the superscript h is the approximate function.

[0171] The introduction of the three stabilization terms is an important improvement in the NMM numerical simulation of the hydro - mechanical coupling problem, which greatly improves the computational stability. However, it must be noted that the ghost - penalty parameters and Their values have a certain impact on numerical calculations. Whether they are too large or too small, they may lead to instability in numerical calculations or difficulties in convergence. Therefore, the optimal value of each ghost penalty parameter should be within a finite range. These parameters are designed to maintain the coercivity and inf-sup conditions of the corresponding terms in the weak form, and their specific values depend on the corresponding materials and other calculation parameters, such as the time step, etc.

[0172] In a specific embodiment, the discretization of the weak form equation with the introduced stabilization term specifically includes: performing time discretization on the weak form formula using the backward Euler method, and its formula is:

[0173]

[0174] where X represents u and p, Δt is the time step, u is the solid skeleton displacement vector, and p is the pore water pressure;

[0175] After performing time discretization, the fully discrete weak form is obtained: for n≥1, given Find the solution such that for any the following equation holds:

[0176]

[0177] where n and n - 1 represent the current time step and the previous time step respectively, and the t in n represents the corresponding value at the nth time step, and p h,n-1 , p h,n represent the finite element approximate solutions of the pore water pressure at the previous time step and the current time step respectively, represents the value of the prescribed displacement on the Dirichlet boundary Γ u at the nth time step.

[0178] In a specific embodiment, the fully discrete weak form is successfully established by adopting the Galerkin method. In order to perform efficient solution on a computer, within the framework of the numerical manifold method, a regular quadrilateral mesh with high quality and uniform distribution is selected as the mathematical cover for accurate numerical solution of the hydro-mechanical coupling problem. In the two-dimensional case, the NMM global approximations of displacement and pore water pressure can be expressed as:

[0179]

[0180] where the displacement shape function matrix N u (x) and the pore water pressure shape function matrix N p (x) are respectively:

[0181]

[0182] Among them, nu and np are the number of interpolation nodes for displacement and water pressure respectively.

[0183] In addition, the expressions for strain and pressure gradient are:

[0184]

[0185] Among them, is the pore water pressure gradient, are its components in the x - direction and y - direction respectively, d p (t) is the nodal vector of pore water pressure, ε h (x,t) is the finite - element approximation function of strain, and its components are expressed as ε xx (x,t), ε yy (x,t), γ xy (x,t), B u (x) is the strain matrix, B p (x) is the pore water pressure gradient matrix,

[0186]

[0187] Volumetric strain The matrix form is:

[0188]

[0189] Among them, is the volumetric strain, B vol is the volumetric strain matrix, d u (t) is the nodal vector of displacement

[0190] B vol = [1 1 0]B u (1.31)

[0191] Furthermore, the following algebraic equation for the fully coupled water - coupled problem is obtained, expressed as:

[0192] Kd = F; (1.32)

[0193] Among them,

[0194]

[0195] In the matrix, the specific expressions of each term are,

[0196]

[0197]

[0198] Among them, K is the system stiffness matrix, K uu, K up , K up,b , K up,I , K pp M pp , is the coefficient matrix, is the nodal displacement degree - of - freedom vector at the nth time step; is the pore - water pressure degree - of - freedom vector at the nth time step; is the value of the force vector acting on the displacement degree of freedom u at the nth time step, is the value of the force vector acting on the pore - water pressure degree of freedom p at the nth time step, is the value of the boundary force vector acting on the pore - water pressure degree of freedom p related, Δt is the time step, Ω is the problem domain, E is the mesh element, represents the discrete form of the Dirichlet boundary, is the discrete form of the interface, is the discrete form of the water - pressure boundary, Γ represents the boundary, e represents the element boundary, F represents the face, T h is the set of all elements, B u is the strain matrix, B p is the pressure - gradient matrix, B vol is the volumetric - strain matrix, D is the elastic matrix, N u is the displacement shape - function matrix, N p is the pore - water pressure shape - function matrix, f h,n is the value of the fluid source f at the nth time step, is defined on the Neumann boundary Γ f the normal flux on it at the nth time step, is the pore pressure defined on the Dirichlet boundary Γ p on it at the nth time step, b h,n is the value of the body force b at the nth time step, is the external - force load defined on the Neumann boundary Γ t on it at the nth time step, is the prescribed displacement defined on the Dirichlet boundary Γ u on it at the nth time step, α is the Biot coefficient, M is the Biot modulus, κ is the permeability tensor, is the user - defined penalty parameter, is the stabilization parameter, and is the ghost penalty parameter, n V and n M is the matrix corresponding to the unit normal vector n, which is used for the calculations between the unit normal vector and vectors as well as second-order tensors, n V and n M is the matrix corresponding to the unit normal vector n, which is used for the calculations between the unit normal vector and vectors as well as second-order tensors, n V is defined as:

[0199] n V =[n x n y (1.35)

[0200] n M is defined as:

[0201]

[0202] Solving the above matrix in MATLAB can obtain the displacement u and pore water pressure p.

[0203] Post-process the finite element calculation results, and draw stress nephograms, displacement vector diagrams, and pore water pressure nephograms through visualization techniques to display the mechanical state of the porous medium inside the hydro-mechanical coupling field. At the same time, according to the optimization objectives and constraints, comprehensively evaluate the optimized structural design scheme, including comparative analysis of indicators such as the maximum stress and maximum displacement.

[0204] Based on the numerical manifold method, the present invention establishes a mismatched finite element method with a mesh not matching the physical domain for a two-field (displacement and pore water pressure as the main variables) poroelastic problem containing weak discontinuities. Compared with the matched finite element method, this method does not need to consider various interfaces of the physical domain. It embeds the physical domain into a relatively simple mesh, greatly simplifying the generation process of the preprocessing. Moreover, when new discontinuous interfaces are generated, the mismatched method can easily handle these new discontinuous interfaces without changing the original mesh structure, greatly simplifying the geometric operations in the simulation process and having greater flexibility and applicability.

[0205] In addition, the mismatched finite element method often has problems such as difficult application of boundary conditions and too large condition number of the coefficient matrix. The present invention uses the Nitsche method to apply boundary conditions, establishes a Nitsche-type symmetric weak form of the hydro-mechanical coupling problem, and then ensures that the bilinear form approximately satisfies the inf-sup condition in cases such as too small permeability or time step by using Taylor-Hood elements for interpolation, and uses the ghost penalty stabilization method to solve the ill-conditioning problem of the coefficient matrix caused by small cuts. Thus, the accuracy and stability of the numerical simulation results obtained by the present invention are ensured.

[0206] The various embodiments in this specification are described in a progressive manner. Each embodiment focuses on the differences from other embodiments. For the same or similar parts among the various embodiments, reference can be made to each other. For the apparatuses disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple. For related parts, reference can be made to the description in the method section.

[0207] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art. The general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to these embodiments shown herein, but rather to the broadest scope consistent with the principles and novel features disclosed herein.

Claims

1. A soil-rock mixture deformation-seepage coupling simulation method, characterized in that: include: According to the actual geometric characteristics of the soil-rock mixture, a polygon matching its shape is established, the polygon information is used as the model boundary, and the corresponding material parameters and control parameters are input to establish a numerical manifold porous media model; Applying boundary conditions to the numerical manifold porous media model using a modified Nitsche method to establish a Nitsche-type symmetric weak form of the hydraulic coupling problem; Retrieve cutting units on the boundary and interface, and perform ghost penalty stabilization on the cutting units; The displacement and pore water pressure of each manifold unit in the numerical manifold porous media model after ghost penalty stabilization are calculated under hydraulic coupling.

2. The soil-rock mixture deformation-seepage coupling simulation method according to claim 1, characterized in that: The modified Nitsche method is used to apply boundary conditions to the numerical manifold porous media model to establish the Nitsche type symmetric weak form of the hydraulic coupling problem, including: The equilibrium equation of the material interface is a Nitsche weak form, and the solution (u,p)∈U is sought. u ×U p , so that for any v∈V u , the following formula holds: a u (v,u)+b(v,p)=l u (v); Nitsche-type weak form of the mass conservation equation, finding solutions (u,p)∈U u ×U p , so that for any q∈V p , the following formula holds: -b1(u,q)+a M (q,p)+a p (p,q)=l p (q); Among them, V u , U u is the displacement test function and trial function space, v p , U p is the test function and trial function space of pore water pressure, α is the Biot coefficient, κ is the hydraulic conductivity tensor, v is the trial function related to displacement, is the divergence of the displacement test function, represents the first-order derivative of the displacement vector u with respect to time t, {q} represents the average value of the pore water pressure test function, M represents the Biot modulus, represents the first-order derivative of pore water pressure p with respect to time t, represents the gradient of pore water pressure, represents the gradient of the pore water pressure test function, Γ f represents the Neumann boundary that specifies the normal flux, Γ p represents the Dirichlet boundary that specifies the pore pressure, is defined on the Dirichlet boundary Γ p , u is the solid skeleton displacement vector, p is the pore water pressure, q is a test function related to the pore water pressure, ε(u) is the strain tensor about u, σ′(v) and σ′(u) are the effective stress tensors about v and u respectively, Ω represents the entire problem domain, Γ u represents the displacement boundary, Γ t is the force boundary condition, Γ I is the material interface, n is the unit external normal vector on the boundary, b is the body force vector, is defined on the boundary Γ t The known surface force vector on is defined on the boundary Γ u The known displacement vector on and are the jump terms about v and u on the interface, {σ′(u)} and {σ′(v)} are the effective stresses on both sides of the interface, {p} is the average term of water pressure, h is the characteristic length of the unit, and is the penalty parameter.

3. The soil-rock mixture deformation-seepage coupling simulation method according to claim 2, characterized in that: The searching of the cutting unit on the boundary and the interface and performing the ghost penalty stabilization process on the cutting unit specifically includes: At the boundary, three ghost penalty stabilization terms are added to the inner edge of each cutting unit; at the interface, three ghost penalty stabilization terms are added to the inner edge of each virtual unit as shown below: in, Represents the displacement test function v i The j-th normal partial derivative of is defined as represents the j-order normal partial derivative of the pore water pressure function p, defined as represents the j-order normal partial derivative of the pore water pressure test function q, defined as h F is the dimension of side F, represents the 2j-1 power of the size of edge F; represents the 2j-1 power of the size of edge F; represents the size of edge F to the power of 2j+1; is a user-defined penalty parameter. The virtual unit is formed by decomposing the cut unit through the interface to form two units on both sides of the interface. Γ is the set of inner edges of the cut unit, defined as: F Γ ={F:F=E1∩E2,E1≠E2,E1∈T Γ or E2∈T Γ }; Among them, d is the spatial dimension, p u is the order of displacement interpolation, is a given penalty parameter, h F is the size of side F, E1, E2 represent two different units, T Γ represents the boundary cutting unit; represents the displacement component u i The j-order normal partial derivative of is defined as: The square brackets [·] represent the difference between the two sides of edge F: [and i ]=and i |F l -and i |F r ; After introducing the ghost penalty stability term, the weak form of the system becomes: Find the solution For any The following holds true: in: v h is the approximate function of the displacement test function, is the jth power of the unit normal vector, F l represents the edge of the left virtual unit, F r represents the edge of the right virtual unit, u h is the displacement function, p h is the pore water pressure function, q h is the pore water pressure test function, and is the displacement test function and the test function space, and is the pore water pressure test function and the test function space; i s (v h ,u h ), with i p (q h ,p h ) is the ghost penalty stabilization term, and the superscript h is the approximate function.

4. The soil-rock mixture deformation-seepage coupling simulation method according to claim 3 is characterized in that: After the cutting unit is subjected to ghost penalty stabilization, the weak form of the ghost penalty stabilization term is also time discretized, and the formula is: Where X represents u and p, Δt is the time step, u is the solid skeleton displacement vector, and p is the pore water pressure; After time discretization, we obtain the fully discrete weak form: for n ≥ 1, Finding a solution So that for any The following holds true: Among them, n and n-1 represent the current time step and the previous time step respectively. and t in n represents the corresponding value at the nth time step, p h,n-1 、p h,n They represent the finite element approximate solutions of pore water pressure in the previous time step and the current time step, respectively. represents the Dirichlet boundary Γ u The value of the prescribed displacement on at the nth time step.

5. The soil-rock mixture deformation-seepage coupling simulation method according to claim 4, characterized in that: The calculation of the displacement and pore water pressure of each manifold unit in the numerical manifold porous medium model after ghost penalty stabilization under hydraulic coupling specifically includes: Global approximation for displacement: [u h (x,t)]=[u x u y ] T =N u (x)d u (t); Global approximation of water pressure: p h (x,t)=N p (x)d p (t); The matrix expression of strain is: [ε h (x,t)]=[ε xx (x,t)ε yy (x,t)γ xy (x,t)] T =B u (x)d u (t); The matrix expression of pore water pressure gradient is: Volumetric strain The expression is: The calculation expressions of displacement and pore water pressure are: Kd = F; Among them, N u (x) and N p (x) are the shape function matrices of displacement and pore water pressure, respectively, u h (x, t) is the finite element approximation function of displacement, p h (x, t) is the finite element approximate function of pore water pressure, is the pore water pressure gradient, are their components in the x and y directions respectively, d u (t) is the node vector of displacement, d p (t) is the node vector of pore water pressure, ε h (x, t) is the finite element approximate function of strain, and its components are represented by ε xx (x,t), ε yy (x,t),γ xy (x,t),B u (x) is the strain matrix, B p (x) is the pore water pressure gradient matrix, B vol is the body strain matrix, K is the system stiffness matrix, d is the node displacement vector, and F is the node force vector.

6. The soil-rock mixture deformation-seepage coupling simulation method according to claim 1, characterized in that: The material parameters include elastic modulus, Poisson's ratio, permeability, seepage rate, permeability coefficient, cross-sectional area, upstream and downstream water pressure difference, seepage length, hydraulic gradient, fluid viscosity, pore pressure gradient, Biot modulus and Biot-Willis constant of each medium.

Images