TECHNISCHE MECHANIK 6(1985)Heft2

Manuskripteingang: 22. 5.1984

Anwendung der Konjugierte-Gradienten-Methode bei der Lösung

elastoplastisdier Aufgaben mit Hilfe der Methode der finiten Elemente

Dobril Christow, Marco Todorow

l. Einleitung

Untersuchungen von plastischen Spannungs- und Verzer-

rungszuständen mit Hilfe der Methode der finiten Ele-

mente benutzen hauptsächlich die Methode der Anfangs-

(Zusatz-)lasten [l], Die Methode der veränderlichen

Steifigkeit [l], [8] wird selten verwendet, weil die Re-

chenzeiten zur Lösung der Gleichungssysteme durch di-

rekte Methoden unvertretbar lang sind. Es ist aber be-

kannt [2], [9], daß sie physikalisch genauer die plasti-

sche Verformung beschreibt und schneller, was die An-

zahl der lterationsschritte betrifft, als die anderen Ver-

fahren konvergiert.

ln der vorliegenden Arbeit wird die Methode der verän-

derlichen Steifigkeit mit der Konjugierte-Gradienten-

Methode (KGM) [3], [4], [6], [10] verbunden, die in

diesem Falle bestimmte Vorteile vor den direkten Ver-

fahren zur Lösung von Gleichungssystemen besitzt.

2. Die Konjugierte-Gradienten-Methode

Der Grundgedanke der Methode der finiten Elemente,

angewendet auf das Funktional der gesamten potentiel-

len Energie, führt bei der linearen Elastizität zur Glei-

chung

H(u)=äuTKu—uTb. (1)

Hier sind:

H(u) — die gesamte potentielle Energie des elastischen

Systems,

u — die unbekannten verallgemeinerten Knotenver-

schiebungen, -

b — die vorgegebenen verallgemeinerten Knoten-

kräfte,

K ~ die Gesamtsteifigkeitsmatrix.

Die Matrix K ist symmetrisch und bei korrekt eingebau-

ten Randbedingungen positiv definit.

Die Forderung nach Extremum von (1) fiihrt zum linea-

ren Gleichungssystem

Ku=b (m

Die Lösung“ von (2) mit Hilfe der KGM erfolgt nach dem

Algorithmus [4]:

‘V\

v1=q=b_Kg, q

m=Kn‚ W

_V. 'l'i

_ T. ‚

Vi Pi;

“i =ui_1+aivi, 5" fiiri=1,2,...,n.

ri+1 =fi—“iPi»

531"“ßi =—T_*~a

v'r-

i

"i+1 = l‘i+1 + 5m

Hier sind

uo — eine beliebige Anfangsnäherung,

ri — die Abweichung,

vi —— der Richtungsvektor,

n — die Anzahl der Gleichungen in

Theoretisch kommt man zur exakten Lösung von (2)

nach n Iterationsschritten. Eine Reihe von Forschungen

[6], [10] zeigen aber, daß man eine fiir die Praxis ausrei-

chende Genauigkeit nach bedeutend weniger Iteratio-

nen m erhalten kann. Mit

m27°n‚ (4)

für die in [6] untersuchten Gleichungssysteme (n =

363 . . . .2421) gilt '7 = 0,1 v. . . 0,23. Die gleichen Verfas-

ser bemerken, daß mit wachsendem n, 7 kleiner wird.

Zur Beendigung der Prozedur (3) kann man verschiedene

Kriterien verwenden. In [10] wird die Befriedigung der

Ungleichung

f’ir'l'i

er: lBuße] <5)

vorgeschlagen. Dabei ist Ei der vorgeschriebene Fehler

für die mittlere quadratische Abweichung vom Gleich-

gewichtszustand der Knoten. Als ein anderes Kriterium

[3] wird die Differenz zwischen zwei aufeinanderfolgen-

den Näherungen benutzt:

Ilui—ui 1" "vi"

: ___;__ 2 . __ < ‚ 6

ex ||ui_1|| “I ||ui_1|| \62 Ü

wobei H. . . ll als Euklid, bzw. als Maximumnorm aufzu-

fassen ist.

Man kann zeigen [4], daß für jedes i der Prozedur (3) die

Ungleichung

"MH)<HM)

erfüllt ist. Dieser Umstand gewährleistet die Konvergenz

der KGM. Man kann weiterhin zeigen, daß gilt

51

Bild l

Vergleich der Konvergenzkriterien

Ä

5% er

5x ex

1.0 0,1 <

QJ - am

0 0‚

1 1o 20 3'0 4'0 5'0 - I Iterationen

fl (“i+1) _ [I : _ ä ai [fir ‚i ‚ (7) cherplätze. Gleichzeitig erhöht sich aber die Anzahl der

Bezeichnet man mit H* eine Näherungsschätzung für die

Arbeit der äußeren Kräfte, so kann man als Kriterium für

die Exaktheit der Lösung auch die Ungleichung

(Xi r? l'i

_ 2 II*

verwenden.

Auf Bild 1 sind die relativen Fehler entsprechend (5),

(6) und (8) über der Anzahl der Iterationen aufgetragen.

Die Kurven beziehen sich auf die Lösung des Glei-

chungssystems für das Beispiel von Bild 6. Man sieht,

daß zwischen den einzelnen Kurven kein großer qualita-

tiver Unterschied existiert Bei der weiteren Arbeit

wurde deshalb das Kriterium (5) angenommen, welches

g E3 i

einen klaren physikalischen Sinn hat.

Die Vorteile der KGM vor der Lösung der linearen FEM-

Gleichungssysteme mit direkten Verfahren sind fol-

gende:

a) Die Steifigkeitsmatrix kann in voll verdichteter Form

dargestellt werden. Dabei werden alle Nullelemente aus-

geschlossen. Das führt zu Verminderung der Anzahl der

Fließkommaoperationen und der erforderlichen Spei-

52

Operationen mit ganzen Zahlen, die aber bedeutend

schneller als die Fließkommaoperationen durchgeführt

werden. Bezeichnet man mit Ä die mittlere Zahl der

Knoten, die einen Beitrag zum Gleichgewicht eines be-

stimmten Knotens haben, und mit K die Anzahl der

Unbekannten pro Knoten, so gilt für die Anzahl N der

von Null verschiedenen Elemente der Matrix K

nv K ..

2 (M + x + 1), <9)

worin nv die Knotenanzahl ist.

Der Koeffizient Ä hängt von den topologischen Eigen-

schaften des Netzes und von der Ordnung der Form-

funktion ab. Für ebene Aufgaben und finite Elemente

mit linearen Ansätzen für die Verschiebungen ist Ä z 7;

fiir quadratische Ansätze ist Ä ä 15. Für räumliche Pro-

bleme und lineare Verschiebungsfunküonen ist Ä z 27.

Die verdichtete Darstellung der Matrix K erfordert zu-

sätzliche Speicherplätze für die auf die Lage der von

Null verschiedenen Elemente hinweisende Informa-

tion. Es kann gezeigt werden, da5 dafür

N:

N1 z 25 (x + 4) (10)

ganze Zahlen notwendig sind. Auf Bild 2 ist die Größe

3h

0,3-

a2.

a1

500

Bild 2

Die erforderlichen Hauptspeicherplätze in Kbyte für die Unter-

bringung des Gleichungssysteme bei:

——— verdichteter Speicherung

—‚—.— sky-line—Speicherung

— über n = nv ' K aufgetragen, wobei C die notwendigenn

Speicherplätze in Kbyte sind.

In den für die verdichtete Speicherung der Matrix K not—

wendigen Speicherplätzen sind auch Plätze für die Ar-

beitsfelder zur Lösung des Gleichungssystems und die

Adressenmatrizen vorgesehen. Auf demselben Bild sind

ebenfalls die notwendigen Speicherplätze aufgetragen,

für den Fall, daß die Matrix K als Band-sky-line-Matrix

gespeichert wird. Der Koeffizient

ß = A + 1 ‚nv

wobei A die größte Knotennummerdifferenz ist, hängt

von der Art der Nummerierung ab. Falls sie günstig ist,

nimmt ß meistens Werte im Intervall 0,07 bis 0,3 an.

Aus Bild 2 ist ersichtlich, daß bei vorhandenem Spei-

cherraum von 64 Kb und verdichteter Matrix, ohne ex-

terne Speicher bis 1000 Gleichungen gelöst werden kön-

nen. Bei Bandmatrix (ß = 0,1) kann man unter gleichen

Bedingungen nicht mehr als 400 Gleichungen lösen.

b) Die Matrix K ändert sich nicht während der Proze-

dur (3) der KGM. Wird es notwendig, die Werkstoff-

konstanten bzw. die geometrischen Kennwerte eines

finiten Elementes zu verändern, trifft dies auf ganzde-

finite Weise leicht bestimmbare Teile der Gesamtsteifig-

keitsmatrix K. Falls man dagegen eine direkte Lösungs-

methode verwendet, inuß man die ganze Matrix K neu

bilden, bzw. von einem externem Speicher neu einlesen.

(11)

c) Die Lösung des Gleichungssystem kann mit einer be-

liebig vorgeschriebenen Genauigkeit durchgeführt wer-

den. Das ist besonders wichtig bei den nichtlinearen Pro-

blemen, da sie zur Lösung einer Reihenfolge linearer

Aufgaben führen. Für die Zwischenschritte ist keine

hohe Genauigkeit notwendig, was insgesamt die Anzahl

der Rechenoperationen verringert.

d) Es verringert sich ebenfalls die Anzahl der Fließkom-

maoperationen. Bei dem Lösen des Gleichungssystems

nach Cholesky ist die Anzahl der Fließkommamultipli-

kationen M1 annähernd

~1

Mpaßgns’ <12)

‘r

und bei verdichteter Matrix und der KGM —

(13)

Auf Bild 3 ist die Größe M/n2 als Funktion der Anzahl

M2 z-ä—ynzax + 1: +11).

der Gleichungen n aufgetragen.

gn2

Ü

,4 af=a_1=27 „120,3 „ y

12

”25. 1=7 r=aa .~ .6

6

4 _ acsg 2-7 [-0.3

___„ /OC%6 1-7 no.7 _ _ 7

2 „aexz 1-7_ [an g _ß'o'

m» ‚000 n

Bild 3

Die erforderliche Anzahl von Fliefikommamultiplikationen bei

der Lösung eines Gleichungssystem nach:

———— dem Cholesky-Verfahren

—.—.— der KGM mit verdichteter Matrix

e) Beim Arbeiten mit verdichteter Matrix K ist die Rei-

henfolge der Knoten- und Elementennummerierung für

die Rechenzeit bedeutimgslos. Dieser Umstand gestat-

tet, bei Notwendigkeit, leichtes Hinzufügen von neuen

finiten Elementen zum ursprünglich untersuchten Ob-

jekt.

Der wesentliche Nachteil der KGM, der sie konkurrenz-

unfähig den direkten Lösungsmethoden gegenüber

macht, ist die Tatsache, daß für jede neue rechte Seite

die Prozedur (3) voll durchzuführen ist. Bei nur einem

Belastungsfall, kann die KGM wie bereits oben gezeigt,

effektiver als die direkten Verfahren sein.

3. Methode der veränderlichen elastischen Para-

meter in Zusammenhang mit der Deforma-

tionstheorie der Plastizitätslehre

Bei der Methode der veränderlichen elastischen Para-

meter [2] (Methode der Sekantensteifigkeit) werden die

physikalischen Gleichungen der Deformationstheorie der

Plastizität auf die Form des ‚Hookeschen Gesetzes zu-

rückgeführt:

1

ex _ [ax—m(oy+az)1,...mzx = a; rzx- Hier sind

E-x- g ai/el

l +1—2“ g ’

3E Ei

> (14)

1 1—2p ‚0|

1 + 1—2“ ‚_l.

3E « Ej

G* - 0i/3ei ‚

53

. .. . . fl 3 2 2 2die veranderlichen elastischen Parameter, ei = —3~— (ex —ey)2 + (<5y —Ez)2 + (#52 ——Ex)2 + E ('Yxy +7yz +7“) (15)

ist die äquivalente Dehnung und 0i ist die entsprechende

Spannung nach der Werkstoffverformungskurve. Die An-

wendung der Methode der veränderlichen elastischen Pa-

rameter ist aus Bild 4 ersichtlich.

4. Verbindung der Konjugierte-Gradienten—Me-

thode der veränderlichen elastischen Parame-

ter

Der Grundgedanke der vorliegenden Arbeit besteht dar-

in, die Iterationsprozesse beider Methoden miteinander

zu verbinden. Zu diesem Zweck wird bei der ursprüng-

lichen Gesamtsteifigkeitsmatrix K und bei voller Bela-

stung b eine Anzahl KGM-lterationen durchgeführt,

ohne hohe Genauigkeit zu fordern (man kann z. B. die

Bedingung (5) bei e] = 0,1 erfüllen). Es ist anzunehmen,

daß diese Anzahl 0,5 bis 0,75 von der gesamten für die

Lösung mit hoher Genauigkeit (61 = 0,01) erforderlichen

Iterationenanzahl sein kann. Es sei danach die äquiva-

lente Spannung

M!

2

*-

(16)

in einigen finiten Elementen höher geworden als die, der

erreichten Dehnung Gil auf der Verformungskurve ent-

sprechende Spannung oil (Bild 4). Die Werkstoffkon-

stanten dieser Elemente werden nach (l4) neuberechnet

und ihre neuen Elementsteifigkeitsmatrizen in die Ge-

samtsteifigkeitsmatrix eingearbeitet. Die Iterationspro-

zedur (3) der KGM wird dann bei der so modifizierten

Matrix fortgesetzt. Als Anfangsnäherung uo wird dabei

die bereits gewonnene Lösung angenommen.

Ein Blockdiagramm des so beschriebenen Algorithmus

ist auf Bild 5 dargestellt.

5. Das Computerprogramm

Nach dem beschriebenen Algorithmus wurde ein Pro-

gramm in FORTRAN für den Minicomputer ISOT 310

zur Lösung von elastoplastischen ebenen Aufgaben‘ge-

schrieben. Die Unterprogramme zur Bildung der Elemen-

tensteifigkeiten sind so aufgebaut, daß sie direkt die Dif-

ferenz zwischen den „neuen” und den „alten” Elemen-

tensteifigkeitsmatrizen berechnen. Die Einarbeitung die-

ser Differenz in die Gesamtsteifigkeitsmatrix führt auto-

matisch zu ihrer Erneuerung. Die kinematischen Rand-

bedingungen werden allgemein in der Form ui = cuj

eingebaut, wobei i und beliebig sind. Ist c = 0, so wird

der entsprechende Freiheitsgrad gleich Null gesetzt.

Diese Randbedingungen sind so in die KGM-Prozedur

eingebaut, dal's sich die Gesamtsteifigkeitsmatrix nicht

ändert. Die Berücksichtigung von Federkonstanten ist

leicht möglich.

Das Programm arbeitet in Dialogbetrieb. Die Anzahl der

durchzuführenden Iterationen wird über die Angabe des

zulässigen Fehlers (5) gesteuert. Man kann die Zwischen—

ergebnisse analysieren und entscheiden, 0b die Lösung

abzubrechen ist.

54

1 2

oil z 7—2— .\/ (ax — (5)2 + (av—"02 + ("Z ‘03)2 + 6(Tx2y +132 +sz)

Berechnung vonu nach KEN

SchLeLfeüber Elem.

Berechnung van 5,— and ÖI”

für das Element

Berechnung von

Effür das Element

l

Berechnung van Ele-

mentsteifigkeitsmatrtx

und Einaröeifung in K

h—w Änderungwnß

Bild 5

Blockdiagramm des Algorithmus

Das Unterprogramm, das die Bedingung

0: (fi) = 0i(€i)

prüft und bei Notwendigkeit die elastischen Parameter

verändert, ist so aufgebaut, daß es sowohl den ideal-

elastoplastischen Werkstoff als auch solchen mit Fließ-

ebene und Verfesügung modellieren kann.

0:670]; Mp07

6. Beispiele

6.1. Kreisscheibe mit Öffnung unter Innendruck

Die analytische Lösung für den elastoplastischen Zustand

einer gelochten Kreisscheibe mit konstanter Dicke unter

Innendruck ist in [2, § 35] angeführt. Bei einem Innen-

durchmesser von 20 mm, Außendurchmesser 40 mm,

/A\

200

AG

/

‚q f7

no

{q I! g 11

Ws/-Ioo f

\4\\

Bild 6

Kreisscheibe mit Öffnung ‘200

unter Innendka

55

Dicke 1 mm und Fliefigrenze (IF = 300 MPa ist der nach

der ‘Gestaltänderungsenergiehypothese berechnete Druck

ps, bei dem das erste Fließen an der Öffnung ansetzt,

gleich 128,6 MPa. Diese Scheibe ist in zehn ringförmige

finite Elemente [5] mit einer Differenz zwischen dem

äußeren und dem inneren Durchmesser jedes Elementes

von 1 mm aufgeteilt. Die Anzahl der Gleichungen ist ll.

Es wurden zwei Lösungen durchgeführt — bei 1,4- ps bzw.

1,7 ps. Die berechneten Radial- und Tangentialspan-

nungen sind im Bild 6 gezeigt, worin die ausgezogenen

Linien die nach l2] analytisch bestimmten Werte dar-

stellen. Bei p = 1,4 ps ist Konvergenz nach 101 Itera—

tionen erreicht, bei p = 1,7 ps — nach 187 Iterationen.

6.2. Celochte Scheibe unter Zugbelastung

Nach der Methode der Anfangsspannungen ist dieses

Problem in [1, S. 363] und nach der Deformations-

theorie in [9] untersucht.

Bei der Vernetzung nach Bild 7 gewinnt man eine ausrei-

chend genaue elastische Lösung (6, < 0,01) nach 63 Ite-

rationen. Es wurden 4 Rechnungen bei nichtelastischen

Werkstoffverhalten durchgeführt. Bei p = 81 N fließt nur

ein einziges Element (Bild 7, 3). Es waren 70 Iterationen

notWendig, davon 3 mit Änderung der elastischen Para-

meter.

Bei p = 143 N und 126 Iterationen, davon l4 mit Plasti-

Zone nac 426 Itera-

Bild?

p 2,, a, 12,, 1p Si°2‘f°i‘5‘i%if:°ä’:'3,‘;‚Ö„‘;“E:;äi‘m„1.....

a) I 5)P-148M Ian!

elast etucher

Werke .Pautische

tic/my.

s l1% l //

"V m w '— ""

V w“ I

c) p-178/V. d) p-fla/V.

Ideal elastoplasti- Werkstoff mit

scher Werkstoff. . linearer Verfestigun .

Nautische Zone nach Plastische lone MC

217 Iterationen. 156 Iterationen.

Keine Konvergenz.

!20mm

WG

Q5

W 8.0 3.0

4,0 '6}

zierung kommt man zu der plastischen Zone nach

(Bild 7, b).

Die Belastung p = 178 N bei ideal elastoplastischem

Werkstoff führt zu Plastizierung über die ganze Breite der

Scheibe (Bild 7, c). In diesem Fall wurde nach 217 Itera-

tionen keine Konvergenz festgestellt. Die gleiche Bela-

stung ergibt, bei Werkstoff mit linearer Verfestigung

(E' = 2250 MPa), eine Fließzone nach Bild 7, d. Sie

wurde nach 156 Iterationen, davon 21 plastische, ge-

wonnen. Auf Bild 8, das dem Bild 2 aus [9] entspricht,

ist die mittlere nominale Zugspannung im Schnitt durch

die Öffnung (Bild 7)

Fz L 16p

“m A W

über die größte Längsdehnung 6 aufgetragen. Man sieht,

daß die hier berechneten Werte näher an den experimen-

tellen Ergebnissen von Theokaris und Marketos [1], [9]

liegen, als die nach anderen Verfahren [l], [9] gewon-

nenen Resultate.

LITERATUR

[ 1 J Zienkiewicz O. C.: Methode der finiten Elemente. Leipzig

19m 1,

[2 ] Mammal! H.: I'Ipnxnaruiaa reopwi macrntmocm

n nonsyqecm. Mocxsa, 1975

[3 ] Mmmoro X., Mnecn T.: Onpenenel-me Koatbnunenra

mrencnsuocm Hanpaxceimü um pacmmaaemoü

macmnbx C Ipeuumoü‚ Bbixozmuiefi Ha HOBerHOCTb.

B c6. Paclier ynpyrnx Koncrpyxunü c ncnonbsona-

HHeM 3BM, TOM l, Hemmrpan, 1974.

Bild 8

Scheibe mit Öffnung: Abhängigkeit der mittleren Zugspannung

von der maximalen Längsdehnung

experimentelle Ergebnisse von Theokaris und

Marketos (nach [9])

ideal elastoplastiseher Werkstoff, Ergebnisse

von Marcal und King (nach [9])

Werkstoff mit Verfestigung. Ergebnisse von

Marcal und King (nach [9])

O Werkstoff mit Verfestigung. Ergebnisse von King [9],

Werkstoff mit Verfestigung Ergebnisse der vorge-

schlagenen Methodik.

O

[4 ] Ralsten A.: A First Course in Numerical Analysis. New

York, 1965.

[5 ] Xpncron IL: I/Iscnensane Ha remneparypara, ripe-

Mecmammra u Hanpexcenfima, cbauauerm or Mm"-

Honeu mmeeI-I snap!)qu MBTO‘IHKK Ha rommna.

Kaummarcxa uneeprauua, Pyce, 1978.

[6 J Yettram A. L., Hirst M. J. S.: The Solution of Structural

Equilibrium Equations by the Conjugate Gradient Method

with Particular Reference to Plane Stress Analysis. Int. J.

Num. Meth. Eng. 7 (1971) p. 349 — 360.

[7] Bathe K. 1., Wilson E. L.: Numerical Methods in Finite

Element Analysis. New Jersey, 1976.

[8 l Argiris J. H., Scharpf D. W., SpoonerJ. 13.: Die elastopla-

stische Berechnung von allgemeinen Tragwerken und

Kontinua. Ing-Archiv 37 (1969), S. 326 — 352.

[9] Kim H. 0.: Finite Element Formulation Based on Na-

dai’s Deformation Theory for Elasto-plastic Analysis:

Int. J. Num. Meth. Eng. 17 (1981), p. 1861 — 1876.

[10] Altenbach J., Sacharov A. S.: Die Methode der finiten

Elemente in der Festkörpermechanik. Leipzig, 1982.

Anschrift der Verfasser:

Doz. Dr.-Ing. Doril Christow

Doz. Dr.-In.g. Marco Todorow

Technische Hochschule

Russe, VR Bulgarien

Lehrstuhl fiir Informatik

(.1

«I