- Email: [email protected]
A conjugate gradient algorithm applied to plane stress-strain problems
COMPUTER METHODS NORTH-HGLLAND
1N APPLIED
A CONSULATE
MECHANICS
AND ENGINEERING
40 (1983) 127-135
GRADIENT ALGORITHM APPLIED STRESS-STRAIN PROBLEMS
Received 7 Aprit f981 Revised manuscript received 11 February
TO PLANE
1982
The numerical methods of optimal control theory have recently been applied effectively to problems arising from the application of finite element methods to transonic and incompressible viscous flows. This paper uses the same approach to show that a conjugate gradient algorithm can solve stress-strain problems with far less computer storage than that required by Conventional methods.
The main methods used for the solution of plane stress and plane strain problems have been those based upon finite differences and finite elements. When finite differences are used a biharmonic equation fl]
is solved and the stresses obtained from the equations
Equation (1) is transformed in matrix form as
to a set of simultaneous,
algebraic equations which may be written
AX=B.
When using a finite element method [2f again a set of simultaneous must be solved. In this case they are given by
(3)
atgebraic equations
where K is the stiffness matrix, S is the dis~~a~eme~t vector and F is the force vector acting on the object. Clearly, in both methods the main work is to solve the high-order systems (3) or (4) which in ~45-~8~5/83/$3.00
@ 1983, Elsevier Science Publishers B.V. (worth-Holland~
128
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
theory is no problem but in practice has the difficulty of computer storage for the matrix A or K. The order of these coefficient matrices is equal to twice the number of nodal points of the mesh being used, and for any realistic problem this will be several hundreds. A number of methods are available for solving large numbers of linear, algebraic equations, and the choice of method depends on the core storage available. If storage is not a problem, then the Gauss elimination or the Cholesky decomposition are better than most methods. If it is not possible to store the whole coefficient matrix, then methods such as the Gauss-Seidel, over-relaxation [2] or incomplete Cholesky-conjugate gradient [33 may be suitable although, even here, all non-zero elements of the coefficient matrix have to be stored. In many situations insufficient core store is available, even to make these methods viable, and so one is forced to look for methods which bypass the solution of large scale systems such as (3) or (4). This paper develops such a method using a conjugate gradient algorithm. The conjugate gradient method and modifications to it has been used extensively in many areas of engineering since Hestenes and Stiefel [4] first described their idea in 1952. The solution of the load-deflection equations arising from finite element discretization has been the subject of very many papers over the last decade. In particular, the use of conjugate gradient techniques in conjunction with finite element discretization has generated a rich literature [s-16]. There have been reports in the literature of convergence problems arising from the use of conjugate gradient techniques when applied to finite element discretization. The algorithm used in this paper is that due to Polak and Rib&e for which convergence has been established, and a rate of convergence can be assessed. The authors have not experienced any convergence difficulties in the present investigations. The Polak-Rib&e algorithm coincides with the well-known Fletcher-Reeves algorithm when applied to a quadratic function [18].
2. The primal problem For conservative, holonomic, mechanical systems under Hamilton’s principle 1171 may be expressed as
the action of potential
forces,
where T and V are the kinetic and potential energies, respectively. In problems of equilibrium equation (5) reduces to Vdt=O.
(6)
Because the system is assumed conservative the total energy T + V is constant and thus the potential energy is independent of time t. Equation (6) then reduces to
sv=o, which is the well-known basic equation for equilibrium.
Q. YE, D.J. Bell, A corzjugate gradient algorithm for plane stress-strain problems
129
In problems of plane strain and plane stress we may write the total potential energy of the system as [2] v=
v,+ vz
(8)
where V, is the stored strain energy and V2 the potential energy of the external forces. The problem is then a variational one: to minimize the functional V with respect to displacements. These displacements of the object in the directions x and y we write as u(x, y) and ~(x, y), (x, y) E R, where R is the domain occupied by the object. The strain e can then be evaluated from du
(
e = (e,,, eyv,e,,>” = -
dv
-,
dx’ay
au
-
ay
dv t . dx >
+-
In the linear, elastic problem the relationship
(9) between stress and strain is
a=De is the stress, and where u = (aXX,a;,, uXxv)f .
(11) Here E* = E and it* = Y for plane stress problems or E* = E/(1 - Y’), Y* = v/(1 - Y) for plane strain problems. E is Young’s modulus, v is Poisson’s ratio, and they are both independent of u and e (and therefore of u and u), determined only by the material properties. Thus, for linear, elastic problems, the stored strain energy in the entire solution domain fz is
(12) The external forces involve the body forces x, t and the forces F., F, in the directions of x and y at the boundary C. The potential energy of the external forces is
v*= - j i, (it%+ %)dx
dy - j. (F’u + F,v)ds .
(13)
Thus, the total potential energy of the system is v=$
II
R
e’De dx dy - 1 I, (% + %)dx dy - 1 (F,u + F,u)ds .
(14)
Here u and v are functionals of x and y so the variational problem inherent in (7) must be solved in the Sobolev space H’(LJ). This means the feasible domain of control variables u and ~1is a Sobolev Space. This is a property which must be considered if we are to use a conjugate gradient method to solve the given extremum problem [3].
130
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
3. The numerical problem
For simplicity we here neglect the body forces. Then v=$
(F,u + F,u)ds
e’De dx dy -
(1%
.
In order to calculate the double integral in equation (15) we must evaluate the strain e through equation (9). This implies the calculation of &/&, dully c~v/& and dully. Many methods are available for this purpose; for example, by the finite element method we take triangular elements and assume the displacements u and v vary linearly over each element. Fig. 1 illustrates the mth eiement of a mesh. We assume the displacements at points (x, y) within the element are given by u(x, y) = G -+ czx -+ GY 3
zJ(X,y) = c4 + csx + csy *
The strain e, is then obtained in terms of the displacements
em= &-
at the vertices i, j, k as
B,&,,
(17)
m
where A, is the area of the element, 8, =
(Ui, Vi,
uj,
vj, Uk,
WCJ
VkY
and B, is the matrix b, 0 ai
0 a,
bi
bj 0 aj
0 aj
bk 0
bj
ak
0
(19)
Y
Ll+----F Fig. 1. Typical element of mesh.
X
131
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
Note that B, is only determined by the mesh and is independent the mesh is given, B, is a constant matrix. The total potential energy (15) now becomes
of the displacements.
When
where s is the number of elements of the mesh and 1 the number of nodes on the boundary. Given the displacements at the nodes of the mesh the total potential energy V can be calculated by equation (20). The variational problem given by (7) therefore has 2n independent variables Ui, L+(i = 1,2, . . . , n), where it is the number of nodes in the mesh. It will be solved in Euclidean space R2”. 4. A conjugate gradient algorithm In order to solve the problem of equation (7) by a conjugate gradient method the gradients of V with respect to the variables Ui, vi must be evaluated. From (20) we obtain %)-Fix, I
i=l,2
,....)
n
(21)
where Fix, Fiy are zero if i is not on the boundary. From the constraints at some nodes the displacements Ui and/or Ui are given. The corresponding derivatives L~v/b’u,and 3V/& therefore vanish. The algorithm used in this paper is the one due to Polak and Ribike [18]. The procedure for this method is as follows. Step 1. Select a normal displacement
u. and v. subject to the boundary condition. Calculate
V(uovo) by (20) and ?f(~ovo>
= K.~(u~v~)V$(U~V~>>
=
@v/au,av/a2))tu=w,o,w
by (21), where V, is the operator (a/aul, . . . , a/au,)’ and similarly for V,. Set i = 0, E > 0. Step 2. If (Vf(uovo)~< E, then stop; else go to Step 3. Step 3. Set h’d = g’6 =
hb = g; = -Vf(uov0) ) Step 4. Compute
-Vf”(uovo)
.
Ai > 0 such that
V(Ui + Aihi, I,+ + Aihy)= min{V(ui + Ah:, Vi+ Ah:))lA 2 0). Step 5. Set ui+l
=Ui+Aih’,
Vi+1
=
Vi+Aihye
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
132
Step 6. Compute
vf(“i+l~ vi+*)= (Vf’(ui+l,
Ui+l)Vf”(U’
t+l,
ui+l))t=
(dV/aU,
~V/dU)'u=~++~,,u=vi+~
Step 7. If (Vf(ui+i, ~i+~)(
gi+l
=
-Vf’(ui+*r
ui+l)
hi+, = g:+1+ y,hl
)
g’i:l = -Vfll(Ui+l, Ui+l),
>
II:+:, = g’i:, + yJ$
with y, = (gl+1-
,
81,gf+1> + w+1- g’i’,S::l> (gT,gi) + (g:‘,gI))
Set i = i + 1 and go to Step 4.
5. A numerical example
To illustrate the method described above we calculate the stress in a square, flat plate of side b and with a small circular hole of radius a at its centre [2]. It is subjected to a uniform tensile stress (T. The thickness of the plate may be assumed sufficiently small that the plane stress approximation is applicable. If b is very large we are able to obtain an analytical solution to this problem. The stress along the line Y = 0 (Fig. 2) is
X
1
2
3
Fig. 2. The mesh.
5
133
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
%Y
oTable 1 Results from Polak-Ribiere
algorithm
S
1.5
2.0
2.5
c?,& 9
2.716 40
2,772 52
2.657 74
I.0 -
Fig. 3. Element point stress; -
and nodal point stresses; analytical solution).
N = 9 c element
stress;
x
nodal
The maximum stress (X = a, Y = 0) is 3~. This solution is accurate near the hole when b % a. The mesh chosen is shown in Fig. 2. The two parameters are N, the number of points in the coordinate directions, and s the constant ratio of radial distances between successive rows of nodes. The numerical values taken for the example are a = 0.05, b = 1.0, E = 1.0, v = 0.3. Table 1 shows the results from the Polak-Rib&e algorithm for the case where N = 9. In Table 1 the ratios of maximum stress to applied stress, &,,,,/a, are the stresses computed for element I shown in Fig. 2. 4 is the number of iterations required for convergence of the conjugate gradient algorithm. Fig. 3 shows the variation of a,,, with distance from the hole along a small part of the X-axis very near to the hole for N = 9 and s = 2. The element stress is shown for elements such as II, III, VI,.VII, . . . of Fig. 2. The nodal point stresses are obtained by averaging; for example, cYY for nodes 2 and 3 are obtained as the averages of the values associated with elements II, III and III, IV, V, VI respectively. These results are close to those obtained in [2]. There the finite element method was used in conjunction with an over-relaxation method to solve the simultaneous, algebraic equations. The accuracy of that method therefore is much the same as the method described in this paper. To improve the accuracy the value of it must be increased. For example, when the values N = 19 and s = 1.3 are used the computed maximum stress for element I (Fig. 2) is 2.976. With the analytical solution of 3 the relative error is only 0.8%. In this case the number of iterations required for the conjugate gradient algorithm was q = 89. Element and nodal point stresses are plotted in Fig. 4. It is clear from this figure that good accuracy has been obtained in this case.
134
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
3’0 5Y -
o-
.I\.
I
ov 0’0s
O-06
0'07
0.08
0.09
O-10 X
Fig. 4. Element
and nodal point stress; N = 19 (. element stress; X nodal point stress; -
analytical solution).
6. Discussion and conclusion
The traditional finite element method involves the solution of a set of 2n simultaneous, linear equations (4), where n is the number of nodal points of the mesh. The elements of the overall stiffness matrix K, of order 2n x 2n, must be computed and stored. In general n is of the order of several hundreds. Fortunately, the matrix K is sparse. Using the rectangular form and Gauss elimination to solve (4) the storage requirement for storing matrix K is 2n x b, when b is the band width of the stiffness matrix and is of the order of 4d/n. On the other hand, by using the Gauss-Seidel or the sequential over-relaxation method to solve (4) only the non-zero elements of matrix K are stored. The number of non-zero coefficients in each row of K is dependent on p, the number of adjacent nodal points. In general [2] p = 9. In this case the storage requirement for storing K is 2n(2p) + n@ + 1) = (5~ + 1)n. The method described in this paper does not involve the solution of equation (4). Hence the elements of the overall stiffness matrix K does not have to be computed and stored. The variational problem (7) has been solved using the Polak-Rib&e method in which we need to store Ui, Ui, g:, gy, g:+l, gy+l, hi, hy (i = 1, 2,. . . , a). The total storage requirement is 8n. Table 2 compares the storage requirements for the various methods. Here u is the storage requirement to store the mesh data. To demonstrate the importance of storage saving we see in Table 2 that, when N = 19, the number of nodes n = 387 and the number of elements M = 698. The storage requirement a to store the mesh data and provide working space is approximately 11 K, i.e. Sn + 12M. On a time-sharing computer such as a PDP-10 the maximum storage available to each user is 40 K,
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
Table 2 Storage requirements
135
for each method Storage requirement
Method Gauss elimination (full matrix used)
4n2+ a
Gauss elimination (rectangular form used)
8n3” + a
Gauss-Seidel or SOR
46n+a
8n + a
Polak-Rib&e
say. This means roughly 20 K when storage for programming is excluded. Gauss-Seidel method, which requires the smallest amount of storage would be insufficient core store. Only a conjugate gradient method such method used in this paper would be suitable in the given circumstances. It should be remarked that if the overall stiffness matrix K of (4) rectangular form, then formulae (20) and (21) can be written as V--2 ‘S’KS-F’S
and
Thus, even with the for matrix K, there as the Polak-Rib&e can be stored in its
g=KG-F.
In this case the computing time can that for the Polak-Ribiere method. estimate the optimal over-relaxation of the matrix K. This is not easy to does not arise.
be significantly reduced and it becomes comparable with However, in the Gauss-Seidel method it is necessary to coefficient which means an estimate of the spectral radius do. By using a conjugate gradient method this difficulty
References [l] H. Ford, Advanced Mechanics of Materials (Longman, London, 1963). [Z] R.T. Fenner, Finite Element Method for Engineers (Macmillan, New York, 1975). [3] M.O. Bristeau, R. Glowinski, J. Periaux, P. Perrier, 0. Pironneau and G. Poirier, Application of optimal control and finite element methods to the calculation of transonic flows and incompressible viscous Rows, in: B. Hunt, ed., Numerical Methods in Applied Ffuids Dynamics (Academic Press, London, 1980) 203-312. [4] M.R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Nat. Bur. Standards 49 (1952) 409-436. [5] P. Concus, G.H. Golub and D.P. O’Leary, A generalized conjugate gradient method for the numerical solution of eihptic partial differential equations, Symp. Sparse Matrix Computations, Argonne National Lab., September 1975. [6] G. Gambolati, Diagonally dominant matrices for the finite element method in hydrology, Internat. J. Numer. Meths. Engrg. 6 (1973) 587-591.
136
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
[7] G. Gambolati, Use of the over-relaxation technique in the simulation of large groundwater basins by the finite element method, Internat. J. Numer. Meths. Engrg. 9 (1975) 219-233. [8] G. Gambolati, Free surface flow to a well by the finite element method and the over-relaxation technique, 2nd Internat. Symp. F.E.M. in Flow Problems, S. Margherita Ligure, Italy, 1976. [9] B.M. Irons, The Conjugate Gradient Method, Calgary Conf. Proc., 1973. [lo] D.S. Kershaw, The incomplete Cholesky-conjugate gradient method for the iterative solution of systems of linear equations, J. Comp. Phys. 26 (1978) 43-65. [ll] S.P. Newman and T.N. Narasimhan, Mixed explicit-implicit iterative finite element scheme for diffusion-type problems: I. Theory, Internat. J. Numer. Meths. Engrg. 11 (1977) 309-323. [12] J.K. Reid, On the method of conjugate gradients for the solution of large sparse systems of linear equations, in: Proc. Conf. Large Sparse Sets of Linear Equations (Academic Press, New York, 1971). [13] J.K. Reid, The use of conjugate gradients for systems of linear equations possessing ‘Property A’, SIAM J. Numer. Anal. 9 (1972) 325-332. [14] J.K. Reid, A discussion on a modified conjugate gradient method, Internat. J. Numer. Meths. Engrg. 8 (1973) 431432. [15] W.L. Wood, Note on a modified conjugate gradient method, Internat. J. Numer. Meths. Engrg. 7 (1973) 228-232. [16] G. Gambolati, Fast solution to finite element flow equations by Newton iteration and modified conjugate gradient method, Internat. J. Numer. Meths. Engrg. 15 (1980) 661-675. [17] R.M. Rosenberg, Analytical Dynamics of Discrete Systems (Plenum Press, New York, 1977). [18] E. Polak, Computational Methods in Optimization (Academic Press, New York, 1971).
1N APPLIED
A CONSULATE
MECHANICS
AND ENGINEERING
40 (1983) 127-135
GRADIENT ALGORITHM APPLIED STRESS-STRAIN PROBLEMS
Received 7 Aprit f981 Revised manuscript received 11 February
TO PLANE
1982
The numerical methods of optimal control theory have recently been applied effectively to problems arising from the application of finite element methods to transonic and incompressible viscous flows. This paper uses the same approach to show that a conjugate gradient algorithm can solve stress-strain problems with far less computer storage than that required by Conventional methods.
The main methods used for the solution of plane stress and plane strain problems have been those based upon finite differences and finite elements. When finite differences are used a biharmonic equation fl]
is solved and the stresses obtained from the equations
Equation (1) is transformed in matrix form as
to a set of simultaneous,
algebraic equations which may be written
AX=B.
When using a finite element method [2f again a set of simultaneous must be solved. In this case they are given by
(3)
atgebraic equations
where K is the stiffness matrix, S is the dis~~a~eme~t vector and F is the force vector acting on the object. Clearly, in both methods the main work is to solve the high-order systems (3) or (4) which in ~45-~8~5/83/$3.00
@ 1983, Elsevier Science Publishers B.V. (worth-Holland~
128
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
theory is no problem but in practice has the difficulty of computer storage for the matrix A or K. The order of these coefficient matrices is equal to twice the number of nodal points of the mesh being used, and for any realistic problem this will be several hundreds. A number of methods are available for solving large numbers of linear, algebraic equations, and the choice of method depends on the core storage available. If storage is not a problem, then the Gauss elimination or the Cholesky decomposition are better than most methods. If it is not possible to store the whole coefficient matrix, then methods such as the Gauss-Seidel, over-relaxation [2] or incomplete Cholesky-conjugate gradient [33 may be suitable although, even here, all non-zero elements of the coefficient matrix have to be stored. In many situations insufficient core store is available, even to make these methods viable, and so one is forced to look for methods which bypass the solution of large scale systems such as (3) or (4). This paper develops such a method using a conjugate gradient algorithm. The conjugate gradient method and modifications to it has been used extensively in many areas of engineering since Hestenes and Stiefel [4] first described their idea in 1952. The solution of the load-deflection equations arising from finite element discretization has been the subject of very many papers over the last decade. In particular, the use of conjugate gradient techniques in conjunction with finite element discretization has generated a rich literature [s-16]. There have been reports in the literature of convergence problems arising from the use of conjugate gradient techniques when applied to finite element discretization. The algorithm used in this paper is that due to Polak and Rib&e for which convergence has been established, and a rate of convergence can be assessed. The authors have not experienced any convergence difficulties in the present investigations. The Polak-Rib&e algorithm coincides with the well-known Fletcher-Reeves algorithm when applied to a quadratic function [18].
2. The primal problem For conservative, holonomic, mechanical systems under Hamilton’s principle 1171 may be expressed as
the action of potential
forces,
where T and V are the kinetic and potential energies, respectively. In problems of equilibrium equation (5) reduces to Vdt=O.
(6)
Because the system is assumed conservative the total energy T + V is constant and thus the potential energy is independent of time t. Equation (6) then reduces to
sv=o, which is the well-known basic equation for equilibrium.
Q. YE, D.J. Bell, A corzjugate gradient algorithm for plane stress-strain problems
129
In problems of plane strain and plane stress we may write the total potential energy of the system as [2] v=
v,+ vz
(8)
where V, is the stored strain energy and V2 the potential energy of the external forces. The problem is then a variational one: to minimize the functional V with respect to displacements. These displacements of the object in the directions x and y we write as u(x, y) and ~(x, y), (x, y) E R, where R is the domain occupied by the object. The strain e can then be evaluated from du
(
e = (e,,, eyv,e,,>” = -
dv
-,
dx’ay
au
-
ay
dv t . dx >
+-
In the linear, elastic problem the relationship
(9) between stress and strain is
a=De is the stress, and where u = (aXX,a;,, uXxv)f .
(11) Here E* = E and it* = Y for plane stress problems or E* = E/(1 - Y’), Y* = v/(1 - Y) for plane strain problems. E is Young’s modulus, v is Poisson’s ratio, and they are both independent of u and e (and therefore of u and u), determined only by the material properties. Thus, for linear, elastic problems, the stored strain energy in the entire solution domain fz is
(12) The external forces involve the body forces x, t and the forces F., F, in the directions of x and y at the boundary C. The potential energy of the external forces is
v*= - j i, (it%+ %)dx
dy - j. (F’u + F,v)ds .
(13)
Thus, the total potential energy of the system is v=$
II
R
e’De dx dy - 1 I, (% + %)dx dy - 1 (F,u + F,u)ds .
(14)
Here u and v are functionals of x and y so the variational problem inherent in (7) must be solved in the Sobolev space H’(LJ). This means the feasible domain of control variables u and ~1is a Sobolev Space. This is a property which must be considered if we are to use a conjugate gradient method to solve the given extremum problem [3].
130
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
3. The numerical problem
For simplicity we here neglect the body forces. Then v=$
(F,u + F,u)ds
e’De dx dy -
(1%
.
In order to calculate the double integral in equation (15) we must evaluate the strain e through equation (9). This implies the calculation of &/&, dully c~v/& and dully. Many methods are available for this purpose; for example, by the finite element method we take triangular elements and assume the displacements u and v vary linearly over each element. Fig. 1 illustrates the mth eiement of a mesh. We assume the displacements at points (x, y) within the element are given by u(x, y) = G -+ czx -+ GY 3
zJ(X,y) = c4 + csx + csy *
The strain e, is then obtained in terms of the displacements
em= &-
at the vertices i, j, k as
B,&,,
(17)
m
where A, is the area of the element, 8, =
(Ui, Vi,
uj,
vj, Uk,
WCJ
VkY
and B, is the matrix b, 0 ai
0 a,
bi
bj 0 aj
0 aj
bk 0
bj
ak
0
(19)
Y
Ll+----F Fig. 1. Typical element of mesh.
X
131
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
Note that B, is only determined by the mesh and is independent the mesh is given, B, is a constant matrix. The total potential energy (15) now becomes
of the displacements.
When
where s is the number of elements of the mesh and 1 the number of nodes on the boundary. Given the displacements at the nodes of the mesh the total potential energy V can be calculated by equation (20). The variational problem given by (7) therefore has 2n independent variables Ui, L+(i = 1,2, . . . , n), where it is the number of nodes in the mesh. It will be solved in Euclidean space R2”. 4. A conjugate gradient algorithm In order to solve the problem of equation (7) by a conjugate gradient method the gradients of V with respect to the variables Ui, vi must be evaluated. From (20) we obtain %)-Fix, I
i=l,2
,....)
n
(21)
where Fix, Fiy are zero if i is not on the boundary. From the constraints at some nodes the displacements Ui and/or Ui are given. The corresponding derivatives L~v/b’u,and 3V/& therefore vanish. The algorithm used in this paper is the one due to Polak and Ribike [18]. The procedure for this method is as follows. Step 1. Select a normal displacement
u. and v. subject to the boundary condition. Calculate
V(uovo) by (20) and ?f(~ovo>
= K.~(u~v~)V$(U~V~>>
=
@v/au,av/a2))tu=w,o,w
by (21), where V, is the operator (a/aul, . . . , a/au,)’ and similarly for V,. Set i = 0, E > 0. Step 2. If (Vf(uovo)~< E, then stop; else go to Step 3. Step 3. Set h’d = g’6 =
hb = g; = -Vf(uov0) ) Step 4. Compute
-Vf”(uovo)
.
Ai > 0 such that
V(Ui + Aihi, I,+ + Aihy)= min{V(ui + Ah:, Vi+ Ah:))lA 2 0). Step 5. Set ui+l
=Ui+Aih’,
Vi+1
=
Vi+Aihye
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
132
Step 6. Compute
vf(“i+l~ vi+*)= (Vf’(ui+l,
Ui+l)Vf”(U’
t+l,
ui+l))t=
(dV/aU,
~V/dU)'u=~++~,,u=vi+~
Step 7. If (Vf(ui+i, ~i+~)(
gi+l
=
-Vf’(ui+*r
ui+l)
hi+, = g:+1+ y,hl
)
g’i:l = -Vfll(Ui+l, Ui+l),
>
II:+:, = g’i:, + yJ$
with y, = (gl+1-
,
81,gf+1> + w+1- g’i’,S::l> (gT,gi) + (g:‘,gI))
Set i = i + 1 and go to Step 4.
5. A numerical example
To illustrate the method described above we calculate the stress in a square, flat plate of side b and with a small circular hole of radius a at its centre [2]. It is subjected to a uniform tensile stress (T. The thickness of the plate may be assumed sufficiently small that the plane stress approximation is applicable. If b is very large we are able to obtain an analytical solution to this problem. The stress along the line Y = 0 (Fig. 2) is
X
1
2
3
Fig. 2. The mesh.
5
133
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
%Y
oTable 1 Results from Polak-Ribiere
algorithm
S
1.5
2.0
2.5
c?,& 9
2.716 40
2,772 52
2.657 74
I.0 -
Fig. 3. Element point stress; -
and nodal point stresses; analytical solution).
N = 9 c element
stress;
x
nodal
The maximum stress (X = a, Y = 0) is 3~. This solution is accurate near the hole when b % a. The mesh chosen is shown in Fig. 2. The two parameters are N, the number of points in the coordinate directions, and s the constant ratio of radial distances between successive rows of nodes. The numerical values taken for the example are a = 0.05, b = 1.0, E = 1.0, v = 0.3. Table 1 shows the results from the Polak-Rib&e algorithm for the case where N = 9. In Table 1 the ratios of maximum stress to applied stress, &,,,,/a, are the stresses computed for element I shown in Fig. 2. 4 is the number of iterations required for convergence of the conjugate gradient algorithm. Fig. 3 shows the variation of a,,, with distance from the hole along a small part of the X-axis very near to the hole for N = 9 and s = 2. The element stress is shown for elements such as II, III, VI,.VII, . . . of Fig. 2. The nodal point stresses are obtained by averaging; for example, cYY for nodes 2 and 3 are obtained as the averages of the values associated with elements II, III and III, IV, V, VI respectively. These results are close to those obtained in [2]. There the finite element method was used in conjunction with an over-relaxation method to solve the simultaneous, algebraic equations. The accuracy of that method therefore is much the same as the method described in this paper. To improve the accuracy the value of it must be increased. For example, when the values N = 19 and s = 1.3 are used the computed maximum stress for element I (Fig. 2) is 2.976. With the analytical solution of 3 the relative error is only 0.8%. In this case the number of iterations required for the conjugate gradient algorithm was q = 89. Element and nodal point stresses are plotted in Fig. 4. It is clear from this figure that good accuracy has been obtained in this case.
134
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
3’0 5Y -
o-
.I\.
I
ov 0’0s
O-06
0'07
0.08
0.09
O-10 X
Fig. 4. Element
and nodal point stress; N = 19 (. element stress; X nodal point stress; -
analytical solution).
6. Discussion and conclusion
The traditional finite element method involves the solution of a set of 2n simultaneous, linear equations (4), where n is the number of nodal points of the mesh. The elements of the overall stiffness matrix K, of order 2n x 2n, must be computed and stored. In general n is of the order of several hundreds. Fortunately, the matrix K is sparse. Using the rectangular form and Gauss elimination to solve (4) the storage requirement for storing matrix K is 2n x b, when b is the band width of the stiffness matrix and is of the order of 4d/n. On the other hand, by using the Gauss-Seidel or the sequential over-relaxation method to solve (4) only the non-zero elements of matrix K are stored. The number of non-zero coefficients in each row of K is dependent on p, the number of adjacent nodal points. In general [2] p = 9. In this case the storage requirement for storing K is 2n(2p) + n@ + 1) = (5~ + 1)n. The method described in this paper does not involve the solution of equation (4). Hence the elements of the overall stiffness matrix K does not have to be computed and stored. The variational problem (7) has been solved using the Polak-Rib&e method in which we need to store Ui, Ui, g:, gy, g:+l, gy+l, hi, hy (i = 1, 2,. . . , a). The total storage requirement is 8n. Table 2 compares the storage requirements for the various methods. Here u is the storage requirement to store the mesh data. To demonstrate the importance of storage saving we see in Table 2 that, when N = 19, the number of nodes n = 387 and the number of elements M = 698. The storage requirement a to store the mesh data and provide working space is approximately 11 K, i.e. Sn + 12M. On a time-sharing computer such as a PDP-10 the maximum storage available to each user is 40 K,
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
Table 2 Storage requirements
135
for each method Storage requirement
Method Gauss elimination (full matrix used)
4n2+ a
Gauss elimination (rectangular form used)
8n3” + a
Gauss-Seidel or SOR
46n+a
8n + a
Polak-Rib&e
say. This means roughly 20 K when storage for programming is excluded. Gauss-Seidel method, which requires the smallest amount of storage would be insufficient core store. Only a conjugate gradient method such method used in this paper would be suitable in the given circumstances. It should be remarked that if the overall stiffness matrix K of (4) rectangular form, then formulae (20) and (21) can be written as V--2 ‘S’KS-F’S
and
Thus, even with the for matrix K, there as the Polak-Rib&e can be stored in its
g=KG-F.
In this case the computing time can that for the Polak-Ribiere method. estimate the optimal over-relaxation of the matrix K. This is not easy to does not arise.
be significantly reduced and it becomes comparable with However, in the Gauss-Seidel method it is necessary to coefficient which means an estimate of the spectral radius do. By using a conjugate gradient method this difficulty
References [l] H. Ford, Advanced Mechanics of Materials (Longman, London, 1963). [Z] R.T. Fenner, Finite Element Method for Engineers (Macmillan, New York, 1975). [3] M.O. Bristeau, R. Glowinski, J. Periaux, P. Perrier, 0. Pironneau and G. Poirier, Application of optimal control and finite element methods to the calculation of transonic flows and incompressible viscous Rows, in: B. Hunt, ed., Numerical Methods in Applied Ffuids Dynamics (Academic Press, London, 1980) 203-312. [4] M.R. Hestenes and E. Stiefel, Methods of conjugate gradients for solving linear systems, J. Res. Nat. Bur. Standards 49 (1952) 409-436. [5] P. Concus, G.H. Golub and D.P. O’Leary, A generalized conjugate gradient method for the numerical solution of eihptic partial differential equations, Symp. Sparse Matrix Computations, Argonne National Lab., September 1975. [6] G. Gambolati, Diagonally dominant matrices for the finite element method in hydrology, Internat. J. Numer. Meths. Engrg. 6 (1973) 587-591.
136
Q. Ye, D.J. Bell, A conjugate gradient algorithm for plane stress-strain problems
[7] G. Gambolati, Use of the over-relaxation technique in the simulation of large groundwater basins by the finite element method, Internat. J. Numer. Meths. Engrg. 9 (1975) 219-233. [8] G. Gambolati, Free surface flow to a well by the finite element method and the over-relaxation technique, 2nd Internat. Symp. F.E.M. in Flow Problems, S. Margherita Ligure, Italy, 1976. [9] B.M. Irons, The Conjugate Gradient Method, Calgary Conf. Proc., 1973. [lo] D.S. Kershaw, The incomplete Cholesky-conjugate gradient method for the iterative solution of systems of linear equations, J. Comp. Phys. 26 (1978) 43-65. [ll] S.P. Newman and T.N. Narasimhan, Mixed explicit-implicit iterative finite element scheme for diffusion-type problems: I. Theory, Internat. J. Numer. Meths. Engrg. 11 (1977) 309-323. [12] J.K. Reid, On the method of conjugate gradients for the solution of large sparse systems of linear equations, in: Proc. Conf. Large Sparse Sets of Linear Equations (Academic Press, New York, 1971). [13] J.K. Reid, The use of conjugate gradients for systems of linear equations possessing ‘Property A’, SIAM J. Numer. Anal. 9 (1972) 325-332. [14] J.K. Reid, A discussion on a modified conjugate gradient method, Internat. J. Numer. Meths. Engrg. 8 (1973) 431432. [15] W.L. Wood, Note on a modified conjugate gradient method, Internat. J. Numer. Meths. Engrg. 7 (1973) 228-232. [16] G. Gambolati, Fast solution to finite element flow equations by Newton iteration and modified conjugate gradient method, Internat. J. Numer. Meths. Engrg. 15 (1980) 661-675. [17] R.M. Rosenberg, Analytical Dynamics of Discrete Systems (Plenum Press, New York, 1977). [18] E. Polak, Computational Methods in Optimization (Academic Press, New York, 1971).