JPH03229156A - Batch simulation system for navier-stokes' equation - Google Patents

Abstract

PURPOSE:To obtain a stable numeric solution and to enable fast calculation by performing numeric simulation while the Navier-Stokes' equation and continu ous equations are made discrete at a time, and making the continuous equations hold at each time for calculation. CONSTITUTION:The Navier-Stokes' equation and continuous equations are made simultaneously discrete in units of a 4X4 block matrix containing a speed and pressure as unknowns to generate a block matrix A. Then block incomplete triangular decomposition is performed as to the block matrix A. An inverse matrix is an approximate matrix. Simultaneous linear equations including the matrix A as a coefficient are solved by the method of reciprocal solution of a conjugate gradient method system for asymmetry which uses said block incom plete triangular decomposition as processing. In this method, the Navier-Stokes' equation and continuous equations are only described as they are, so the easiness of use is improved. Further, the calculation is carried out while the continuous equations are made to hold at each time, so the calculation becomes fast and stability to time step-size is increased.

Description

[Detailed Description of the Invention] [Industrial Application Field] The present invention is applicable to numerical simulations of incompressible viscous fluids. It is concerned with numerical calculations using computers to find numerical solutions to problems.

[Conventional technology]

Regarding conventional incompressible viscous flow, see "Numerical Simulation" edited by Japan Society of Mechanical Engineers and published by Corona Publishing, pages 197 to 19.
8, the SMAC method (Simplified m
Arker and cell method) were used. Since this document does not perform two-dimensional scaling using the Re number, here we will show equations that are three-dimensional scaling based on the Re number, and the SMAC method corresponding to equations (1) to (4) of the embodiment of the present application. The superscript n indicates the calculated value of t=nΔ in time, and n+1 indicates the calculated value of t=(n+1)Δt. Here, Δt is the time step width.

W=W −Δt・(VV'grad(Wn)+−±2
wnaz Re Second step: Solve the following Poisson equation for φ using the auxiliary variable φ.

Third step: Using φ obtained above, U°゛1 Find P0+1.

, ■ n+1 , W n+1 KO:TEVV0= (U', V', Wo) LV'' indicates the Laplacian.

The drawback of this method is that the equation of continuity is not solved as is, but is modified using the auxiliary variable φ so that the equation of continuity is satisfied at the next time. Therefore, since the equation of continuity is not satisfied at the same time, the solution may become unstable, and there is a drawback that the time step width Δt cannot be made very large.

[Problem to be solved by the invention] The conventional technology has the above-mentioned problems.

The purpose of the present invention is to simultaneously discretize the Navier-Stokes equations and continuity equations, which are the governing equations of incompressible fluids, and perform numerical simulations, thereby achieving faster speeds than conventional methods and being less sensitive to time-step convergence and domain division methods. The objective is to provide a stable numerical solution method that is not influenced by In addition, in the calculation method of the block conjugate gradient system with incomplete triangular decomposition used at this time, we have devised a method for handling cases where there is an all-zero row in the diagonal block generated from five consecutive equations, and we The purpose is to find iterative solutions to linear equations.

[Means to solve the problem]

In order to achieve the above objective, for the Navier-Stokes equation and the continuity equation, which are the basic equations of flow, we created a simultaneous discretization method that kept the block length the same and a system of simultaneous linear equations using the same block length as coefficients. This allows numerical solutions to be obtained using a solution method for conjugate gradient groups with complete triangular decomposition.

In order to use a simultaneous discretization method that keeps the block length the same, set the diagonal elements on the diagonal block corresponding to the fixed variable to 1 at fixed points (points where the value is fixed), Set all values of matrix elements in the same row to zero. In addition, when using a staggered mesh (the position defined for each variable is different; in the case of a three-dimensional flow, velocity is placed on a surface and pressure is placed at the center point), the velocity U in the x, y, and z directions is
.

The number of definition points is different for v, W and pressure P. In this case, we create dummy points to keep the number of variables the same and the block length the same.

Incomplete triangular decomposition uses the following method to obtain an approximate inverse matrix of a diagonal block after triangular decomposition. Therefore, when σ is the row norm of the diagonal block.

If σ=0, this is generated from the continuity equation and the value depends only on the mesh decomposition, so σ is changed to the row norm including off-diagonals. Furthermore, when calculating the inverse matrix of the diagonal block after decomposition, if the pivot value P is P/σ<ε, then P=
Obtain an approximate inverse matrix that eliminates singularity by replacing it with σ×εx (sign of p). Furthermore, in order to make the diagonal block more dominant than the off-diagonal block, the diagonal block x (1+α) is used in the incomplete triangular decomposition, and the parameters ε and α are determined by numerical experiments. However, in many cases i=0.1. α; A value of 1.0 is sufficient.

[Effect]

The reason for using the simultaneous m-patterning method that maintains the same block length is that block incomplete triangular decomposition can be easily performed and a calculation method for supercomputers can be adopted.

Obtaining an approximate inverse matrix of a diagonal block matrix after decomposition using incomplete triangular decomposition is because the diagonal block matrix is singular and an accurate inverse matrix cannot be obtained, and the inverse matrix used for preprocessing of the conjugate gradient group This is because the complete triangular decomposition itself is already an approximate decomposition, and it is only necessary to exhibit the effect of accelerating the convergence of the solution method for conjugate gradient groups.

〔Example〕

The incompressible fluid is standardized by the Re number (Reynolds number), and the velocity in each axial direction is U in the Cartesian coordinate system.

When written as v and W, if physical strength does not work, the Navier-Stokes equation and the equation of continuity can be expressed as follows. Ko,
::1'P is pressure, VV= (U, V, W), and the continuity equation is equation (4).

Jinju B+B=O(4) ax ay az An embodiment of the present invention will be described below with reference to FIG. 1 for this equation. FIG. 1 is a PAD diagram showing a batch solution method using a pretreated block conjugate gradient group for an incompressible viscous fluid. 1
05 indicates that flow analysis is performed until time tmax. 1
01 indicates that equations (1) to (4) above are simultaneously discretized. As a result, a block matrix A of 4×4 units is created.

102 and 1σ3 are iterative solutions for finding the solution X of simultaneous linear equations Ax=b with matrix A as a coefficient.

At 102, the matrix A is subjected to incomplete block decomposition to prepare for iterative calculation. The inverse matrix of the diagonal block after block incomplete triangular decomposition is approximated as an inverse matrix using the procedure shown in FIG. 103
Now, using this incomplete triangular decomposition in preprocessing, we can calculate U, V, W at the relevant time by iterative solution of the conjugate gradient group for block asymmetry.
, P at the same time. Here, the BCG method (Bi Conjugate Gradie method) is used as the conjugate gradient method for asymmetry.
nt Method) and CGS method (Conjugate
Gradient 5 squared Method)
etc. At 104, the time is advanced by Δt (time step width).

FIG. 2 shows the form of a matrix created by simultaneously discretizing the Navier-Stokes equations and continuity equations shown in equations (1) to (4) using the three-dimensional difference method. 201 is a diagonal block matrix, and 202 and 203 are block matrices generated in relation to the preceding and following nodes. These occur at a position offset from the diagonal by one block. 204 and 205
is a block matrix generated in relation to nodes on one side away. This block matrix is shifted downward and upward by ml from the diagonal in units of blocks. 206 and 207
is a block matrix generated in relation to nodes on one plane apart.

This block matrix is m2 in blocks from the diagonal.
It is shifted downward and upward. Here, m2 is a multiple of ml. 208 indicates that the mark appearing in the matrix is a 4×4 block matrix. This block matrix is 1,2
, 3.4 columns correspond to variables U, V, W and P, respectively,
Lines 1, 2, 3.4 are (1), (2), (3
), which corresponds to the equation (4).

Figure 3 shows the concrete form of the matrix created by simultaneously discretizing the two-dimensional Navier-Stokes equation and the continuity equation using the two-dimensional difference method, and the analysis target area with numbered nodes. . Reference numeral 301 denotes an analysis target area, in which water flows in from the left side at a constant speed and flows from the right side. 30
2 indicates a node where the flow velocities U and V are fixed, and 303 indicates a node where the flow velocity U is fixed.
, ■ are interior points that are unknowns. 304 indicates the node number of the upper region. Due to space limitations, only a portion of the matrix is shown.

In 305, since the flow velocity U and ■ are fixed, the diagonal elements of the corresponding diagonal block are set to 1. Also, for 306, all corresponding row elements other than diagonal elements are set to O for the same reason. Here, * enters a value other than O, and blanks indicate all zeros. 307 is a value generated by discretizing the equation of continuity of diagonal blocks corresponding to interior points. This part is all zero due to the nature of the equation. This is the biggest reason why block incomplete triangular decomposition is difficult.

FIG. 4 is a PAD diagram for performing incomplete triangular decomposition of the simultaneously discretized block matrix A. 406 indicates that i is increased by 1 from 1 until the number n of diagonal blocks is reached, and the processes 401 to 405 are repeated. Step 401 determines whether the row norm of the j-th row of the i-th diagonal block is zero. When this row norm is not zero, σLJ sets the row norm of the j-th row of the i-th diagonal block. On the other hand, this row norm becomes zero when the continuity equation is discretized, so σ1J sets the row norm of the entire matrix A corresponding to the j-th row of the i-th diagonal block. Furthermore, the parameters ε and α are obtained from numerical experiments and are set as follows.

In many cases, it is sufficient to set ε to 0.1 and α=1.0. 404 is a part that calculates a diagonal block matrix as a result of incomplete block decomposition. Here, in order to save storage capacity, the incompletely decomposed lower triangular matrix and upper triangular matrix are the same as the original matrix A except for the diagonal blocks. As a result, only the diagonal blocks need to be calculated. The calculation formula for Wi shown here is intended for simultaneous discretization in the form of a matrix shown in FIG. Step 405 calculates an approximate inverse matrix Di of the just calculated Wi. Since Wl is often a singular matrix, the pivot value PiJ is obtained by exchanging the axes, and Pij/σ1. If <ε, Pij is replaced by σ, xεX (sign of Pij) to eliminate singularity.

Figure 5 shows a three-dimensional incompressible viscous fluid using the simultaneous I LUCGS method (I ncol) I
plete L U decompot, 1ti
on Conjugat, eGradient 5q
The procedure for solving the problem using the UARED Method is shown below.

501 indicates that flow analysis is to be performed until time tmax.

502 is the Navier-Stokes equation and continuity equation corresponding to equations (1) to (4), but the difference is that the nonlinearity of the Navier-Stokes equation is removed by using VVO and the previous value instead of vv. The point is that

503 is a simultaneous discretization of the equation of 502, and 4×4 at each node.
The process of creating a block matrix of is shown below. In 504, the fourth
Using the incomplete triangular decomposition shown in the figure, create a system of simultaneous linear equations whose coefficients are the matrices created in step 503.
Solve using the S method to obtain U, V, W, and P(7) values simultaneously.

In step 505, time t is advanced by △t and the newly obtained value of vv is assigned to vv0. Here, vv is (U, V,
W), and VVO is a vector (UOVO, WO).

FIG. 6 shows the calculation procedure of the simultaneous ILUCGS method. Uppercase letters indicate matrices, lowercase letters indicate vectors, and Greek letters indicate scalar values. Also, (q, s) indicates the calculation of the inner product of vectors q and S.

601 shows the block incomplete triangular decomposition shown in FIG. 4, and shows the result as a matrix product in the form of LDU.

602 is the forward direction of the LDU for the residual vector b-Ax.
Apply backward substitution and calculate (LDU) (b-Ax). 603 is preparation for iterative calculation. 604 is a calculation procedure for iterative calculation. Vector X gradually approaches the solution vector through repeated calculations. Since the vector r is used to calculate the residual, iterative calculations can be aborted using the norm of r.

Although this calculation example was shown for a matrix discretized by the finite difference method in an orthogonal coordinate system, it can also be applied to the boundary fit method (the finite difference method in a general coordinate system) and the finite element method. In the boundary fit method, only the positions of non-zero blocks are different.

In the finite element method, the positions of the non-zero glocks are different and the velocities U, V.

Since W and pressure P generally use elements of different orders (P is one order lower), dummy elements are created for the missing nodes. This calculation example can also be applied to staggered meshes.

In this case, for missing variables, dummy elements are used to create a block matrix with block lengths of the same length. For incomplete triangular decomposition, which is the most problematic problem, the measures taken using the finite difference method can be used as is.

〔Effect of the invention〕

According to the present invention, in ultra-high-level languages that can describe partial differential equations almost as they are, such as Hitachi's DEQSOL (simulation language for partial differential equations), the Navier-Stokes equations and continuity equations can be written as they are. The ease of use in fluid analysis is greatly improved compared to the conventional SMAC method.

Furthermore, according to the present invention, since the calculation is performed while establishing the equation of continuity for each time, the time step width is S M
It can be 10 to 100 times larger than the AC method, enabling high-speed calculation, and changes in the time step width can be absorbed by changes in the number of iterations in the simultaneous ILUCGS method, making it possible to calculate in 1 hour steps. Increased stability relative to width.

[Brief explanation of drawings]

FIG. 1 is a PAD diagram showing an example of the processing procedure of the present invention;
Figure 2 shows an example of the overall form of a matrix when the three-dimensional Navier-Stokes and continuity equations are simultaneously discretized, and Figure 3 shows the partial form of a two-dimensional matrix. A diagram showing an example,
FIG. 4 is a diagram showing an example of a block incomplete triangular decomposition method, which is a key point of the present invention. FIG.
Figure 6 shows an example of applying the simultaneous ILUCGS method to the present invention using
The figure shows an example of the calculation procedure of the simultaneous I LUCGS method. 101... Creation of block matrix by simultaneous discretization, 1
02 Block incomplete triangular decomposition, 103 Iterative solution of simultaneous linear equations with block matrix A as a coefficient, +04
・Advance the time by the time step width △t, 105t m
Unsteady calculation for ax, 201...diagonal block matrix,
203 to 207...off-diagonal block matrix, 208...
- Shape of one block constituting the block matrix, 301 - Area to be analyzed, 302 - Nodes with fixed values, 303 -
・Internal nodes, 304 Node number, 305・Diagonal elements of nodes whose values are fixed, 306・Off-diagonal elements for nodes in 305, 307・Diagonal blocks generated by discretizing the continuity equation of internal nodes Value, 401... Zero determination of row norm of 1 diagonal block, 402... Setting of row norm of diagonal block, 403... Processing when the row norm of diagonal block is zero, 404... Invalid Diagonal block calculation method using complete triangular decomposition, 405...Calculation of approximate inverse matrix of the diagonal block calculated in 404, 406.Repetition for calculation of all blocks, 501...j Unsteady calculation up to 1QaX, 502...Navier-Stokes equation and continuity equation standardized by Re number, 503...Simultaneous discretization and the form of one block at that time, 504...Simultaneous linear equation Ax=
Iterative solution of b, 505... Update of time t and update of flow velocity vector, incomplete block decomposition of 60 block matrix A,
602 Calculation of residual vector, 603... Preparation for iterative calculation, 604, ILUCGS method (Conjugate GradientSqua with incomplete triangular decomposition)
red Method) An internal calculation method for iterative calculations. The first episode of the first episode of the first episode / The first episode of the first episode

Claims (1)

[Claims] 1. In the numerical simulation method of fluid phenomena, the Navier-Stokes equations and continuity equations, which are the governing equations, are written as they are, and the flow velocity and pressure in three directions are discretized using the finite element method and the finite element method in Cartesian and curved coordinate systems. A step of creating simultaneous linear equations that are discretized as a set of blocks by keeping the block length of the same length including fixed points, and using a calculation method of a conjugate gradient method system with preprocessing for the simultaneous linear equations. The Navier-Stokes equation is characterized by a step of replacing the inverse matrix of the diagonal block that occurs during incomplete triangular decomposition with its approximate inverse matrix even if there is an all-zero row in the diagonal block in order to solve it numerically. Batch simulation method.