JP2002328913A - Numerical analysis device by discrete nabla operator and non-storage type coefficient matrix - Google Patents

Abstract

PROBLEM TO BE SOLVED: To attain a high-speed calculation while maintaining precision of calculation with a small storage capacity in a numerical analysis device which utilizes a finite volume method or a finite element method of a non-structured grating. SOLUTION: A discrete nabla operator is defined, taking note of the fact that a nabla operator is a vector volume and is a volume not dependent on a coordinate system. A coefficient matrix element described by the discrete nabla operator is disintegrated into the product of a part dependent on a shape of the element and a part not dependent on a shape of the element. Only the part dependent on a shape of the element is stored while the other part is found analytically. Each time a need arises in matrix calculation, the coefficient matrix element is created from this. When this discrete nabla operator is used in the finite volume method or the finite element method of a non-structured grating, a non-storage type coefficient matrix can be created. Therefore, the coefficient matrix can be found without having to store all the coefficient matrix elements, thus offering low-capacity storage unit and high speed calculation.

Description

DETAILED DESCRIPTION OF THE INVENTION

[0001]

BACKGROUND OF THE INVENTION 1. Field of the Invention The present invention relates to a numerical analysis apparatus, and more particularly, to a numerical analysis apparatus which performs calculations by a finite volume method or a finite element method of an unstructured grid using a discrete nabla operator and a non-memory type coefficient matrix. About.

[0002]

2. Description of the Related Art A finite element method is one of the methods for solving the boundary value problem of a partial differential equation which appears in problems in fluid dynamics, structural dynamics, and the like by numerical calculation using a computer. The finite element method will be briefly described. For details, see Reference 1.
[Takahiko Tanahashi, Finite Element Analysis of Flows I and II] (199
7), published by Asakura Shoten].

The partial differential equation L (φ) = f (on the region Ω) of the field function φ is converted to the boundary condition φ = φ ₀ (on the boundary Γ ₁ ) (∂φ / ∂n) = q (boundary Γ _This will be explained by taking the case of solving with ₂ ) as an example. Here, n is a normal vector of the boundary. Boundary gamma ₁ is a boundary where the value phi ₀ of the function phi is given (Dirichlet boundary). The boundary Γ ₂ is a function φ
Is a boundary (Neumann boundary) to which the differential value q of is given.
A specific partial differential equation is, for example, L (φ) = ∇ ² φ = f.

An approximate solution of this partial differential equation is defined as φ ^* .
The approximate solution is represented by a linear combination of appropriate basis functions. In the exact solution, the residual with respect to the value on the left side of the partial differential equation r = L (φ ^* ) − L (φ) = L (φ ^* ) − f is assigned a weight function w _j (j = 1, 2,..., n; n is as the integral over the number) number i.e. undetermined coefficients of the weighting function becomes 0, as _{(w j, r) ≡∫w j} rdΩ = 0 ( integration range region Omega), the approximate Determine the solution φ ^* . For any weighting function, r that satisfies this integral equation is 0, so the approximate solution φ ^* matches the exact solution at that time.

Similarly, the residual of the boundary condition is also integrated by multiplying by a weighting function. Then, ∫ {L (φ ^* ) − f} w _j dΩ (the integration range is the area Ω) + ∫ {φ
^* −φ ₀ } w _j dΓ ₁ (the integration range is the boundary Γ ₁ ) −∫ {(∂φ ^* /
{N) -q} w _j dΓ ₂ (the integration range is the boundary Γ ₂ ) = 0. Also in this case, if this integral equation is satisfied for an arbitrary weight function, the approximate solution φ ^* matches the exact solution.

Here, to simplify the integral equation,
Imposing conditions on the weight function, the approximate solution φ ^*Relax the restrictions.
Weight function w_jAnd the boundary Γ₁Imposes the condition 0 above
And ∫ {L (φ^*) -F} w_jdΩ (integration range is area Ω) = ∫
{(∂φ^*/ {N) -q} w_jdΓ (the range of integration is the boundary Γ_Two). This is a basic expression of the finite element method.

[0007] When the basic formula is partially integrated, a weak form is obtained. For example, when L (φ) = ∇ ² φ = f, ∫ {∇ ² φ ^* −f} w _j dΩ (the integration range is the area Ω)
= ∫ {(∂φ ^* / ∂n) -q} w _j dΓ (the integration range is the boundary Γ ₂ ), so using Green's first formula, ∫ (∇φ ^* · ∇w _j ) dΩ ( Integration range is area Ω) + ∫∇ ² φ
^* w _j dΩ (integration range is area Ω) = ∫w _j (∂φ ^* / ∂n) d
Assuming that Γ (integration range is boundary Γ ₂ ), the weak form is ∫ {(∇φ ^* · ∇w _j ) + fw _j } dΩ (integration range is area Ω)
= {Qw _j d} (the integration range is the boundary Γ ₂ ). According to the variation principle, a weak form can be obtained from the condition that the variation of the functional J becomes zero. There are various solutions depending on how the weight function is selected. For example, the finite element method Ga
In the lerkin method, the weight function is a shape function, and the mapping function is the same function. The shape function and the mapping function will be described later.

As an example, a method of solving ∇ ² φ = f by the Galerkin method will be described. The region Ω is divided into tetrahedral elements having n nodes as vertices. Let the shape function of node i be N _i
And Let the value of the approximate solution φ ^* at node i be φ _i (unknown)
And When the value of the weight function w _j at the node j and _{^{w j *, φ * = Σ}} i N i φ i + Σ i N i h i ( sum, i = 1 to n) and w _j = N _j w _j ^* Become. However, h _i is a value that is known at the boundary gamma _1, in a region inside is zero.

[0009] By substituting this to the weak _{form, ∫ {(∇ (Σ i} N i φ i + Σ i N i h i) · ∇ (N j w j *)) ( sum, i = 1~n) + f ^{_{(N j w j *)}}} dΩ ( integration range region _{Ω) = ∫q (N j w} j *) dΓ ( integral range boundary gamma ₂₎ becomes. Since w _j ^* is arbitrary, its coefficient is 0.
_{Thus, ∫ {(∇ (Σ i} N i φ i + Σ i N i h i) · ∇N j) + fN j} dΩ
(Integration range is area Ω, sum is i = 1 to n) = ∫ (qN _j ) d
Γ (the range of integration is the boundary Γ ₂ ). To summarize, Σ _i φ _i ∫ (∇N _i · ∇N _j ) dΩ (the integration range is the area Ω, the sum is i = 1 to n) = ∫ (qN _j ) dΓ (the integration range is the boundary Γ
₂ ) -∫ (fN _j ) dΩ (integration range is area Ω) -Σ _i hi _i ∫
(∇N _i · ∇N _j ) dΩ (the integration range is the area Ω, and the sum is i =
1 to n).

Here, A _ij = ∫ (∇N _i · ∇N _j ) dΩ (integration range is area Ω) B _j = ∫ (qN _j ) dΓ (integration range is boundary Γ ₂ ) −∫ (f
N _j ) dΩ (the integration range is the area Ω) −Σ _i h _i ∫ (∇N _i · ∇
N _j) dΩ (integration range region Omega, sum, when i = 1 to n) _{to, Σ i A ij φ i =} B j ( sum, i = 1 to n) it becomes. That is, (φ _i ) = (A _ij ) ⁻¹ (B _j ). Here, (φ _i ) is a vertical vector having φ _i as the i-th component, (B _j ) is a vertical vector having B _j as the j-th component, (A
_ij ) ^-1 is an inverse matrix of a matrix having A _ij as i-th and j-th components. Thus, the boundary value problem of the partial differential equation can be solved by solving the simultaneous linear equations.

In this case, it is important to calculate the coefficient matrix element A _ij = ∫ (∇N _i · jN _j ) dΩ (the integration range is an area Ω). To calculate this easily, we map the physical space to the computation space and discretize the nabla operator ∇. Solve the boundary value problem of the PDE by describing the coefficient matrix elements with the discrete nabla operator and solving a system of linear equations. The size of the coefficient matrix element is (number of nodes) ×
(Number of elements), so to perform precise calculations,
It requires enormous storage capacity and requires a long calculation time.

By the way, a system equation of a nonlinear dynamical system can be described as a time evolution equation. In the example above,
We have explained PDEs that do not include time variables.
The finite element method can also be applied to partial differential equations including time variables. In the time evolution equation, there are two types of differential operators,
∂ / ∂t and ∇ are included. The time integration of ∂ / ∂t yields various time-marching methods. On the other hand, the spatial integration of the nabla operator ∇ produces various discretization methods.

In the computer simulation of nonlinear system equations, the time progress method and the discretization of space are major issues in solving a large-scale problem of a complex system. In particular, in a finite volume method or a finite element method using an unstructured grid, a discretization procedure of a space is generally complicated.

In the “Numerical Analysis Apparatus, Numerical Analysis Method, and Integrated Circuit” disclosed in Japanese Patent Application Laid-Open No. H10-177591, the present inventors use a discrete nabla operator to realize a simple configuration and a short A numerical analysis device capable of numerical analysis in time was proposed. In this numerical analysis device, based on the element information of each element obtained by dividing the field to be analyzed,
A discrete nabla operator calculating means for calculating a discrete nabla operator obtained by discretizing the nabla operator is provided. Using the calculated discrete nabla operator, the physical quantity is iteratively modified so as to satisfy the equation of the continuity of the field to be analyzed, and the physical quantity that satisfies the equation of the continuity is calculated by the physical quantity calculating means.

[0015]

However, in the conventional finite element method, there is a problem that it takes a long time to calculate the coefficient matrix as it is even if the calculation is performed using the discrete nabla operator. Further, once the coefficient matrix calculated once is stored, there is a problem that an enormous storage capacity is required in a large-scale problem.

The present invention solves the above-mentioned conventional problems and solves the above problem by using a finite element method using a discrete nabla operator.
An object of the present invention is to speed up the calculation while reducing the storage capacity.

[0017]

According to the present invention, there is provided a numerical analysis apparatus for calculating a solution of a boundary value problem of a partial differential equation by a finite volume method or a finite element method of an unstructured grid. , A coefficient matrix base storage means for storing a factor depending on the element shape of the coefficient matrix element represented by the discrete nabla operator in association with the element shape, and an element of the coefficient matrix element represented by the discrete nabla operator A partial coefficient factor generating means for non-storing factors independent of shape for each type to generate a partial coefficient factor, and a coefficient matrix base and a partial coefficient factor obtained by referring to a coefficient matrix base storing means each time a matrix calculation is required And a calculating means for generating a coefficient matrix element from the partial coefficient factor generated by the generating means. With this configuration, it is possible to speed up the calculation while reducing the storage capacity.

That is, noting that the nabla operator ∇ is a vector quantity and a quantity independent of a coordinate system,
Define the nabla operator. Based on this basic idea, the coefficient matrix element expressed by the discrete nabla operator is decomposed into two parts. One is a part that depends on the element shape, and the other is a part that does not depend on the element shape. Only the part that depends on the element shape is stored, and the other parts are calculated analytically. Each time it is needed inside the computer, a coefficient matrix element is generated from the two parts and stored. In the unstructured grid finite volume or finite element method, this discrete nabla
Using operators is superior in the following points.

(1) When the discrete nabla operator is used, discretization of a space for an arbitrary element can be extremely simplified. (2) The discrete nabla operator is suitable for the finite element method and the finite volume method of an unstructured grid. (3) The discrete nabla operator is extremely easy to turn into an object. (4) A non-memory type coefficient matrix can be created using the discrete nabla operator. Therefore, it contributes to lower capacity and higher speed. (5) The use of the discrete nabla operator facilitates the development of new simulation codes for the next generation. (6) With a special-purpose computer that converts discrete nabla operators into LSI,
Large-scale calculations of complex systems can be performed with low cost and high accuracy. (7) Using the discrete nabla operator can speed up the Poisson Solver.

[0020]

DESCRIPTION OF THE PREFERRED EMBODIMENTS Embodiments of the present invention will be described below in detail with reference to FIGS.

(Embodiment) An embodiment of the present invention stores a factor depending on the element shape of the discrete nabla operator when calculating by the finite volume method or the finite element method of an unstructured grid, A numerical analysis device that analytically calculates factors that do not depend on the element shape of the discrete nabla operator, generates a coefficient matrix element each time it is needed, and solves the coefficient matrix without storing all the coefficient matrix elements It is.

FIG. 1 is an explanatory diagram of the gradient, divergence, and rotation of a one-dimensional single element used in the numerical analysis device according to the embodiment of the present invention. FIG. 2 is a schematic diagram of the 2
It is explanatory drawing of gradient, divergence, and rotation of a dimensional single element. FIG. 3 is an explanatory diagram of gradient, divergence, and rotation of a three-dimensional single element used in the numerical analysis device.

FIG. 4 shows the gradient, divergence, and divergence of one-dimensional symmetric elements used in the numerical analysis device according to the embodiment of the present invention.
It is explanatory drawing of rotation. FIG. 5 is an explanatory diagram of gradient, divergence, and rotation of a two-dimensional symmetric element used in the numerical analysis device. FIG. 6 is an explanatory diagram of the gradient, divergence, and rotation of a three-dimensional symmetric element used in the numerical analysis device.

FIG. 7 is a conceptual diagram of a Dell chip of a numerical analysis device according to an embodiment of the present invention. In FIG.
The element node coordinate generating means 71 is means for outputting the node coordinates of the element to be calculated. Del chip 72 discrete nabla
LS to calculate element volume and contravariant basis vector by operator
It is an I circuit. The data use unit 73 is a calculation unit that uses the output data g _e ⁱ and Ω _e .

FIG. 8 is a flowchart showing a processing procedure of the numerical analysis device according to the embodiment of the present invention. FIG. 9 is a functional block diagram of the numerical analysis device according to the embodiment of the present invention. In FIG. 9, the input device 1 is a keyboard or the like for inputting data such as parameters. The parameter boundary condition storage device 2 is a memory for storing data of the parameters of the partial differential equations and the boundary conditions. The node position vector storage device 3 is a memory for storing coordinate values of all nodes. The known vector calculation means 4
This is a means for calculating a known vector from parameters and boundary conditions. The known vector storage device 5 is a memory for storing known vectors. The partial differential equation information storage device 6 is a memory for storing information on a computable partial differential equation. The display device 7 is a CRT display or the like for displaying the analysis result. The coefficient matrix basis calculation means 8 is a means for calculating an invariant basis vector, an element volume, and the like from the nodal coordinates. The coefficient matrix base storage device 9 is a memory for storing coefficient matrix bases. The shape function storage device 10 is a memory for storing data related to a shape function. The partial coefficient factor generating means 11 is a means for generating a partial coefficient factor from the data of the shape function calculated analytically. The coefficient matrix element operation means 12 is means for calculating a coefficient matrix element from a coefficient matrix base and a partial coefficient factor.
The matrix calculation means 13 is a means for obtaining a function vector from a coefficient matrix element and a known vector. The matrix element storage work area 14 is a memory for temporarily storing coefficient matrix elements.

The operation of the numerical analysis apparatus according to the embodiment of the present invention configured as described above will be described. First, the definition of the discrete nabla operator will be described. Consider how the nabla operator ∇≡i (∂ / ∂x) + j (∂ / ∂y) + k (∂ / ∂z) defined in the continuous space should be defined in the discrete space.
Here, i, j, and k are unit vectors in the axial direction of x, y, and z, respectively.

Let φ be an arbitrary tensor quantity. The scalar is the zeroth-order tensor, the vector is the first-order tensor, and the stress tensor is the second-order tensor. For example, φ is the pressure p
(A scalar or zero-order tensor), a velocity v (vector or first-order tensor), or a stress tensor T (tensor or second-order tensor).

The star (*) product represents gradient, divergence, and rotation. For example, ∇p (gradient), ∇ · v
(Divergence (dot product)), ∇ × v (rotation (cross product)), ∇v
(Tensor product (open product)), ∇ · T (divergence of tensor), ∇ × T (rotation of tensor), ∇T (gradient of tensor). These operations are all represented by ∇ * φ.

The element mean value of the finite element Ω _e (r) of ∇ * φ is expressed as <∇ * φ> _e ≡ (1 / Ω _e (r)) ∫ {∇ * φdΩ (r)} (the integration range is Ω _e ). Here, Ω _e (r) means the volume of a tetrahedron / 5 / pentahedron / 6 element in the case of three dimensions, and the area of a triangle / 4 element in the case of two dimensions. Means the length of the line element. Ω _e is an arbitrary closed region, and if its elements are made sufficiently small, naturally, ∇ * φ = lim <∇ * φ> _e (the limit is Ω _e → 0).

Element Ω_eIs the element Ω_eNode number of
Is expressed as a, φ = Σ_aN_aφ_a (The sum is a = 1 to m, m is the node of the element
Number). Node a is also a node of an adjacent element
Therefore, the shape function should also be distinguished by adding the element number
However, since it is complicated and there is no risk of being lost, _aWhen
I decide to write. Where N_aIs a shape function. Confused
If there is no danger, omit the symbol Σ and_aφ_a Sometimes it is written. In this case, the subscript a follows the rules of the sum
That. By substituting this φ into the definition formula of the element average, <∇ * φ>_e= Σ_a(1 / Ω_e(r)) ∫ {∇N_adΩ (r)} * φ_a (Integration range is Ω_e) = Σ_a∇_a* Φ_a(The sum is a = 1 to m).

[0031] Here, a discrete nabla operator ∇ _a corresponding to the node a region _{Ω e (r), ∇ a} ≡ (1 / Ω e (r)) ∫ {∇N a dΩ (r)} ( integration The range is defined by Ω _e ). That is, the discrete nabla operator is an element average value of the gradient of the shape functions N _a. This definition is the basic idea. Then, using this definition equation, using element average _{value, <∇p> e = Σ a} ∇ a p a, <∇ · v> e = Σ a ∇ a · v a, <∇ × v> e _{= Σ a ∇ a × v a} , <∇v> e = ∇ a v a, <∇ · T> e = Σ a ∇ a · T a, <∇ × T> e = Σ a ∇ a × T a, _{<∇T> e = Σ a ∇} a T a and the like can be calculated immediately.

[0032] In the same manner, element average value of the viscous _{stress, τ e = <μ [∇v} + (∇v) T]> e = μ e Σ a (∇ a v a + v a ∇ a) and can be displayed. Here, μ is a viscosity coefficient. μ _e is
Element mean value of μ. The remaining problem is how to calculate the discrete nabla operator. By defining the discrete nabla operator in this way, the gradient, divergence, and rotation
Can be described using operators.

Second, a method of calculating the discrete nabla operator will be described. There are two methods of calculating the discrete nabla operator: a field calculation method that calculates by volume integral, and a surface calculation method that calculates by area. Elements distorted in the physical space can be normalized by mapping them to the computation space. Therefore, the calculation of the discrete nabla operator is preferably performed in the calculation space. For example, r = (x,
y, z) = (x ¹ , x ² , x ³ ) space, ξ = (ξ, η, ζ) = (ξ ¹ ,
ξ ² , ξ ³ ) Map to space. The mapping is performed using the shape function N _a as follows: r (ξ) = Σ _a N _a (ξ) r _a (the sum is a = 1 to m). r _a is the position vector of the node a.

The nabla operator ∇ is a vector, and uses the fact that it does not depend on the coordinate system, and finds the nabla operator in the calculation space. That is, the nabla operator, in the physical space and the computational space, in three dimensions, ∇ = ∇x ¹ (∂ / ∂x ¹ ) + ∇x ² (∂ / ∂x ² ) + ∇x ³ (∂ /
∂x ³ ) Physical space ∇ = ∇ξ ¹ (∂ / ∂ξ ¹ ) + ∇ξ ² (∂ / ∂ξ ² ) + ∇ξ ³ (∂ /
∂ξ ³ ) Computation space can be displayed, and both are equal. Here, if the contravariant basis vector g ⁱ ≡∇ξ ⁱ (i = 1, 2, 3) is defined, the nabla operator simply describes ∇ = Σ _i g ⁱ (∂ / ∂ξ ⁱ ) = g ¹ It can be described as (＋ / ∂ξ ¹ ) + g ² (∂ / ∂ξ ² ) + g ³ (∂ / ∂ξ ³ ).

If Jacobian of the variable transformation is represented by J, then dΩ (r) = JdΩ (ξ) J = {∂ (x, y, z) / ∂ (ξ, η, ζ)} holds. Here, if the element mean value is defined as J _e ≡ <J> _e ≡Ω _e (r) / Ω _e (ξ) g _e ⁱ ≡ <g ⁱ > _e , the discrete nabla operator becomes ∇ _a ≡ _{(1 / Ω e (r)} ) ∫ {∇N a dΩ (r)} ( integration range _{Ω e) = (1 / Ω} e (r)) ∫ {∇N a JdΩ (ξ)} ≒ Σ i ( J _e / Ω _e (r)) ∫ {g ⁱ (∂N _a / ∂ξ ⁱ ) dΩ (ξ)} ≒ Σ _i g _e ⁱ (1 / Ω _e (ξ)) ∫ {(∂N _a / ∂ ^{ξ i) dΩ (ξ)}} = Σ i g e i <∂N a / ∂ξ i> e ( sum, i = 1 to 3) can be approximated. Therefore, it can be described ^{_{∇ a = Σ i g e i}} <∂N a / ∂ξ i> e ( body calculus) and. This is called the field calculation method of the discrete nabla operator. That is, the discrete nabla operator can be represented by the product sum of the element average value of the invariant base vector depending on the element shape and the first derivative of the shape function independent of the element shape.

Another calculation method is possible. (∇ *
The element mean of φ) is calculated by the boundary integral using Green's theorem.
It can also be represented. That is, <∇ * φ>_e= (1 / Ω_e(r)) ∫ {n * φdΓ (r)} (integral range
The enclosure is Γ_e). Where Γ_e= ∂Ω_eIs Ω_eAt the boundary of
This is the unit normal vector of the surface. This is expressed as the sum of each surface.
In other words, <∇ * φ>_e≒ (1 / Ω_e(r)) Σ_sΔΓ^s(r) * φ_s = Σ_sΔΓ^s(r) * φ_s/ Ω_e(r) (Surface calculation method). △ Γ^sIs the area vector of the microelement. Σ_s
Is the boundary element ∂Ω_eMeans the sum of the divided elements. to s
Σ may be omitted in the sense of the sum of φ _sIs
This is the surface average of φ.

As described above, the discrete nabla operator can be calculated by the field calculation method and the surface calculation method. The field calculation method is simple to display, and is convenient for performing systematic calculations on a computer.
The surface calculation method is suitable for the physical interpretation of the discrete nabla operator. Tables 1, 2, and 3 summarize the above results.
As is clear from Table 2, the discrete nabla operator can be calculated if the difference Δr between the position vectors of the nodes of the elements is obtained. And the discrete nabla operator plays the most fundamental role in the finite volume or finite element method of unstructured grids.
Furthermore, the concept of mixing the advantages of both the finite volume method and the finite element method has the potential to create a new calculation method for the next generation.

[0038] [Table 1] discrete nabla operator _{<∇ * φ> e = (} 1 / Ω e (r)) {∫∇N a dΩ (r)} * φ a ( integration range Ω _e) = Σ _{a ^{Σ i g e i <∂N a}} / ∂ξ i> e * φ a <∇ * φ> e = (1 / Ω e (r)) {∫φdΓ} * φ ( integration range Γ _e) = Σ _s ΔΓ ^s (r) * φ _s / Ω _e (r) <∇ * φ> _e = Σ _a ∇ _a * φ _a ∇ _a ≡ (1 / Ω _e (r)) {∫∇N _a dΩ (r)} = Σ _i g _e ⁱ <∂N _a / ∂ξ ⁱ > _e (definition formula)

[0039] [Table 2] discrete nabla operator display type (hexahedral ^{_{element) ∇ a ≡Σ i g e i}} <∂N a / ∂ξ i> e ( sum, i = 1 to 3) (definition formula ) ∇ _a = (1 / Ω _e (r)) (A ^* ₁ ξ _a + A ^* ₂ η _a + A ^* ₃ ζ _a ) ∇ _a = (1/4 Ω _e (r)) (S ^* ₁ ξ _a + S ^* ₂ η _a + S ^* ₃ ζ _a ) where ξ ¹ = ξ, ξ ² = η, ξ ³ = ζ g _e ¹ ≡ <g ¹ > _e = <g ₂ × g ₃ > _e / J _e (inverse basis Vector 1) g _e ² ≡ <g ² > _e = <g ₃ × g ₁ > _e / J _e (contravariant basis vector 2) g _e ³ ≡ <g ³ > _e = <g ₁ × g ₂ > _e / J _e (invariant basis vector 3) g _i ≡∂r / ∂ξ ⁱ (i = 1, 2, 3) (covariant basis vector) A ^* ₁ ≡ <g ₂ × g ₃ > _e = <(∂r ^{/ ∂ξ 2) × (∂r /} ∂ξ 3)> e A * 2 ≡ <g 3 × g 1> e = <(∂r / ∂ξ 3) × (∂r / ∂ξ 1)> e A ^* ₃ ≡ <g ₁ × g ₂ > _e = <(∂r / ∂ξ ¹ ) × (∂r / ∂ξ ² )> _e S ^* ₁ ≡ (Δr) ^* ₂ × (Δr) ^* ₃ S ^* ₂ ≡ (Δr) ^* ₃ × (Δr) ^* ₁ S ^* ₃ ≡ (Δr) ^* ₁ × (Δr) ^* ₂ (A hexahedron is defined as the vertex with 8 points where ξ, η and ζ are ± 1 respectively. If cubes were normalized to that, xi] _a is .Ita _a is xi] coordinates of the node a is .Zeta _a is η coordinates of the node a is the ζ coordinate of the node a.) (Δr) ^* ₁ = (1/4) {(Δr) ₂₁ + (Δr) ₃₄ + (Δr) ₆₅ + (Δ
r) ₇₈ } (Δr) ^* ₂ = (1/4) {(Δr) ₃₂ + (Δr) ₄₁ + (Δr) ₇₆ + (Δ
r) ₈₅ } (Δr) ^* ₃ = (1/4) {(Δr) ₅₁ + (Δr) ₆₂ + (Δr) ₇₃ + (Δ
r) ₈₄ } (Δr) _ab = r _a −r _b (a and b are node numbers)

[0040] [Table 3] Display type discrete nabla operator _{∇ a ≡ (1 / Ω e} (r)) ∫ {∇N a dΩ (r)} = Σ i g e i <∂N a / ∂ξ i > in the case of _e symmetry elements _{∇ a = Σ i (1 /} Ω e (r)) g e i ξ a i ∇ a * φ a = Σ i S * i * (δφ) * i / Ω e (r) alone In the case of an element ∇ _a ≡ (1 / Ω _e (ξ)) g ^a ∇ _a * φ _a = (1 / Ω _e (r)) g ^a * φ _a where i is the coordinate axis number and a is the element node ^{_{numbers, g e i = <g i}} > e (i = 1,2,3) ( element average contravariant basis vectors)

Third, the gradient, divergence, and rotation of a single element will be described. One-dimensional line segments, two-dimensional triangles, and three-dimensional tetrahedrons are collectively referred to as single elements. Calculate gradient, divergence, and rotation for a single element. Since the contravariant basis vector g ⁱ = ∇ξ ⁱ (i = 1, 2, 3) is constant for a single element, the element average value g _e ⁱ
It is equal to g ^i. As a result, the field calculation method of the discrete nabla operator method matches the surface calculation method, and the discrete nabla operator matches not the approximate expression but the analysis display. Then, Σ _i ξ ⁱ = 1, Σ _i g ⁱ = 0 according to each dimension.

Further, since it is N ^_a = ξ _a, equal to (when _{i ≠ a) <∂N a /} ∂ξ i> e = δ ia = 1 (i = when a) = 0. At this time, ∇ _a = g ^a . Thus, discretization of the general formula becomes ^{_{∇ a * φ a = g a}} * φ a. This is specifically shown in FIG. 1 for a one-dimensional case, as shown in FIG. 2 for a two-dimensional case, and as shown in FIG. 3 for a three-dimensional case.

The one-dimensional case will be described with reference to FIG. FIG. 1 shows a node 1 and a node 2.
_Let a line segment of length Le from to be an element. The x-coordinate and x _1, the x-coordinate and x ₂ of. _{Thus, x 2 -x 1 = L e} L e = | (1,1) T (x 1, x 2) T | become. Let the unit vector of the x coordinate be i. i is a unit vector pointing right. Cities origin of xi] ¹ coordinates, and the left is a positive direction, normalized by the length L _e. City origin of xi] ² coordinates, and the right is a positive direction, normalized by the length L _e. Denoting xi] ¹ coordinates and xi] ² coordinates x-coordinate, xi] ¹ = a _{(x 2 -x) / L e} ξ 2 = (x-x 1) / L e. Therefore, ξ ¹ + ξ ² = 1. The invariant basis vectors g ¹ and g ² are g ¹ = ∇ξ ¹ = -i / L _e g ² = ∇ξ ² = i / L _e . Therefore, g ¹ + g ² = 0. Shape function N ₁ is in is 1, in from a linear function becomes 0, N ₁ = a _{(x 2 -x) / L e} = ξ 1. Since the shape function N ₂ is a linear function that becomes 1 and becomes 0 when, N ₂ = (x−x ₁ ) / L _e = な る² . Thus, written as ^{_{N a = ξ a (a =}} 1,2).

The two-dimensional case will be described with reference to FIG. ,, are node numbers. A triangle is used as an element. Let the (x, y) coordinates of the node be (x ₁ , y ₁ ).
The (x, y) coordinates of the node are set to (x ₂ , y ₂ ). Nodal
The (x, y) coordinates are (x ₃ , y ₃ ). The area of the triangle is S
When _{e, 2S e = | (1,1,1} ) T (x 1, x 2, x 3) T (y 1, y 2, y 3) T | become. Δr is the amount of change in the position vector r of the element. (Δr) ₁ is a vector from the node to the node. (Δr) ₂ is a vector from the node to the node. (Δr) ₃ is a vector from the node to the node. The vector k is a unit vector that is perpendicular to the paper surface and that faces upward. ξ The origin of ^one coordinate is on the side, and the direction perpendicular to the side and toward the node is the positive direction. ξ
The origin of the ^two coordinates is on the side, and the direction perpendicular to the side and toward the node is the positive direction.と り The origin of the ^three coordinates is on the side, and the direction perpendicular to the side and toward the node is the positive direction. Denoting xi] ¹ at (x, y) ^{coordinates, ξ 1 = ax + by +} c 1 = ax 1 + by 1 + c 0 = ax 2 + by 2 + c 0 = ax 3 + by 3 + c a = - (y 3 -y 2) / 2S _eb = (x ₃ −x ₂ ) / 2S _e c = (x ₂ y ₃ −x ₃ y ₂ ) / 2S _e Therefore, the contravariant basis vector g ¹ is g ¹ = ∇ξ ¹ = −i (y ₃ −y ₂ ) / 2S _e + j (x ₃ −x ₂ ) / 2S _e .

By the way, (Δr) ₁ = i (x ₃ −x ₂ ) + j (y ₃ −y ₂ ) i = −k × j, j = k × ik × (Δr) ₁ = j (x ₃ − since x ₂₎ is -i (y ₃ -y _2), the ^{g 1 = ∇ξ 1 = k ×} (Δr) 1 / 2S e. Similarly for xi] ² and xi] ^3, the ^{g 2 = ∇ξ 2 = k ×} (Δr) 2 / 2S e g 3 = ∇ξ 3 = k × (Δr) 3 / 2S e. ξ ¹ , ξ ^2, and ξ ³ are not linearly independent, and (x ₃ −x ₂ ) + (x ₂ −x ₁ ) + (x ₁ −x ₃ ) = 0 (y ₃ −y ₂ ) + (y ₂ −y ₁ ) + (y ₁ −y ₃ ) = 0 (x ₂ y ₃ −x ₃ y ₂ ) + (x ₃ y ₁ −x ₁ y ₃ ) + (x ₁ y ₂ −x ₂ y
₁₎ = from a 2S _e, clearly, a ^{ξ 1 + ξ 2 + ξ 3} = 1. g ¹ and g ² and g ³ is also not linearly _{independent, (Δr) 1 + (Δr} ) 2 + (Δr) from a ₃ = 0, obviously, is g ¹ + g ² + g ³ = 0.

The shape function N ₁ node will become 1 at node, because a linear function becomes 0 at node obviously coincides with xi] ^1. Since same is true for xi] ² and xi] ^3, written as ^{_{N a = ξ a (a =}} 1,2,3).

Let the value at the node of φ be φ₁And the nodes
The value at φ_TwoAnd the value at node 3 is φ_ThreeAnd
The general formula is_a* Φ_a= G^a* Φ_a = G¹* Φ₁+ G^Two* Φ_Two+ G^Three* Φ_Three =-(G^Two+ G^Three) * Φ₁− (G^Three+ G¹) * Φ_Two− (G¹+ G^Two) * Φ_Three = -G¹* (Φ_Two+ Φ_Three) -G^Two* (Φ_Three+ Φ₁) -G^Three* (Φ₁+ Φ_Two) = − K × (Δr)₁/ S_e* (Φ_Two+ Φ_Three) / 2 −k × (Δr)_Two/ S_e* (Φ_Three+ Φ₁) / 2 −k × (Δr)_Three/ S_e* (Φ₁+ Φ_Two) / 2 = (1 / S_e) {S₁* Φ₁ ^-+ S_Two* Φ_Two ^-+ S_Three* Φ_Three ^-}. Where S₁= (Δr)₁× k, φ₁ ^-= (Φ_Two+ Φ_Three) / 2 S_Two= (Δr)_Two× k, φ_Two ^-= (Φ_Three+ Φ₁) / 2 S_Three= (Δr)_Three× k, φ_Three ^-= (Φ₁+ Φ_Two) / 2. Also, f₁ ^-= (F_Two+ F_Three) / 2 f_Two ^-= (F_Three+ F₁) / 2 f_Three ^-= (F₁+ F_Two) / 2, ∇_af_a= (1 / S_e) {S₁f₁ ^-+ S_Twof_Two ^-+ S_Threef_Three ^-} Can be written. Similarly, ∇_a・ V_a= (1 / S_e) {S₁・ V₁ ^-+ S_Two・ V_Two ^-+ S_Three・ V_Three
^-} (∇_a× v_a) · K = (1 / S_e) {(Δr)₁・ V₁ ^-+ (Δr)_Two・
v_Two ^-+ (Δr)_Three・ V_Three ^-} Can be written.

The three-dimensional case will be described with reference to FIG.
You. ,,, are node numbers. Tetrahedron
Is an element. The (x, y, z) coordinates of the node are expressed as (x₁, y₁,
z₁). The (x, y, z) coordinates of the node are expressed as (x_Two, y_Two,
z_Two). The (x, y, z) coordinates of the node are expressed as (x_Three, y_Three,
z_Three). The (x, y, z) coordinates of the node are expressed as (x_Four, y_Four,
z _Four). Let the volume of the tetrahedron be V_eThen, 6V_e= | (1,1,1,1)^T (x₁, x_Two, x_Three, x_Four)^T (y₁, y_Two,
y_Three, y_Four)^T (z₁, z_Two, z_Three, z_Four)^T| Let the position vector of the node be r₁And the position of the node
The vector_TwoAnd the position vector of the node is r_ThreeWhen
And the position vector of the node_FourAnd ξ¹To (x, y,
z) When expressed in coordinates, ξ¹= Ax + by + cz + d 1 = ax₁+ By₁+ Cz₁+ D 0 = ax_Two+ By_Two+ Cz_Two+ D 0 = ax_Three+ By_Three+ Cz_Three+ D 0 = ax_Three+ By_Three+ Cz_Three+ D a =-(S₁・ I) / 3V_e b =-(S₁・ J) / 3V_e c =-(S₁・ K) / 3V_e d = | S₁| / 3V_e Becomes Where S₁= (R_Three-R_Two) × (r_Four-R_Two) / 2. Therefore, the contravariant basis vector g¹Is g¹= ∇ξ¹ = -I (S₁・ I) / 3V_e−j (S₁・ J) / 3V_e-K (S₁
・ K) / 3V_e = -S₁/ 3V_e Becomes ξ^TwoAnd ξ^ThreeAnd ξ^FourSimilarly, for g^Two= ∇ξ^Two= -S_Two/ 3V_e g^Three= ∇ξ^Three= -S_Three/ 3V_e g^Four= ∇ξ^Four= -S_Four/ 3V_e Becomes ξ¹And ξ^TwoAnd ξ^ThreeAnd ξ^FourIs not linear independence, apparent
In, ξ¹+ Ξ^Two+ Ξ^Three+ Ξ^Four= 1. g¹And g^TwoAnd g^ThreeAnd g^FourAlso not linear independence, obvious
G¹+ G^Two+ G^Three+ G^Four= 0.

The shape function N ₁ node will become 1 at node, because a linear function becomes 0 at node obviously coincides with xi] ^1. Since same is true for xi] ² and xi] ³ and xi] ^4, written as ^{_{N a = ξ a (a =}} 1,2,3,4).

The value at the node of φ is φ ₁ , the value at the node is φ ₂ , the value at node 3 is φ _3, and the value at node 4 is φ ₄ . General _{formula, Σ a ∇ a * φ a} = Σ a g a * φ a = g 1 * φ 1 + g 2 * φ 2 + g 3 * φ 3 + g 4 * φ 4 = - (g 2 + g 3 + g 4) * Φ ₁ − (g ³ + g ⁴ + g ¹ ) * φ ₂ −
(g ⁴ + g ¹ + g ² ) * φ _3- (g ¹ + g ² + g ³ ) * φ ₄ = -g ¹ * (φ ₂ + φ ₃ + φ ₄ ) -g ² * (φ ₃ + φ ₄ + φ ₁ ) -g
³ * (φ ₄ + φ ₁ + φ ₂ ) -g ⁴ * (φ ₁ + φ ₂ + φ ₃ ) = S ₁ / V _e * (φ ₂ + φ ₃ + φ ₄ ) / 3 + S ₂ / V _e * (φ ₃ +
φ ₄ + φ ₁ ) / 3 + S ₃ / V _e * (φ ₄ + φ ₁ + φ ₂ ) / 3 + S _{4
/ V e * (φ 1 +} φ 2 + φ 3) / 3 = (1 / V e) {S 1 * φ 1 - + S 2 * φ 2 - + S 3 * φ 3 - + S 4 *
a} ^- phi _4. Where S ₁ = (r ₃ −r ₂ ) × (r ₄ −r ₂ ) / 2, φ ₁ ⁻ = (φ ₂ + φ ₃ +
φ ₄ ) / 3 S ₂ = (r ₃ −r ₁ ) × (r ₄ −r ₁ ) / 2, φ ₂ ⁻ = (φ ₃ + φ ₄ +
φ ₁ ) / 3 S ₃ = (r ₄ −r ₂ ) × (r ₁ −r ₂ ) / 2, φ ₃ ⁻ = (φ ₄ + φ ₁ +
φ ₂ ) / 3 S ₄ = (r ₃ −r ₁ ) × (r ₂ −r ₁ ) / 2, φ ₄ ⁻ = (φ ₁ + φ ₂ +
φ ₃ ) / 3. F ₁ ⁻ = (f ₂ + f ₃ + f ₄ ) / 3 f ₂ ⁻ = (f ₃ + f ₄ + f ₁ ) / 3 f ₃ ⁻ = (f ₄ + f ₁ + f ₂ ) / 3 f ₄ ⁻ = (f When _{1 + f 2 + f 3)} / 3, Σ a ∇ a f a = (1 / V e) {S 1 f 1 - + S 2 f 2 - + S 3 f 3 - + S
₄ f ₄ ^-} and write. _{Similarly, Σ a ∇ a · v a} = (1 / V e) {S 1 · v 1 - + S 2 · v 2 - + S 3 ·
^{_{v 3 - + S 4 · v}} 4 -} Σ a ∇ a × v a = (1 / V e) {S 1 × v 1 - + S 2 × v 2 - + S 3 ×
v ₃ ⁻ + S ₄ × v ₄ ⁻ }.

Fourth, the gradient, divergence, and rotation of the symmetric element will be described. One-dimensional line segments, two-dimensional tetragons, and three-dimensional hexahedrons are examined. These elements, when mapped to calculation space by mapping function _{r = Σ a N a r a} , since the symmetry with respect to each coordinate axis, referred to as a symmetrical element. Shape function of symmetry elements, for each _{dimension, = N a (1/2) (} 1 + ξ a ξ) (1 _{-dimensional) N a = (1/4) (} 1 + ξ a ξ) (1 + η a η) ( Two-dimensional) N _a = (1/8) (1 + ξ _a ξ) (1 + η _a η) (1 + ζ _a ζ) (three-dimensional). ξ _a , η _a , ζ _a is +1 or −1 determined by the node. In the case of one dimension, N ₁ = (1/2) (1-ξ); ξ _a = −1 N ₂ = (1/2) (1 + ξ); ξ _a = + 1.

[0052] Thus, the element average of each ^{_{dimension, <∂N a / ∂ξ i>}} e = (1/2) ξ a i (1 -D (i = 1)) <∂N a / ∂ξ i> e = (1/4) ξ _a ⁱ (2 dimensions (i = 1,2)) <∂N _a / ∂ξ ⁱ > _e = (1/8) ξ _a ⁱ (3 dimensions (i = 1,2,3 )) Can be displayed. Than _{this, ∇ a * φ a = Σ} i g e i * (1/2) ξ a i φ (1 _{-dimensional) ∇ a * φ a = Σ} i g e i * (1/4) ξ a i φ ( _{2D) ∇ a * φ a = Σ} i g e i * (1/8) to obtain xi] _a ⁱ phi (the 3-dimensional). One-dimensional, two-dimensional, and three-dimensional are combined to obtain a generalized formula, ∇ _a * φ _a = Σ _i (1 / Ω _e (ξ)) g _e ⁱ * (ξ _a ⁱ φ _a ). Here, Ω _e (ξ) is the volume of each dimension of the computation space. The result for the one-dimensional case is shown in FIG. The result in the case of two dimensions is shown in FIG.
As shown in FIG.

In particular, for a two-dimensional square element, the area S _i (i = 1, 2, 3, 4) for each side is defined as follows. (Δr) ₁ × k = S ₁ k × (Δr) ₃ = S ₃ (Δr) ₂ × k = S ₂ k × (Δr) ₄ = S ₄ At this time, the area formula <∇ * φ> _e = Σ with _{a ∇ a * φ a = Σ} s ΔΓ s (x) * φ s / Ω e (x), divergence and _{rotation, Σ a ∇ a * v a} = (1 / S e) [S 1 · ^{_{v 1 * + S 2 · v}} 2 * + S 3
^{_{· V 3 *] Σ a (}} ∇ a × v a) · k = (1 / S e) [(Δr) 1 · v 1 * + (Δ
r) ₂ · v ₂ ^* − (Δr) ₃ · v ₃ ^* − (Δr) ₄ · v ₄ ^* ].

[0054] In the case of parallelogram, this result is, ∇ results sigma _a calculated volume fraction Official _{a * v a = (1 /} S e) [S 1 * · δv 1 * + S 2 * · δv 2
^{_{*] Σ a (∇ a ×}} v a) · k = (1 / S e) [S 1 * × δv 1 * + S 2 *
× δv ₂ ^* ] · k. The same applies to the case of three dimensions. In other words, in the case of a parallelepiped, the results of the area formula and the volume integral formula match. And for the distorted element, the calculation is a little complicated, but if precision is important, the area formula should be used.

Fifth, the analytical display and non-storage of the coefficient matrix will be described. For finite element analysis of magneto-thermal fluid,
Many coefficient matrices appear. These coefficient matrices are listed in Tables 4 and 5. Table 4 shows discrete nabla operator representations of the Open matrix, Cross matrix, and Dot matrix, which are coefficient matrices of the second derivative. Also, discrete n of advection matrix, diffusion matrix and BTD matrix
Table 5 shows the abla operator. Here, in the case of a Cross matrix resulting from the curl curl operator, care must be taken in the order of the products. The Dot matrix arising from Laplacian is
Also called a diffusion matrix. The BTD matrix is a matrix resulting from the upstreamization of the advection term.

However, the product expression of the discrete nabla operator of the coefficient matrix is defined including the hourglass term. For example, in the diffusion matrix _{D ab = Ω e ∇ a ·} ∇ b = J e g e i · g e j <(∂N a / ∂ξ i) (∂N b / ∂ξ j)> e = J e g _e ⁱ · g _e ^j <∂N _a / ∂ξ ⁱ > _e <∂N _b / ∂ξ ^j > _e + (ho
urglass section). In other words, when ∇ _a and ∇ _b are independently dot-products, the hourglass term (secondary part)
Need to be added. Since the coefficient matrix is well known, further detailed description will be omitted, but if necessary, refer to Document 1 and the like.

[0057] [Table 4] coefficient matrix and discrete nabla operator (Part 1 2 derivative) Open matrix asymmetric matrix _{O ab = ∫∇N a ∇N b dΩ} (r) = Ω e (r) ∇ a ∇ b = Σ _i Σ _j Ω _e (r) g _e ⁱ g _e ^j <(∂N _a / ∂ξ ⁱ ) (∂N _b / ∂ξ
^j)> Cross matrix antisymmetric matrix _{R ab = ∫∇N a × ∇N b} dΩ (r) = Ω e (r) ∇ a × ∇ b = Σ i Σ j Ω e (r) g e i × g e ^j <(∂N _a / ∂ξ ⁱ ) (∂N _b / ∂
ξ ^j)> Dot matrix symmetric matrix _{D ab = ∫∇N a · ∇N b} dΩ (r) = Ω e (r) ∇ a · ∇ b = Σ i Σ j Ω e (r) g e i · g e ^j <(∂N _a / ∂ξ ⁱ ) (∂N _b / ∂
ξ ^j )>

[Table 5] Coefficient matrix and discrete nabla operator (No. 2) Mass matrix (symmetric matrix) M _ab = ∫N _a N _b dΩ (r) Concentrated mass matrix (n is the number of nodes of elements) M ^* _{ab = J e δ ab J e =} Ω e (r) / n advection matrix (non-symmetric _{matrices) A ab = ∫N a v ·} ∇N b dΩ (r) = Σ j Ω e (r) v e · g e j _{<N a (∂N b / ∂ξ} j)> diffusion matrix (Dot _{matrix) D ab = ∫∇N a · ∇N} b dΩ (r) = Ω e (r) ∇ a · ∇ b = Σ i Σ j Ω _e (r) g _e ⁱ · g _e ^j <(∂N _a / ∂ξ ⁱ ) (∂N _b / ∂
xi] ^j)> BTD matrix (symmetric _{matrix) B ab = ∫v · ∇N a} v · ∇N b dΩ (r) = Ω e (r) v e · ∇ a v e · ∇ b = Σ i Σ j Ω _{e (r) v e · g} e i v e · g e j <(∂N a / ∂ξ i)
(∂N _b / ∂ξ ^j )> What is needed for the calculation is <∂r / ∂ξ ⁱ > (i = 1, 2, 3)
And _{<N a (∂N b / ∂ξ} j)>, <(∂N a / ∂ξ i) (∂N b / ∂
ξ ^j )> only.

Storing this coefficient matrix requires an enormous amount of memory. If all calculations are performed each time without storing, the calculation time becomes longer. Therefore, non-memory and a high-speed calculation method are required. For this purpose, in the matrix calculation, the matrix elements are decomposed into two parts and displayed as a sum of products. That is, (coefficient matrix element) = Σ (coefficient matrix base) × <partial coefficient factor>.

<Partial factor> is obtained analytically in a portion that does not depend on the shape of each element. Can be calculated by the calculation space <partial coefficient factor> ^{_{basically, <(∂N a / ∂ξ i}} ) (∂
_{^{N b / ∂ξ j)> e}} , only two _{<N a (∂N b / ∂ξ} i)> e. These are generated and unstored each time they are needed. On the other hand, the remaining part is called (coefficient matrix base). (Coefficient matrix basis) is a part that depends on the element shape. This can be easily calculated using the covariant basis vector g _i (≡∂r / ∂ξ ⁱ (i = 1, 2, 3)). (Coefficient matrix basis) calculates only the types of elements. For example, if the whole is composed of eight types of elements, it is sufficient to calculate eight elements. 100 million of even elements, everyone if an element of the same shape, only information about one of the elements <∂r / ∂ξ ^i> _e
(i = 1, 2, 3) should just be stored. This part is stored in the coefficient matrix basis storage device. When coefficient matrix elements are needed for matrix calculation, read (coefficient matrix basis)
A <partial coefficient factor> is generated and combined by a product-sum operation to generate a coefficient matrix element, thereby making it non-memory.

Sixth, the conversion of the contravariant basis vector into an object will be described. Discrete nabla operator is defined by ^{_{∇ e ≡Σ i g e i <}} ∂N a / ∂ξ i> e. Here, the calculation of the contravariant basis vector g _e ⁱ is important. And the invariant basis vector g _e ⁱ is common to the finite volume method and the finite element method. As the types of elements increase, the amount of calculation also increases. In addition, the storage capacity also increases. Then, given the type of element and its coordinates, the routine for calculating this contravariant basis vector g _e ⁱ is converted into an object and an LSI, and if node coordinates are sent as a message, g _e ⁱ can be obtained. Make it available as a built-in function. If a three-dimensional 4/5 / 6-hedron element is formed into chips, analysis using a mixed element becomes possible, and grid generation becomes extremely easy. This is a Dell chip (see FIG. 7). The name del is an alias for the nabla operator.

Referring to FIGS. 8 and 9, a non-memory type calculation procedure of the finite element method using the discrete nabla operator will be described. In step 1, the type of partial differential equation is selected, and parameters and boundary conditions are input. Information about the computable partial differential equations stored in the partial differential equation information storage device 6 is displayed on the display device 7, and necessary information is input from the input device 1. The input parameters and boundary condition data of the PDE are stored in the parameter boundary condition storage device 2. In step 2, the vector of the node position is input to the input device 1
Enter from. This includes not only inputting individual data but also generating a grid automatically. The coordinate values of all nodes are stored in the node position vector storage device 3.

In step 3, the known vector calculation means 4 calculates a known vector from the parameters and the boundary conditions, and stores the known vector in the known vector storage device 5. In step 4, the coefficient matrix basis calculation means 8 calculates the invariant basis vector and the element volume from the nodal coordinates, and stores the coefficient matrix basis in the coefficient matrix basis storage device 9. As the coefficient matrix base calculation means 8, the Delchip 72 of FIG. 7 can be used.

In step 5, the coefficient matrix element calculating means is obtained from the partial coefficient factors generated by the partial coefficient factor generating means 11 from the analytically calculated shape function data and the coefficient matrix base read from the coefficient matrix base storage device 9. According to 12, the coefficient matrix element is calculated. The coefficient matrix elements are temporarily stored in the matrix element storage work area 14. In step 6, the matrix operation means 13 calculates a solution vector from the coefficient matrix element and the known vector. Steps 5 and 6 are repeated to temporarily store the coefficient matrix elements calculated from the coefficient matrix bases and the partial coefficient factors using the work area and continue the matrix calculation to obtain a solution. In step 7, the analysis result is displayed on a display device 7 such as a CRT display. As described above, since the coefficient matrix element is calculated each time it becomes necessary in the matrix calculation, a memory for storing the entire coefficient matrix is not required. Although there are only (number of nodes) × (number of nodes) coefficient matrix elements, non-zero values take about several tens per row.
The work area is small because it is concentrated near the diagonal components.

As described above, in the embodiment of the present invention,
A numerical analysis device that calculates by the finite volume method or the finite element method of unstructured grids stores the factors that depend on the element shape of the coefficient matrix element described by the discrete nabla operator, and analyzes the factors that do not depend on the element shape Since the coefficient matrix element is generated from now on whenever it is needed, the calculation can be speeded up while reducing the storage capacity.

[0066]

As is apparent from the above description, according to the present invention, a numerical analysis apparatus for calculating a solution of a boundary value problem of a partial differential equation by a finite volume method or a finite element method of an unstructured grid is provided by a discrete nabla operation. Coefficient matrix base storage means for storing factors depending on element shapes of coefficient matrix elements represented by children in correspondence with element shapes, and independent of element shapes of coefficient matrix elements represented by discrete nabla operators A partial coefficient factor generating means for non-storing factors for each type to generate a partial coefficient factor, and a coefficient matrix base and a partial coefficient factor generating means obtained by referring to a coefficient matrix base storing means each time a matrix calculation is required Since it is configured to include a calculation means for generating a coefficient matrix element from the obtained partial coefficient factors, an effect is obtained that the calculation can be sped up while the storage capacity is reduced.

[Brief description of the drawings]

FIG. 1 is an explanatory diagram of gradient, divergence, and rotation of a one-dimensional single element used in a numerical analysis device according to an embodiment of the present invention;

FIG. 2 is an explanatory diagram of gradient, divergence, and rotation of a two-dimensional single element used in the numerical analysis device according to the embodiment of the present invention;

FIG. 3 is an explanatory diagram of gradient, divergence, and rotation of a three-dimensional single element used in the numerical analysis device according to the embodiment of the present invention;

FIG. 4 is an explanatory diagram of gradient, divergence, and rotation of a one-dimensional symmetric element used in the numerical analysis device according to the embodiment of the present invention;

FIG. 5 is an explanatory diagram of gradient, divergence, and rotation of a two-dimensional symmetric element used in the numerical analysis device according to the embodiment of the present invention;

FIG. 6 is an explanatory diagram of gradient, divergence, and rotation of a three-dimensional symmetric element used in the numerical analysis device according to the embodiment of the present invention;

FIG. 7 is a conceptual diagram of a Dell chip of the numerical analysis device according to the embodiment of the present invention.

FIG. 8 is a flowchart showing a processing procedure of the numerical analysis device according to the embodiment of the present invention;

FIG. 9 is a functional block diagram of the numerical analysis device according to the embodiment of the present invention.

[Explanation of symbols]

 REFERENCE SIGNS LIST 1 input device 2 parameter boundary condition storage device 3 node position vector storage device 4 known vector calculation device 5 known vector storage device 6 partial differential equation information storage device 7 display device 8 coefficient matrix base calculation device 9 coefficient matrix base storage device 10 shape function Storage device 11 Partial coefficient factor generation means 12 Coefficient matrix element calculation means 13 Matrix calculation means 14 Matrix element storage work area 71 Element node coordinate generation means 72 Delchip 73 Data utilization means

Claims (4)

[Claims] 1. A numerical analysis apparatus for calculating a solution of a boundary value problem of a partial differential equation by a finite volume method or a finite element method of an unstructured grid, wherein the element shape of a coefficient matrix element represented by a discrete nabla operator is used. Coefficient matrix basis storage means for storing dependent factors in association with element shapes, discrete nabla
A partial coefficient factor generating means for storing a factor independent of the element shape among coefficient matrix elements represented by an operator for each type to generate a partial coefficient factor; and a coefficient matrix base storing means each time a matrix calculation is required. And a calculating means for generating a coefficient matrix element from a coefficient matrix base obtained by referring to (1) and the partial coefficient factor generated by the partial coefficient factor generating means. 2. A numerical analysis method for calculating a solution of a boundary value problem of a partial differential equation by a finite volume method or a finite element method of an unstructured grid, wherein the element shape of a coefficient matrix element represented by a discrete nabla operator is used. Dependent factors are stored as coefficient matrix bases corresponding to element shapes, and factors that do not depend on element shape among the coefficient matrix elements expressed by the discrete nabla operator are generated for each type as partial coefficient factors, which are necessary for matrix calculation A numerical analysis method, wherein a coefficient matrix element is generated from the coefficient matrix base and the partial coefficient factor each time the matrix matrix becomes, and a matrix calculation is executed without storing all the coefficient matrix elements. Discrete nabla operator ∇ _a corresponding to the node a of 3. Element _{e, ∇ a = Σ i g} e i <∂N a / ∂ξ i> e ( sum, i = 1 to k,
k is the number of dimensions (where a is the node number of the element, i is the number of the coordinate axis, N _a is the shape function, ξ ⁱ is the calculated coordinate, and g _e ⁱ is g _e ⁱ ≡ <g ⁱ > _e = <∇ξ ⁱ > _e = (1 / Ω _e (x)) ∫ {∇ξ ⁱ dΩ (x)} (the integration range is Ω _e ). The numerical analysis method according to claim 2, wherein: The partial coefficient factor that is independent of the element shapes of the discrete nabla operator corresponding to the node a as claimed in claim 4 Element e, calculated in ^{_{<∂N a / ∂ξ i> e}} (i = 1~k) Discrete nabla corresponding to node a of element e
The coefficient matrix base depending on the element shape of the operator is obtained by calculating g _e ⁱ = <g ⁱ > _e (i = 1 to k), and (coefficient matrix element) = (coefficient matrix base) × 3. The numerical analysis method according to claim 2, wherein a coefficient matrix element is calculated by a product-sum operation called <partial coefficient factor>.