KR20090069604A - Method for 2d finite element numerical analysis of fluid flow with cdg method - Google Patents

Abstract

A 2D finite element numerical analysis method of fluid flow using a CDG method is provided to exactly simulate fluid flow of various types generated under various conditions. Basic data including topographical coordinates, an initial value, an illumination coefficient, a boundary condition, a boundary value, a time unit, and a calculation time interval of nodal points configuring an element are inputted(S10). A direction of a normal vector of a waterway sidewall required when an action boundary condition of the waterway sidewall is set is calculated(S20). A matrix induced by applying a finite element scheme based on each element is configured(S40). A matrix list configuring a column and a row of an element matrix is written(S50).

Description

Method for 2D Finite Element Numerical Analysis of Fluid Flow with CDG Method

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to the field of hydraulic engineering, and more particularly, to a numerical analysis method for analyzing a flow of a fluid such as a river based on a two-dimensional finite element.

Digitization method that can be applied for numerical analysis of fluid flow is largely finite difference method, finite volume method, finite element method.

In the finite-difference and finite volume methods, artificial viscous terms have been added to the governing equations to solve the problem of numerical vibration of discontinuous fluid flow, and numerous numerical techniques have been developed to suppress the numerical instability.

Nevertheless, the finite difference method or finite volume method does not have smooth lattice structure when applied to natural rivers, and there are disadvantages such as the need to test variables according to the simulation object.

Attempts have been made to apply coordinate transformations to complex shaped channels, but there are disadvantages that the discretization equations become complicated and time consuming.

The finite element method has an advantage that a triangular element or a square element can be used to form a smoother grid than the finite difference method or the finite volume method.

Moreover, when comparing simulation results with the same object, it is known to be superior in numerical accuracy and calculation time compared with the difference method and the volume method.

When the traditional finite element method is applied using the two-dimensional flow governing equations, the problems of numerical vibration can occur like the finite difference method or the finite volume method.

The SMS-RMA2 model of the US Corps of Engineers, a two-dimensional finite element flow model applied with the elementary and traditional finite element technique, which is currently used in the field of hydraulic engineering in Korea, causes numerical vibration problems when applied to point currents, rapid currents, and currents. Yields a numerical solution that is not obvious.

Thus, there is a need for a new two-dimensional finite element numerical method that can simulate complex and discontinuous fluid flows.

The present invention has been made to solve the above problems, a discontinuous flow that occurs in the river, that is, a discontinuous fluid flow that can accurately simulate the various types of fluid flow that occurs in a variety of conditions that were not possible to simulate so far Its purpose is to provide a two-dimensional finite element numerical analysis method of.

In order to achieve the object described above, the present invention inputs basic data including topographic coordinates, initial values, roughness coefficients, boundary conditions, boundary values, time units, and calculation time intervals of a plurality of nodes constituting each element. Receiving basic data input step; A normal direction calculation step of calculating a direction of a normal vector of the channel side wall required when setting an active boundary condition of the channel side wall; Β calculation step of discretizing partial derivatives of time and updating the solution calculated at the previous time level when applying the negative solution method; An element matrix constructing step of constructing a matrix derived by applying a finite element method based on each element; It proposes a two-dimensional finite element numerical analysis method of the fluid flow, including; matrix list creation step of preparing a matrix list constituting the rows and columns of the element matrix.

The normal direction calculation step is preferably calculated by the equation (1).

The β calculation step is preferably calculated by the equations (2) to (4).

(Where n is the known time level, n + 1 is the unknown time level, θ is the parameter representing implicitness, U is the depth and flow rate component, and Δt is the time level)

When n = 0, it is preferable to apply β ¹ = αU ⁰ .

The step of constructing the element matrix includes calculating interpolation of variables such as spatial coordinates, depth, and flow rate and interpolation of the spatial derivatives of depth, flow rate, and shape function with respect to Gauss points represented by normalized coordinate axes ξ and η. desirable.

Preferably, the element matrix constructing step includes calculating Δx and Δy by Equations 8 to 9.

Preferably, the element matrix constructing step includes calculating an up-weighted matrix by Equations 10 to 13.

The element matrix constructing step may include calculating E ₁ , E ₂ , and E ₃ by Equations 14 to 16.

Preferably, the element matrix constructing step includes calculating f _i by Equation 7.

The element matrix constructing step may include constructing a matrix of Equation 6 using the fi.

The step of constructing the element matrix includes the step of performing the coordinate transformation by Equations 17 to 18 at the node designated as the activity boundary after performing all the steps for all Gaussian points in the one element. desirable.

In order to use the element criterion rather than the nodal criterion in the process of collecting the element matrix into the entire matrix, the final position of the element configuration node is evaluated, a matrix list is prepared, and then the total matrix according to the matrix list is collected. By creating pivot rows, it is desirable to reduce computation time and reduce storage such as memory.

Evaluating arrangement information of nodes constituting each element; A pivot row creating step of evaluating and storing pivot rows after each of the element matrices is collected as a whole matrix; Depth and flow rate calculation step of calculating the depth and flow rate component through the post-substitution process after completing the pivot row for all the above elements; It is preferable to further include.

After the step of calculating the depth and flow rate component, compared to the previous step, the depth and flow rate component at a new time level when the convergence range is calculated and ends when the given time range is reached; preferably includes.

As another means for achieving the above object, the present invention provides a basic data including topographic coordinates, initial values, roughness coefficients, boundary conditions, boundary values, time units, calculation time intervals of a plurality of nodes constituting each element. Receiving input and evaluating arrangement information of nodes constituting each element; A normal direction calculation step of calculating a direction of a normal vector of the channel side wall required when setting an active boundary condition of the channel side wall; Β calculation step of discretizing partial derivatives of time and updating the solution calculated at the previous time level when applying the negative solution method; An element matrix constructing step of constructing a matrix derived by applying a finite element method based on each element; A matrix list generating step of preparing a matrix list constituting rows and columns of the element matrix; A pivot row creating step of evaluating and storing pivot rows after each of the element matrices is collected as a whole matrix; Depth and flow rate calculation step of calculating the depth and flow rate component through post-substitution process after completing the pivot row for all elements; Comparing with the previous step, we propose a two-dimensional finite element numerical analysis method of the fluid flow, including; the step of calculating the depth and flow rate components at the new time level and ending when the given time range is reached.

The normal direction calculation step is preferably calculated by the equation (1).

The β calculation step is preferably calculated by the equations (2) to (4).

When n = 0, it is preferable to apply β ¹ = αU ⁰ .

The step of constructing the element matrix includes calculating interpolation of variables such as spatial coordinates, depth, and flow rate and interpolation of depth, flow rate, and shape function with respect to Gauss points represented by normalized coordinate axes ξ and η; Calculating Δx and Δy by Equations 8 to 9; Calculating an upweighting matrix according to Equations 10 to 13; Calculating E ₁ , E ₂ , E ₃ by Equations 14 to 16; Calculating f _i by equation (7); Constructing a matrix of Equation 6 using the fi; After performing all the steps for all the Gaussian points in the one element, and performing the coordinate transformation by the equation (17) to the node designated as the activity boundary; preferably.

The present invention provides a method for two-dimensional finite element numerical analysis of a fluid flow that can stably and accurately simulate the flow of a stream in a stream as well as in upstream and quadrature conditions as well as in upstream and quadrature conditions. For maintenance and maintenance of rivers, evaluation of stability when designing river structures, prediction of results according to river restoration, results of flood calculations, comparison of indoor experimental channel results, etc. do.

Hereinafter, embodiments of the present invention will be described in detail with reference to the accompanying drawings.

1 is a flowchart illustrating a two-dimensional finite element numerical analysis method flow of a fluid flow according to the present invention, Figure 2 is a flow chart illustrating a process of creating an element matrix using the CDG technique. 3 to 6 are views illustrating a matrix list in a simple finite element lattice construction and numerical analysis process for better understanding, and the present invention will be described with reference to the following.

Basic data input and node array information evaluation step (S10) includes topographical data such as topographic coordinates ( x , y , z ) of nodes and information about nodes constituting one element, and initial values, roughness coefficients, boundary conditions and boundaries. This is a process of receiving basic values such as values, time units, and calculation time intervals, and performing basic calculations on the nodes that make up elements.

3 is an embodiment of an element grid composed of a total of 12 nodes and 8 elements.

The upstream flow boundary conditions are set at nodes 1, 2, and 3, and the downstream water level or depth boundary conditions are set at nodes 10, 11, and 12.

One element can have a rectangular or triangular shape, and depending on whether you use a one-dimensional or two-dimensional grid, you can use four or eight nodes for the rectangular element and three or six nodes for the triangle element.

For element-based matrix calculations, it is necessary to find and mark the last node position that appears among the element construct nodes.

4 is a diagram illustrating negative display at a position where each node appears last after receiving element and node arrangement information.

The normal direction calculation step (S20) is a step for calculating the normal direction necessary for ignoring the flow component in the normal vector direction when the active boundary condition is set on the sidewall of the simulation area, and required for setting the upstream inflow flow boundary condition.

4, 6, 7, and 9 nodes in Fig. 3 correspond to the channel sidewalls.

An active boundary condition in which no fluid inflow from the outside along the channel sidewall or an outflow of fluid from the interior of the channel occurs, but only a flow in a direction parallel to the sidewall, or an inactive boundary condition in which no fluid flow exists on the sidewall. Can be set.

When the inactive boundary condition is set, the flow component at the side wall is excluded from the numerical calculation. When the active boundary condition is set, only the flow component in the direction of the tangent vector parallel to the side wall should be calculated.

In addition, even when there is an inflow of the flow upstream, as in nodes 1, 2, and 3 of FIG.

Equation 1 below calculates the angle of the normal to the outward direction of the boundary at the node corresponding to the boundary.

Here, Φ is the normal angle in the counterclockwise direction and the positive direction of the x-axis, N _k is the shape function of the node k in one element, A is the area integral in the element. However, for nodes that consist of only one element, nodes that make up edges, and nodes that can produce inadequate results, you must set the appropriate normal angle.

β calculation step (S30) is a new β value based on the calculation result of the previous stage time level from the derivative term at the unknown time level represented by Equations 2 to 4 below, which is derived by discretizing the partial derivative of time. The process of calculating it.

Where n is a known time level, n + 1 represents an unknown time level, θ is a parameter representing implicitness, U represents depth and flow rate components, and Δt represents time level .

However, when n = 0, β ¹ = αU ⁰ is applied. Using Equation 3, it is the process of updating the β values for the x and y components of the depth and flow rate before entering a new time level.

The element matrix constructing step (S40) is formed by applying the Newton-Raphson method represented by Equation 5 to solve the nonlinear equation system derived by applying the weighted residual method to the depth-integrated two-dimensional flow equation. The process of constructing a matrix by elements.

The matrix generated at the i th node constituting one element by Equation 5 is expressed as Equation 6.

Here, the matrix on the left side represents the Jacobian matrix, and each component on the right side is represented by Equation (7).

Where _i represents the integral in the entire simulation region, N _i is a shape function, ω is an upwinding coefficient, and Δx and Δy are represented by Equations 8 and 9, respectively.

Here, ξ and η are coordinate axes in normalized local coordinates, each having a value of [-1, 1] in a rectangular element and a value of [0, 1] in a triangular element.

In Equation 7, W _x and W _y are up-weighted matrices in the x and y directions, and are represented by Equations 10 to 13 below, and an up-weighted matrix calculation technique represented by Equations 10 to 11 is called CDG. It is called.

Where g is gravity acceleration, h is depth, and p and q are the x and y flow rate components, respectively. In Equations 10 to 11, the inverse of the square root may be calculated using Cayley-Hamilton's theorem.

In Equation 7, the vector E = [E ₁ E ₂ E ₃ ] ^T is a flow governing equation in which Equation 2 is applied to a partial derivative of time, and is represented by Equations 14 to 16 below.

(Where ν _t is the turbulent kinematic coefficient and z _b is the riverbed)

The element matrix constructing step S40 will now be described with reference to FIG. 2.

For the Gaussian points represented by the normalized coordinate axes ξ and η, interpolation of variables such as spatial coordinates, depth, and flow rate, and interpolation about the spatial derivative of the depth, flow rate, and shape function is calculated (S41). Δx and Δy expressed are calculated (S42).

Calculate the up-weighted matrix represented by Equations 10 to 13 (S43), calculate E ₁ , E ₂ , E ₃ represented by Equations 14 to 16 (S44), and Equation 6 represented by Equation 7 The right side f _i is calculated (S45), and the left-side Jacobian matrix of Equation 6 is formed (S46).

After performing all the above steps for all the Gaussian points in one element, the coordinate transformation represented by Equations 17 to 18 below is performed at the node designated as the activity boundary (S47).

The matrix list preparing step S50 is a process of creating a list of rows and columns of the matrix created through the element matrix constructing step S40.

5 is an example of creating a matrix list for the element of FIG.

The variables that need to be calculated for each node are the depth (h), the x, y components (p, q) of the flow rate, or the depth, the normal and the tangential components (q _n , q _s ) of the flow rate, so the matrix per node The number of lists is also increased by a factor of three, arranged in the order h, p, q or h, q _n , q _s .

Pivot row creation step (S60) collects the entire matrix using the element matrix calculated in the element matrix construction step (S40) and the matrix list created in the matrix list creation step (S50), and then creates a pivot row through the pivot process. That's the process.

In the process of collecting each element matrix as an entire matrix, the rows and columns occupied by the variables of the nodes given as boundary conditions must be erased. The pivot row is created by setting the largest of the negative component values in the matrix list to pivot. After that, remove the pivot rows and columns from the entire matrix.

6 is an example illustrating a process of combining the elements of FIG. 3 into an entire matrix.

When factor 1 is combined into the entire matrix, the matrix list represented by -2 and -3 corresponds to the variable specified as the boundary condition at node 1, so the corresponding rows and columns are deleted from the entire matrix.

Set the components of the matrix corresponding to matrix list -1 to pivot, save the pivot rows separately, and then delete them from the entire matrix.

When element 2 is combined into an entire matrix, the matrix list corresponding to nodes 2 and 5 already exists in the entire matrix, so all you have to do is add an element matrix at that position.

Since -5 and -6 variables correspond to the boundary conditions, they are erased from the entire matrix, and the -4 component is set to pivot to store the pivot rows separately and then to delete them from the entire matrix.

The same process is performed for all element collections, and eventually pivots on all variables that are not set as boundary conditions.

In this way, the processes from steps S40 to S60 are repeated for all elements to obtain a pivot row.

The depth and flow rate component calculation step S70 is a process of calculating the x and y components of the depth and flow rate at each node using the pivot rows obtained from the pivot row creation step S60.

Since the pivot rows are the upper triangular matrices that are generated by using the Gaussian elimination method after collecting the element matrices, the solution can be obtained through the backward substitution method from the last pivot row.

If the maximum value of the difference between the newly calculated solution and the solution in the previous step is larger than the convergence value given as the input value, the process from steps S40 to S70 is repeated to calculate a new solution.

If the difference between the newly calculated solution and the solution in the previous step is smaller than the convergence value, the calculation at the corresponding time is terminated, the process moves to the next time level and the steps S30 to S70 are repeated, and when the final time is reached, The calculation is finished (S80).

In the above description of the present invention, the two-dimensional finite element numerical analysis method of the fluid flow using the CDG technique has been described in detail and specifically, but the present invention can be variously modified and changed by those skilled in the art. Should be construed as falling within the protection scope of the present invention.

1 to 6 are views for explaining an embodiment of the present invention,

1 is a flow diagram of a two-dimensional finite element numerical method procedure of a fluid flow.

2 is a flowchart illustrating a process of creating an element matrix using a CDG technique.

3 to 6 are diagrams showing matrix lists in a finite element lattice construction and numerical analysis process.

Claims (19)

A basic data input step of receiving basic data including topographic coordinates, initial values, roughness coefficients, boundary conditions, boundary values, time units, and calculation time intervals of a plurality of nodes constituting each element; A normal direction calculation step of calculating a direction of a normal vector of the channel side wall required when setting an active boundary condition of the channel side wall; Β calculation step of discretizing partial derivatives of time and updating the solution calculated at the previous time level when applying the negative solution method; An element matrix constructing step of constructing a matrix derived by applying a finite element method based on each element; A matrix list creation step of creating a matrix list constituting the rows and columns of the element matrix; 2D finite element numerical analysis method of a fluid flow comprising a. The method of claim 1, The normal direction calculation step is a two-dimensional finite element numerical method of fluid flow, characterized in that calculated by Equation 1. The method of claim 1, The β calculation step A two-dimensional finite element numerical analysis method of a fluid flow, characterized in that calculated by equations (2) to (4). The method of claim 3, 2 n finite element numerical analysis method of fluid flow, characterized in that when n = 0, β ¹ = αU ⁰ . The method of claim 1, The element matrix construction step Calculating interpolation of variables such as spatial coordinates, depth, and flow rate and interpolation of depth, flow rate, and shape function with respect to Gauss points represented by normalized coordinate axes ξ and η; A two-dimensional finite element numerical method of fluid flow comprising a. The method of claim 5, The element matrix construction step Computing Δx, Δy by the equations (8) to (9); two-dimensional finite element numerical analysis method of a fluid flow comprising a. The method of claim 6, The element matrix construction step Computing the up-weighted matrix by the equation (10 to 13); two-dimensional finite element numerical analysis method of the fluid flow comprising a. The method of claim 7, wherein The element matrix construction step Computing E ₁ , E ₂ , E ₃ by Equation 14 to 16; a two-dimensional finite element numerical analysis method of a fluid flow comprising a. The method of claim 8, The element matrix construction step Computing f _i by the equation (7); two-dimensional finite element numerical analysis method of a fluid flow comprising a. The method of claim 9, The element matrix construction step Comprising the step of constructing the matrix of the equation (6) using the fi; 2D finite element numerical analysis method of the fluid flow comprising a. The method of claim 10, The element matrix construction step After performing all the steps for all Gaussian points in the one element, performing coordinate transformation by Equations 17 to 18 at the node designated as an activity boundary; A two-dimensional finite element numerical method of fluid flow comprising a. The method of claim 1, In order to use the element criterion rather than the nodal criterion in the process of collecting the element matrix into the entire matrix, the final position of the element configuration node is evaluated, a matrix list is prepared, and then the total matrix according to the matrix list is collected. A two-dimensional finite element numerical analysis method of a fluid flow characterized by reducing a calculation time and a storage space such as a memory by creating a pivot row. The method of claim 12, Evaluating arrangement information of nodes constituting each element; A pivot row creating step of evaluating and storing pivot rows after each of the element matrices is collected as a whole matrix; Depth and flow rate calculation step of calculating the depth and flow rate component through post-substitution process after completing the pivot row for all the elements; The two-dimensional finite element numerical analysis method of the fluid flow further comprising. The method of claim 13, After the depth and flow rate calculation step, Comparing with the previous step in the convergence range, calculating the depth and flow rate components at a new time level and ending when the given time range; comprising a two-dimensional finite element numerical analysis method of a fluid flow comprising a. It receives basic data including topographic coordinates, initial value, roughness coefficient, boundary condition, boundary value, time unit, and calculation time interval of a plurality of nodes constituting each element, and obtains array information of nodes constituting each element. Evaluating; A normal direction calculation step of calculating a direction of a normal vector of the channel side wall required when setting an active boundary condition of the channel side wall; Β calculation step of discretizing partial derivatives of time and updating the solution calculated at the previous time level when applying the negative solution method; An element matrix constructing step of constructing a matrix derived by applying a finite element method based on each element; A matrix list generating step of preparing a matrix list constituting rows and columns of the element matrix; A pivot row creating step of evaluating and storing pivot rows after each of the element matrices is collected as a whole matrix; Depth and flow rate calculation step of calculating the depth and flow rate component through post-substitution process after completing the pivot row for all elements; Comparing with the previous step in the convergence range, calculating the depth and flow rate component at the new time level and ending when the given time range; comprising a two-dimensional finite element numerical analysis method of the fluid flow. The method of claim 15, The normal direction calculation step is a two-dimensional finite element numerical method of fluid flow, characterized in that calculated by Equation 1. The method of claim 15, The β calculation step A two-dimensional finite element numerical analysis method of a fluid flow, characterized in that calculated by equations (2) to (4). The method of claim 17, 2 n finite element numerical analysis method of fluid flow, characterized in that when n = 0, β ¹ = αU ⁰ . The method of claim 15, The element matrix construction step Calculating interpolation of variables such as spatial coordinates, depth, and flow rate and interpolation of depth, flow rate, and shape function with respect to Gauss points represented by normalized coordinate axes ξ and η; Calculating Δx and Δy by Equations 8 to 9; Calculating an upweighting matrix according to Equations 10 to 13; Calculating E ₁ , E ₂ , E ₃ by Equations 14 to 16; Calculating f _i by equation (7); Constructing a matrix of Equation 6 using the fi; After performing all the steps for all Gaussian points in the one element, performing coordinate transformation by Equations 17 to 18 at the node designated as an activity boundary; A two-dimensional finite element numerical method of fluid flow comprising a.

Images