Saturated ground analysis method and program

Abstract

To provide an analysis method and a program therefor capable of improving accuracy by enhancing resolution related to a spatial distribution of a pressure gradient (or a hydraulic gradient) in finite element analysis of a saturated ground.SOLUTION: In a saturated ground analysis method for discretizing and solving a continuity equation described by a u-p formulation in which a displacement u of a soil skeleton and a pore water pressure p are unknowns by application of a finite volume method, a permeability in an element is discretized into nodes by an equivalent permeability matrix shown in Equation 1 so that an outflow amount of pore water in the element can be expressed using a pressure gradient vector (or a hydraulic gradient vector) at a node represented by the pore water pressure p.SELECTED DRAWING: None

Description

[Technical Field]

[0001] The present invention relates to a method for analyzing saturated ground and a program using the same. [Background technology]

[0002] When it comes to the behavior of saturated ground during civil engineering and construction work, earthquakes, heavy rain, etc., the governing equations discretized using methods such as the finite element method and finite volume method are numerically analyzed to obtain approximate solutions for ground displacement and pore water pressure, and their effects are sometimes examined. The governing equations for saturated ground can be derived using Biot's saturated porous media theory. In this theory, saturated soil is modeled as a two-phase system consisting of the soil skeleton (solid phase) and pore water (liquid phase), and quantities such as the density and stress of the saturated soil are expressed as the superposition of the quantities of each phase.

[0003] The governing equations based on the above theory can be formulated in various ways, depending on how unknowns are handled and the approximation method used. For example, the UP formulation, in which the saturated soil balance equation shown in Equation 1, the pore water balance equation shown in Equation 2, and the continuity equation shown in Equation 3 are treated as conservative equations, and the unknowns are the displacement u of the soil skeleton and the pore water pressure p, is one such formulation. Soil-water coupled finite element analysis based on the discretized forms of Equation 1 and Equation 3 is widely applied (see, for example, Non-Patent Documents 1 to 3).

[0004]

number

[0005] Regarding the up formulation, when the displacement of the soil skeleton is discretized into nodes using first-order finite elements and the pore water pressure is discretized into elements using zero-order shape functions, the finite element method has been applied to the spatial discretization of the equilibrium equation for saturated soil, and the difference method and finite volume method have been applied to the spatial discretization of the continuity equation. [Prior art documents] [Non-patent literature]

[0006] [Non-Patent Document 1] Akai and Tamura: Numerical Analysis of Multidimensional Consolidation Using Elastic-Plastic Constitutive Equations, Journal of the Japan Society of Civil Engineers, No. 269, pp. 95-104, 1978. [Non-patent document 2] Asaoka, Noda and Fernando:Effects of changes in geometry on the linear elastic consolidation deformation, Soils and Foundations, Vol.37, No.1, pp.29-39, 1997. [Non-patent document 3] Oka and Kimoto: Computational modeling of multiphase geomaterials, pp.171-210, CRC press, Taylor & Francis group, 2012. [Non-patent document 4] Takeyama, Iizuka, and Ohta: Spatial Discretization of Hydraulic Head Using Linear Function Approximation, 41st Geotechnical Engineering Research Conference, pp.321-322, 2006. Summary of the Invention [Problem to be solved by the invention]

[0007] In the conventional discretization format of the continuum type, the flow associated with the pressure gradient (or hydraulic gradient) on the element boundary surface is determined only by the water pressure difference (or head difference) between two adjacent elements at the boundary surface and geometric conditions, so the amount of pore water outflow in the element and other analysis results are overly dependent on the mesh division, which may result in a decrease in analysis accuracy (see, for example, non-patent document 4).

[0008] The objective of this invention is to propose an analytical method and a program for running this method on a computer that enables improved accuracy by increasing the resolution of the spatial distribution of pressure gradients (or hydraulic gradients) in finite element analysis of saturated ground based on the up formulation, when soil skeleton displacements are discretized into nodes and pore water pressures are discretized into elements. [Means for solving the problem]

[0009] In order to solve the above problem, the method for analyzing saturated ground of the present invention discretizes and solves the local form of the continuity equation described by the up-p formulation by applying the finite volume method, and discretizes the permeability coefficient within an element into nodes using the equivalent permeability coefficient matrix shown in Equation 4.

[0010]

number

[0011] The program of the present invention is a program that causes a computer to execute finite element analysis based on the up formulation using the saturated ground analysis method to obtain a numerical solution for the behavior of a ground model.

[0012] This method and program for analyzing saturated ground makes it possible to spatially discretize the hydraulic conductivity within an element with equivalent values ​​at each node. Furthermore, this allows the pressure gradient vector (or hydraulic gradient vector) at the node, which is represented by the pore water pressure of the element, to be used to more precisely represent the flow velocity vector within the ground model at each node. This in turn allows the pressure gradient vector (or hydraulic gradient vector) to be used to more precisely represent the amount of pore water outflow in the element, thereby improving the accuracy of analysis methods based on the UP formulation.

[0013] In order to express the flow velocity vector at a node using the pressure gradient vector at that node, which is expressed by the pore water pressure of the element, it is desirable to discretize the local form of the pore water balance equation shown in Equation 2 into the global stiffness equation shown in Equation 5 using the finite element method. In addition, it is desirable to find the pore water pressure to be applied to satisfy the flow boundary condition by adding the flow constraint equation (Equation 6) to the element boundary surface borne by the node that forms the flow boundary. Furthermore, it is desirable to express the flow velocity vector at the node N of interest using the equivalent permeability matrix and pressure gradient vector as in Equation 7.

[0014]

number

[0015]

number

[0016] In Equation 7, the flow associated with the pressure gradient is expressed by the second, third, and fifth terms on the right-hand side. In other words, the pressure gradient vector and hydraulic gradient vector at the node N of interest are expressed as Equations 8 and 9.

[0017]

number

[0018] The saturated ground analysis method and program of the present invention can improve the accuracy of finite element analysis based on the up formulation, where soil skeleton displacement is discretized into nodes and pore water pressure is discretized into elements. In other words, when predicting the behavior of saturated ground due to construction work, earthquakes, etc., through numerical simulation, this contributes to improving the reliability of evaluations of the impact of ground displacement, pore water pressure, and other changes on surrounding structures. Furthermore, it enables studies aimed at realizing effective countermeasures and rational foundation structures. [Brief explanation of the drawings]

[0019] [Figure 1] FIG. 10 is a diagram showing an example of calculation of an equivalent hydraulic conductivity matrix. [Figure 2] FIG. 10 is a diagram showing an example of setting a flow velocity vector at a node, expressed using an equivalent permeability matrix and a pressure gradient vector. [Figure 3] FIG. 1 shows the ground model used for verification by one-dimensional vertical seismic analysis of saturated ground. [Figure 4]10 is a graph showing the results of verification by one-dimensional vertical seismic analysis of saturated ground, where (a) is the analysis result of a comparative example, and (b) is the analysis result of an example. [Figure 5] 1 shows the ground models used for verification by static seepage analysis, where (a) is an I-type model and (b) is an L-type model. [Figure 6] 10 is a graph showing the verification results of static penetration analysis, where (a) is the analysis result obtained using an I-type model, and (b) is the analysis result obtained using an L-type model. [Figure 7] This figure shows the ground model consisting of quadrilateral elements and triangular elements used for verification by quasi-static consolidation analysis. (a) is a 4m wide ground model (Case 1), and (b) is a 1m wide ground model (Case 2). [Figure 8] 10 is a graph showing the verification results of quasi-static consolidation analysis, where (a) is the analysis result of a comparative example, and (b) is the analysis result of an example. DETAILED DESCRIPTION OF THE INVENTION

[0020] Below, we will explain a method for analyzing saturated ground according to an embodiment of the present invention. Here, the entire analysis domain is V, and we will target a soil-water coupled finite element analysis in which the displacement u of the soil skeleton is discretized into nodes using first-order finite elements, and the pore water pressure p is discretized into elements using zeroth-order shape functions. We will also assume that strain is infinitesimal, soil particles are incompressible, and temperature changes are negligible. The signs of stress and strain are positive for tensile components, and the sign of pore water pressure is positive for compression. <Introduction of equivalent hydraulic conductivity matrix> First, the hydraulic conductivity within an element is spatially discretized with equivalent values ​​at the nodes. To do this, an "equivalent hydraulic conductivity matrix" is introduced. The equivalent hydraulic conductivity matrix is ​​a square matrix whose order is the number of degrees of freedom for node displacement, and is also a diagonal matrix. The components of the equivalent hydraulic conductivity matrix corresponding to node N are calculated as shown in Equation 4.

[0021]

number

[0022] Figure 1 shows an example of the calculation of the equivalent hydraulic conductivity matrix. The calculation example in Figure 1 is for a case where two elements of two dimensions, rectangle, and equal volume (volume: BH) are adjacent, and the hydraulic conductivity (k E1 , k E2 ) is assumed to be homogeneous and isotropic. In the calculation example of Figure 1, for each node in the upper, middle, and lower rows, the "total volume of elements discretized at the node in the discretization format of the pore water balance equation," "total volume of elements weighted by the inverse of the hydraulic conductivity (= 1 / k)," and "components of the equivalent hydraulic conductivity matrix corresponding to each node" are shown in Table 1.

[0023] [Table 1]

[0024] <Spatial discretization of flow velocity vectors and expression using equivalent permeability matrix> Assuming that the average relative acceleration of pore water with respect to the soil skeleton is sufficiently small compared with the acceleration of the soil skeleton, the movement of pore water is described by the pore water equilibrium equation shown in Equation 2.

[0025]

number

[0026] When the local form of the pore water balance equation shown in Equation 2 is discretized using the finite element method, the global stiffness equation shown in Equation 5 is obtained.

[0027]

number

[0028] Here, the pore water mass matrix and hydraulic resistance matrix are lumped matrices with off-diagonal elements of 0. Furthermore, if the element boundary surface that makes up the element is a flow boundary, then in order to satisfy this condition, the pore water pressure is applied as an undetermined multiplier to each range borne by the constituent nodes of said boundary surface. This pore water pressure can be found using Lagrange's undetermined multiplier method, with Equation 6 as the constraint equation for flow rate.

[0029]

number

[0030] From the discretized form of the pore water balance equation shown in Equation 5, the flow velocity vector at node N can be formally calculated as Equation 10. Here, since the hydraulic resistance matrix is ​​a lumped matrix, its inverse matrix is ​​also a lumped matrix. Furthermore, since the pore water mass matrix is ​​also a lumped matrix, the flow velocity vector at the node N of interest is expressed only by the quantities related to node N using the equivalent hydraulic conductivity matrix and pressure gradient vector, as shown in Equation 7, and is uncoupled from the quantities of other nodes.

[0031]

number

[0032] <Setting constraint equations for flow boundary> Furthermore, the pore water pressure that each node bears on the element boundary surface to satisfy the flow rate boundary condition is found by adding a constraint equation for the flow rate shown in Equation 6 to the element boundary surface using Lagrange's undetermined multiplier method. In this embodiment, Equation 11, which is rewritten by substituting Equation 7 into Equation 6, is used as the constraint equation.

[0033]

number

[0034] Equation 7 is the general form of the velocity vector at a node, and indicates that the flow associated with the pressure gradient is determined by the second, third, and fifth terms on the right-hand side. Here, as shown in the example setting in Figure 2, if the node N of interest constitutes only the element boundary surface where elements are adjacent to each other (hereinafter referred to as the adjacent surface), the flow associated with the pressure gradient will be determined only by the second term on the right-hand side, and the expression of the velocity vector at node N will be simplified to Equation 12.

[0035]

number

[0036] <Discretization of continuity equations using the finite volume method> The local form of the continuity equation is a form that can take into account the effect of the spatial gradient of the hydraulic conductivity, as shown in Equation 3. Here, the left side of Equation 3 represents the "pore water outflow amount." The right side, using two terms, represents the "volume compression amount per unit time of saturated soil equivalent to the pore water outflow amount," which is expressed by subtracting the volume compression amount per unit time of pore water due to an increase in water pressure from the volume compression amount per unit time of the soil skeleton.

[0037]

number

[0038] Multiplying the continuity equation in the local form shown in Equation 3 by the weight w, integrating it over the region V, and applying Gauss's divergence theorem, we obtain the weak form shown in Equation 13. Furthermore, the area V of Equation 13 is the area V of element E. E and the weight w is replaced by the area V E If we set the value of the continuity equation to 1 inside and 0 outside, we obtain the continuity equation (14) discretized into element E by the finite volume method.

[0039]

number

[0040] Here, the left side of Equation 14 is expressed as a flow velocity vector at the element constituent node. The left side of Equation 14 is the pore water outflow term, and for the element E of interest, it can be written as Equation 15.

[0041]

number

[0042] The right side of Equation 15, i.e., the left side of Equation 14, can be expressed as Equation 16 using Equation 7.

[0043]

number

[0044] The right side of equation 14 is expressed by equations 17 and 18, as in the prior art.

[0045]

number

[0046] By replacing the left side of equation 14 with the right side of equation 16 and substituting equations 17 and 18 into the right side of equation 14 and rearranging, the continuity equation for element E, which has been discretized using the finite volume method, is obtained as shown in equation 19, and by superimposing the continuity equations for each element, the continuity equation for the entire analysis domain is constructed as shown in equation 20. In equation 20, a tilde (~) is added to each coefficient matrix to indicate the superposition for each element. According to Equation 20, the first term on the right-hand side appropriately discretizes the effect of earthquake-induced body forces on pore water flow, improving the accuracy of analysis of earthquake behavior near stratum boundaries and flow rate boundaries where permeability changes significantly.

[0047]

number

[0048] <Discretization of equilibrium equations for saturated soil using the finite element method> As with the prior art, the equilibrium equation for saturated soil is obtained by discretizing the local form shown in Equation 1 using the finite element method and assembling and applying the global stiffness equation of Equation 21.

[0049]

number

[0050] <Calculating the numerical solution of the ground model by running the program> The numerical solution of the ground model is calculated using a computer, which performs sequential analysis using a program that executes soil-water coupled finite element analysis based on the saturated ground analysis method of this embodiment. When the program is executed via a computer, the equilibrium equation for saturated soil discretized using the finite element method (Equation 21), the continuity equation discretized using the finite volume method (Equation 20), and the constraint equation for the flow rate at the element boundary surface borne by the nodes forming the flow rate boundary (Equation 11) are solved as simultaneous equations, thereby calculating numerical solutions for various behaviors of the ground model.

[0051] As described above, the saturated ground analysis method of this embodiment can improve the accuracy of finite element analysis based on the up formulation, where soil skeleton displacements are discretized into nodes and pore water pressures are discretized into elements. In other words, when predicting the behavior of saturated ground due to construction work, earthquakes, etc., through numerical simulation, this contributes to improving the reliability of evaluations of the impact of ground displacement, pore water pressure, and other changes on surrounding structures. Furthermore, it becomes possible to consider the realization of effective countermeasures and rational foundation structures.

[0052] In conventional analysis methods, the effect of earthquake-induced body forces on porewater flow is ignored in the process of deriving the local form of the continuity equation, based on the assumption that the spatial gradient of the permeability coefficient is sufficiently small, resulting in a problem of reduced analytical accuracy of earthquake behavior near stratum boundaries or flow rate boundaries where permeability changes significantly.However, with the analysis method for saturated ground of this embodiment, the effect of earthquake-induced body forces on porewater flow is appropriately taken into account in the discretized form of the continuity equation, thereby improving the analytical accuracy of earthquake behavior.

[0053] Furthermore, in the discretized continuous equation format used in conventional analysis methods, the flow associated with the pressure gradient (or hydraulic gradient) on the element boundary surface is determined only by the water pressure difference (or head difference) between two adjacent elements at that boundary surface and geometric conditions, which can result in the amount of pore water outflow in the element and other analysis results being overly dependent on mesh division, resulting in reduced analysis accuracy.However, the saturated ground analysis method of this embodiment can suppress mesh dependency and the resulting reduction in analysis accuracy.

[0054] On the other hand, according to the method for analyzing saturated ground of this embodiment, the effect of earthquake-induced body forces on the flow of pore water is appropriately taken into account in the discretized continuous form, thereby improving the accuracy of analyzing behavior during an earthquake near stratum boundaries and flow rate boundaries where permeability changes significantly.

[0055] Furthermore, according to the method for analyzing saturated ground of this embodiment, the pressure gradient is defined at the nodes based on the discretization form of the pore water balance equation using the finite element method, and the resolution regarding the spatial distribution of the pressure gradient is increased, thereby improving the accuracy of analysis when the analysis mesh contains distorted quadrangles that are not rectangular or square, or when the permeability direction is distributed in a complex manner.

[0056] Furthermore, according to the method for analyzing saturated ground of this embodiment, the discretization form of each governing equation is formulated in vector matrix representation, so the phenomena under consideration are not limited to one- and two-dimensional ones, and numerical solutions can be obtained for phenomena with three-dimensional geometric shapes, improving the accuracy of the analysis regardless of the dimension. [Example]

[0057] Next, three analytical examples performed to confirm the effects of this embodiment are shown. Note that in the following analytical examples, an analysis using conventional technology was also performed for comparison, and the pressure gradient on adjacent surfaces in that case was expressed by dividing the water pressure difference between two adjacent elements by the distance between the centers of gravity of the elements, with reference to Non-Patent Document 3. First, we present the analytical results when vertical earthquake motion is input to a saturated poroelastic soil column model consisting of two layers with different permeabilities. Figure 3 shows the soil column model (layer thickness: 10 m) used in the analysis, and Table 2 shows the conditions such as the seismic conditions and soil properties. As shown in Table 2, the soil column model is composed of two strata with different permeabilities, with a depth of 5.0 m as the strata boundary. The analysis method used was dynamic analysis that took into account the inertial term. Furthermore, since there is no theoretical solution for the one-dimensional vertical vibration problem in this example, the numerical solution obtained by the UWP formulation, which is more rigorous than the UP formulation in this example, was used as the reference solution.

[0058] [Table 2]

[0059] Figure 4 shows the time history of pore water pressure obtained in an element at a central depth of 5.1 m when the seismic frequency was 0.2 Hz, as a representative analysis result near a layer boundary. Figure 4(a) shows the analysis result using conventional technology (comparison example), and Figure 4(b) shows the analysis result using this embodiment (example). Each figure also shows the time history of the reference solution. Table 3 shows the error from the reference solution of pore water pressure at a depth of 5.1 m obtained in the comparative example and the example. Here, the "error from the reference solution" is the "maximum difference from the reference solution" divided by the "maximum amplitude value of the reference solution." The numerical solution according to the comparative example has a larger pore water pressure amplitude than the reference solution, as shown in Figure 4(a), but the numerical solution according to the embodiment has a much smaller difference from the reference solution than the comparative example, as shown in Figure 4(b) and Table 3, confirming the improvement in the analytical accuracy of the seismic behavior according to this embodiment.

[0060] [Table 3] [Example]

[0061] Next, we present the steady-state results of static seepage analysis carried out using two soil models consisting of distorted quadrilateral elements. The ground models used in the analysis were the I-type model shown in Figure 5(a) and the L-type model shown in Figure 5(b), and certain element boundaries were made undrained (flow boundaries where no pore water flows in or out) so that the macroscopic flow of pore water occurs from the upstream water pressure boundary (water pressure: 49.05 kPa) to the downstream water pressure boundary (water pressure: 9.81 kPa), in the I-type model as an I-type, and in the L-type model as an L-type. Furthermore, the conditions for the soil properties were set as shown in Table 4, and an analysis was also performed using a conventional analysis method as a comparative example.

[0062] [Table 4]

[0063] Figure 6(a) shows the steady-state pore water pressure distribution obtained for the I-type model, and Figure 6(b) shows the L-type model. The theoretical solution shown in each figure is calculated as the permeability distance, which is the length of the model's central axis (= 6m) connecting three points: the center of the upstream hydraulic boundary surface, the center of the adjacent surface between elements, and the center of the downstream hydraulic boundary surface. As shown in Figures 6(a) and 6(b), the numerical solution of the comparative example (prior art) has a large error from the theoretical solution in the case of the L-shaped model. On the other hand, the numerical solution of the example (present embodiment) is in good agreement with the theoretical solution for both the I-shaped and L-shaped ground models. In other words, in Example 2, the effect of improving the resolution of the spatial distribution of the pressure gradient was confirmed in the case where the element shape is a distorted quadrangle. [Example]

[0064] Next, the results of a quasi-static consolidation analysis conducted using a ground model consisting of diamond-shaped quadrilateral elements and isosceles triangular elements are presented. The ground model used in the analysis was set up with two analysis meshes (Case 1 and Case 2) with different horizontal division widths, as shown in Figure 7(a) and Figure 7(b). The boundary conditions were a consolidation load of 100 kN / m 2 The drainage conditions were set to one-sided drainage, and one-dimensional consolidation conditions were used. Furthermore, the conditions for the soil skeleton constitutive law and ground properties were set as shown in Table 5, and an analysis using conventional technology was also performed as a comparative example.

[0065] [Table 5]

[0066] Figure 8 shows the time history of the average degree of compaction as a representative analysis result. Here, Figure 8(a) shows the analysis result using conventional technology, and Figure 8(b) shows the analysis result using this embodiment. Each figure also shows the theoretical solution based on Terzachi's one-dimensional consolidation theory. As shown in Figure 8(a), in the analysis using the conventional technology, the error from the theoretical solution is large in both cases. In addition, the degree of consolidation progresses faster in Case 2 than in Case 1, and differences due to the analytical mesh are recognized. On the other hand, in the analysis using this embodiment, as shown in Figure 8(b), the mesh dependency of the numerical solution as in the case of the conventional technology is not observed, and both Case 1 and Case 2 accurately reproduce the theoretical solution. Example 3 confirmed that the accuracy of analysis when the analysis mesh contains distorted quadrilaterals can be improved by specifying the pressure gradient at the nodes, even when the deformation of the soil skeleton and the flow of pore water are in an unsteady state.

[0067] <Modification> The present invention is not limited to the above-described embodiment, and the above-described components can be appropriately modified without departing from the spirit and scope of the present invention. For example, in the embodiment, to satisfy the flow boundary condition, Lagrange's undetermined multiplier method is applied, and the pore water pressure acting on the element boundary surface borne by each node is set as an undetermined multiplier, but the penalty method may be applied instead of Lagrange's undetermined multiplier method. By setting the penalty with respect to the flow in the normal direction of the flow boundary, the pore water pressure acts on the element boundary surface borne by each node forming the flow boundary, and the flow boundary condition can be approximately satisfied.

[0068] When applied to finite element analysis of unsaturated soil, for example, the soil skeleton displacement, pore water pressure, and pore air pressure are formulated as unknowns, and the unsaturated soil equilibrium equation, pore water continuity equation, and pore air continuity equation are conservative equations. The soil skeleton displacement and pore water pressure are discretized as in the previous embodiment, and the pore air pressure is discretized into elements using a zeroth-order shape function. Furthermore, the unsaturated soil equilibrium equation is discretized using the finite element method, and the pore water and pore air continuity equations are discretized using the finite volume method. In the discretization process of the pore air continuity equation using the finite volume method, a matrix is ​​introduced that discretizes the air permeability coefficient within the element into nodes so that the pore air outflow rate in the element can be expressed using the pore air pressure gradient vector at the node.

[0069] The shape of the finite element can be set arbitrarily within the category of first-order elements, and the placement of integration points and the method of numerical integration can also be selected appropriately.

[0070] In the above-described embodiment, the permeability coefficient is assumed to be homogeneous and isotropic within the element, but the saturated ground analysis method of the present invention is not limited to this. It can also be applied to cases where the permeability coefficient within the element is not homogeneous, and the anisotropy of the permeability coefficient can also be taken into account.

[0071] In the above-described embodiment, Newmark's β method was applied to dynamic problems and the backward difference method was applied to quasi-static problems for time integration, but the saturated ground analysis method of the present invention is not limited to these. Any method can be used for time integration, such as Wilson's θ method or forward difference method. Furthermore, the constitutive model that defines the relationship between the effective stress increment and the strain increment of the soil skeleton is not limited to a linear elastic model, but can be any model depending on the purpose and conditions, such as an elastoplastic constitutive model or other nonlinear models. Furthermore, in the above-described embodiment, it is assumed that the soil particles are incompressible and that temperature changes are negligible, but the saturated ground analysis method of the present invention is not limited to this. It can also be applied to cases where the compressibility of soil particles or temperature changes are taken into account.

[0072] The nonlinearity of the soil-water coupled finite element analysis according to the embodiment of the present invention can take into account not only material nonlinearity due to an elasto-plastic constitutive model, but also geometric nonlinearity due to the total Lagrange method, updated Lagrange method, or the like.

Claims

1. This is a method for analyzing saturated ground based on the up formulation, which has the soil skeleton displacement u and pore water pressure p as unknowns. When the soil skeleton displacement u is discretized into nodes and the pore water pressure p is discretized into elements, the continuity equation, which is one of the governing equations in the up formulation, is discretized and solved by applying the finite volume method. The method is characterized in that the flow velocity vector at a node is expressed using a pressure gradient vector or hydraulic gradient vector at that node, which is represented by the pore water pressure of the element, and the permeability coefficient within the element is discretized into nodes using the equivalent permeability matrix shown in Equation 1, so that the amount of pore water outflow in the element composed of that node can be expressed using the pressure gradient vector or hydraulic gradient vector via the flow velocity vector. [Equation 1]

2. A method for analyzing saturated ground as described in claim 1, characterized in that the local form of the pore water balance equation is discretized into the overall stiffness equation of Equation 2 using the finite element method in order to express the flow velocity vector at the node using the pressure gradient vector at the node represented by the pore water pressure of the element. [Equation 2]

3. The method for analyzing saturated ground described in claim 2, characterized in that the pore water pressure to be applied to satisfy the flow boundary condition is set as an undetermined multiplier, and this is obtained by adding a constraint equation regarding the flow rate shown in Equation 3. [Equation 3]

4. 2. The method for analyzing saturated ground according to claim 1, wherein the flow velocity vector at the node is expressed by Equation 4 using the equivalent permeability matrix and the pressure gradient vector. [Equation 4]

5. A program for causing the saturated ground analysis method according to any one of claims 1 to 4 to function on a computer.

Images