A Picard-PCG iteration based accelerating algorithm for solving unsaturated seepage

Abstract

The application discloses a non-saturated seepage solving acceleration algorithm based on a Picard-PCG iteration method, and comprises the following steps: (1) constructing a nonlinear control equation and a soil moisture characteristic curve model of non-saturated seepage; (2) obtaining an expression of phase mobility M from the nonlinear control equation, and obtaining expressions of capillary water capacity C(h) and effective saturation θ from the soil moisture characteristic curve model; (3) linearizing the nonlinear control equation of non-saturated seepage into a linear control equation by using an iteration-accelerated Picard method, and performing internal iteration and external iteration on the linear control equation by using a preconditioned conjugate gradient method and the iteration-accelerated Picard method respectively; and (4) substituting a non-saturated seepage value obtained by solving into a visualization software to perform seepage analysis, and obtaining a distribution of rainfall infiltration and seepage flow. θ e The application realizes second-order iteration acceleration solving of the non-saturated seepage value by combining the preconditioned conjugate gradient method and the iteration-accelerated Picard method.​

Description

Technical Field

[0001] This application relates to an accelerated algorithm for solving unsaturated seepage based on the Picard-PCG iterative method, belonging to the technical field of research on the relative permeability coefficient of unsaturated soil in civil and geotechnical engineering. Background Technology

[0002] Unsaturated seepage properties play a crucial role in improving the accuracy of numerical solutions, and the study of the permeability characteristics of unsaturated soils is gradually gaining momentum. Seepage problems in unsaturated soils require numerical calculation methods for solution, enabling the prediction and analysis of soil moisture content distribution and changes. In the analysis of archaeological site deterioration, rainwater has a significant impact. The application of unsaturated seepage analysis mainly includes rainfall infiltration and slope stability analysis. Therefore, how to effectively and accurately simulate and analyze unsaturated seepage problems has always been a popular research topic.

[0003] The fundamental equation for unsaturated seepage is the Richards equation, which is a strongly nonlinear and degenerate partial differential equation. Therefore, solving this equation is a key issue in unsaturated seepage analysis. In the calculation of the Richards equation, the backward Euler method is currently the primary method for time discretization, while spatial discretization methods include the finite element method, the finite element volume method, and the finite element difference method. However, the performance of both spatial and temporal discretization strategies largely depends on nonlinear iterative algorithms. The Picard iteration method is currently the most widely used algorithm for solving unsaturated seepage equations, but its convergence is difficult to guarantee, and oscillating convergence may occur during the convergence process, leading to slow convergence or even non-convergence. Summary of the Invention

[0004] In view of this, an accelerated algorithm for solving unsaturated seepage problems based on the Picard-PCG iterative method is proposed. The proposed algorithm is a nonlinear iterative acceleration algorithm based on Picard-PCG with second-order convergence. First, the equation system is linearized by the Picard method with iterative acceleration, and then the linear equation system is solved by the preprocessed conjugate gradient method PCG. It can solve unsaturated seepage problems accurately and efficiently.

[0005] The specific plan is as follows:

[0006] An accelerated algorithm for solving unsaturated seepage based on the Picard-PCG iterative method includes:

[0007] (1) Construct the nonlinear control equations and soil moisture characteristic curve model for unsaturated seepage;

[0008] (2) The phase mobility M is calculated from the nonlinear control equation. θThe capillary water capacity C(h) and effective saturation θ are calculated from the soil moisture characteristic curve model using the given expressions. e The expression;

[0009] (3) The nonlinear control equations of unsaturated seepage are linearized by Picard method to obtain linear control equations. The preprocessed conjugate gradient method PCG and the iteratively accelerated Picard method are used to perform internal and external iterations on the linear control equations, respectively, thereby achieving second-order convergence and accelerated solution of the linear control equations.

[0010] (4) Substitute the unsaturated seepage value obtained from the solution into the visualization software to perform seepage analysis, thereby obtaining the distribution of water content and seepage flow of rainfall infiltration on the slope site.

[0011] Preferably, the nonlinear equation for the unsaturated seepage is expressed based on the Richards equation for pressure head as follows:

[0012]

[0013] In equation (1), C(h) is the capillary water capacity, h is the water head pressure in meters, K(h) is the hydraulic conductivity of the unsaturated porous medium, and z is the vertical distance.

[0014] Preferably, the solution obtained from (1) is characterized by...

[0015]

[0016] In equation (2) ρ θ Let M be the phase density, ||g||² represent the amplitude of the gravitational field, and M be the gravitational field amplitude. θ Let M be the phase mobility; thus, the phase mobility M is derived. θ The expression is:

[0017]

[0018] In equation (3), K is the inherent permeability of the porous medium, μ θ k is the viscosity of the liquid. r,θ Let be the relative permeability, which is a function of saturation θ.

[0019] Preferably, the soil moisture characteristic curve is represented using the Van Genuchten model, where the relationship between pressure head and saturation is as follows:

[0020]

[0021] In equation (4) θ s and θ r Here, saturated water content and residual water content are respectively, and α and m are parameters of the empirical model.

[0022] Preferably, the capillary water capacity C(h) is expressed using the Van Genuchten model as follows:

[0023]

[0024] The effective saturation θ is obtained from equation (5). e The expression is:

[0025]

[0026] In equations (5) and (6), α and m are parameters of the empirical model, and θ s and θ r These represent saturated moisture content and residual moisture content, respectively, θ (h) , calculated by equation (4), represents the saturation degree when the soil suction is h.

[0027] Preferably, the Picard-PCG iterative method is as follows:

[0028] First, the nonlinear governing equations of the unsaturated seepage are normalized using Picard to obtain the linear governing equations.

[0029] Then, the preprocessed conjugate gradient method (PCG) is used to perform internal iteration on the linear control equations, and the Picard method with accelerated iteration is used to perform external iteration on the linear control equations, thereby achieving second-order convergence and accelerating the solution of the linear control equations.

[0030] Preferably, the specific process of the Picard-PCG iterative method is as follows:

[0031] a. The linear governing equations obtained by linearizing the nonlinear governing equations of the unsaturated seepage using the Picard method with iterative acceleration are as follows:

[0032]

[0033] In equation (7), h j ρ represents the pressure head value at time j. θ Let M be the phase density, ||g||² represent the amplitude of the gravitational field, and M be the gravitational field amplitude. θ Phase mobility;

[0034] b. The pre-processed conjugate gradient method with diagonal incomplete Cholesky decomposition is used to perform the first calculation and solve equation (7) to obtain the pressure head value h for the (k+1)th time. j+1,k+1 ;

[0035] The pressure head value h at the (k+1)th time j+1,k+1 Substituting into equation (2) to update the relative permeability k r,θAnd the capillary water capacity C(h), then the pressure head value h of k+1 times. j+1,k+1 Substitute into the linear control equation (7) to calculate the first residual r k+1 ;

[0036] c. Using the pre-processed conjugate gradient method with diagonal incomplete Cholesky decomposition, equation (7) is calculated a second time to obtain the pressure head value h for the (k+2)th time. j+1,k+2 The pressure head value h of the (k+2)th time j+1,k+2 Substituting into equation (2) to update the relative permeability k r,θ and capillary water capacity C(h), pressure head h j+1,k+2 Substitute into equation (7) to calculate the second residual r k+2 ,

[0037] d. Based on the two iteration values ​​of b and c, recalculate the pressure head value h of the (k+1)th iteration using equation (7). j +1,k+1 The recalculated pressure head value is:

[0038]

[0039] e. The recalculated pressure head value h j+1,k+1 Substituting into equation (2) to update the relative permeability k r,θ The capillary water capacity C(h) is calculated and substituted into the linear control equation to calculate the final residual r. If the final residual r is less than the user-set ε, the loop exits and jumps to the next time step, outputting the current head pressure h and the current water pressure h. Moisture content θ r From this place Water head pressure h and Water content θ r As the final calculated value; if r is not less than ε, then jump to a and proceed to the next iteration.

[0040] Preferably, for equation (7), the first residual r obtained in the first solution is... k+1 The second residual r obtained from the second solution k+2 If the calculated first residual r k+1 Second residual r k+2 If the residual is less than the user-defined residual ε, exit the loop, jump to the next time step, and output the current pressure head and water content. If the calculated first residual r is less than the user-defined residual ε, then exit the loop, jump to the next time step, and output the current pressure head and water content. k+1 Second residual r k+2 If the residual is not less than the user-defined residual ε, proceed to the next step.

[0041] The beneficial effects that this application can produce include:

[0042] This application employs a combination of the pre-processed conjugate gradient method (PCG) and the iteratively accelerated Picard method to achieve second-order convergence and accelerate the solution of the nonlinear equations for unsaturated seepage, thereby obtaining the unsaturated seepage value. This improves the accuracy, efficiency, and stability of the algorithm for unsaturated permeability. Specifically, the PCG method is used for internal iteration of the linear equation system, while the Picard method is used to control the external iteration of the linear equation system. This achieves second-order convergence and accelerated solution of the linear control equation system through a double-loop system, accurately and efficiently solving for the unsaturated seepage value. This allows for the prediction and analysis of soil moisture content distribution and changes, providing a reliable guarantee for slope stability analysis. Attached Figure Description

[0043] Figure 1 This is a flowchart of an accelerated algorithm for solving unsaturated seepage based on the Picard-PCG iterative method, as described in this application.

[0044] Figure 2 This figure compares the convergence speed of the accelerated solution algorithm for unsaturated seepage based on the Picard-PCG iterative method of this application with that of the general Picard method. Detailed Implementation

[0045] The present application is described in detail below with reference to the embodiments, but the present application is not limited to these embodiments.

[0046] Example

[0047] As attached Figure 1 This invention discloses an accelerated algorithm for solving unsaturated seepage flow based on the Picard-PCG iterative method. It includes the following steps:

[0048] First, the Richards equation for unsaturated seepage based on head pressure is constructed as shown in equation (1), and the Richards equation for unsaturated seepage is discretized based on, but not limited to, the preferred element volume method.

[0049]

[0050] In equation (1), C(h) is the capillary water capacity, h is the head pressure in meters, K(h) is the hydraulic conductivity of the unsaturated porous medium, and z is the vertical distance. Equation (2) is obtained by solving equation (1), and equation (2) is expressed as:

[0051]

[0052] In equation (2) ρ θ Let M be the phase density, ||g||² represent the amplitude of the gravitational field, and M be the gravitational field amplitude. θ The phase mobility M is obtained from equation (2). θThe expression for , as shown in equation (3), can be expressed as:

[0053]

[0054] In equation (3), K is the inherent permeability of the porous medium, μ θ k is the viscosity of the liquid. r,θ Let k be the relative permeability, a function of the saturation θ, where the relative permeability k is... r,θ It can be represented using a water-holding model such as the Van Genuchten model.

[0055] Meanwhile, the relationship between pressure head and saturation is represented using the Van Genuchten model, as shown in equation (4).

[0056]

[0057] Equation (4) yields the capillary water capacity C(h) and effective saturation θ, respectively. e The expression for the capillary water capacity C(h) is shown in equation (5).

[0058]

[0059] Effective saturation θ e The expression for is shown in equation (6).

[0060]

[0061] In equations (4), (5), and (6), θ s and θ r Here, α and m are the saturated water content and residual water content, respectively; θ is the parameter of the empirical model. (h) The value obtained from equation (4) represents the saturation of the soil when the suction force is h.

[0062] Then, the Picard-PCG iterative method is used to accelerate the solution of the unsaturated seepage equation. The specific solution process is as follows:

[0063] The Picard-PCG iterative method is used to accelerate the solution of the unsaturated seepage equation;

[0064] Accelerated Iteration Algorithm:

[0065] Start loop

[0066] a. Using the Picard method to linearize the nonlinear governing equations of unsaturated seepage, the pressure head h at the (k+1)th time is... j +1,k+1 The calculation formula is:

[0067]

[0068] Where h j ρ represents the pressure head value at time j. θ Let M be the phase density, ||g||² represent the amplitude of the gravitational field, and M be the gravitational field amplitude. θ This represents the phase mobility.

[0069] b. The pressure head h at the (k+1)th time is obtained by calculating the preprocessed conjugate gradient method (PCG) with diagonal incomplete Cholesky decomposition. j+1,k+1 .

[0070] According to the new pressure head value h j+1,k+1 Update the relative permeability and capillary capacity, and substitute them into the governing equation to calculate the first residual r. k+1 ,

[0071] If the first residual r k+1 If the value is less than ε, where ε is the residual set by the user, then exit the loop, jump to the next time step, and output the current pressure head and water content.

[0072] Calculate the pressure head h for the (k+2)th time according to formula a. j+1,k+2 The results were obtained through DIC-PCG calculations, the relative permeability and capillary capacity were updated, and the results were substituted into the governing equations to calculate the second residual r. k+2 If the second residual r k+2 If the value is less than ε, where ε is the residual set by the user, then exit the loop, jump to the next time step, and output the current pressure head and water content.

[0073] Based on the values ​​from the previous two iterations, recalculate the pressure head h for the (k+1)th iteration. j+1,k+1 :

[0074]

[0075] Update the relative permeability and capillary capacity, and input them into the governing equation to calculate the final residual r. If the final residual r is less than ε (where ε is the user-defined residual), exit the loop, jump to the next time step, and output the current pressure head and water content. If r is not less than ε, jump to step a for the next iteration.

[0076] Wherein, the design residual ε is the minimum allowable residual set by the user, and the residual r is the value obtained by substituting the calculated head pressure h into equation (2). The first residual r k+1 Second residual r k+2 The comparison with the design residual ε is necessary to ensure that the algorithm terminates promptly once the user-defined computational precision is achieved, thereby reducing unnecessary loop steps.

[0077] The new pressure head h is recalculated using equation (8). j+1,k+1 The value is the actual iteration value.

[0078] Finally, the obtained unsaturated seepage values ​​are substituted into visualization software such as Paraview for seepage analysis to obtain the distribution and changes in rainfall or precipitation moisture content.

[0079] In this application:

[0080] By using the PCG algorithm to achieve second-order convergence and by solving equation (8) through data processing, the Picard method iteration, which is accelerated by iteration, also achieves second-order convergence. Thus, both internal and external iterations achieve second-order convergence accuracy, achieving an acceleration effect.

[0081] Second-order convergence can also be achieved using the Newton algorithm to solve the unsaturated seepage equation. However, while the Newton algorithm can achieve second-order convergence accuracy, it requires additional calculation of the Jacobian matrix, which is too time-consuming.

[0082] To verify the algorithm of Example 1, this application used the free drainage experiment at Los Alamos National Laboratory for numerical simulation. This application chose this example because it demonstrates the transition from saturated to unsaturated conditions due to free drainage at the bottom of the soil column.

[0083] This experiment used a 6m soil column. To ensure that no water entered the soil column, zero flux was applied at the top boundary of the soil column, and zero pressure head was maintained at the bottom boundary of the soil column so that water could be continuously discharged.

[0084] The experimental procedure is as follows:

[0085] The soil column was discretized into 600 regular elements, and each element was initialized with zero pressure head to establish saturation conditions within the soil column. Over time, the soil column gradually transitioned from a saturated state to an unsaturated state. The Van Genuchten model parameters were set to α = 2.5, m = 1.506, and θ... s θ is 0.33. r It is 0.

[0086] Depend on Figure 2 It can be seen that for transient problems, under given boundary conditions, the first time step deviates the most from the final stable state, and the number of iterations is also the highest. The choice of time step size also significantly affects the number of iterations; the larger the time step, the more iterations. We selected a time step of 2 seconds. The relationship between the number of iterations and the residuals for the two algorithms within the first time step is shown in the figure. The general Picard method exhibits oscillatory convergence with a very slow convergence speed, while our accelerated Picard method does not exhibit oscillatory convergence, has higher stability, and its convergence accuracy and speed are far faster than the general algorithm.

[0087] It should be noted that:

[0088] In this application, the Van Genuchten model is an empirical formula describing the characteristic curve of soil moisture. This curve reflects the relationship between soil moisture energy and quantity, and the relationship between pressure head and saturation.

[0089] The above description is merely an embodiment of this application and is not intended to limit this application in any way. Although this application discloses preferred embodiments as described above, it is not intended to limit this application. Any changes or modifications made by those skilled in the art without departing from the scope of the technical solution of this application using the disclosed technical content are equivalent to equivalent implementation cases and fall within the scope of the technical solution.

Claims

1. An accelerated algorithm for solving unsaturated seepage based on the Picard-PCG iterative method, characterized in that, include: (1) Construct the nonlinear governing equations and soil moisture characteristic curve model for unsaturated seepage; (2) Obtain the phase mobility from the nonlinear control equation. The expression for capillary water capacity is obtained from the soil moisture characteristic curve model. and effective saturation The expression; (3) The nonlinear control equation of unsaturated seepage is obtained by Picard method. The linear control equation is then obtained by using the preprocessed conjugate gradient method PCG and the iteratively accelerated Picard method to perform internal and external iterations on the linear control equation, thereby achieving second-order convergence and accelerated solution of the linear control equation. (4) Substitute the obtained unsaturated seepage values ​​into the visualization software for seepage analysis to obtain the distribution and changes of rainfall or precipitation moisture content; The nonlinear equation for the unsaturated seepage, based on the Richards equation for pressure head, is expressed as follows: Equation (1); In equation (1), For capillary water capacity, This refers to the head pressure, measured in meters (m). The hydraulic conductivity of unsaturated porous media The distance is in the vertical direction; Solving equation (1) yields: Equation (2); In formula (2) For phase density, This represents the amplitude of the gravitational field. Phase mobility; The Picard-PCG iterative method is as follows: First, the nonlinear governing equations of unsaturated seepage are linearized using the Picard method with iterative acceleration to obtain linear governing equations. Then, the preprocessed conjugate gradient method PCG is used to perform internal iteration on the linear control equations, and the Picard method with iterative acceleration is used to perform external iteration on the linear control equations, thereby achieving second-order convergence and accelerating the solution of the linear control equations. The specific process of the Picard-PCG iterative method is as follows: a. The nonlinear governing equations of the unsaturated seepage flow are linearized using the Picard method with iterative acceleration to obtain the linear governing equations: Equation (7); In equation (7), This represents the pressure head value at time j. For phase density, This represents the amplitude of the gravitational field. Phase mobility; b. The pre-processed conjugate gradient method with diagonal incomplete Cholesky decomposition is used to solve equation (7) for the first time to obtain the pressure head value of the (k+1)th time. ; The pressure head value at the (k+1)th time Substituting into equation (2) to update the relative penetration rate and capillary water capacity Then, the pressure head value of k+1 times. Substitute into the linear control equation (7) to calculate the first residual r k+1 ; c. The pre-processed conjugate gradient method with diagonal incomplete Cholesky decomposition is used to perform a second calculation on equation (7) to obtain the pressure head value for the (k+2)th time. The pressure head value of the (k+2)th time Substituting into equation (2) to update the relative penetration rate and capillary water capacity Pressure head value Substitute into equation (7) to calculate the second residual r k+2 ; d. Based on the two iteration values ​​of b and c, recalculate the pressure head value for the (k+1)th iteration using equation (7). The recalculated pressure head value is: Equation (8); e. The obtained pressure head value will be recalculated. Substituting into equation (2) to update the relative penetration rate and capillary water capacity Then, substitute the values ​​into the linear control equation to calculate the final residual r. If the final residual r is less than the user-set value... If the loop terminates, the loop will exit and proceed to the next time step, outputting the current head pressure h and water content. Based on the water head pressure h and water content at this location As the final calculated value; If the final residual r is not less than If the condition is met, then jump to a and proceed to the next iteration.

2. The accelerated algorithm for solving unsaturated seepage based on the Picard-PCG iterative method according to claim 1, characterized in that, The phase mobility can be obtained by solving equation (2). The expression is: Equation (3); In formula (3) The inherent permeability of porous media, The viscosity of the liquid. Relative permeability is relative to saturation. The function.

3. The accelerated algorithm for solving unsaturated seepage based on the Picard-PCG iterative method according to claim 1, characterized in that, The soil moisture characteristic curves are represented using the Van Genuchten model, where the relationship between pressure head and saturation is as follows: Equation (4); In formula (4) and These are saturated moisture content and residual moisture content, respectively. and These are the parameters of the empirical model.

4. The accelerated algorithm for solving unsaturated seepage based on the Picard-PCG iterative method according to claim 3, characterized in that, Capillary capacity Represented using the Van Genuchten model: Equation (5); The effective saturation is obtained from equation (5). The expression is: Equation (6); In equations (5) and (6) and For the parameters of the empirical model, and These are saturated moisture content and residual moisture content, respectively. θ （h） The value is calculated from equation (4), which represents the saturation of the soil when the suction force is h.

5. The accelerated algorithm for solving unsaturated seepage based on the Picard-PCG iterative method as described in claim 1, characterized in that, For equation (7), the first residual r obtained in the first solution is respectively k+1 The second residual r obtained from the second solution k+2 If the calculated first residual r k+1 Second residual r k+2 Less than the user-set residual If the first residual r is obtained, the loop exits, jumps to the next time step, and outputs the current pressure head and water content. k+1 Second residual r k+2 Not less than the user-defined residual Then proceed to the next step.

Images