[ Article ]

Transactions of the Korean Society for Noise and Vibration Engineering - Vol. 29, No. 3, pp.327-346

ISSN: 1598-2785 (Print) 2287-5476 (Online)

Print publication date 20 Jun 2019

Received 14 Feb 2019 Revised 23 Apr 2019 Accepted 23 Apr 2019

DOI: https://doi.org/10.5050/KSNVE.2019.29.3.327

Development of Structural Displacement Finite Element Codes Based on Physical Coordinates and High Performance Computing Methods Including PARDISO, PCG, GMRES and MATLAB or OCTAVE Intrinsic Operator ‘\’

Seok-Tae Park^†

물리적 좌표기반 구조 변위해석 유한 요소 코드 및 PARDISO, PCG, GMRES와 MATLAB 또는 OCTAVE 연산자 ‘\’ 를 포함하는 여러 고성능 계산기법 개발

박석태^†

Correspondence to: ^†Member, Dept. of Academic Affairs, Chungbuk Health and Science University E-mail : stpark@chsu.ac.kr # A part of this paper was presented at the KSNVE 2019 Annual Spring Conference.

‡ Recommended by Editor Won Ju Jeon

Abstract

This paper presents physical coordinate-based 3-node triangular, 4-node rectangular, and 4-node quadrilateral FE codes on the Kirchhoff thin plate theory. Isoparametric 4-node quadrilateral and 4-node plate shell FE codes based on natural coordinates are also presented. The calculation time needed to solve a system equation was compared using direct, iterative, and full matrix and sparse matrix solvers including Gauss elimination method, CG, DGESV, DSGESV, SOR, PCG, GMRES, FGMRES, MGMRES, PARDISO, etc. The full matrix and the sparse matrix were applied to compare the computing time of the system equation using OCTAVE or the MATLAB intrinsic operator ‘\’. The computation times for PCG methods using four types of preconditioners were also compared. The computation time of the finite element method is compared with that of the PCG method after processing the system matrix with the full matrix and sparse matrix forms. We showed that PARDISO and sparse MATLAB or OCTAVE operator ‘\’ could be useful and effective methods for conducting large-scale finite element analyses using personal computers.

초록

이 논문에서는 키르히호프 평판 이론에 기초하여 물리적 좌표기반 3절점 삼각형 유한요소, 4절점 직사각형 유한요소, 4절점 사각형 유한요소 코드들을 제시하였다. Natural coordinates에 기반한 isoparametric 4절점 사각형 유한요소 코드와 4절점 plate shell 유한요소 코드도 제시하였다. 직접법, 반복법 및 full matrix 해법들과 희박행렬 해법들을 이용하여 시스템 방정식을 푸는데 필요한 계산시간을 비교하였다. Gauss 소거법, CG, DGESV, DSGESV, SOR, PCG, GMRES, FGMRES, MGMRES, PARDISO 법 등을 비교 검토하였다. OCTAVE와 MATLAB에 내장된 연산자 ‘\’를 full matrix 형태와 희박행렬 형태 시스템 방정식에 적용하여 계산시간을 비교하였다. 4가지 종류의 preconditioners를 채택한 PCG 방법들에 대한 계산시간도 비교 검토하였다. 또한, 유한요소법에서 시스템 행렬을 full matrix와 희박행렬로 처리한 후에 각각의 경우에 PCG 방법을 적용했을 때의 계산시간도 비교하였다. PC에서도 PARDISO와 sparse MATLAB 연산자 ‘\’를 사용하여 대규모의 유한요소 해석을 수행할 수 있음을 보였다.

Keywords:

Analytical Integration, CG, Conjugate Gradient, Direct Solver, GMRES, Generalized Minimal RESidual, FGMRES, Flexible Generalized Minimal RESidual, Iterative Solver , Kirchhoff-Love Plate, Krylov Subspace Method, MATLAB or OCTAVE Intrinsic Operator ‘\’, MGMRES, Multiple (Restarted) Generalized Minimal RESidual, Mindlin-Reissner Plate, PARDISO, PARallel DIrect Sparse sOlver, PCG, Preconditioned Conjugate Gradient, Physical Coordinates Based FEM, Sparse Matrix Solver

키워드:

해석적 적분, 켤레 기울기법, 직접 해법, 일반화된 최소 잔차법, 유연한 일반화된 최소 잔차법, 반복 해법, 키르히호프-러브 평판, 크릴로프 부 공간 법, 매트랩 또는 옥타브 내재 연산자 ‘\’, 다중 일반화된 최소잔차법, 민들린-라이즈너 평판, 병렬 직접 희박 해법, 선조건 적용 켤레 기울기법, 물리적 좌표기반 유한요소법, 희박행렬 해법

1. History and Research Motivation of Thin Plate Analysis

In this paper, the development of structural displacement finite element method (FEM) codes and the computational time comparisons of system equation using several solvers were described. For the first topic, we dealt with the Kirchhoff flat plate theory in Section 2, and described the formulations of the finite element (FE) in Section 3. We verified the validity of the developed FE codes in Section 4. As a second theme, in Section 5, we compared the computation times obtained by various solving methods of the system equation within structural static displacement FE codes, and proposed an appropriate analysis methods. Acoustic BEM (boundary element method) code and structural vibration FEM code were needed to perform structural-acoustical coupling analysis by developing high-speed computations of large degrees of freedom without using commercial software. The reason for developing the structural displacement analysis FEM program was that Holmström⁽¹⁾ omitted the FEM program for vibration analysis with blackbox in the program presented in the vibration-acoustic coupled analysis paper. (1) In this paper, the first developed FEM codes were the same as rectangular FE codes presented by Irvine^{(2 ~ 4)}. It was a physical coordinate based rectangular FE codes that could be applied only when the four sides were all parallel to the x or y axis⁽²⁾, developed according to the assumptions used in the Kirchhoff plate theory (thin plate theory). (2) Semie suggested a physical coordinate based 3-node triangular plate FE codes that could model arbitrary two-dimensional shapes and analyzed the displacement of a thin plate⁽⁵⁾. He investigated two cases where four sides of the square plate were simply supported or fixed (clamped). However, the analytical results were significantly deviated from the values of the commercial software ANSYS^® and Kirchhoff thin plate’s theory. In this paper, we developed a physical coordinate based triangular FE codes according to the methodology suggested by Semie, and found that there would be some errors in the FE codes presented by Semie. (3) We have also developed a physical coordinate based 4-node quadrilateral FE codes that can model any two-dimensional shape in order to complement the disadvantages of FEM using rectangular FEs parallel to the x- or y-axis. (4) On the other hand, the above three FE codes were assumed to ignore transverse shear strain effects because of the Kirchhoff thin plate theory. Ferreira proposed an isoparametric quadrilateral FE codes based on the Mindlin-Reissner plate theory, which considered transverse shear deformation⁽⁶⁾. We also developed the natural coordinates based isoparametric quadrilateral FE codes which were based on the thick plate theory, and compared with the displacement analysis results of the Kirchhoff plate theory. (5) The displacement analysis FEM codes based on the Kirchhoff plate and Mindlin-Reissner plate theories were able to only analyze the z-direction displacement in the normal direction of the 2-D coplanar plate. Kwon proposed a 3-dimensional isoparametric plate shell FE codes for static displacement analysis in 3-D space using MATLAB^® language⁽⁷⁾. We have translated this codes into FORTRAN ones to improve memory scalability and computational speed. Consequently, four 2-D coplanar structural FE codes and one 3-D structural displacement FE codes were described in this paper. The differences between Kirchhoff theory, Mindlin-Reissner plate theory and three dimensional plate shell element were described for two dimensional plate. In order to evaluate the performance of FEM system equation, which was a system matrix equation Ax=b, the computation times of MATLAB^®, OCTAVE^® and FORTRAN codes were compared in both direct and iterative methods. We also compared the computation time for full matrix solvers that processed the coefficient matrix A of system equation as a full matrix form and sparse matrix solvers for the sparse matrix one. The computation times for PCG methods using four kinds of preconditioners were also compared in sparse system matrix. The finite element method was compared with the computation time of PCG method after processing the system matrix with full matrix and sparse matrix form, respectively.

2. Kirchhoff Thin Plate Theory

According to Semie⁽⁵⁾, in 1776 Euler first attempted to solve the free vibration problem of the plate. Germain developed differential equations for flat plates without considering warping term, and Lagrange proposed general plate differential equations by considering warping term. Poisson proposed the Germain-Lagrange plate equation assuming a constant flexural rigidity for the plate subjected to static loads. Navier has established a flat bending theory that uses flexural rigidity as a function of plate thickness and proposed a method to obtain exact solution using Fourier trigonometric series. Kirchhoff proposed the theory of plate bending using basic assumptions known as ‘Kirchhoff's hypothesis.’ Kirchhoff's thin plate’s hypothesis is expressed: the cross sections of the plate are remained as straight and unstretched and also normal to midsurface even after the plate is deformed. That is, γ_xz=0, γ_yz=0, σ_zz=0 and ε_zz=0. Let the reference in the z direction be the middle of the plate thickness, and the displacements in the x and y directions in the plate u, v respectively, and the displacement in the z direction w (Fig. 1).

Fig. 1

Slice cut of plate normal to y axis

u = - z θ x = - z ∂ w ∂ x

(1)

v = - z θ y = - z ∂ w ∂ y

(2)

For inplane strains, strain-displacement relationships:

ɛ x = ∂ u ∂ x = - z ∂ 2 w ∂ x 2

(3)

ɛ y = ∂ v ∂ y = - z ∂ 2 w ∂ y 2

(4)

γ x y = ∂ u ∂ y + ∂ v ∂ x = - 2 z ∂ 2 w ∂ x ∂ y

(5)

The relationship between plane stress and plane strain for isotropic materials according to Kirchhoff’s assumptions is as follows. Stress-strain relationships:

σ x = E 1 - υ 2 ɛ x + υ ɛ y

(6)

σ y = E 1 - υ 2 ɛ y + υ ɛ x

(7)

τ x y = G γ x y

(8)

Where, shear modulus $G = E 2 1 + υ$ , Young’s modulus E, Poisson’s ratio υ

Bending moments⁽⁷⁾ along the x, y, and xy edge (Fig. 2):

Fig. 2

Free body diagram of the plate element

M x = ∫ - t / 2 t / 2 z σ x d z = - D ∂ 2 w ∂ x 2 + υ ∂ 2 w ∂ y 2

(9)

M y = ∫ - t / 2 t / 2 z σ y d z = - D ∂ 2 w ∂ y 2 + υ ∂ 2 w ∂ x 2

(10)

M x y = ∫ - t / 2 t / 2 z τ x y d z = - D 1 - υ ∂ 2 w ∂ x ∂ y

(11)

Where, bending rigidity of the plate $D = E t 3 12 1 - υ 2$ and t is thickness of the plate.

By the moment equilibrium,

∂ Q x ∂ x + ∂ Q y ∂ y + q = 0

(12)

∂ M x ∂ x + ∂ M x y ∂ y - Q x = 0

(13)

∂ M y ∂ y + ∂ M x y ∂ x - Q y = 0

(14)

Inserting Eq. (9) ~ Eq. (11) into Eq. (13) and Eq. (14), after $∂ E q . 13 ∂ x + ∂ E q . 14 ∂ y + ∂ E q . 12$

D ∂ 4 w ∂ x 4 + 2 ∂ 4 w ∂ x 2 ∂ y 2 + ∂ 4 w ∂ y 4 = q

(15)

This is the governing equation for Kirchhoff thin plate bending theory with transversely loaded.

2.1 Analytical Solution of Thin Plate

Boundary conditions for a simple supported rectangular plate: width and length of the plate are a and b, respectively.

W = 0, ∂ 2 w ∂ x 2 = 0, a t x = 0, a

W = 0, ∂ 2 w ∂ y 2 = 0, a t y = 0, b

Using the infinite double Fourier trigonometric series satisfying the boundary condition in the method proposed by Navier^(5,8), out of plate displacement w(x, y) can be expressed as Eq. (16).

w x, y = ∑ m = 1 ∞ ∑ n = 1 ∞ w m n s i n m π x a s i n m π y b

(16)

q = p x, y = ∑ m = 1 ∞ ∑ n = 1 ∞ p m n s i n m π x a s i n m π y b

(17)

p m n = 4 a b ∫ 0 a ∫ 0 b p x, y s i n m π x a s i n n π y b d x d y

(18)

Eq. (16) and Eq. (17) can be inserted into Eq. (15) and solved w(x, y), but Levy⁽⁸⁾ proposed a simple series method and got the maximum displacement at the center of the square plate (a=b).

w m a x = 0.00406 p 0 a 4 D

(19)

The displacement due to concentrated load P at the center of the square plate is followed⁽⁸⁾.

w m a x = 0.01160 P a 2 D

(20)

The maximum displacements at the center of the square plate with respect to the distributed loads p₀ and concentrated loads P for the fixed boundary conditions are given by Eq. (21) and Eq. (22)⁽⁸⁾.

w m a x = 0.00126 p 0 a 4 D f o r p 0

(21)

w m a x = 0.00560 P a 2 D f o r P

(22)

3. Formulations of Various Structural Displacement Finite Elements

In this paper, three-node triangular, four-node rectangular and quadrilateral FE codes based on physical coordinates and Kirchhoff thin plate theory were presented and also natural coordinates based isoparametric four-node quadrilateral FE codes according to Reissner-Mindlin thick plate theory and 3-D plate shell FE codes. According to Kirchhoff’s thin plate theory, the transverse shear strain energy is negligible compared to the bending energy. Therefore, the total strain energy of the plate⁽²⁾ is as Eq. (23).

U = 12 ∫ v σ t ɛ d v = D 2 ∬ ∂ 2 w ∂ x 2 2 + ∂ 2 w ∂ y 2 2 + 2 υ ∂ 2 w ∂ x 2 ∂ 2 w ∂ y 2 + 2 1 - υ ∂ 2 w ∂ x ∂ y 2 d x d y

(23)

σ t = σ x σ x τ x y

(24)

ɛ t = ɛ x ε x γ x y

(25)

Total potential energy Π, potential energy of external forces V

Π = U - V

(26)

Where, $V = ∫ υ w t f d υ$

3.1 Triangular Finite Element Based on Physical Coordinates

Batoz proposed a DKT (discrete Kirchhoff theory) element⁽⁹⁾ with the origin of the local coordinate of the geometric centroid of the triangle. But, in this paper, the origin of the local coordinates was defined as in Fig. 3. In Fig. 3, let the longest side of the triangular element coincide with the x-axis of the local coordinates. The local triangle with new coordinate values was constructed by translational and rotational movement so that the vertex facing the longest side was in the first quadrant. FE codes using 9 degrees of freedom (DOF) triangular element with 3 DOF per the triangular node were developed. The node numbering order was counterclockwise. Local FEs constructed from local coordinates were converted to their original global coordinates by rotational movement. The global total stiffness matrix was constructed by assembling each global FEs. W was out of plane displacement in z coordinate and was expressed by physical coordinates like as Eq. (27)⁽⁵⁾.

Fig. 3

Local and global coordinates

w x, y, t = a 1 t + a 2 t x + a 3 t y + a 4 t x 2 + a 5 t x y + a 6 t y 2 + a 7 t x 3 + a 8 t x 2 y + x y 2 + a 9 t y 3 = A t t Z

(27)

A t t = a 1 t ⋯ a 9 t

(28)

Z t = 1 x y x 2 x y y 2 x 3 x 2 y + x y 2 y 3

(29)

For node i,

w 0,0 = w i = a 1 θ x i = ∂ w / ∂ x = a 2, θ y i = ∂ w / ∂ y = a 3

After applying the same processes for node j and node k,

w e t = w i θ x i θ y i w j θ x j θ y j w k θ x k θ y k

(30)

w e = B A

(31)

That is,

w i θ x i θ y i w j θ x j θ y j w k θ x k θ y k = 100000000010000000001000000 1 x j 0 x j 2 00 x j 3 00 010 2 x j 00 3 x j 2 00 0010 x j 00 x j 2 0 1 x k y k x k 2 x k y k y k 2 x k 3 x k 2 y k + x k y k 2 y k 3 010 2 x k y k 0 3 x k 2 2 x k y k + y k 2 0 0010 x k 2 y k 0 x k 2 + 2 x k y k 3 y k 2 a 1 a 2 a 3 a 4 a 5 a 6 a 7 a 8 a 9

(32)

A = B - 1 w e = C w e

(33)

Where, [C]=[B]^-1

Equation (23) can be rewritten as Eq. (34).

U = D 2 w e t C t ∬ S x, y d x d y C w e = 12 w e t k e w e

(34)

Local element stiffness matrix [k]e,

k e = D C t ∬ S (x, y) d x d y [C]

(35)

S x, y 9 x 9 = ∂ 2 Z ∂ x 2 ∂ 2 Z ∂ x 2 t + ∂ 2 Z ∂ y 2 ∂ 2 Z ∂ y 2 t + υ ∂ 2 Z ∂ x 2 ∂ 2 Z ∂ y 2 t + ∂ 2 Z ∂ y 2 ∂ 2 Z ∂ x 2 t + 2 1 - υ ∂ 2 Z ∂ x ∂ y ∂ 2 Z ∂ y ∂ x t

(36)

Where, β=2(1-υ)

Analytical area integration (Fig. 4) for S₄₇ = 12x,

Fig. 4

Analytical integration

∬ S 47 d x d y = ∫ 0 y k ∫ y / p 1 y - q / p 2 12 x d x d y = 2 y k 3 q 2 - 3 q y k = y k 2 / p 22 - 2 y k 3 / p 12

Where, p₁=y_k/x_k, p₂=y_k/(x_k-x_j), q=-x_jp₂

Global element stiffness matrix [k]g was obtained like as follows.

k g = T t k e T

(37)

Where, rotation matrix $T = T 1 00 0 T 1 0 00 T 1$

T 1 = 100 0 cos ⁡ - ϕ - sin ⁡ - ϕ 0 sin ⁡ - ϕ cos ⁡ - ϕ

By variational approaches for Eq. (26),

k g w g = f g

(38)

Total global stiffness matrix, displacement vector and force vector were followed, respectively.

K = ∑ k g

(39)

W = ∑ w g

(40)

F = ∑ f g

(41)

System equation was expressed as Eq. (42).

K W = F

(42)

3.2 Rectangular Finite Element Based on Physical Coordinates

According to Tom Irvine⁽²⁾, the out of displacement W can be expressed as physical coordinates as 4 DOF at each of four vertices of a rectangle.

w x, y, t = a 1 t + a 2 t x + a 3 t y + a 4 t x 2 + a 5 t x y + a 6 t y 2 + a 7 t x 3 + a 8 t x 2 y + a 9 t x y 2 + a 10 t y 3 + a 11 t x 3 y + a 12 t x y 3

(43)

The FE could be obtained in the same way as Sec 3.1. For S₄₇ = 12x,

∬ S 47 d x d y = ∫ 0 b ∫ 0 a 12 x d x d y = 6 a 2 b

Where, a and b was side length of rectangular element.

Analytical area integration was simple compared to triangular FE, but rectangular FE suggested by Irvine could be applied only if all sides were parallel to x or y axis. However, in this paper, even though the local rectangular element was not parallel to the x or y axis, we have modified the FE codes so that it could be applicable to any orientation of rectangle FE by applying the rotational transformation as shown in Eq. (37).

3.3 Quadrilateral Finite Element Based on Physical Coordinates

Batoz developed the DKQ (discrete Kirchhoff quadrilateral) element by extending the DKT concept. However, the results of the convergence rates during displacement analysis were not good compared to McNeal or LORA by Robinson and Haggenmacher⁽¹⁰⁾. In this paper, quadrilateral FE were developed by the same physical coordinate based methodology as the triangular FE mentioned in Sec 3.1. The FE formulation process was the same as the triangular FE in Sec 3.1, but the analytical area integration $∬ S x, y d x d y$ was very lengthy and was obtained using software.

3.4 Isoparametric Quadrilateral Finite Element (Mindlin-Reissner Thick Plate Theory) Based on Natural Coordinates

Mindlin-Reissner, Naghdi plate theory^{(6,7,11 ~ 13)} used a generalization of the Kirchhoff hypothesis: points of the plate which are on the normal to the undeformed middle surface remain on a straight line even after deformations are proceeded but, which are not necessarily normal to the deformed middle surface⁽¹¹⁾. That means γ_xz≠0, γ_xz≠0, σ_zz=0 and ε_zz=0. Therefore, shear deformation term was added to traditional Kirchhoff plate theory, Eq. (23) as a second term.

U = 12 ∫ v σ t ɛ d v + α 2 ∫ v σ c t ɛ c d v

(44)

σ c t = τ x z τ y z

(45)

ɛ c t = γ x z γ y z

(46)

Where, α is shear correction factor⁽⁶⁾, 5/6, {σ_c} and {ε_c} are stress and strain due to shear deformations. Transverse shear strain-displacement relationship⁽⁶⁾:

γ x z = ∂ w ∂ x + ∂ u ∂ z = ∂ w ∂ x + θ x

(47)

γ y z = ∂ w ∂ y + ∂ v ∂ z = ∂ w ∂ y + θ y

(48)

σ = D ɛ

(49)

σ c = D c ɛ c

(50)

Where, $D = E 1 - υ 2 1 υ 0 υ 10 00 1 - υ 2, D c = G 0 0 G$

With the shape function,

w = ∑ i = 1 4 N i ζ, η w i, θ x = ∑ i = 1 4 N i ζ, η θ x i, θ y = ∑ i = 1 4 N i ζ, η θ y i

(51)

ɛ = - z B w e

(52)

ɛ c = B c w e

(53)

B = 0 ∂ N 1 ∂ x 00 ∂ N 2 ∂ x 00 ∂ N 3 ∂ x 00 ∂ N 4 ∂ x 0 00 ∂ N 1 ∂ y 00 ∂ N 2 ∂ y 00 ∂ N 3 ∂ y 00 ∂ N 4 ∂ y 0 ∂ N 1 ∂ y ∂ N 1 ∂ x 0 ∂ N 2 ∂ y ∂ N 2 ∂ x 0 ∂ N 3 ∂ y ∂ N 3 ∂ x 0 ∂ N 4 ∂ y ∂ N 4 ∂ x

(54)

B c = ∂ N 1 ∂ x N 1 0 ∂ N 2 ∂ x N 2 0 ∂ N 3 ∂ x N 3 0 ∂ N 4 ∂ x N 4 0 ∂ N 1 ∂ y 0 N 1 ∂ N 2 ∂ y 0 N 2 ∂ N 3 ∂ y 0 N 3 ∂ N 4 ∂ y 0 N 4

(55)

w e t = w i θ x i θ y i w j θ x j θ y j w k θ x k θ y k w l θ x l θ y l

(56)

U = 12 w e t ∫ A B t t 3 12 D B d A w e + 12 w e t ∫ A α t B c t D c B c d A w e = 12 w e t k e w e

(57)

k e = t 3 12 ∫ A B t t 3 12 D B d A + ∫ A α t B c t D c B c d A

(58)

f e = ∫ - 1 1 ∫ - 1 1 N P J d ξ d η

(59)

k e w e = f e

(60)

3.5 Plate Shell Finite Element Based on Natural Coordinates

Kwon proposed MATLAB^® FE codes using plate shell finite element, 5 DOFs per node, which combined membrane element and plate bending element using natural coordinates based isoparametric elements⁽⁷⁾. Since the inplane rotational term was not included in plate shell element, it left null or zero values in the stiffness matrix. So-called drilling DOF incurred singularity in structural stiffness matrix if all the elements were coplanar⁽⁷⁾. To solve this problem, Knight⁽¹⁴⁾ suggested that a very small value be specified for the stiffness of the drilling DOF so that the contribution to the strain energy equation from this term would be zero⁽⁷⁾. The local stiffness matrix in the natural local coordinates was transformed into the global coordinates and assembled to obtain the total global stiffness matrix.

4. Validation

The FE codes suggested in Section 3 were validated using a square plate with four sides fixed. The plate specifications and the analysis conditions used in the analysis were as followed. Young’s modulus E=2×10¹¹Pa, Poisson ratio υ=0.285, density ρ= 7874kg/m³ and width×length, 1m×1m with thickness 0.1m. A square plate with a uniform distribution load of 50000Pa was modeled on a 10×10, 20×20, 40×40 grids on a flat plate. Theoretical calculations, results calculated from the Semie model⁽⁵⁾, ANSYS^® shell 63, ANSYS^® shell 163 and LISA^® software are shown in Table 1. The results using the five FE codes suggested in Section 3 are also presented. The displacement at the center of the square plate calculated by the theoretical formula (Eq. (21)) was ‒3.47μm. Semie model using triangular FE showed a displacement of ‒4.81μm when modeling a square plate at 100×100. It was 38.6% error compared with the value predicted by the theoretical Eq. (21), and it was also reported that the error was ‒42.4% and 384.15%, respectively, when using ANSYS^® shell 63 and ANSYS^® shell 163 element ⁽⁵⁾. However, in Sec 3.1, the result using the triangle FE codes formulated by the methodology proposed by Semie was ‒3.50 μm, which was within 0.9% of the theoretical value. The FE codes presented in Sec 3.1 ~ Sec 3.3 which were formulated according to the Kirchhoff thin plate theory showed a maximum error of 0.9% with the theoretically predicted value. Therefore, there would be some errors in the FE codes programmed by Semie, and it was assured that the analys using ANSYS^® software was also applied incorrectly. Consequently, it showed that the FEM codes formulated according to the Kirchhoff thin plate theory were suitable for use. The analysis result by the commercial software LISA^® version 8.1.0 beta version was 16.14%. In case of Mindlin-Reissner thick plate theory considering shear deformation, it showed 8.6% error with Kirchhoff thin plate theory. On the other hand, 3 dimensional plate shell FE codes showed 19.02% error with Kirchhoff thin plate theory.

Table 1

Comparison of displacement at center of square plate under distributed load 50000Pa with fixed boundary condition by several methods(Unit: μm)

5. Comparison of Computation Time for System Matrix According to Solving Methods

It can be classified by two regimes to obtain the solution of the system equation, Ax=b. One is full matrix solving methods with A(i, j) form of square coefficient matrix A and the other sparse matrix solving methods with coefficient matrix A as sparse form.

5.1 Full Matrix Solving Methods

The FE codes developed in Section 3 initially were coded in OCTAVE^® language because of necessity of graphic postprocessing in the coding stage. OCTAVE^® intrinsic operator ‘\’ was adopted to obtain the vector x in the system equation Ax=b. But, it needed a long computation time, and codes were translated to FORTRAN language in order to improve its computing performance. In order to find the solution vector in the system equation Ax=b, which was the form [K]{W}={F} used in the structural displacement FE codes in Section 3, and compare computing performance, the following eight full matrix solving methods including four methods described in reference⁽¹⁵⁾ were examined.

(1) Gauss elimination method with partial pivoting algorithm, direct solving method.
(2) Successive over-relaxation (SOR(ω=1.95)) algorithm, iterative solving method.
(3) Conjugate gradient (CG) Algorithm, iterative solving methods, tolerance=10^-11.
(4) Preconditioned conjugate gradient algorithm with diagonal preconditioner (PCG-D), iterative solving method, tolerance=10^-11.
(5) Direct Gauss elimination method with partial pivoting (DGESV), FORTRAN MKL^® library.
(6) Iterative Gauss elimination method with partial pivoting (DSGESV), FORTRAN MKL^® library.
(7) OCTAVE^® ‘\’, full matrix solver (OCT-F).
(8) MATLAB^® ‘\’, full matrix solver (MAT-F).

Methods(1) to (4) showed that the same algorithmic codes were simply translated in OCTAVE^®, MATLAB^®, and FORTRAN languages, respectively, rather than using specialized built-in functions in each language. There were two classes of iteration methods⁽¹⁶⁾. One was stationary iteration methods or fixed-point iteration methods like as Jacobi, Gauss-Seidel, SOR and the other was projection methods, which was called Krylov subspace method, like as CG, PCG, generalized minimum residual method(GMRES). Projection methods were required that approximation solution x extracted from the subspace K should make residual x=b‒Ax orthogonal to the subspace L. In the first stage, the system coefficient matrix A used in the simulation was a symmetric positive definite (SPD) matrix, a tridiagonal matrix with a bandwidth of 1+2 × nx, with one diagonal with four off-diagonal terms added⁽¹⁷⁾. Nx was the distance between diagonal term and element locations and n was the size of matrix A.

n x = n A i, i = 4 × 1 + n x 2 A i, i + 1 = A i + 1, i = A i, i + n x = A i + n x, i = - 1 + n x 2, i = 1,2 …, n

The vector b is obtained by assuming x(1:n) = 1 in Ax=b. The personal computer (PC) used was Intel^® Core i7 4790K@4.0GHz processor, 32GB RAM, and operating on Windows 7 64bit. FORTRAN codes were compiled using Intel^® Parallel Studio XE 2015 MKL compiler 64bit. In the case of 2000DOF, the calculation time ratio was 27 883 sec: 144.24 sec:32.26 sec = 864 : 4 : 1 when the Gauss elimination method(solid line) was excuted with OCTAVE^®, MATLAB^®, and FORTRAN language, respectively (Fig. 5). In the case of the PCG-D algorithm, it showed OCTAVE^® : MATLAB^®:FORTRAN=21436sec:70.31sec:1.12 sec=19139:63:1. When the Gauss method and the SOR method were used, the calculation time ratio of OCTAVE^®/MATLAB^® was 193 and 197, CG method and PCG-D was 355 and 305, respectively. That is, it showed that the computing speed was 193 to 355 times faster to operate in MATLAB^® than in OCTAVE^® for same algorithmic codes. In addition, executing similar codes in FORTRAN was 4 to 63 times faster than running in MATLAB^®, and 831 to 2269 times faster than in OCTAVE^® between 200 DOF and 10000 DOF. The PCG-D method in FORTRAN showed the fastest calculation time.

Fig. 5

Comparison of computing time using full matrix solvers for OCTAVE®, MATLAB® and FORTRAN codes for 5 model cases, marker: Ｏ: OCTAVE®, ◇: MATLAB®, □: FORTRAN, linestyle: -: Gauss, --: SOR, -.: CG, ▫▫: PCG-D, range: 200 DOF ~ 10 000 DOF (Unit: sec)

Fig. 6 showed the computing time when the operator ‘\’ built in OCTAVE^® and MATLAB^®, the FORTRAN MKL library direct solver DGESV and iterative solver DSGESV were used. The shortest computation time for 2000 DOF was 0.051 sec for the direct method DGESV and 21436sec (≒10hours) for the slowest case, OCTAVE^®-PCG-D (marker: blue ●). It showed OCTAVE^®-PCG-D:MATLAB^®-PCG-D: FORTRAN-PCG-D:OCTAVE ‘\’:MATLAB ‘\’:DGESV =21436:70.31:1.12:0.637:0.065:0.051=420314: 1379:22:12.5:1.27:1. In the case of MKL FORTRAN library, iteration method was faster than direct method at 5000 DOF or more. In view of above discussions, it showed that the program language and solving method to be used were very important in computing time of system equation. The intrinsic operator ‘\’ used in MATLAB^® was shown to be an effective method. Using the MATLAB^® ‘\’ (-◆-) showed 1080 times faster than the computing time of the PCG-D programmed in MATLAB^® language(◆) in case of 2000DOF. For 10000DOF system equation, FORTRAN-PCG-D was 69.3sec (Fig. 5), DGESV 4.5sec, OCTAVE^® ‘\’ 104.4sec, and MATLAB^® ‘\’ 3.3sec(Fig. 6), respectively.

Fig. 6

Comparison of computing time using full matrix solvers, marker-line: ●: OCTAVE®-PCG-D, Ｏ-: OCTAVE® ‘\’, ◆: MATLAB®-PCG-D, ◇-: MATLAB® ‘\’, ■: FORTRAN-PCG-D, □-: DSGESV, □--: DGESV, range: 200 DOF ~ 10 000 DOF (Unit: sec)

So far we have compared the computation times for the cases where the coefficient matrix A was ideally diagonally banded SPD matrix. Gauss and SOR method, which was stationary iteration method with slow computation speed, were excluded from comparison. Using the data obtained from the FEM codes in Section 3.2, the computational speeds for the system equations were compared. In the formal FEM, the stiffness matrix [K] typically exhibited the characteristics of SPD matrix⁽¹⁸⁾. The stiffness matrix and the force vector were obtained using FEM codes using the rectangular FE codes in Sec. 3.2. The square plate were modeled to 20×20, 40×40, 80×80, 100×100, 140×140 FE meshes, respectively, which was called Ex1, Ex2, Ex3, Ex4 and Ex5. Since the DOF per node was 3, each model size was 1323, 5043, 19683, 30603, 59643 DOF, respectively. The specifications, boundary conditions and force conditions of the square plate were described in Section 4. The stiffness matrices and force vectors were obtained by applying fixed boundary conditions to the four sides of the square plate. The stiffness matrices and force vectors for the five data sizes were output to data files and then used in each solving method programs. Fig. 7 and Fig. 8 showed the stiffness matrices before and after applying the boundary condition to the system equation of the case of square plate (4×4) (75DOF). The modified stiffness matrix pattern (Fig. 8) was a general matrix form and obtained by applying fixed boundary conditions to the stiffness matrix pattern in Fig. 7.

Fig. 7

Pattern of stiffness matrix K of 75×75 before applied boundary condition, white part: ‘0’ value

Fig. 8

Pattern of stiffness matrix K of 75×75 after applied boundary condition, white part: ‘0’ value

Table 2 showed the calculation times obtained using the following seven solving methods for the five matrix size models. Methods(1) ~ (4) were compiled with Intel^® Parallel Studio XE 2015 FORTRAN compiler. Fig. 9 showed the computation time results for each solver.

Table 2

Comparison of computing time using full matrix solvers for 5 models(Unit: sec)

Fig. 9

Comparison of computing time using full matrix solvers for 7 methods, marker-line: □-: DSGESV, □--: DGESV, □-.: CG, □▫▫: PCG-D, Ｏ-: OCTAVE® ‘\’, ◆- : MATLAB® ‘\’, ▲-: LISA®, range: 1323 DOF ~ 59 643 DOF (Unit: sec)

(1) CG codes⁽¹⁵⁾, iterative method.
(2) PCG-D codes⁽¹⁵⁾, diagonal preconditioner, iterative method.
(3) DSGESV, iterative Gauss elimination method with partial pivoting.
(4) DGESV, direct Gauss elimination method with partial pivoting.
(5) OCTAVE^® operator ‘\’, OCT-F.
(6) MATLAB^® operator ‘\’, MAT-F.
(7) LISA^® software, version 8.1.0 Beta version.

In the direct solution method DGESV (□ dashed line), it showed that the calculation time was nearly proportional to DOF³. In MATLAB^® operator ‘\’, it was 0.065sec at 1323DOF and 3443sec at 59643 DOF, the computation time was also proportional to DOF³. In Fig. 6, iterative method DSGESV showed faster calculation speed than direct method DGESV, but in Fig. 9, iterative method showed slower calculation speed. This was because the coefficient matrix A applied in Fig. 6 was a sparse SPD matrix and the matrix A used in Fig. 9 was a non-SPD matrix and a relatively dense matrix. In OCTAVE^® and MATLAB^®, the same operator ‘\’ was used. As the DOF increases, the calculation time difference between MATLAB^® and OCTAVE^® becomes larger. When calculating the FE model of 19683DOF (Ex3), the computation time ratio was 55 times as 47 seconds in MATLAB^® and 2581 seconds in OCTAVE^®. In the case of 59643 DOF (Ex5), it was 3443sec for MATLAB^® ‘\’ and 1218 sec for DGESV routine. There was a rapid increment in computing time at 59643DOF for DSGESV. The reason why the computing time using DSGESV, an iterative method, surged compared with the direct method DGESV, was that it was caused by a reduction in efficiency caused by using virtual memory due to lack of CPU memory. Unlike DGESV, DSGESV, an iterative method, required additional workspace memory equal in size to system matrix A. To calculate an Ex5 model using DSGESV, we needed at least 56GB of memory, but used PC had hardware specifications with 32GB of RAM, so it needed virtual memory and so, it utilized hard disk as a virtual memory source, which resulted in a slower overall memory accessing time than using only CPU memory in direct solver DGESV. That is, by using hard disk memory rather than CPU memory as swap memory, the DOF was increased 1.9times in Ex5 compared to Ex4, but the calculation time was 31236sec which was 155 times larger than 202sec in Ex4. Therefore, in the case of Ex5, 59643DOF, DSGESV was needed 8hours 41minutes (31236sec), MATLAB^® solver ‘\’ 57minutes (3443sec) and DGESV took 20min (1218sec). In all cases, the direct solver, the DGESV routine showed the fastest calculation speed. For the commercial software LISA^®, the solving time was 106sec. A comparative analysis of Table 2, Fig. 5 and Fig. 6 revealed interesting facts. It was summarized the calculation time for 5000 DOF in Table 3.

Table 3

Comparison of computing time using full matrix solvers for general matrix (Ex2) and sparse SPD matrix(Unit: sec)

Direct solver DGESV showed almost the same calculation time in both cases. The iterative solvers, CG and DSGESV, showed 4-8 times faster computing time with the SPD coefficient matrix in system equation. This was because the DOF was almost the same, but the SPD matrix was a sparser matrix than Ex2. Even when the MATLAB^® operator ‘\’ was used, the calculation time was half shorter in the sparse matrix SPD. As an interesting result, the direct method DGESV in Ex2 showed faster calculation speed than the iterative method DSGESV, but the iterative method showed faster calculation time than the direct method in sparse matrix SPD. That is, when the number of non-zero coefficients was small, the iteration method showed faster calculation speed than the direct method.

5.2 Sparse Matrix Solving Method

In Section 5.1, full-matrix solving methods were used to examine the computing times in system equation for the SPD coefficient matrix (Figs. 5 ~ 6, 200 DOF ~ 10000 DOF) and for the stiffness matrix obtained in the structural displacement FEM (Table 2, 1323 ~ 59643 DOF). Practically, there were limitations on the size of matrix that could be solved due to full matrix memory space. Also, since the matrix operation was performed on the ‘0’ values in the coefficient matrix A (Figs. 7 ~ 8), lots of calculation time was required for full matrix solvers. If there were many ‘0’ values in coefficient matrix A, using a sparse solver would be benefit from memory savings and computation speed. In order to use the sparse matrix solving methods, the system matrix A must be converted from full matrix form to a sparse matrix one. In OCTAVE^® and MATLAB^®, the sparse matrix solver used the ‘\’ operator, which was the same form as the full matrix solver. However, the system matrix A must be transformed into a sparse matrix form. In this paper, it was adopted in order to make the system coefficient matrix A a sparse matrix form, one-based indexing CSR (compressed sparse row) format^{(19 ~ 21)} in FORTRAN codes and CSC (compressed sparse column) format in MATLAB^® and OCTAVE^® codes, respectively. NNZ (number of non-zero) was defined the number of non-zero coefficients in the coefficient matrix A whose size was n×n is. Sparsity (%) is defined as sparsity (%)=NNZ/(n×n)×100. Table 4 and Table 5 showed the DOF, NNZ, and sparsity(%) of the data to be used in this section and the amount of memory required to store matrix A in full matrix and sparse matrix. Table 4 showed that the 50000 DOF (#7) model required 20GB of memory to process a full matrix of coefficients A, but a sparse matrix required only 2MB, which was 0.01% of the original memory requirement. 300000DOF (#10) required 720GB of memory for the full matrix, but it could not be handled by general PCs. However, if sparse matrix form is used, only 12MB of 0.002% of raw data is needed. Table 6 compared computing times of full matrix solvers and sparse matrix solvers for from #1 to #5 data models.

Table 4

Comparison of memory sizes of full and sparse for SPD matrix, NNZ, sparsity (%)

Table 5

Comparison of memory sizes of full and sparse matrix for FEM stiffness matrix K, NNZ, sparsity (%)

Table 6

Comparison of computing time of Ax=b using full matrix solvers and sparse matrix solvers, (Unit: sec)

When using the full matrix solver, operator ‘\’ in models with a range of 200 DOF ~ 10000 DOF, it showed computing time ratio of OCTAVE^®/MATLAB^® was 2 ~ 32. However, using sparse solver showed OCTAVE^®/MATLAB^®≈1. The method using sparse solver ‘\’ showed the same performance of MATLAB^® and OCTAVE^®. In the 5000DOF model (#4), the sparse solver (0.006sec) of MATLAB^® and OCTAVE^® was 122 times faster than the full matrix direct solver DGESV (0.732sec). Sparse PCG-D (0.016sec) showed 664 times faster than full matrix PCG-D (10.62sec). Sparse solvers were used to calculate system equation with the SPD coefficient matrix from 200 to 20000000 DOF using eight sparse matrix methods.

(1) PARDISO^{(19,22 ~ 24)}: PARallel sparse Direct and multi-recursive Iterative linear SOlvers
(2) OCTAVE^® ‘\’, OCT-S.
(3) MATLAB^® ‘\’, MAT-S.
(4) CG^{(15,16,21,24)}
(5) PCG-D^{(15,16,21,24)}.
(6) PCG-split-D^{(15,16,21,24)}, PCG with diagonal split preconditioner.
(7) PCG-IC(0)^{(15,16,21,24,25)}: PCG with incomplete Cholesky decomposition with no fillin preconditioner
(8) PCG-ILU(0)^{(15,16,21,24)}: PCG with incomplete LU decomposition with no fillin preconditioner

In the case of the 2000DOF mentioned in Section 5.1, the PCG-D (marker: black ●) implemented with OCTAVE^® language was 21436sec, and the OCTAVE^® full matrix solver ‘\’ (marker: red ●), 0.637sec and OCTAVE^® sparse solver ‘\’ (-●-), 0.0018sec, respectively(Fig. 10). The calculation time ratio was 12000000: 425:1. OCTAVE^® sparse solver ‘\’ showed twelve million times faster calculation time than PCG-D translated by OCTAVE^® language. The computing time using the sparse solver PCG-D compiled with FORTRAN was 3.3×10^-3sec, which was 1.8 times slower than the OCTAVE^® sparse solver ‘\’. PCG-D adopted the diagonal terms of the coefficient matrix A as the preconditioner M and PCG-split-D adopted the preconditioner as the (diagonal terms)^1/2 of the coefficient matrix A, respectively. Fig. 10 and Fig. 11 showed that PCG-split-D was more efficient than PCG-D because the iteration number in PCG-split-D is 40% smaller than in PCG-D above 500,000 DOF model. PCG-IC(0) showed similar calculation times as PCG-D and PCG-split-D, even though the number of iterations was 1/2 to 1/3 times that of PCG-D or PCG-split-D (Fig. 11). Thus, it meant that the computation amount per iteration of PCG-D or PCG-split-D was smaller than PCG-IC(0). In the SPD coefficient matrix data, iteration number in PCG-D showed 140% greater than the number of iteration in PCG-split-D and needed 140% of the calculation time. In the case of 2000000 DOF (#13), the calculation time in OCTAVE^® ‘\’ was 47.8sec, but the time to read the data file took 6301sec. It took 114 sec for MATLAB^® and 5.2sec for FORTRAN PCG to read the data file #13. Saad showed that when the coefficient matrix A was the SPD matrix and the preconditioner M was incomplete Cholesky product, A≈M = LL^t, it was equal to the number of iterations when using M as the preconditioner or when splitting it into L and L^t for PCG method⁽²¹⁾.

Fig. 10

Comparison of computing time according to sparse solvers for SPD matrix, marker-line: Δ-: PARDISO, Ｏ-: OCT-S ‘\’, ◇-: MAT-S ‘\’, +-: CG, □-: PCG-D, □--: PCG-split-D., □-.: PCG-IC( 0), □▫▫: PCG-ILU(0), ●: OCTAVE®-PCG-D, full matrix solver, ●: OCTAVE® ‘\’ -full matrix solver, ◆: MATLAB®-PCG-D, full matrix solver, ◆: MATLAB® ‘\’ -full matrix solver, range: 200 DOF ~ 20 000 000 DOF (Unit: sec)

Fig. 11

Comparison of iterations according to sparse solvers for SPD matrix, marker-line: Ｏ-: CG , ◇-: PCG-D, □-: PCG-split-D, +-: PCG-IC(0), range: 200 DOF ~ 20 000 000 DOF, (Unit: number)

On the other hand, when the coefficient matrix A was a general matrix, a preconditioner was obtained with incomplete LU factorization instead of incomplete Cholesky decomposition. Split GMRES using L and U as preconditioners in GMRES instead of PCG was most effective when the coefficient matrix A was nearly symmetric and was not effective from using preconditioner M = LU in other cases⁽²¹⁾. In this paper, GMRES was applied to FEM data, which showed that it was inferior to PCG-D or PCG-split-D in terms of calculation speed. Especially, when using FEM data, it showed a tendency to diverge above 1323 DOF (Ex1 model). FGMRES (RCI Flexible Generalized Minimal RESidual method with ILUT Preconditioner)⁽²⁶⁾ and MGMRES(restarted GMRES) suggested by Burkardt⁽²⁷⁾ were less likely to diverge than GMRES. However, according to the characteristics of the coefficient matrix A, the user had to determine the inner loop and the outer loop number in advance, and the divergence or convergence was determined according to these values. Also, if the number of inner loop was large, the solution converged rapidly but the memory to be memorized increased and the memory efficiency decreased. Here, interestingly, Saad⁽²¹⁾ described the difference in computational complexity in PCG when choosing preconditioner using incomplete Cholesky decomposition for SPD matrix. In this paper, we showed that the PCG could also be applied to non-SPD matrices and the same description applied to the SPD matrix when we adopted diagonal preconditioner or split diagonal preconditioner and it could be especially emphasized to non-SPD coefficient matrix. The PCG method with incomplete Cholesky decomposition preconditioner was known to apply only to SPD matrix⁽²⁸⁾. If the coefficient matrix was non-SPD matrix, Cholesky decomposition could not be applied; instead, LU decomposition applicable to a general matrix could be applied. The PCG-IC(0) could be applied only when the coefficient matrix was SPD, so it could not be applied to the FEM system matrix equation applied the actual boundary condition as in Fig. 8.

In Table 5 (1323 DOF ~ 59 643 DOF), the PCD-IC(0) could not be applied because the matrix was non- SPD. In order to solve the general matrix in the PCG method, the preconditioner could be obtained by incomplete LU decomposition instead of incomplete Cholesky decomposition. We could use PCG-ILU(0), which used L and U obtained by incomplete LU decomposition with no fillin as the preconditioner, to perform the same size LU decomposition using the coefficient matrix A. The preconditioner ILU(0) was efficient in terms of memory because the size and position of the matrix were the same as the coefficient matrix A, but there was a restriction on its use because it tended to diverge. In case of non-SPD coefficient matrix as in Fig. 12, PCG-ILU(0) showed divergence when the size of the matrix was larger than 1323 DOF. As an alternative, we could use PCGILU(τ) by obtaining L and U with Crout incomplete LU decomposition. However, the size of L and U obtained by Crout incomplete LU decomposition was large, which was not useful in terms of memory. The Crout method could be used to solve the ILU when using the Crout ILU method. However, when the Crout method was used to obtain the preconditioner LU, it showed NNZ (ILU(τ)) >> NNZ(A) (Table 7). Therefore, it could not be applied to a large size non- SPD coefficient matrix. On the other hand, as shown in the appendix, PCG-D or PCG-split-D could perform very fast convergences if coefficient matrix A was not asymmetric. In Fig. 12, the full matrix direct solver DGESV (marker: red ▼) was 1218 sec and the MATLAB^® ‘\’ -full matrix solver (marker: blue ◆) was 3443 sec in the case of 59 643 DOF.

Table 7

Comparison of number of non-zero of coefficient matrix A, incomplete LU with no fillin (ILU(0)) and incomplete LU with Crout ILU(τ), τ: drop tolerance, τ=10-9

Fig. 12

Comparison of computing time using sparse solvers for non-SPD matrix, marker-line: Δ-: PARDISO, Ｏ-: OCTAVE® ‘\’, ◇-: MATLAB® ‘\’, +-: CG, □-: PCG-D, □--: PCG-split-D., □▫▫: PCG-ILU(0), ◆: MATLAB® ‘\’ - full matrix solver, ▼: DGESV-full matrix solver, range: 75 DOF ~ 59 643 DOF (Unit: sec)

The MATLAB^® sparse matrix solver ‘\’ (marker: blue -◇-) was 0.592 sec. Therefore, when using the MATLAB^® operator ‘\’, computing time ratio of the full matrix solver / sparse matrix solver = 5816 times. That is, the sparse MATLAB^® solver was about six thousand times faster than the full matrix MATLAB^® solver.

Ratio of number of iterations in the sparse PCG method was PCD-D / PCD-split-D = 1.4 in the SPD matrix (Fig. 11). When the coefficient of the matrix was non-SPD (Fig. 13), it showed ratio of number of iterations of PCD-D / PCD-split-D = 4 and the convergence speed of PCG-split-D was 4 times faster than PCG-D. That is, in the case of the asymmetric matrix, PCG-split-D showed shorter computation time than PCG-D. The PCG-split-D, 4.805 sec, was 4.7 times faster than the PCG-D, 22.792 sec, in the case of the 59 643 DOF, where the coefficient matrix was a little asymmetric matrix (Fig. 12). When the coefficient matrix was SPD (Fig. 10), PCG-split- D was compared with 40% faster than PCG-D. In sparse solvers, the speed of calculation depended on the non-zero number of the coefficient matrix and the non-SPD matrix, not on the size of the matrix. In Fig. 14, it was compared the computing times according to NNZ of SPD and non-SPD coefficient matrices. In SPD coefficient matrix 20 000 000 DOF, PARDISO (mark: ▼) and MATLAB^® sparse ‘\’ (mark: ▲) solvers showed 193 sec and 98 sec, respectively. Using the FEM stiffness matrix, non-SPD matrix, for example, Ex5 (59 643 DOF, NNZ = 1 564 440), system equation was solved by PCG-split-D, resulting in 4.805 sec, which was nearly equivalent time to solve 200,000 DOF SPD matrix (#9) problem. In other words, the sparse solver was not dependent on the size of the coefficient matrix A but on the number of non-zero of the coefficient matrix A. When NNZ was the same, SPD matrix converged faster than the non-SPD matrix in CG method and PCGD converges equally to SPD or non- SPD matrix. However, in case of equal NNZs of coefficient matrices, it showed that PCG-split-D showed faster convergence in non-SPD matrices than SPD ones.

Fig. 13

Comparison of iterations using sparse solvers for non-SPD matrix, marker-line: Ｏ-: CG, ◇-: PCG-D, □-: PCG-split-D, range: 75 DOF ~ 59 643 DOF, (Unit: number)

Fig. 14

Comparison of computing time according to NNZ of sparse solvers for SPD and non-SPD matrix, marker-line: Ｏ-: CG-SPD, ◇-: PCG-DSPD, □-: PCG-split-D-SPD, Ｏ--: CG-FEM, ◇--: PCG-D-FEM, □--: PCG-split-D-FEM, range: 200 DOF ~ 20 000 000 DOF (Unit: sec)

The NNZ of the non-SPD matrix Ex5 (59 643 DOF, NNZ: 1 564 440) was equal to the NNZ of SPD matrix #10 (300 000 DOF, NNZ: 1 498 902). However, the time required for calculating using Ex5 using PCG-split-D was equal to the time for solving #9 (200 000 DOF, NNZ: 999 104). The reason why was that the number of iterations when using PCG-split-D for SPD matrix #9 and #10 was 1533 and 1846, respectively. However, when applying PCG-split-D to non-SPD matrix Ex5, the number of iterations was 2929. The non-SPD matrix Ex5 showed 48 828 and 14 475 iterations, respectively, when CG or PCG-D was used to solve system equation. Consequently, the computing time was the shortest since the number of iterations was smaller than CG or PCG-D when solving non- SPD matrix Ex5 with PCG-split-D.

6. Conclusions

(1) In this paper, we proposed physical coordinated based 3-node triangular FE, 4-node rectangular FE, and 4-node quadrilateral FE codes according to the Kirchhoff thin plate theory.

(2) FE codes using isoparametric quadrilateral FE and 3-D plate shell FE elements based on the Mindlin- Reissner thick plate theory were also presented.

(3) We found that there would be some errors in the FE codes presented by Semie⁽⁵⁾ through validation test using a three-node triangular FE codes developed in this paper.

(4) The displacements were compared at the center point of a square plate with four sides fixed and under constant pressure. Three FE codes, which were physical coordinated based 3-node triangular FE, 4- node rectangular FE, and 4-node quadrilateral FE program developed based on the Kirchhoff plate theory, showed within a deviation of 0.9% from the theoretical value obtained by the Levy method. Therefore, three developed FEM codes were proved to be adequate to use.

(5) We also presented results using natural coordinates based isoparametric quadrilateral FE and 3-D plate shell FE codes based on the Mindlin-Reissner thick plate theory. The displacement was greater than 8 % of the value obtained by the Kirchhoff thin plate theory. The commercial software LISA analysis results were also compared and it was 16% higher than Kirchhoff’s theory without taking shear strain into account.

(6) The coefficient matrix from the actual FEM analysis was a non-SPD matrix. In case of 59 643 DOF non-SPD coefficient matrix, it showed that the calculation time was 1218 sec of the MKL library full- matrix direct solver DGESV, 4.805 sec in sparse PCG-split-D and 0.590 sec in sparse MATLAB ‘\’, respectively, indicating that the sparse solvers were superior to the full matrix solvers. Especially, sparse PARDISO took 0.345 sec.

(7)Sparse solvers showed excellent memory saving and computation speed when coefficient matrix was sparse when solving large system equations.

(8)FGMRES and MGMRES were less likely to diverge than GMRES when solving system equations with non-SPD coefficient matrix. However, according to the characteristics of the coefficient matrix, the number of inner loop and the number of outer loop must be predetermined by the user, and divergence or convergence was determined according to these values. Also, if the iteration number of inner loop was large, the solution vector converged rapidly but the memory to be memorized increased and the memory efficiency decreased. Calculation time was needed more than sparse PCG-D or PCG-split-D.

(9)In the case of the SPD coefficient matrix, the calculation performance of PCG-D and PCG-split-D was similar, but in the case of non-SPD coefficient matrix, PCG-split-D showed faster computation time than PCG-D. In the non-SPD matrix, 59643DOF, PCG-split-D was five times faster than PCG-D due to faster convergence.

(10)OCTAVE^® showed that reading data files was very slower than in MATLAB^® or FORTRAN.

(11)It showed that PCG-D and PCG-split-D almost were required the same calculation amount per iteration to solve system equations with SPD coefficient matrix.

(12)If the coefficient matrix was asymmetric, it showed that it may be happened to inaccurate results or divergency in PCG-D, PCG-split-D and PCGILU(0). The PCG-ILU(τ) using Crout ILU was accurate, but required more memory than ILU(0), which may not be suitable for large-size problem.

(13)Sparse solver was not dependent on the size of the coefficient matrix A but on the number of non-zero of the coefficient matrix A.

(14)In case of system equation with large non-SPD coefficient matrices, it showed that computing times were ordered with PARDISO, MATLAB^® sparse operator ‘\’ and sparse PCG-split-D in short time order. But, in case of large SPD matrix, for example, 20000000DOF, sparse MATLAB^® ‘\’, PARDISO, and sparse PCG-split-D exhibited 98sec, 193sec, and 5632sec, respectively.

References

Holmström, F., (2001), Structure-acoustic Analysis Using BEM-FEM implementation in MATLAB, M.S. Thesis, Lund University, Lund, Sweden.
Irvine, T., (2011), Plate Bending Frequencies via the Finite Element Method with Rectangular Elements, Revision A.
Irvine, T., (2011), The Natural Frequency of a Rectangular Plate with Fixed-free-fixed-free Boundary Conditions.
Irvine, T., (2011), The Natural Frequency of a Rectangular Plate with Fixed-fixed-fixed-fixed-fixed Boundary Conditions Revision B.
Semie, A. G., (2010), Numerical Modelling of Thin Plates Using the Finite Element Method, M.S. Thesis, Addis Ababa University, Addis Ababa, Ethiopia.
Ferreira, A. J. M., (2009), MATLAB Codes for Finite Element Analysis, New York, NY, Springer.
Kwon, Y. W., and Bang, H., (2000), The Finite Element Method Using MATLAB, 2nd edition, Boca Raton, FL, CRC Press.
Timoshenko, S. P., and Woinowsky-Krieger, S., (1989), Theory of Plates and Shells, New York, NY, McGraw-Hill.
Batoz, J.-L., Bathe, K.-J., and Ho, L.-W., (1980), A Study of Three-node Triangular Plate Bending Elements, International Journal for Numerical Methods in Engineering, 15(12), p1771-1812. [https://doi.org/10.1002/nme.1620151205]
Batoz, J. L., and Tahar, M. B., (1982), Evaluation of a New Quadrilateral Thin Plate Bending Element, International Journal for Numerical Methods in Engineering, 18(11), p1655-1677. [https://doi.org/10.1002/nme.1620181106]
Bischoff, M., (2008), Modeling of Shells with Three-dimensional Finite Elements, Proceedings of the 6th International Conference on Computation of Shell and. Spatial Structures, p12-15.
Petyt, M., (1990), Introduction to Finite Element Vibration Analysis, Cambridge, England, Cambridge University Press. [https://doi.org/10.1017/cbo9780511761195]
Reddy, J. N., (1997), Mechanics of Laminated Composite Plates, New York, NY, CRC Press.
Knight Jr., N. F., (1997), Raasch Challenge for Shell Elements, AIAA Journal, 35(2), p375-381.
Burden, R. L., and Fairs, J. D., (2011), Numerical Analysis, 9th edition, Pacific Grove, CA, Brooks/Cole.
Du, L., and Liu, S., (2015), Study on Preconditioned Conjugate Gradient Methods for Solving Large Sparse Matrix in CSR Format.
Burgerscentrum, J. M., (2012), Iterative Solution Methods.
Nikishkov, G. P., (2004), Introduction to the Finite Element Method (Lecture Notes), Aizu-Wakamatsu, Japan, University of Aizu.
Intel, (2014), Intel^® Math Kernel Library (Version 11.2) Reference Manual.
Frydendall, J., (2009), A Parallel Implementation of a Finite Element Solver, IMM-Technical Report-2009-10, Lyngby, Denmark, DTU Informatics.
Saad, Y., (2003), Iterative Methods for Sparse Linear Systems, 2nd edition, Philadelphia, PA, Society for Industrial and Applied Mathematics. [https://doi.org/10.1137/1.9780898718003]
Schenk, O., and Gärtner, K., (2018), Parallel Sparse Direct and Multi-recursive Iterative Linear Solvers User’s Guide Version 6.1.0.
Shenk, O., and Gärtner, K., (2003), Solving Unsymmetric Sparse Systems of Linear Equations with PARDISO, Preprint Submitted to Elsevier preprint.
W. Layton, and M. Sussman, (2014), Numerical Linear Algebra, Morrisville, NC, Lulu Press.
Stahel, A., (2014), Numerical Methods.
http://sep.stanford.edu/sep/claudio/Research/Prst_ExpRefl/ShtPSPI/intel/mkl/10.0.3.020/examples/solver/source/dcsrilut_exampl2.f, .
Burkardt, J., (2019), MGMRES: Restarted GMRES Solver for Sparse Linear Systems, https://people.sc.fsu.edu/~jburkardt/m_src/mgmres/mgmres.html.
Mathworks, (2018), pcg (Preconditioned Conjugate Gradients Method), https://kr.mathworks.com/help/matlab/ref/pcg.html#f93-998415.

Appendix

A study on the applicability of PCG method by the asymmetry of coefficient matrix in system equation

We applied PCG method to SPD coefficient matrix #1~#16 (Table 4, Fig. 10) and non-SPD matrix Ex1~Ex5 (Table 5, Fig. 12) similar to symmetric matrix as a whole. And we got solutions within an acceptable accuracy. However, CG or PCG methods should be used with caution when solving system equations with non-SPD matrices other than SPD matrices. We investigated two examples like as followed.

Example 1: In the system equation, Ax=b, if the coefficient matrix A is a general matrix with strong asymmetry.

A = 5317626004027300004510307, x = 11111 t, b = 221212911 t,

CG, PCG with ILU(0) or PCG with diagonal preconditioner (PCG-D) and PCG with split diagonal preconditioner (PCG-split-D) showed incorrect results but only PCG with Crout version showed correct results (Table A1).

Example 2: If the coefficient matrix is not asymmetric in the system equation Ax=b,

A = 5317636204027337034264327, x = 11111 t, b = 2215151611 t,

CG only showed the wrong results (TableA1).

Table A1

Analysis results of PCG method according to characteristics of coefficient matrix A

Seok-Tae Park, B.S. in Hanyang University, Mechanical Eng., 1984., M.S. in KAIST, Mechanical Eng., 1986., Ph.D. in Ajou University, Systems Eng., 1999. 1986 ~ 1989., KAIST (KIST) Mechanical Eng., researcher. 1989 ~ 1992., Ssangyong Motor Co., senior researcher. 1993 ~ 1999., IAE Automotive Technical Lab. principal researcher. 2000 ~ present, in Chungbuk Health & Science University Assistant Professor. His principal interest is the fields of sound quality index codes programming, noise control, and subjective evaluation, loudspeaker, horn and microspeaker system analysis and design, boundary integral analysis in Acoustics.

	Elem×elem	10×10	20×20	40×40
1	Theory, Eq.(21)	‒3.47	‒3.47	‒3.47
2	Semi⁽⁵⁾	-	-	‒4.81
3	ANSYS 63⁽⁵⁾	-	-	‒2.00
4	ANSYS 163⁽⁵⁾	-	-	‒16.80
5	LISA	‒4.01	‒4.02	‒4.03
6	Triangular FE	‒3.65	‒3.54	‒3.50
7	Rectangular FE	‒3.56	‒3.51	‒3.49
8	Quadri. FE	‒3.56	‒3.51	‒3.49
9	Mindlin-Reissner	‒3.74	‒3.76	‒3.77
10	Plate shell	‒4.12	‒4.13	‒4.13

DOF	Ex1	Ex2	Ex3	Ex4	Ex5
DOF	1323	5043	19683	30603	59643
Blank spaces mean omitted because of expecting excessive computing times. LISA®, 4 node quadrilateral element
(1) CG	1.5	133.9	5657.6	13543.5
(2) PCG-D	0.32	12.06	836.8	3078.7
(3) DSGESV	0.11	1.65	55.8	202.0	31235.7
(4) DGESV	0.02	0.75	31.5	120.8	1218.4
(5) OCTAVE	0.27	14.84	2580.9
(6) MATLAB	0.06	1.12	47.1	168.2	3442.8
(7) LISA	-	1.0	13.0	26.0	106.0

	DOF	Full	NNZ	Sparse	Spar.(%)
#1	200	.32MB	970	0.01MB	2.43
#2	1000	8MB	4934	0.04MB	0.49
#3	2000	32MB	9908	0.08MB	0.25
#4	5000	200MB	24856	0.20MB	0.10
#5	10000	800MB	49798	0.40MB	0.05
#6	30000	7.2GB	149652	1.20MB	0.02
#7	50000	20GB	249550	2.00MB	0.01
#8	0.1M	80GB	499366	3.99MB	5×10^-3
#9	0.2M	320GB	999104	7.99MB	3×10^-3
#10	0.3M	720GB	1498902	11.99MB	1×10^-3
#11	0.5M	2TB	2498584	19.99MB	5×10^-4
#12	1M	8TB	4997998	39.98MB	5×10^-4
#13	2M	32TB	9997170	79.98MB	3×10^-4
#14	5M	200TB	24995526	0.2GB	1×10^-4
#15	10M	800TB	49993674	0.4GB	5×10^-5
#16	20M	3.2PB	99991054	0.8GB	3×10^-5

	DOF	Full	NNZ	Sparse	Spar.(%)
Ex1	1323	14MB	29453	0.24MB	1.68
Ex2	5043	203MB	123381	0.99MB	0.49
Ex3	19683	3.10GB	505325	4.04MB	0.13
Ex4	30603	7.49GB	790713	6.33MB	0.08
Ex5	59643	28.46GB	1564440	12.52MB	0.04

	DOF		OCT.	MAT.	PCG-D	DSG.	DGE.
F and S in column 3 mean full matrix solver and sparse matrix solver, each. PCG-D: PCG with diagonal preconditioner, full matrix codes are from Ref.(26) and sparse matrix codes are modified from Ref.(25). *DSGESV and DGESV are full matrix solver.
#1	200	F	2.0×10^-3	0.	2.0×10^-3	5.0×10^-3	0.002
#1	200	S	1.0×10^-3	1.0×10^-3	0.	5.0×10^-3	0.002
#2	1000	F	6.3×10^-2	3.0×10^-2	1.6×10^-1	1.2×10^-2	0.01
#2	1000	S	1.0×10^-3	1.0×10^-3	0.	1.2×10^-2	0.01
#3	2000	F	6.4×10^-3	6.5×10^-2	1.1	6.9×10^-2	0.05
#3	2000	S	2.0×10^-3	2.0×10^-3	1.6×10^-2	6.9×10^-2	0.05
#4	5000	F	1.3×10¹	5.4×10^-1	10.6	4.1×10^-1	0.73
#4	5000	S	6.0×10^-3	6.0×10^-3	1.6×10^-2	4.1×10^-1	0.73
#5	10000	F	1.0×10²	3.3	63.4	2.5	4.52
#5	10000	S	1.4×10^-2	1.1×10^-2	4.7×10^-2	2.5	4.52

DOF	Ex2	SPD	Ratio, Ex2/SPD
DOF	5043	5000	Ratio, Ex2/SPD
CG	133.85	15.90	8.42
PCG-D	12.06	10.62	1.14
DSGESV	1.65	0.41	3.99
DGESV	0.75	0.73	1.02
MATLAB^® ‘\’	1.12	0.54	2.07

		X(1)	X(2)	X(3)	X(4)	X(5)
Exact sol.		1.	1.	1.	1.	1.
Ex1	CG	3473	-772	1143	-575	-1044
	PCG-D	10565	-3604	3729	-13531	3484
	PCG-split-D	10565	-3604	3729	-13531	3484
	PCG-ILU(0)	1.012	0.95	0.97	1.12	0.93
	PCG-Crout	1.000	1.000	1.000	1.000	1.000
Ex2	CG	0.998	1.001	1.00	1.01	1.000
	PCG-D	1.000	1.000	1.000	1.000	1.000
	PCG-split-D	1.000	1.000	1.000	1.000	1.000
	PCG-ILU(0)	1.000	1.000	1.000	1.000	1.000
	PCG-Crout	1.000	1.000	1.000	1.000	1.000