Gauss-Seidel Iteration Applied to Nonlocal Averaging

For the purpose of finite element analysis, we will now assume that subscripts p and v label the numerical integration points of finite elements, rather than the individual microcracks. This means that the micro – cracks are represented by their mean statistical characteristics sampled only at the numerical integration points.

In finite element programs, nonlinearity is typically handled by iterations of the loading steps. Let us, therefore, examine the iterative solution of (13.3.8) or (13.3.13), which represents a system of N linear algebraic equations for N unknowns Дif AS^P are given. The matrix of is, in general, nonsymmetric (because the influence of a large crack on a small crack is not the same as the influence of a small crack on a large crack). This nonsymmetry seems disturbing until one realizes that this is so only because of our choice of variables AS»^ and Дsj^ which do not represent thermodynamically

conjugate pairs of generalized forces and generalized displacements. If ASjP were expressed in terms of the crack opening volumes, then the matrix of the equation system resulting from (13.3.8) or (13.3.13) would have to be symmetric (because of Betti’s theorem) and also positive definite (if the body is stable). These are the attributes mathematical ly required for convergence of the iterative solution by Gauss-Seidel method (e. g., Reklorys 1969; Collate 1960; Korn and Korn 1968; Varga 1962; Fox 1965; Strang 1980). Aside from that, convergence of the iterative solution of (13.3.8) or (13.3.13) must also be expected on physical grounds (because it is mechanically equivalent to the relaxation method, which always converges for stable clastic systems).

In the r-th iteration, the new, improved values of the unknowns, labeled by superscripts [r – f – 1], are calculated from the previous values, labeled by superscript [r], either according to the recursive relations;


Д;3;І+1) – Д^,І + Л^Др!;1 (13.3.27)

12= I

_______ N

д^ОМг-м] = д5.0) +^л(і„Д5(|)М (р = 1,…,Л’) (13.3.28)


or according to the recursive relations:

fi-f N

Др|[+,1 = Д^ + ^Л;шДр!)-+,1+ X. Л^ДрМ (p – І,…,Л0 (13.3.29)

и -~1 12 = //+ 1

Д5(1)Н-1] = Д50) + ^Л^Д^+Ч + ]Г Л^Д^М (р I (13.3.30)

12=1 12=// + ]

Equation (13.3.28), also known as the Gauss method or Jacobi method, is less, but normally only slightly less efficient than (13.3.30), in which the latest approximations are always used. The values of AS^ may be used as the initial values of Дб’/j1^ in the first iteration.

It is possible to derive Eqs. (13.3.27) and (13.3.28) more directly, rather than from (13.3.5). To this end, we note that the sequence of iterations is identical to a solution by the relaxation method in which one crack after another is relaxed (i. e., its pressure reduced to zero) while all the other cracks are frozen (so in each relaxation step, one has a problem with one crack only), as illustrated in Fig. 13.3.1b. Each relaxation produces pressure on the previously relaxed cracks. After relaxing all the cracks one by one, the cycle through all the cracks is repeated again and again. This kind of relaxation is known in mechanics to converge in general (this was numerically demonstrated for a system of cracks and inclusions by Pijaudier-Cabot and Bazant 1991). The solution to which the relaxation process converges is obviously that defined by Eq. (13.3.8). (Note also that this relaxation argument in fact represents a simple way to prove the superposition equation (13.3.5).)

For structural engineers, it is interesting to note the similarity with the Cross method (moment distri­bution method) for elastic frames. Relaxing the pressure at one crack while all the other microcracks are frozen (glued) is analogous to relaxing one joint in a frame while all the other joints are held fixed.

Repeating this for each joint, ami then repeating the cycles of such relaxations of all the joints, eventually converges to the exact solution of the frame.

The macro-continuum counterpart of the Gauss-Seidel iterative method, which converges to the solution of the Fredholm integral equation (13.3.9), is analogous to (13.3.28) and is given by the following relation for successive approximations (iterations):

Д5(|)1г+|1(х) = ДУОІ’Й f f A(x.£)Ab’(l)W(0<W(£) (13.3.31)


The discrete approximation of the last relation is the equation that ought to be used in finite element programs with iterations in each step. Wc sec that the form of averaging is different from that assumed in the phenomenological models we described. There are now two additive spatial integrals: one for close-range averaging of (he inelastic stresses from the local constitutive relation, and one for long-range crack interactions based on the latest iterates of the inelastic stresses.

In programming, the old iterates need not be stored in the computer memory. So the superscripts [r] and [r 4-1] may be dropped and equations (13.3.9) and (13.3.31) may be replaced by (he following assignment statements:

_____ x

AS’;p <- AS^ ]ГЛ„„Д§l]) (ц -1,2, …N) (13.3.32)

V – I

Д5’(|,(х) г-Д^х) i – [ A(x,£)AS(l)(£)dV(£) (13.3.33)


A strict implementation of Gauss-Seidel iterations suggests programming each iteration loop for F. q. (13.3.32) to be contained within another loop for the iterations of the loading step in which the displacement and strain increments in the structure are solved. However, it is computationally more efficient to use one common iteration loop serving both purposes. Then, of course, the iteration solution is not exactly the Gauss-Seidel method because the strains are also being updated during each iteration. There is already some computational experience (Ozbolt and Bazant 1996) showing that convergence is still achieved.

The common iteration loop has the advantage that it permits using the explicit load-step algorithm for structural analysis. In each loading step of this algorithm, one evaluates in each iteration at each integration point the elastic stress increments EAe and the local inelastic stress increments AS from fixed strains Де; then one uses (13.3.32) to calculate from AS the nonlocal inelastic stress increments AS for all the integration points, solves new nodal displacements by clastic structural analysis, and, finally, updates the strains.