1
00:00:00,040 --> 00:00:02,470
The following content is
provided under a Creative
2
00:00:02,470 --> 00:00:03,880
Commons license.
3
00:00:03,880 --> 00:00:06,920
Your support will help MIT
OpenCourseWare continue to
4
00:00:06,920 --> 00:00:10,570
offer high-quality educational
resources for free.
5
00:00:10,570 --> 00:00:13,470
To make a donation, or view
additional materials from
6
00:00:13,470 --> 00:00:17,872
hundreds of MIT courses, visit
MIT OpenCourseWare at
7
00:00:17,872 --> 00:00:19,122
ocw.mit.edu.
8
00:00:21,230 --> 00:00:22,980
PROFESSOR: Ladies and gentlemen,
welcome to this
9
00:00:22,980 --> 00:00:25,550
lecture on nonlinear finite
element analysis of solids and
10
00:00:25,550 --> 00:00:26,820
structures.
11
00:00:26,820 --> 00:00:29,390
In this lecture, I like to
discuss with you solution
12
00:00:29,390 --> 00:00:32,580
methods that we use to solve the
finite element equations
13
00:00:32,580 --> 00:00:34,790
in nonlinear static analysis.
14
00:00:34,790 --> 00:00:38,520
In the previous lectures, we
derived this set of equations,
15
00:00:38,520 --> 00:00:42,020
where tK is a tangent
stiffness matrix.
16
00:00:42,020 --> 00:00:46,360
Delta Ui is a nodal point
vector of incremental
17
00:00:46,360 --> 00:00:49,760
displacements corresponding
to iteration i.
18
00:00:49,760 --> 00:00:54,210
t plus delta t R is the
externally applied load vector
19
00:00:54,210 --> 00:00:58,260
corresponding to time
t plus delta t.
20
00:00:58,260 --> 00:01:02,920
And t plus delta t F i minus 1
is equal to the nodal point
21
00:01:02,920 --> 00:01:06,530
force vector corresponding to
the internal elements stresses
22
00:01:06,530 --> 00:01:09,120
at time t plus delta t,
and at the end of
23
00:01:09,120 --> 00:01:12,060
iteration i minus 1.
24
00:01:12,060 --> 00:01:15,770
The displacements are updated
as shown here.
25
00:01:15,770 --> 00:01:17,620
Delta Ui, of course,
is calculated
26
00:01:17,620 --> 00:01:19,450
from this set of equations.
27
00:01:19,450 --> 00:01:22,850
We add delta Ui to the previous
displacements that
28
00:01:22,850 --> 00:01:25,210
corresponded to iteration
t plus delta --
29
00:01:25,210 --> 00:01:28,720
to time t plus delta t and
iteration i minus 1.
30
00:01:28,720 --> 00:01:31,880
It's the end of iteration
i minus 1.
31
00:01:31,880 --> 00:01:37,110
And this right-hand side gives
us the new displacement vector
32
00:01:37,110 --> 00:01:42,340
corresponding to iteration
i, end of iteration i.
33
00:01:42,340 --> 00:01:45,680
Now this is one scheme that it's
used to solve the finite
34
00:01:45,680 --> 00:01:47,770
element equations.
35
00:01:47,770 --> 00:01:50,220
We refer to this as
the mortified
36
00:01:50,220 --> 00:01:51,940
Newton-Raphson method.
37
00:01:51,940 --> 00:01:53,790
However, there are other
schemes as well.
38
00:01:53,790 --> 00:01:55,600
And schemes that are in certain
39
00:01:55,600 --> 00:01:58,550
problems much more effective.
40
00:01:58,550 --> 00:02:01,770
Since it is so important to use
an effective method for
41
00:02:01,770 --> 00:02:04,430
the solution of the finite
element equations because the
42
00:02:04,430 --> 00:02:07,470
costs can otherwise be very
high, we should be familiar
43
00:02:07,470 --> 00:02:09,360
with those other techniques.
44
00:02:09,360 --> 00:02:13,980
And I'd like to share
those now with you.
45
00:02:13,980 --> 00:02:16,400
Therefore, we will look in
this lecture at other
46
00:02:16,400 --> 00:02:18,810
techniques to solve the finite
element equations.
47
00:02:18,810 --> 00:02:22,240
And we need, of course, to
discuss convergence criteria.
48
00:02:22,240 --> 00:02:23,980
It's very important to use
49
00:02:23,980 --> 00:02:25,880
appropriate convergence criteria.
50
00:02:25,880 --> 00:02:30,850
Otherwise, you iterate too long
and/or on the other side,
51
00:02:30,850 --> 00:02:33,930
you may actually stop iterating
before you have
52
00:02:33,930 --> 00:02:36,250
achieved proper convergence.
53
00:02:36,250 --> 00:02:41,310
So let me then go over to my
view graphs and let us start
54
00:02:41,310 --> 00:02:42,580
the discussion here.
55
00:02:42,580 --> 00:02:47,316
The basic set of equations that
we would like to solve
56
00:02:47,316 --> 00:02:48,530
are given here.
57
00:02:48,530 --> 00:02:52,510
Notice R is the vector of
externally applied nodal point
58
00:02:52,510 --> 00:02:55,390
forces at time t plus delta t.
59
00:02:55,390 --> 00:03:00,580
And F is the vector of nodal
point forces corresponding to
60
00:03:00,580 --> 00:03:04,560
the internal element stresses
at time t plus delta t.
61
00:03:04,560 --> 00:03:07,870
Of course, this vector is
unknown, and we want to
62
00:03:07,870 --> 00:03:13,560
iterate somehow to find it and
make sure, of course, that at
63
00:03:13,560 --> 00:03:17,150
the end of that iteration, if we
have the proper vector F, R
64
00:03:17,150 --> 00:03:21,670
is equal to F, or R minus
F is equal to 0.
65
00:03:21,670 --> 00:03:24,570
We assume that the loading is
deformation-independent.
66
00:03:24,570 --> 00:03:29,820
If say, loading is
deformation-dependent, we can
67
00:03:29,820 --> 00:03:31,250
also deal with that situation.
68
00:03:31,250 --> 00:03:34,890
But we will have to add
additional terms to our
69
00:03:34,890 --> 00:03:35,980
iterative scheme.
70
00:03:35,980 --> 00:03:38,650
And I will get back to
that a little later.
71
00:03:38,650 --> 00:03:41,940
Notice that F in the total
Lagrangian formulation the way
72
00:03:41,940 --> 00:03:47,030
we have discussed it earlier is
evaluated by this product
73
00:03:47,030 --> 00:03:50,290
here integrated over the
original volume for a single
74
00:03:50,290 --> 00:03:54,270
element and for the UL
formulation, F is calculated
75
00:03:54,270 --> 00:03:55,310
as shown here.
76
00:03:55,310 --> 00:03:58,500
Here, of course, we integrate
over the current volume, or
77
00:03:58,500 --> 00:04:02,580
the volume that we actually
want to calculate.
78
00:04:02,580 --> 00:04:05,360
Notice in the iteration, of
course, this volume here would
79
00:04:05,360 --> 00:04:12,020
be t plus delta t V i minus 1
as we will discuss just now.
80
00:04:12,020 --> 00:04:15,980
Because in the iteration, we
always update this integral
81
00:04:15,980 --> 00:04:20,089
and that integral with the
iteration counter that we are
82
00:04:20,089 --> 00:04:24,950
having on the right-hand side
of the system of equations.
83
00:04:24,950 --> 00:04:31,170
The methods that we use are
based on the Newton-Raphson
84
00:04:31,170 --> 00:04:38,020
method, which is really used
very abundantly to find roots
85
00:04:38,020 --> 00:04:40,260
of an equation.
86
00:04:40,260 --> 00:04:41,960
A small historical note:
87
00:04:41,960 --> 00:04:46,360
Newton gave a first version
of the method in 1669.
88
00:04:46,360 --> 00:04:49,560
Raphson then generalized and
simplified the method
89
00:04:49,560 --> 00:04:52,860
actually, in 1690.
90
00:04:52,860 --> 00:04:55,940
Both mathematicians used
the same concept.
91
00:04:55,940 --> 00:04:59,290
And both algorithms gave really,
the same results.
92
00:04:59,290 --> 00:05:03,090
But it is very appropriate to
refer to the methods that
93
00:05:03,090 --> 00:05:06,490
we're using as the
Newton-Raphson method because
94
00:05:06,490 --> 00:05:11,010
Raphson really contributed quite
a bit to that method.
95
00:05:11,010 --> 00:05:14,110
So let us now consider a single
96
00:05:14,110 --> 00:05:16,220
Newton-Raphson iteration.
97
00:05:16,220 --> 00:05:19,030
The way you may have actually
encountered it already earlier
98
00:05:19,030 --> 00:05:23,220
in your studies of solution
methods for the roots of an
99
00:05:23,220 --> 00:05:24,730
algebraic equation.
100
00:05:24,730 --> 00:05:27,070
Basically, what you're
doing is this.
101
00:05:27,070 --> 00:05:32,510
You're saying if you have xi
minus 1, you calculate
102
00:05:32,510 --> 00:05:34,660
f at xi minus 1.
103
00:05:34,660 --> 00:05:39,900
You divide this value by the
f prime at xi minus 1.
104
00:05:39,900 --> 00:05:43,910
And then this right-hand side
gives you better value, a
105
00:05:43,910 --> 00:05:47,340
better approximation to the
root of the equation.
106
00:05:47,340 --> 00:05:50,540
Once you have xi, you repeat
the process and
107
00:05:50,540 --> 00:05:53,090
you get xi plus 1.
108
00:05:53,090 --> 00:05:58,630
You keep on repeating the
process until, basically f xi
109
00:05:58,630 --> 00:06:01,070
here is close to 0.
110
00:06:01,070 --> 00:06:03,980
Because then you have a root.
111
00:06:03,980 --> 00:06:17,470
Well, if we use that
Newton-Raphson formula, it is
112
00:06:17,470 --> 00:06:21,010
quite interesting to see how
it has been derived.
113
00:06:21,010 --> 00:06:25,130
We can write for any point xi,
a neighboring point xi minus
114
00:06:25,130 --> 00:06:29,970
1, directly this equation here
by a Taylor series expansion.
115
00:06:29,970 --> 00:06:33,450
And if we neglect the higher
order terms, we get that f of
116
00:06:33,450 --> 00:06:37,180
xi is approximately equal to
this relationship on the
117
00:06:37,180 --> 00:06:39,220
right-hand side.
118
00:06:39,220 --> 00:06:44,335
Well, since f of xi is supposed
to be 0, because we
119
00:06:44,335 --> 00:06:47,060
are looking for a 0 of
the equation, we
120
00:06:47,060 --> 00:06:48,430
set it equal to 0.
121
00:06:48,430 --> 00:06:50,860
And here you see directly
the formula that
122
00:06:50,860 --> 00:06:52,120
I showed you earlier.
123
00:06:52,120 --> 00:06:55,360
So this is a very quick
derivation of the
124
00:06:55,360 --> 00:06:59,020
Newton-Raphson procedure.
125
00:06:59,020 --> 00:07:03,040
Let us look at a mathematical
example to see how the method
126
00:07:03,040 --> 00:07:07,370
works, and to get some insight
into the technique.
127
00:07:07,370 --> 00:07:10,960
Let us choose a very simple
example, not much to do with
128
00:07:10,960 --> 00:07:12,680
finite element analysis.
129
00:07:12,680 --> 00:07:16,280
Let us say that f of x equals
sine x and that our starting
130
00:07:16,280 --> 00:07:19,080
value is equal to 2
for the iteration.
131
00:07:19,080 --> 00:07:21,770
You always have to choose, of
course, a starting value for
132
00:07:21,770 --> 00:07:23,030
the iteration.
133
00:07:23,030 --> 00:07:28,750
Then, in this column you see the
values that are calculated
134
00:07:28,750 --> 00:07:30,830
in the successive iterations.
135
00:07:30,830 --> 00:07:34,020
And we also show
here the error.
136
00:07:34,020 --> 00:07:37,700
It is interesting to observe
that when you are close to the
137
00:07:37,700 --> 00:07:40,210
root, you have quadratic
convergence.
138
00:07:40,210 --> 00:07:45,410
Meaning that the error here,
epsilon becomes an error.
139
00:07:45,410 --> 00:07:48,570
Epsilon squared here.
140
00:07:48,570 --> 00:07:51,880
This is summarized on
this view graph.
141
00:07:51,880 --> 00:07:58,790
Mathematically, if we have that
the error, Ei minus 1, is
142
00:07:58,790 --> 00:08:00,220
given by this equation here.
143
00:08:00,220 --> 00:08:02,200
Of course, there's an
approximate sign here.
144
00:08:02,200 --> 00:08:04,340
Then the error in the
next iteration
145
00:08:04,340 --> 00:08:07,906
will be of this order.
146
00:08:07,906 --> 00:08:12,590
So the convergence is really
very rapid once we are close
147
00:08:12,590 --> 00:08:15,660
to the root.
148
00:08:15,660 --> 00:08:19,170
Here we can see, however,
that if we are not
149
00:08:19,170 --> 00:08:21,300
close to the root--
150
00:08:21,300 --> 00:08:26,690
in other words, if we are too
far from the root, and we pick
151
00:08:26,690 --> 00:08:31,390
a bad value, then we don't
get to the desired
152
00:08:31,390 --> 00:08:34,080
root, which is pi.
153
00:08:34,080 --> 00:08:38,919
But we get a value that
is far away from pi.
154
00:08:38,919 --> 00:08:42,159
It is also root, but certainly
not the root that we are
155
00:08:42,159 --> 00:08:43,635
interested in.
156
00:08:43,635 --> 00:08:47,550
So Newton-Raphson iteration
is not always convergent
157
00:08:47,550 --> 00:08:51,420
directly, converging directly
to the result that we would
158
00:08:51,420 --> 00:08:52,960
like to obtain.
159
00:08:52,960 --> 00:08:57,090
There has to be some care when
we use that iterative scheme.
160
00:08:57,090 --> 00:09:02,720
And in fact, we have to be close
enough just to basically
161
00:09:02,720 --> 00:09:04,190
say what has to happen.
162
00:09:04,190 --> 00:09:07,110
We have to be close enough to
the root in order to have
163
00:09:07,110 --> 00:09:09,930
convergence to the root that
we're looking for.
164
00:09:09,930 --> 00:09:12,750
And ultimately, if we do
converge to that root, we will
165
00:09:12,750 --> 00:09:15,060
get quadratic convergence
to the root.
166
00:09:15,060 --> 00:09:17,420
And that is, of course,
very desirable.
167
00:09:17,420 --> 00:09:21,040
But the starting value is
critical, as you can see here
168
00:09:21,040 --> 00:09:22,400
from this view graph.
169
00:09:22,400 --> 00:09:27,050
Let us look pictorially at
the solution process.
170
00:09:27,050 --> 00:09:30,020
Here we have our
sine function.
171
00:09:30,020 --> 00:09:31,770
And we are looking
for this root.
172
00:09:34,280 --> 00:09:41,170
In iteration 1, we basically set
up a tangent to that sine
173
00:09:41,170 --> 00:09:44,380
function at the starting
value, x0.
174
00:09:44,380 --> 00:09:52,140
And we calculate x1 using this
tangent as another estimate to
175
00:09:52,140 --> 00:09:55,820
the value that we're
interested in.
176
00:09:55,820 --> 00:10:03,540
Now in the iteration 2, we lay
a tangent to the f function.
177
00:10:03,540 --> 00:10:06,520
Which, of course, is sine
x, at the new value
178
00:10:06,520 --> 00:10:09,050
corresponding to x1.
179
00:10:09,050 --> 00:10:12,520
And with this tangent,
we get to this value.
180
00:10:12,520 --> 00:10:15,150
And you can see already
that we got closer
181
00:10:15,150 --> 00:10:16,220
to the desired value.
182
00:10:16,220 --> 00:10:21,070
We were this far away
originally.
183
00:10:21,070 --> 00:10:23,460
Then that far away.
184
00:10:23,460 --> 00:10:26,910
And now we are only
that far away.
185
00:10:26,910 --> 00:10:29,950
Notice our sine function
has this value here.
186
00:10:29,950 --> 00:10:32,680
And certainly, it's not 0 yet.
187
00:10:32,680 --> 00:10:37,540
So we lay another tangent
at that point.
188
00:10:37,540 --> 00:10:42,120
And that's being done on
this view graph here.
189
00:10:42,120 --> 00:10:45,890
Another tangent, and now
we get much closer
190
00:10:45,890 --> 00:10:46,990
to the desired value.
191
00:10:46,990 --> 00:10:50,000
Right here, x3 is already
quite close.
192
00:10:50,000 --> 00:10:54,340
But not close enough as
shown by this value.
193
00:10:54,340 --> 00:10:59,380
It's not yet close enough
to 0, this length here.
194
00:10:59,380 --> 00:11:03,830
And we lay another tangent
to that point now.
195
00:11:03,830 --> 00:11:08,500
And that's shown on this
view graph now.
196
00:11:08,500 --> 00:11:12,110
And we can see with the blue
line being the tangent at this
197
00:11:12,110 --> 00:11:18,470
point, we get x4 very close
to the 0 value that we are
198
00:11:18,470 --> 00:11:20,770
actually interested in.
199
00:11:20,770 --> 00:11:23,840
Certainly on this view graph,
you can't see any difference
200
00:11:23,840 --> 00:11:28,290
anymore between x4 and
the desired value.
201
00:11:28,290 --> 00:11:31,220
So this is basically the process
that we are following
202
00:11:31,220 --> 00:11:34,510
through when we go from
iteration 1 to iteration 2 to
203
00:11:34,510 --> 00:11:38,540
iteration 3 and to
iteration 4.
204
00:11:38,540 --> 00:11:41,580
Beautiful convergence in this
particular example.
205
00:11:41,580 --> 00:11:47,290
And the reason being that x0 is
close enough to the desired
206
00:11:47,290 --> 00:11:50,730
value, the desired root.
207
00:11:50,730 --> 00:11:52,650
Which, itself, of course,
is close to x4.
208
00:11:52,650 --> 00:11:53,840
Four
209
00:11:53,840 --> 00:12:01,420
Well, if however, we were to
take a value that is not close
210
00:12:01,420 --> 00:12:07,420
enough and then we, as I showed
earlier, we do not
211
00:12:07,420 --> 00:12:11,050
converge to the actual root.
212
00:12:11,050 --> 00:12:13,020
And this is shown here.
213
00:12:13,020 --> 00:12:18,500
For example, if we take a
tangent at a point where f of
214
00:12:18,500 --> 00:12:24,920
prime is almost 0, then of
course, this tangent throws us
215
00:12:24,920 --> 00:12:27,080
far away from this root.
216
00:12:27,080 --> 00:12:31,220
And we will converge
to some other root.
217
00:12:31,220 --> 00:12:34,120
Of course, we don't want to use
an exact 0 here because we
218
00:12:34,120 --> 00:12:37,310
have f of x divided
by f prime of x.
219
00:12:37,310 --> 00:12:40,080
So with an exact 0, we could
not [UNINTELLIGIBLE].
220
00:12:40,080 --> 00:12:45,330
But a tangent that is close to
0, an f prime value that is
221
00:12:45,330 --> 00:12:49,220
close to 0, would throw us far
away from the desired root.
222
00:12:49,220 --> 00:12:53,180
And that shows some of the
difficulties using--
223
00:12:53,180 --> 00:12:59,040
in fact, almost all iterative
schemes that we may not get
224
00:12:59,040 --> 00:13:00,750
fast enough to our root.
225
00:13:00,750 --> 00:13:03,520
And in fact, we may never get
to the root that we are
226
00:13:03,520 --> 00:13:05,010
looking for.
227
00:13:05,010 --> 00:13:09,380
Well, let us then look at the
Newton-Raphson iteration for
228
00:13:09,380 --> 00:13:12,130
multiple degrees of freedom.
229
00:13:12,130 --> 00:13:15,350
Here we have our R minus
F that we want to
230
00:13:15,350 --> 00:13:17,700
be setting to 0.
231
00:13:17,700 --> 00:13:21,130
Or we want to rather, solve for
the displacements t plus
232
00:13:21,130 --> 00:13:27,400
delta t U, which are the
solution because this t plus
233
00:13:27,400 --> 00:13:32,060
delta t u displacements would
give us this force vector.
234
00:13:32,060 --> 00:13:36,610
And this force vector
equilibrates the R vector.
235
00:13:36,610 --> 00:13:41,580
Meaning f of U, this little
f of U, is equal to 0.
236
00:13:41,580 --> 00:13:44,990
And that, of course, is what
we want to achieve.
237
00:13:44,990 --> 00:13:47,130
Notice that we have here,
of course, now
238
00:13:47,130 --> 00:13:48,720
n degrees of freedom.
239
00:13:48,720 --> 00:13:51,660
Therefore, n equations
to solve.
240
00:13:51,660 --> 00:13:55,290
Previously, we only looked
at one equation.
241
00:13:55,290 --> 00:14:01,070
If we apply the same principle
as before to this little f,
242
00:14:01,070 --> 00:14:06,840
namely we expand f as a Taylor
series about t plus delta t U
243
00:14:06,840 --> 00:14:15,000
i minus 1, then we get
this equation on
244
00:14:15,000 --> 00:14:17,770
the right-hand side.
245
00:14:17,770 --> 00:14:20,300
And of course, on the left-hand
side, we have f of t
246
00:14:20,300 --> 00:14:22,120
plus delta t Ui.
247
00:14:22,120 --> 00:14:26,660
Notice they are higher terms,
which we will neglect just as
248
00:14:26,660 --> 00:14:28,630
we have done earlier for
the single equation.
249
00:14:28,630 --> 00:14:31,650
Because we are looking for
Taylor series expansion about
250
00:14:31,650 --> 00:14:34,790
t plus delta t U i minus 1.
251
00:14:34,790 --> 00:14:41,560
If we neglect this part here,
we obtain directly this
252
00:14:41,560 --> 00:14:44,010
equation here.
253
00:14:44,010 --> 00:14:47,020
0 on the left-hand side, because
once again, we are
254
00:14:47,020 --> 00:14:54,440
looking for the displacements
values for which f is 0.
255
00:14:54,440 --> 00:14:56,770
So therefore, we set this
deliberately equal to 0.
256
00:14:56,770 --> 00:15:00,750
And on the right-hand side, we
have f t plus delta t U i
257
00:15:00,750 --> 00:15:04,540
minus 1 plus partial f with
respect to U times the delta
258
00:15:04,540 --> 00:15:07,190
U. This delta U is unknown,
of course,
259
00:15:07,190 --> 00:15:08,420
that we want to calculate.
260
00:15:08,420 --> 00:15:12,910
And this delta U is nothing
else then the difference
261
00:15:12,910 --> 00:15:16,610
between the U value in iteration
i and the U value in
262
00:15:16,610 --> 00:15:19,140
iteration i minus 1.
263
00:15:19,140 --> 00:15:22,670
Notice that this f is
evaluated at t plus
264
00:15:22,670 --> 00:15:24,350
delta t U i minus 1.
265
00:15:27,510 --> 00:15:29,530
Let us look now at
this equation.
266
00:15:29,530 --> 00:15:36,620
And if we write it a little bit
more out, we find that on
267
00:15:36,620 --> 00:15:40,300
the left-hand side, we have,
of course, a vector of 0's.
268
00:15:40,300 --> 00:15:43,990
On the right-hand side, we
have a vector of f i
269
00:15:43,990 --> 00:15:47,070
components, f1 to fn.
270
00:15:47,070 --> 00:15:52,600
All of which are evaluated at
t plus delta t U i minus 1.
271
00:15:52,600 --> 00:15:55,940
And then we have here a matrix,
which is a square
272
00:15:55,940 --> 00:16:02,560
matrix in which the individual
elements are partials of f i
273
00:16:02,560 --> 00:16:04,240
with respect to Uj.
274
00:16:04,240 --> 00:16:09,020
For example, f1, partial f1
with respect to U1 here.
275
00:16:09,020 --> 00:16:13,990
Partial fn with respect to
U n here, et cetera.
276
00:16:13,990 --> 00:16:18,130
This matrix here will give us
a tangent stiffness matrix.
277
00:16:18,130 --> 00:16:21,120
Notice this matrix is multiplied
by the vector of
278
00:16:21,120 --> 00:16:24,090
incremental nodal point
displacements corresponding to
279
00:16:24,090 --> 00:16:25,340
iteration i.
280
00:16:28,840 --> 00:16:37,040
If we remember that f, the
little f, at t plus delta t U
281
00:16:37,040 --> 00:16:41,450
i minus 1 is equal to R minus
2 plus delta t F i minus 1.
282
00:16:41,450 --> 00:16:44,980
Capital F here, little
f there, remember?
283
00:16:44,980 --> 00:16:50,230
Then, the partial of little f
with respect to U at t plus
284
00:16:50,230 --> 00:16:53,733
delta t U i minus 1 is given
via this equation.
285
00:16:53,733 --> 00:16:58,720
Now, please recognize that
partial of R with respect to U
286
00:16:58,720 --> 00:17:03,870
is equal to 0 if the loads are
deformation-independent.
287
00:17:03,870 --> 00:17:06,770
Of course, if the loads depend
on the deformations, in other
288
00:17:06,770 --> 00:17:10,670
words, if the loads depend on
the U displacements, then this
289
00:17:10,670 --> 00:17:15,010
would not be 0, and we would
have to actually put a term in
290
00:17:15,010 --> 00:17:17,979
here, carry that term along
in the solution of
291
00:17:17,979 --> 00:17:20,069
the nonlinear equations.
292
00:17:20,069 --> 00:17:24,050
But for the moment,
we neglect this.
293
00:17:24,050 --> 00:17:27,210
We assume that the loads are
deformation-independent.
294
00:17:27,210 --> 00:17:30,230
And then, this term
is equal to 0.
295
00:17:30,230 --> 00:17:33,770
This here is now giving us the
tangent stiffness matrix.
296
00:17:33,770 --> 00:17:37,430
And it is written down here.
297
00:17:37,430 --> 00:17:39,870
The tangent stiffness matrix
is nothing else than the
298
00:17:39,870 --> 00:17:45,440
partials of capital F with
respect to the U's, to the
299
00:17:45,440 --> 00:17:46,690
displacements.
300
00:17:48,520 --> 00:17:51,990
The final result then is given
on this view graph.
301
00:17:51,990 --> 00:17:55,540
We substitute all the
information that I shared with
302
00:17:55,540 --> 00:18:00,790
you into the Taylor series
expansion around t plus delta
303
00:18:00,790 --> 00:18:02,470
t U i minus 1.
304
00:18:02,470 --> 00:18:05,110
We get directly this
equation here.
305
00:18:05,110 --> 00:18:10,300
Notice tangent stiffness matrix
corresponding to time t
306
00:18:10,300 --> 00:18:14,360
plus delta t and iteration
i minus 1.
307
00:18:14,360 --> 00:18:20,580
Delta U i corresponding to
iteration i we have to be very
308
00:18:20,580 --> 00:18:23,330
clear about this, that this is
an increment in displacement
309
00:18:23,330 --> 00:18:28,460
from time t plus delta t
iteration i minus 1 to
310
00:18:28,460 --> 00:18:30,220
iteration i.
311
00:18:30,220 --> 00:18:33,380
And of course, the nodal
point force vector.
312
00:18:33,380 --> 00:18:36,020
Externally applied nodal point
forces go in here.
313
00:18:36,020 --> 00:18:38,470
And this is the nodal point
force vector corresponding to
314
00:18:38,470 --> 00:18:43,090
the internal element stresses at
time t plus delta t and at
315
00:18:43,090 --> 00:18:45,090
the end of iteration
i minus 1.
316
00:18:45,090 --> 00:18:49,100
Evaluated differently in the
total Lagrangian formulation
317
00:18:49,100 --> 00:18:52,480
from the way we're evaluating
it in the updated Lagrangian
318
00:18:52,480 --> 00:18:53,070
formulation.
319
00:18:53,070 --> 00:18:55,690
We talked about this vector
abundantly in
320
00:18:55,690 --> 00:18:57,330
the previous lectures.
321
00:18:57,330 --> 00:19:00,510
We talked about this vector,
of course, also in the
322
00:19:00,510 --> 00:19:01,630
previous lectures.
323
00:19:01,630 --> 00:19:05,970
We had in the previous lectures
here at tK, a
324
00:19:05,970 --> 00:19:10,260
constant tangent stiffness
matrix, which was set up at
325
00:19:10,260 --> 00:19:14,210
the beginning of this whole
interactive process, and was
326
00:19:14,210 --> 00:19:16,030
never updated.
327
00:19:16,030 --> 00:19:19,850
Well, here we now update
that stiffness matrix.
328
00:19:19,850 --> 00:19:23,650
And I think if you look through
the information that I
329
00:19:23,650 --> 00:19:26,080
just discussed with you,
you recognize why
330
00:19:26,080 --> 00:19:27,360
we're updating it.
331
00:19:27,360 --> 00:19:29,840
We are doing so because we are
always starting with the new
332
00:19:29,840 --> 00:19:37,120
Taylor series expansion about
the point i minus 1.
333
00:19:37,120 --> 00:19:41,990
This set of equations if of
course solved for delta Ui.
334
00:19:41,990 --> 00:19:44,070
We add delta Ui to what we
had already in terms of
335
00:19:44,070 --> 00:19:48,020
displacements to get our new
estimate for the nodal point
336
00:19:48,020 --> 00:19:50,160
displacements.
337
00:19:50,160 --> 00:19:54,450
It is important to realize that
the K matrix, which we
338
00:19:54,450 --> 00:19:57,440
are using in the solution
process, is symmetric.
339
00:19:57,440 --> 00:20:01,190
Because first of all, we use
symmetric stress and strain
340
00:20:01,190 --> 00:20:03,900
measures in our governing
equation.
341
00:20:03,900 --> 00:20:06,990
Remember that when we apply the
principle of virtual work,
342
00:20:06,990 --> 00:20:10,760
we use Cauchy stresses
and infinitesimally
343
00:20:10,760 --> 00:20:12,250
small virtual strains.
344
00:20:12,250 --> 00:20:16,850
Now notice that both of these
tenses are symmetric tenses.
345
00:20:16,850 --> 00:20:19,980
Of course, these are the tenses
we're using in the UL,
346
00:20:19,980 --> 00:20:22,280
in the updated Lagrangian
formulation.
347
00:20:22,280 --> 00:20:25,180
In the total Lagrangian
formulation, we use the second
348
00:20:25,180 --> 00:20:28,340
Piola-Kirchhoff stress in
Green-Lagrange strain.
349
00:20:28,340 --> 00:20:32,350
Both of these measures are
again, symmetric measures.
350
00:20:32,350 --> 00:20:34,650
Both tenses are symmetric
tenses.
351
00:20:34,650 --> 00:20:38,880
If we had introduced for a
formulation non-symmetric
352
00:20:38,880 --> 00:20:42,210
tenses, for example a
non-symmetric stress tenser,
353
00:20:42,210 --> 00:20:44,960
and of course, an energy
conjugate non-symmetric strain
354
00:20:44,960 --> 00:20:48,780
tenser, we would have obtained
a non-symmetric K matrix,
355
00:20:48,780 --> 00:20:53,500
which would be much more
difficult to deal with.
356
00:20:53,500 --> 00:20:57,680
Much less cost effective in
the solution process.
357
00:20:57,680 --> 00:21:01,470
Also, please recognize that
we interpolated the real
358
00:21:01,470 --> 00:21:05,810
displacements and the virtual
displacements with exactly the
359
00:21:05,810 --> 00:21:07,820
same functions.
360
00:21:07,820 --> 00:21:11,020
Whereas this point here is a
continuum mechanics point,
361
00:21:11,020 --> 00:21:13,680
this is really a finite
element point.
362
00:21:13,680 --> 00:21:17,590
If we had used different
interpolation functions for
363
00:21:17,590 --> 00:21:21,930
the real displacements, then for
the virtual displacements,
364
00:21:21,930 --> 00:21:26,300
or vice-versa, then we would
have not necessarily obtained
365
00:21:26,300 --> 00:21:27,860
a symmetric K matrix.
366
00:21:27,860 --> 00:21:32,110
Of course, the most natural
procedure is to use the same
367
00:21:32,110 --> 00:21:35,460
kind of interpolation functions
for the virtual
368
00:21:35,460 --> 00:21:37,940
displacements as we used for
the real displacements, and
369
00:21:37,940 --> 00:21:39,340
that's what we did.
370
00:21:39,340 --> 00:21:42,460
Finally, we also assumed
that the loading was
371
00:21:42,460 --> 00:21:43,710
deformation-independent.
372
00:21:45,880 --> 00:21:49,000
If we have deformation-dependent
loading,
373
00:21:49,000 --> 00:21:52,680
then if you go more back to the
earlier view graph, then
374
00:21:52,680 --> 00:21:57,020
of course, this right-hand
side vector R here would
375
00:21:57,020 --> 00:21:59,750
depend on the displacements.
376
00:21:59,750 --> 00:22:02,630
And there are basically two
different approaches
377
00:22:02,630 --> 00:22:04,310
that one can take.
378
00:22:04,310 --> 00:22:08,480
In the first approach, one
simply updates this vector
379
00:22:08,480 --> 00:22:11,770
with the iteration, or with
taking into account the
380
00:22:11,770 --> 00:22:15,910
current or last calculated
volume and surface areas for
381
00:22:15,910 --> 00:22:16,920
the elements.
382
00:22:16,920 --> 00:22:21,430
So we would simply put here
an i minus 1 up there.
383
00:22:21,430 --> 00:22:24,660
The left-hand side matrix,
we leave unchanged.
384
00:22:24,660 --> 00:22:27,630
And if that converges fast,
certainly it's a very
385
00:22:27,630 --> 00:22:29,540
effective approach to use.
386
00:22:29,540 --> 00:22:32,950
We use that approach
abundantly.
387
00:22:32,950 --> 00:22:36,490
However, another approach
is to actually take the
388
00:22:36,490 --> 00:22:39,980
differentiation of R with
respect to the U's, the way we
389
00:22:39,980 --> 00:22:41,890
have been looking at it earlier
on an earlier view
390
00:22:41,890 --> 00:22:45,260
graph, and then you get some
components that you
391
00:22:45,260 --> 00:22:47,110
add to the K matrix.
392
00:22:47,110 --> 00:22:51,200
And that K matrix may then
out to be non-symmetric.
393
00:22:51,200 --> 00:22:54,400
However, if the loading is
deformation-independent, then
394
00:22:54,400 --> 00:22:57,910
the differentiation of R with
respect to the displacements U
395
00:22:57,910 --> 00:23:01,770
is equal to 0, and there's
no component from that
396
00:23:01,770 --> 00:23:04,250
differentiation coming
into the K matrix.
397
00:23:04,250 --> 00:23:07,550
And our K, of course, provided
these are also
398
00:23:07,550 --> 00:23:10,850
satisfied, is symmetric.
399
00:23:10,850 --> 00:23:13,840
The iterative scheme that we
just discussed is referred to
400
00:23:13,840 --> 00:23:16,780
as the full Newton-Raphson
method.
401
00:23:16,780 --> 00:23:19,630
Full because we are setting
up a new K matrix at the
402
00:23:19,630 --> 00:23:22,340
beginning of each iteration.
403
00:23:22,340 --> 00:23:26,120
The full Newton-Raphson method
shows mathematically quadratic
404
00:23:26,120 --> 00:23:28,650
convergence the way we discussed
it a bit earlier in
405
00:23:28,650 --> 00:23:31,640
the lecture.
406
00:23:31,640 --> 00:23:36,400
And that, of course, is always
the case provided you are
407
00:23:36,400 --> 00:23:38,050
close enough to the root.
408
00:23:38,050 --> 00:23:41,550
This quadratic convergence only
holds provided you are
409
00:23:41,550 --> 00:23:45,780
close enough to the root when
you solve your equations.
410
00:23:45,780 --> 00:23:49,420
In finite element analysis, it
is also important to recognize
411
00:23:49,420 --> 00:23:52,410
that a number of requirements
must be fulfilled.
412
00:23:52,410 --> 00:23:57,410
For example, in elasto-plastic
analysis, the stresses and
413
00:23:57,410 --> 00:23:59,020
strains must be properly--
414
00:23:59,020 --> 00:24:03,520
plastic strains must be
properly updated.
415
00:24:03,520 --> 00:24:07,050
And similarly, the rotations
in a shell analysis must be
416
00:24:07,050 --> 00:24:07,870
properly updated.
417
00:24:07,870 --> 00:24:11,380
So it is not necessarily the
case that you automatically
418
00:24:11,380 --> 00:24:14,780
get quadratic convergence when
you do finite element
419
00:24:14,780 --> 00:24:18,250
announces with a full
Newton-Raphson method.
420
00:24:18,250 --> 00:24:22,120
It is very important to also,
on the level of updating the
421
00:24:22,120 --> 00:24:27,800
stresses and the rotations,
really do things properly in
422
00:24:27,800 --> 00:24:33,040
quotes in order to obtain the
full quadratic convergence of
423
00:24:33,040 --> 00:24:36,900
the Newton-Raphson method.
424
00:24:36,900 --> 00:24:39,580
We can depict the iteration
process in
425
00:24:39,580 --> 00:24:41,350
two equivalent ways.
426
00:24:41,350 --> 00:24:45,300
The first way is shown
here, the left.
427
00:24:45,300 --> 00:24:48,520
And it's really the way we've
been discussing just now the
428
00:24:48,520 --> 00:24:54,280
solution of f equals R minus
capital F. Where we want to
429
00:24:54,280 --> 00:24:58,150
solve for the root, the
0 of this equation.
430
00:24:58,150 --> 00:25:02,980
Notice here we have in red
the little f depicted.
431
00:25:02,980 --> 00:25:07,280
Notice that at this point here,
we have t plus delta t
432
00:25:07,280 --> 00:25:10,650
capital F i minus 1.
433
00:25:10,650 --> 00:25:14,180
And t plus delta t
R is this value.
434
00:25:14,180 --> 00:25:18,130
Now, as we get closer to the
root, which is of course, the
435
00:25:18,130 --> 00:25:23,800
point where the red curve
crosses this U-axis, as we get
436
00:25:23,800 --> 00:25:28,620
closer to the root, this capital
F will get closer and
437
00:25:28,620 --> 00:25:33,140
closer to the R. And that is,
of course, when the little f
438
00:25:33,140 --> 00:25:34,310
is equal to 0.
439
00:25:34,310 --> 00:25:36,080
That's all we are
showing here.
440
00:25:36,080 --> 00:25:41,830
Notice that we are setting up
a slope f prime at the point
441
00:25:41,830 --> 00:25:44,330
corresponding to the i minus
first iteration.
442
00:25:44,330 --> 00:25:49,180
And that slope brings us into
this point, which will be the
443
00:25:49,180 --> 00:25:51,240
next point for--
444
00:25:51,240 --> 00:25:55,330
or the next starting point
for the next iteration.
445
00:25:55,330 --> 00:25:57,930
I should say, the point
of starting
446
00:25:57,930 --> 00:25:59,600
with the next iteration.
447
00:25:59,600 --> 00:26:02,890
Now, this is one way
of looking at
448
00:26:02,890 --> 00:26:04,420
the iterative process.
449
00:26:04,420 --> 00:26:06,960
We can also look at
it as shown here.
450
00:26:06,960 --> 00:26:10,830
Notice we have here the R
plotted now horizontally.
451
00:26:10,830 --> 00:26:13,240
The displacement vertically.
452
00:26:13,240 --> 00:26:19,680
Notice the R at a particular
time t plus delta t is shown
453
00:26:19,680 --> 00:26:22,160
by this dashed line.
454
00:26:22,160 --> 00:26:27,110
And we really want to find this
particular point here and
455
00:26:27,110 --> 00:26:30,970
the corresponding U displacement
of course.
456
00:26:30,970 --> 00:26:34,380
At this point, the little f is
0 and the corresponding U
457
00:26:34,380 --> 00:26:36,250
displacement is down here.
458
00:26:36,250 --> 00:26:40,240
Notice that at the iteration
i minus 1, it's the end of
459
00:26:40,240 --> 00:26:44,570
iteration i minus 1, we have
obtained this point here.
460
00:26:44,570 --> 00:26:47,850
The U displacement corresponding
to that point is
461
00:26:47,850 --> 00:26:52,180
obtained by projecting down
on the displacement axes.
462
00:26:52,180 --> 00:26:57,800
And this slope here, the blue
line, gives us a tangent
463
00:26:57,800 --> 00:27:00,350
stiffness matrix slope.
464
00:27:00,350 --> 00:27:02,700
These are two equivalent ways.
465
00:27:02,700 --> 00:27:06,070
This is here more like a force
deflection curve, and we will
466
00:27:06,070 --> 00:27:10,270
use this representation
now abundantly.
467
00:27:10,270 --> 00:27:13,500
But keep in mind, it's really
the same thing if you look at
468
00:27:13,500 --> 00:27:15,300
it closely.
469
00:27:15,300 --> 00:27:17,890
Well, modifications
to the iterative
470
00:27:17,890 --> 00:27:21,410
scheme are the following.
471
00:27:21,410 --> 00:27:25,100
If we leave the stiffness matrix
constant throughout a
472
00:27:25,100 --> 00:27:26,670
complete solution process--
473
00:27:26,670 --> 00:27:30,070
in other words, tau is
equal to 0 here.
474
00:27:30,070 --> 00:27:31,390
It's never updated.
475
00:27:31,390 --> 00:27:34,570
We talk about the initial
stress method.
476
00:27:34,570 --> 00:27:39,610
If tau is equal to t, where t,
of course, corresponds to a
477
00:27:39,610 --> 00:27:44,670
particular solution step and
meaning that the K matrix is
478
00:27:44,670 --> 00:27:48,620
constant, but it is updated
as the beginning
479
00:27:48,620 --> 00:27:50,010
of a solution step.
480
00:27:50,010 --> 00:27:53,160
Then we talk about the modified
Newton method.
481
00:27:53,160 --> 00:27:55,860
Modified Newton-Raphson
method.
482
00:27:55,860 --> 00:27:59,670
Or, we may also find it
effective to update the K
483
00:27:59,670 --> 00:28:04,040
matrix only at particular
solution steps, at certain
484
00:28:04,040 --> 00:28:06,570
times only.
485
00:28:06,570 --> 00:28:09,590
We note that the initial stress
method and the modified
486
00:28:09,590 --> 00:28:13,810
Newton methods are, of course,
much less expensive than the
487
00:28:13,810 --> 00:28:16,646
full Newton method per
iteration, however.
488
00:28:16,646 --> 00:28:21,070
We should add that many more
iterations are necessary to
489
00:28:21,070 --> 00:28:27,150
achieve the same accuracy if we
don't set up a new K matrix
490
00:28:27,150 --> 00:28:29,210
in every iteration.
491
00:28:29,210 --> 00:28:32,030
The initial stress method and
modified Newton iteration
492
00:28:32,030 --> 00:28:36,240
technique do not exhibit
quadratic convergence because
493
00:28:36,240 --> 00:28:38,770
to obtain quadratic convergence,
you need to set
494
00:28:38,770 --> 00:28:41,750
up a new K matrix in
each iteration.
495
00:28:41,750 --> 00:28:44,920
And of course, you have to be
close enough to the root, to
496
00:28:44,920 --> 00:28:46,860
the solution that you're
looking for.
497
00:28:46,860 --> 00:28:51,200
Let us now look at
a simple example.
498
00:28:51,200 --> 00:28:55,250
An example where we have just
one degree of freedom.
499
00:28:55,250 --> 00:28:57,790
And where we deal with
two load steps.
500
00:28:57,790 --> 00:29:01,580
Here we show the force
applied to the
501
00:29:01,580 --> 00:29:04,210
problem, or to the structure.
502
00:29:04,210 --> 00:29:07,460
In the first load step,
we apply 1R.
503
00:29:07,460 --> 00:29:10,530
And then in the second
load step, 2R.
504
00:29:10,530 --> 00:29:14,570
Notice horizontally here, I'm
plotting displacements.
505
00:29:14,570 --> 00:29:17,460
And the force displacement
curve is
506
00:29:17,460 --> 00:29:20,130
shown by the red line.
507
00:29:20,130 --> 00:29:26,020
Now, to obtain the solution for
this displacement, U1, and
508
00:29:26,020 --> 00:29:31,340
that displacement U2, we, as
we discussed, set up a K
509
00:29:31,340 --> 00:29:34,600
matrix, which corresponds
to the slope, a
510
00:29:34,600 --> 00:29:36,940
slope on the red curve.
511
00:29:36,940 --> 00:29:40,490
And we iterate towards
the sought solutions.
512
00:29:40,490 --> 00:29:44,680
Here, we have the iterative
process using the initial
513
00:29:44,680 --> 00:29:45,950
stress method.
514
00:29:45,950 --> 00:29:49,540
In other words, where tau is
equal to 0 in our governing
515
00:29:49,540 --> 00:29:52,680
finite element equations.
516
00:29:52,680 --> 00:29:58,380
This means that we are setting
up a K matrix initially, and
517
00:29:58,380 --> 00:30:04,960
that K matrix is depicted here
by the slope, the blue slope
518
00:30:04,960 --> 00:30:06,280
that you're seeing.
519
00:30:06,280 --> 00:30:12,380
Now with this slope and this
load applied, we get the
520
00:30:12,380 --> 00:30:14,390
displacement shown down here.
521
00:30:14,390 --> 00:30:19,125
Notice the left superscript
means load step 1.
522
00:30:19,125 --> 00:30:24,200
The right superscript means
solution after iteration 1.
523
00:30:24,200 --> 00:30:28,980
And this is the solution
obtained from this triangle.
524
00:30:28,980 --> 00:30:33,990
Notice there's a blue line up
there, going up there, and
525
00:30:33,990 --> 00:30:39,650
where that blue line crosses the
dashed line, we pick up a
526
00:30:39,650 --> 00:30:41,730
vertical fine, black line.
527
00:30:41,730 --> 00:30:45,890
And that vertical black line
cuts through the displacement
528
00:30:45,890 --> 00:30:48,570
axes at this particular value.
529
00:30:48,570 --> 00:30:51,600
So this is the solution of
our first iteration.
530
00:30:51,600 --> 00:30:55,740
Now, having calculated this
displacement value, we go up
531
00:30:55,740 --> 00:31:00,145
vertically and get to the
red line, the red curve.
532
00:31:00,145 --> 00:31:04,710
And that red curve corresponds,
of course, to the
533
00:31:04,710 --> 00:31:06,560
internal forces.
534
00:31:06,560 --> 00:31:09,340
To the internal force I should
say, corresponding to this
535
00:31:09,340 --> 00:31:11,320
displacement.
536
00:31:11,320 --> 00:31:16,200
There's a red dot right
on that red curve.
537
00:31:16,200 --> 00:31:20,690
And at that red dot now, we've
set up a K matrix again.
538
00:31:20,690 --> 00:31:24,440
Now notice, in the initial
stress method, we keep the K
539
00:31:24,440 --> 00:31:26,000
matrix constant.
540
00:31:26,000 --> 00:31:32,410
This means that the slope here
of this blue line that goes
541
00:31:32,410 --> 00:31:37,830
through this red point is the
same as the slope that we had
542
00:31:37,830 --> 00:31:40,060
earlier right here.
543
00:31:40,060 --> 00:31:44,670
With that slope and the out of
balance load-- and let's look
544
00:31:44,670 --> 00:31:46,330
now very closely here--
545
00:31:46,330 --> 00:31:50,190
that corresponds to the distance
between this red
546
00:31:50,190 --> 00:31:56,330
point and that dashed line with
this out of balance load
547
00:31:56,330 --> 00:32:00,200
and this blue slope, we get an
increment in displacement
548
00:32:00,200 --> 00:32:02,730
shown right here.
549
00:32:02,730 --> 00:32:05,260
And that increment in
displacement added to the
550
00:32:05,260 --> 00:32:11,780
previous displacement gives
us this value down here.
551
00:32:11,780 --> 00:32:16,240
And in fact, this value is
very close to the correct
552
00:32:16,240 --> 00:32:20,380
solution, the exact solution,
which is the solution
553
00:32:20,380 --> 00:32:24,750
corresponding to the dashed
vertical line here.
554
00:32:24,750 --> 00:32:27,180
Notice that the dashed vertical
line is the exact
555
00:32:27,180 --> 00:32:30,160
solution to that dashed
horizontal line.
556
00:32:30,160 --> 00:32:34,920
And our vertical black line,
which we are calculating, is
557
00:32:34,920 --> 00:32:37,970
virtually on that dashed
black line.
558
00:32:37,970 --> 00:32:40,900
So we can accept now,
convergence.
559
00:32:40,900 --> 00:32:42,960
Let's now go into the
next load step.
560
00:32:42,960 --> 00:32:46,330
In the next load step,
we want to, again,
561
00:32:46,330 --> 00:32:47,620
deal with a K matrix.
562
00:32:47,620 --> 00:32:50,640
And in the initial stress
method, we use the same K
563
00:32:50,640 --> 00:32:52,970
matrix throughout
the solution.
564
00:32:52,970 --> 00:32:56,980
Meaning that the slope now,
which we are laying onto the
565
00:32:56,980 --> 00:33:02,860
red curve, that slope being this
blue line, is the same
566
00:33:02,860 --> 00:33:08,470
slope that we used before in
the first two iterations.
567
00:33:08,470 --> 00:33:13,750
With this slope and the out of
balance load corresponding to
568
00:33:13,750 --> 00:33:17,720
the distance between this
horizontal line, the dashed
569
00:33:17,720 --> 00:33:21,780
horizontal line, and that dashed
horizontal line because
570
00:33:21,780 --> 00:33:25,040
we conversed in the
first load step.
571
00:33:25,040 --> 00:33:27,450
With that out of balance load,
we get an incremental
572
00:33:27,450 --> 00:33:34,610
displacement shown from
here to there.
573
00:33:34,610 --> 00:33:37,220
And that incremental
displacement then, added to
574
00:33:37,220 --> 00:33:39,990
the previous displacement,
gives us this value in
575
00:33:39,990 --> 00:33:40,720
displacement.
576
00:33:40,720 --> 00:33:46,390
Notice time t plus delta t or
load step t plus delta t,
577
00:33:46,390 --> 00:33:47,600
which is here.
578
00:33:47,600 --> 00:33:52,620
Load step 2 and iteration
1, end of iteration 1.
579
00:33:52,620 --> 00:33:54,960
This is now our current
displacement.
580
00:33:54,960 --> 00:33:58,670
With this current displacement,
once again,
581
00:33:58,670 --> 00:34:05,760
which comes from this triangle
here, with this current
582
00:34:05,760 --> 00:34:10,610
displacement, we go vertically
up and get to that red point,
583
00:34:10,610 --> 00:34:13,940
which lies on the force
displacement curve--
584
00:34:13,940 --> 00:34:17,500
internal force displacement
curve more accurately.
585
00:34:17,500 --> 00:34:24,510
And now we have this red point,
and we set up, or we
586
00:34:24,510 --> 00:34:27,170
use now, again, a K matrix.
587
00:34:27,170 --> 00:34:29,110
Of course, in the initial
stress method we use the
588
00:34:29,110 --> 00:34:30,110
constant K matrix.
589
00:34:30,110 --> 00:34:32,949
So this slope is the
same as that slope.
590
00:34:32,949 --> 00:34:35,320
And now watch closely.
591
00:34:35,320 --> 00:34:42,960
From this triangle here, which
corresponds to this out of
592
00:34:42,960 --> 00:34:45,130
balance load, which is
the same as that
593
00:34:45,130 --> 00:34:47,820
out of balance load.
594
00:34:47,820 --> 00:34:50,580
Out of balance load obtained
by putting here
595
00:34:50,580 --> 00:34:53,469
a horizontal line.
596
00:34:53,469 --> 00:34:57,100
Notice and the difference
between the dashed line up
597
00:34:57,100 --> 00:35:00,965
there and the horizontal line
given by my pointer is the out
598
00:35:00,965 --> 00:35:02,270
of balance load.
599
00:35:02,270 --> 00:35:05,820
Which of course, is the same as
the distance between this
600
00:35:05,820 --> 00:35:08,430
red point and that
dashed line.
601
00:35:08,430 --> 00:35:12,580
This out of balance load with
this slope gives us another
602
00:35:12,580 --> 00:35:15,190
increment in displacement,
which brings us
603
00:35:15,190 --> 00:35:17,930
to this value here.
604
00:35:17,930 --> 00:35:22,660
We repeat that process and
iteratively, we get closer and
605
00:35:22,660 --> 00:35:26,340
closer to the correct solution,
the exact solution,
606
00:35:26,340 --> 00:35:30,050
which is the dashed vertical
line here for the
607
00:35:30,050 --> 00:35:31,700
displacement.
608
00:35:31,700 --> 00:35:36,040
Because at that dashed vertical
line, we are getting
609
00:35:36,040 --> 00:35:38,660
to that red curve.
610
00:35:38,660 --> 00:35:41,070
Now, this is the initial
stress method.
611
00:35:41,070 --> 00:35:46,030
And once again, the point here
is that we kept the same K
612
00:35:46,030 --> 00:35:48,900
matrix throughout the
solution process.
613
00:35:48,900 --> 00:35:55,710
Surely, if you now were to look
at this iterative scheme
614
00:35:55,710 --> 00:35:59,120
to obtain a faster solution,
you would want to
615
00:35:59,120 --> 00:36:00,770
set up a new K matrix.
616
00:36:00,770 --> 00:36:04,490
In other words, once you have
this displacement value
617
00:36:04,490 --> 00:36:12,670
calculated, you want to update
the slope corresponding to the
618
00:36:12,670 --> 00:36:16,180
slope of the red curve.
619
00:36:16,180 --> 00:36:18,980
And if you do so--
620
00:36:18,980 --> 00:36:22,470
in other words, you update this
blue slope corresponding
621
00:36:22,470 --> 00:36:26,980
to the red dots that you
have obtained in
622
00:36:26,980 --> 00:36:28,640
your iterative scheme.
623
00:36:28,640 --> 00:36:31,300
These red dots, of course,
being different from what
624
00:36:31,300 --> 00:36:32,740
you're seeing here now.
625
00:36:32,740 --> 00:36:37,440
Then you would have the full
Newton-Raphson method.
626
00:36:37,440 --> 00:36:42,750
If you only update the K matrix
or the blue tangent
627
00:36:42,750 --> 00:36:47,230
that we are seeing here whenever
you have converged to
628
00:36:47,230 --> 00:36:50,680
an equilibrium point.
629
00:36:50,680 --> 00:36:54,870
Meaning when you have converged
corresponding for
630
00:36:54,870 --> 00:36:58,710
the first load step in
this particular case.
631
00:36:58,710 --> 00:37:01,120
If you only update the stiffness
matrix once you have
632
00:37:01,120 --> 00:37:03,560
converged to this load--
633
00:37:03,560 --> 00:37:06,690
for this load step, then you
would have the modified
634
00:37:06,690 --> 00:37:09,060
Newton-Raphson method.
635
00:37:09,060 --> 00:37:13,420
I think if you look at this
picture with these
636
00:37:13,420 --> 00:37:16,700
possibilities, you will realize
that surely, the full
637
00:37:16,700 --> 00:37:21,080
Newton-Raphson method with
updating the slope after each
638
00:37:21,080 --> 00:37:23,970
iteration will converge
fastest.
639
00:37:23,970 --> 00:37:26,690
And the modified Newton-Raphson
method will
640
00:37:26,690 --> 00:37:29,150
converge a little slower than
the full Newton-Raphson
641
00:37:29,150 --> 00:37:34,260
method, but still faster than
the initial stress method.
642
00:37:34,260 --> 00:37:38,210
Well, this then brings us to the
end of the discussion of
643
00:37:38,210 --> 00:37:39,080
this example.
644
00:37:39,080 --> 00:37:42,800
I like to now introduce you to
another scheme that can be
645
00:37:42,800 --> 00:37:46,430
very important when we solve
nonlinear equations, finite
646
00:37:46,430 --> 00:37:47,770
element equations.
647
00:37:47,770 --> 00:37:52,290
And that scheme is the scheme
of line searches.
648
00:37:52,290 --> 00:37:55,480
The basic equations that we are
looking at that we want to
649
00:37:55,480 --> 00:37:57,380
solve are once more
depicted here.
650
00:37:57,380 --> 00:37:59,990
Tangent stiffness matrix,
incremental displacement
651
00:37:59,990 --> 00:38:03,350
vector, externally applied
loads, nodal point forces
652
00:38:03,350 --> 00:38:05,790
corresponding to the internal
element stresses at time t
653
00:38:05,790 --> 00:38:10,140
plus delta t, and at the end
of iteration i minus 1,
654
00:38:10,140 --> 00:38:12,640
beginning of iteration
i, of course.
655
00:38:12,640 --> 00:38:15,070
Notice I don't have an i on
here, I just want to denote
656
00:38:15,070 --> 00:38:18,220
this for the moment as an
incremental nodal point
657
00:38:18,220 --> 00:38:22,200
displacement vector, which
carries, however, a bar.
658
00:38:22,200 --> 00:38:29,560
Let us consider forming Fi using
this equation where beta
659
00:38:29,560 --> 00:38:31,800
is an unknown.
660
00:38:31,800 --> 00:38:35,460
We want to choose beta
now judiciously.
661
00:38:35,460 --> 00:38:40,380
We want to choose beta such
that, in some sense, R minus
662
00:38:40,380 --> 00:38:44,580
Fi is small.
663
00:38:44,580 --> 00:38:48,515
And we have to, of course, look
at the mechanism that we
664
00:38:48,515 --> 00:38:53,560
use to make R minus F
small in some sense.
665
00:38:53,560 --> 00:38:59,040
Well, as a side note, we
recognize that when this
666
00:38:59,040 --> 00:39:06,180
equation here is 0 for all
possible U's, then clearly
667
00:39:06,180 --> 00:39:08,640
this must be 0.
668
00:39:08,640 --> 00:39:14,380
In other words, if we were to
say to ourselves, let's make
669
00:39:14,380 --> 00:39:17,410
this equation 0 for
all possible U's.
670
00:39:17,410 --> 00:39:22,770
Then really, we would have
solved this equal to 0.
671
00:39:22,770 --> 00:39:25,640
The reason, of course, for it
is that we could choose here
672
00:39:25,640 --> 00:39:31,310
for U a vector, which carries
0's everywhere except a 1 in
673
00:39:31,310 --> 00:39:32,610
one location.
674
00:39:32,610 --> 00:39:36,880
And if that is a particular
row, then corresponding to
675
00:39:36,880 --> 00:39:44,170
that row, surely this vector
here, this R minus F
676
00:39:44,170 --> 00:39:46,900
corresponding to that row
must be equal to 0.
677
00:39:46,900 --> 00:39:50,390
So if we were simply to use here
the identity matrix, or
678
00:39:50,390 --> 00:39:54,660
if we were to construct the
identity matrix by allowing n
679
00:39:54,660 --> 00:39:57,730
U's that are linearly
independent, of course,
680
00:39:57,730 --> 00:40:03,480
because they all make up such
vector but with a 1 in a
681
00:40:03,480 --> 00:40:04,500
different location.
682
00:40:04,500 --> 00:40:08,950
As a matter of fact, we want to
use U1 with a 1 here, 0's
683
00:40:08,950 --> 00:40:09,850
everywhere.
684
00:40:09,850 --> 00:40:14,650
U2, 1 here, 0 here, everywhere
else 0, and so on.
685
00:40:14,650 --> 00:40:16,870
To construct an identity
matrix here.
686
00:40:16,870 --> 00:40:20,500
Then with that identity matrix,
clearly R minus F
687
00:40:20,500 --> 00:40:23,090
would be 0 everywhere.
688
00:40:23,090 --> 00:40:27,250
Of course, that is not a
viable scheme really.
689
00:40:27,250 --> 00:40:34,170
But what we want to do is to
choose a particular vector
690
00:40:34,170 --> 00:40:40,430
here so as to make R minus
F still 0 in some sense.
691
00:40:40,430 --> 00:40:44,760
And the vector that is quite
effectively chosen is the
692
00:40:44,760 --> 00:40:49,460
vector that we obtain from the
solution of K delta U bar
693
00:40:49,460 --> 00:40:51,430
equals R minus F.
694
00:40:51,430 --> 00:40:58,070
If we use that vector here,
basically projecting the R,
695
00:40:58,070 --> 00:41:01,650
capital R, minus F
onto that vector.
696
00:41:01,650 --> 00:41:06,770
And we look for this to be 0.
697
00:41:06,770 --> 00:41:12,350
Then we have so to say, searched
in the direction of
698
00:41:12,350 --> 00:41:20,210
delta U bar and made delta U
bar transposed times this
699
00:41:20,210 --> 00:41:22,360
vector here equal to 0.
700
00:41:22,360 --> 00:41:23,970
How does this scheme work?
701
00:41:23,970 --> 00:41:29,240
Well, we add beta times delta U
bar to the value that we had
702
00:41:29,240 --> 00:41:31,600
already from the previous
iteration.
703
00:41:31,600 --> 00:41:39,340
And we search along beta such
that this top here is much
704
00:41:39,340 --> 00:41:41,360
smaller than the bottom.
705
00:41:41,360 --> 00:41:45,650
In other words, STOL shall be
a convergence tolerance.
706
00:41:45,650 --> 00:41:47,940
Notice the way it works.
707
00:41:47,940 --> 00:41:49,630
We choose a value of beta.
708
00:41:49,630 --> 00:41:53,550
With beta known, we add this
quantity to this one.
709
00:41:53,550 --> 00:41:55,470
We get a value here.
710
00:41:55,470 --> 00:42:00,470
We take this value, calculate
Fi corresponding to that
711
00:42:00,470 --> 00:42:02,510
value, see whether
the convergence
712
00:42:02,510 --> 00:42:03,520
tolerance is satisfied.
713
00:42:03,520 --> 00:42:06,090
If not, we have to choose
a new beta.
714
00:42:06,090 --> 00:42:13,590
And like that, we choose new
betas until we have satisfied
715
00:42:13,590 --> 00:42:15,360
this convergence tolerance.
716
00:42:15,360 --> 00:42:17,320
That's what we call
line search.
717
00:42:17,320 --> 00:42:21,700
We are searching along the
line given by delta U bar
718
00:42:21,700 --> 00:42:25,110
until this is satisfied.
719
00:42:25,110 --> 00:42:29,210
It's a very important scheme,
and this scheme can be
720
00:42:29,210 --> 00:42:32,590
combined with modified Newton
iteration, with the initial
721
00:42:32,590 --> 00:42:35,590
stress method, and with
full Newton iteration.
722
00:42:35,590 --> 00:42:40,680
And it adds a lot to the
stability of the solution, of
723
00:42:40,680 --> 00:42:43,660
the overall solution or the
nonlinear equations.
724
00:42:43,660 --> 00:42:47,470
We will find, or discuss
applications just now.
725
00:42:47,470 --> 00:42:52,890
But at this point, we need to
stop for a minute because we
726
00:42:52,890 --> 00:42:54,390
need to go onto a new tape.
727
00:42:58,160 --> 00:43:00,970
Finally, I would like to discuss
with you a method that
728
00:43:00,970 --> 00:43:03,540
we also find to be very
effective in nonlinear finite
729
00:43:03,540 --> 00:43:05,050
element computations.
730
00:43:05,050 --> 00:43:07,590
Namely, the BFGS method.
731
00:43:07,590 --> 00:43:08,835
BFGS stands for Broyden-Fletcher
732
00:43:08,835 --> 00:43:10,085
-Goldfarb-Shanno.
733
00:43:12,200 --> 00:43:15,930
And in this method, we
proceed as follows.
734
00:43:15,930 --> 00:43:20,850
We calculate or define
a delta i as shown
735
00:43:20,850 --> 00:43:23,070
here in this equation.
736
00:43:23,070 --> 00:43:28,560
And we define a gamma i as
shown in this equation.
737
00:43:28,560 --> 00:43:32,200
We want to calculate the
coefficient matrix such that
738
00:43:32,200 --> 00:43:34,110
this equation is
here satisfied.
739
00:43:34,110 --> 00:43:37,720
In other words, given delta i
and gamma i, we want to use
740
00:43:37,720 --> 00:43:41,040
the coefficient matrix that
satisfies this equation.
741
00:43:41,040 --> 00:43:45,680
And that's basically what's
being done in the BFGS scheme.
742
00:43:45,680 --> 00:43:50,610
Pictorially then, we can show
for a single degree of freedom
743
00:43:50,610 --> 00:43:55,260
system, say the F plotted
along this axis.
744
00:43:55,260 --> 00:43:58,730
The U displacements plotted
along this axis.
745
00:43:58,730 --> 00:44:01,770
Notice here in green
we show delta i.
746
00:44:01,770 --> 00:44:04,520
And we show here in
green gamma i.
747
00:44:04,520 --> 00:44:08,890
And notice that the matrix
that we want to use is
748
00:44:08,890 --> 00:44:11,760
indicated by this blue slope.
749
00:44:11,760 --> 00:44:16,560
Which is, in other words, given
by delta i and gamma i.
750
00:44:16,560 --> 00:44:21,300
Well, the BFGS method is an
iterative algorithm, which
751
00:44:21,300 --> 00:44:26,090
produces successive
approximations to efficient
752
00:44:26,090 --> 00:44:27,980
stiffness matrix.
753
00:44:27,980 --> 00:44:30,410
Actually, we actually work
with the inverse of that
754
00:44:30,410 --> 00:44:34,450
stiffness matrix as I will
elaborate upon just now.
755
00:44:34,450 --> 00:44:36,770
It's really a compromise
between the full Newton
756
00:44:36,770 --> 00:44:38,600
iteration method and
the modified
757
00:44:38,600 --> 00:44:40,590
Newton iteration method.
758
00:44:40,590 --> 00:44:43,350
It shows very excellent
convergence characteristics in
759
00:44:43,350 --> 00:44:45,950
the solution of many
types of problems.
760
00:44:45,950 --> 00:44:48,680
Let's look at the individual
steps that we use when we
761
00:44:48,680 --> 00:44:50,900
apply the BFGS method.
762
00:44:50,900 --> 00:44:54,530
We, first of all, calculate
the direction of the
763
00:44:54,530 --> 00:44:56,160
displacement increment.
764
00:44:56,160 --> 00:45:03,090
Here we have delta U bar i is
equal to K inverse, the
765
00:45:03,090 --> 00:45:06,120
inverse of the K matrix,
corresponding to
766
00:45:06,120 --> 00:45:08,050
iteration i minus 1.
767
00:45:08,050 --> 00:45:12,150
And here we have the out
of balance load vector.
768
00:45:12,150 --> 00:45:15,170
Notice once again we do not
calculate actually an inverse
769
00:45:15,170 --> 00:45:16,080
of a matrix.
770
00:45:16,080 --> 00:45:20,730
We calculate the LDL transpose
factorization of a matrix, and
771
00:45:20,730 --> 00:45:25,690
then we, as I will just now
show, update that inverse in
772
00:45:25,690 --> 00:45:31,445
the BFGS iteration by specific
matrices of rank 2.
773
00:45:31,445 --> 00:45:33,610
We get to that just now.
774
00:45:33,610 --> 00:45:39,360
Well, having calculated this
delta U bar i, we use that
775
00:45:39,360 --> 00:45:44,850
delta U bar i in the line search
to obtain a better
776
00:45:44,850 --> 00:45:49,650
displacement vector, a
displacement vector
777
00:45:49,650 --> 00:45:51,550
corresponding to iteration i.
778
00:45:51,550 --> 00:45:56,450
A better value then just adding
the delta U bar i to it
779
00:45:56,450 --> 00:45:57,910
by adjusting beta.
780
00:45:57,910 --> 00:46:02,960
In other words, we choose beta
in such way as to satisfy this
781
00:46:02,960 --> 00:46:04,250
equation here.
782
00:46:04,250 --> 00:46:06,060
This is, of course, the equation
that we are now
783
00:46:06,060 --> 00:46:11,680
familiar with as being equation
that shows that
784
00:46:11,680 --> 00:46:16,300
convergence has been reached
in a line search.
785
00:46:16,300 --> 00:46:22,650
Having now calculated the t
plus delta t Ui and the
786
00:46:22,650 --> 00:46:28,390
corresponding t plus delta t Fi,
we can calculate the delta
787
00:46:28,390 --> 00:46:30,660
i and gamma i.
788
00:46:30,660 --> 00:46:33,930
And therefore, apply
the BFGS method.
789
00:46:33,930 --> 00:46:36,630
And that's now step three,
calculation of
790
00:46:36,630 --> 00:46:38,760
the new secant matrix.
791
00:46:38,760 --> 00:46:46,420
K inverse of iteration i is
equal to Ai transposed k
792
00:46:46,420 --> 00:46:50,100
inverse of iteration
i minus 1 times Ai.
793
00:46:50,100 --> 00:46:55,393
This A matrix is obtained
as shown here.
794
00:46:55,393 --> 00:46:57,910
v wi transposed.
795
00:46:57,910 --> 00:47:04,390
The vi and wi's are given as a
function of the delta i, gamma
796
00:47:04,390 --> 00:47:10,450
i, and K i minus 1 value.
797
00:47:10,450 --> 00:47:15,120
You should look at the textbook
to see the details of
798
00:47:15,120 --> 00:47:18,590
actually evaluating
the vi and wi.
799
00:47:18,590 --> 00:47:23,430
But once you have vi, wi, you
can calculate A. And clearly,
800
00:47:23,430 --> 00:47:29,100
by this equation then, you get
a new stiffness matrix.
801
00:47:29,100 --> 00:47:33,670
It is important to note that in
this iterative scheme, only
802
00:47:33,670 --> 00:47:38,550
vector products are needed
to obtain vi and wi.
803
00:47:38,550 --> 00:47:41,680
And it's also then important
to note that only vector
804
00:47:41,680 --> 00:47:45,520
products are used to calculate
delta U bar i.
805
00:47:45,520 --> 00:47:50,040
If this would not be the case,
particularly here.
806
00:47:50,040 --> 00:47:51,540
Then of course, the
iterative scheme
807
00:47:51,540 --> 00:47:53,470
would be very expensive.
808
00:47:53,470 --> 00:47:59,960
Let me show you a bit of the
details why we only need to
809
00:47:59,960 --> 00:48:03,230
use vector products to
get the solution.
810
00:48:03,230 --> 00:48:09,230
Here, we have one typical
step of the iteration.
811
00:48:09,230 --> 00:48:13,730
Notice on the left-hand side,
we have delta U bar i.
812
00:48:13,730 --> 00:48:15,890
The value of displacement
increment
813
00:48:15,890 --> 00:48:17,880
that we want to calculate.
814
00:48:17,880 --> 00:48:22,460
On the right-hand side, we
have i plus wi minus 1 vi
815
00:48:22,460 --> 00:48:23,990
minus 1 transposed.
816
00:48:23,990 --> 00:48:28,890
In other words, the Ai
minus 1 transposed.
817
00:48:28,890 --> 00:48:33,860
And we would have such matrices
going on here until
818
00:48:33,860 --> 00:48:37,670
we are coming to the A1
transposed matrix.
819
00:48:37,670 --> 00:48:40,830
Then we have the inverse
of a stiffness matrix.
820
00:48:40,830 --> 00:48:43,860
Here, corresponding
to time tau.
821
00:48:43,860 --> 00:48:50,590
And then we have the Ai matrices
as shown here.
822
00:48:50,590 --> 00:48:57,590
Now, notice that this total
amount then is multiplied by
823
00:48:57,590 --> 00:49:02,690
the R minus F i minus 1, the
out of balance load vector.
824
00:49:02,690 --> 00:49:06,370
Let's now go through the
computations that are required
825
00:49:06,370 --> 00:49:10,500
to get delta U bar i.
826
00:49:10,500 --> 00:49:14,850
Here we have delta R. This delta
R, the out of balance
827
00:49:14,850 --> 00:49:17,030
load vector, which I
just call delta R,
828
00:49:17,030 --> 00:49:19,110
multiplies this matrix.
829
00:49:19,110 --> 00:49:22,010
But if you look at this
multiplication, we find that
830
00:49:22,010 --> 00:49:29,430
this delta R times w i minus 1
transposed is simply a number.
831
00:49:29,430 --> 00:49:31,460
So it's one vector
multiplication
832
00:49:31,460 --> 00:49:33,020
that gives us a number.
833
00:49:33,020 --> 00:49:39,110
This number multiplied by
v i minus 1 is just the
834
00:49:39,110 --> 00:49:41,100
multiplication of a number
times a vector.
835
00:49:41,100 --> 00:49:42,930
It's very simple.
836
00:49:42,930 --> 00:49:48,230
And notice, of course, this
delta R i minus 1 is also
837
00:49:48,230 --> 00:49:49,850
multiplying the identity
matrix.
838
00:49:49,850 --> 00:49:53,910
So in this step here, we have
really nothing else then a
839
00:49:53,910 --> 00:49:55,080
vector multiplication.
840
00:49:55,080 --> 00:49:58,040
The only vector multiplication
that, in fact, is there, is
841
00:49:58,040 --> 00:50:02,030
this value, this vector here,
times that vector there,
842
00:50:02,030 --> 00:50:03,510
transposed.
843
00:50:03,510 --> 00:50:04,610
That gives us a number.
844
00:50:04,610 --> 00:50:07,120
This number times a vector
[UNINTELLIGIBLE]
845
00:50:07,120 --> 00:50:08,040
operation.
846
00:50:08,040 --> 00:50:11,780
And of course, this one here
times the i matrix is just
847
00:50:11,780 --> 00:50:13,500
this vector by itself.
848
00:50:13,500 --> 00:50:16,690
So in this step, what
we obtain is
849
00:50:16,690 --> 00:50:19,950
simply another vector.
850
00:50:19,950 --> 00:50:24,240
We take this vector and
repeat the process.
851
00:50:24,240 --> 00:50:28,500
And since the repetition is
just what we discussed
852
00:50:28,500 --> 00:50:34,470
already, we find that we only
need vector multiplications to
853
00:50:34,470 --> 00:50:38,790
roll up all of these
multiplications up to here.
854
00:50:38,790 --> 00:50:43,045
Then we have a vector multiplied
by K inverse.
855
00:50:43,045 --> 00:50:47,380
Well, that involves only vector
multiplications again,
856
00:50:47,380 --> 00:50:49,810
because we have already
the LDL transposed
857
00:50:49,810 --> 00:50:51,460
factorization of tau k.
858
00:50:54,140 --> 00:50:57,080
These vector multiplications
result again, in a vector.
859
00:50:57,080 --> 00:51:00,880
And that vector multiplied by
these matrices here in the
860
00:51:00,880 --> 00:51:04,330
same way as we have proceeded
here, means again, only vector
861
00:51:04,330 --> 00:51:05,930
multiplications.
862
00:51:05,930 --> 00:51:10,440
So this very important step
here only requires vector
863
00:51:10,440 --> 00:51:11,140
multiplications.
864
00:51:11,140 --> 00:51:16,520
And that makes this whole
process really efficient.
865
00:51:16,520 --> 00:51:20,740
In summary then, we have the
following solution procedures
866
00:51:20,740 --> 00:51:23,130
that we feel are
very effective.
867
00:51:23,130 --> 00:51:27,020
The mortified Newton-Raphson
iteration with line searches.
868
00:51:27,020 --> 00:51:29,540
Here we have the basic modified
869
00:51:29,540 --> 00:51:31,190
Newton iteration equation.
870
00:51:31,190 --> 00:51:34,110
And here we have the line
search equation.
871
00:51:34,110 --> 00:51:38,930
We, of course, after each such
solution, evaluate an
872
00:51:38,930 --> 00:51:39,810
efficient beta.
873
00:51:39,810 --> 00:51:43,600
A beta that makes U--
874
00:51:43,600 --> 00:51:48,660
t plus delta U i a good vector
when measured on the line
875
00:51:48,660 --> 00:51:52,010
search criteria the way
I discussed it.
876
00:51:52,010 --> 00:51:56,210
And these two equations then,
summarize really the modified
877
00:51:56,210 --> 00:51:58,475
Newton iteration with
the line search.
878
00:51:58,475 --> 00:52:01,750
The BFGS method is another
effective scheme.
879
00:52:01,750 --> 00:52:05,210
We use the BFGS method always
with line searches.
880
00:52:05,210 --> 00:52:07,950
And then, as a third option,
the full Newton-Raphson
881
00:52:07,950 --> 00:52:10,740
iteration with or without
line searches.
882
00:52:10,740 --> 00:52:13,220
The full Newton-Raphson
iteration with line searches
883
00:52:13,220 --> 00:52:14,890
is, of course, most powerful.
884
00:52:14,890 --> 00:52:18,210
But it is also most expensive
per iteration.
885
00:52:18,210 --> 00:52:21,020
We should point out that these
techniques that I discussed so
886
00:52:21,020 --> 00:52:24,760
far cannot be applied for
post-buckling analysis.
887
00:52:24,760 --> 00:52:29,160
We will consider the solution of
the response after buckling
888
00:52:29,160 --> 00:52:32,520
has been occurring, after the
item load level has been
889
00:52:32,520 --> 00:52:36,020
reached in a forthcoming
lecture.
890
00:52:36,020 --> 00:52:41,610
The next view graphs summarize
these schemes: the modified
891
00:52:41,610 --> 00:52:44,450
Newton-Raphson, the BFGS
method, and the full
892
00:52:44,450 --> 00:52:46,840
Newton-Raphson method.
893
00:52:46,840 --> 00:52:50,030
And I'd like to just very
quickly go through the whole
894
00:52:50,030 --> 00:52:54,070
solution process for each
of these methods.
895
00:52:54,070 --> 00:52:57,700
In the modified Newton iteration
technique, we
896
00:52:57,700 --> 00:53:03,580
initialize our displacements
corresponding to the beginning
897
00:53:03,580 --> 00:53:09,610
of the first iteration, the end
of the 0 iteration as tU.
898
00:53:09,610 --> 00:53:12,920
And similarly, we initialize our
force vector corresponding
899
00:53:12,920 --> 00:53:14,470
to the internal element
stresses.
900
00:53:14,470 --> 00:53:19,010
We set i is equal to 1, the
iteration counter equal to 1.
901
00:53:19,010 --> 00:53:22,610
We calculate a K matrix at the
beginning of the load step,
902
00:53:22,610 --> 00:53:25,510
and we keep that K
matrix constant.
903
00:53:25,510 --> 00:53:29,450
We then go into this step here
where we calculate the out of
904
00:53:29,450 --> 00:53:30,700
balance loads.
905
00:53:33,210 --> 00:53:36,790
Calculate a new displacement
vector.
906
00:53:36,790 --> 00:53:41,240
And in this box here, we
measure whether we have
907
00:53:41,240 --> 00:53:42,160
converged already.
908
00:53:42,160 --> 00:53:45,390
Of course, this convergence
would only be measured after
909
00:53:45,390 --> 00:53:49,560
we have gone say, at
least twice through
910
00:53:49,560 --> 00:53:51,070
the iteration process.
911
00:53:51,070 --> 00:53:54,650
Or rather, we have calculated
this increment in displacement
912
00:53:54,650 --> 00:53:56,020
vector twice.
913
00:53:56,020 --> 00:53:59,560
Because at the beginning
corresponding to i equal to 1,
914
00:53:59,560 --> 00:54:03,400
we have now this fairly large
out of balance load.
915
00:54:03,400 --> 00:54:05,870
We have just incremented
the load vector.
916
00:54:05,870 --> 00:54:10,880
Therefore, this right inside
should be non-zero.
917
00:54:10,880 --> 00:54:13,580
And we would get an increment
in displacements.
918
00:54:13,580 --> 00:54:17,500
And to measure how large that
increment in displacement is,
919
00:54:17,500 --> 00:54:20,350
we really want to go through
this whole cycle at
920
00:54:20,350 --> 00:54:22,170
least one more time.
921
00:54:22,170 --> 00:54:25,610
We then, having calculated this
incremental displacement
922
00:54:25,610 --> 00:54:30,620
vector with the bar on there,
we perform a line search the
923
00:54:30,620 --> 00:54:34,230
way we discussed it to update
our displacements
924
00:54:34,230 --> 00:54:36,720
corresponding to the
iteration i.
925
00:54:39,960 --> 00:54:42,300
Of course, we only have come so
far once through it, so i
926
00:54:42,300 --> 00:54:43,520
is equal to 1 still.
927
00:54:43,520 --> 00:54:48,360
And this gives us the accepted
value of displacements
928
00:54:48,360 --> 00:54:50,936
corresponding to the
first iteration.
929
00:54:50,936 --> 00:54:55,560
We now increment our iteration
counter, and go in to the
930
00:54:55,560 --> 00:54:57,370
second iteration.
931
00:54:57,370 --> 00:55:01,060
Notice now in the second
iteration, we have R minus t
932
00:55:01,060 --> 00:55:03,670
plus delta t F1 here.
933
00:55:03,670 --> 00:55:06,250
We calculate delta U bar 2.
934
00:55:06,250 --> 00:55:10,040
And now we would certainly
measure convergence.
935
00:55:10,040 --> 00:55:13,360
And I will just now talk about
how we measure convergence.
936
00:55:13,360 --> 00:55:17,240
If we have not converged yet,
we keep on cycling through
937
00:55:17,240 --> 00:55:22,080
here until we do actually
get convergence.
938
00:55:22,080 --> 00:55:25,840
The distinguishing feature of
the modified Newton iteration
939
00:55:25,840 --> 00:55:29,642
is that we are having
a constant K matrix.
940
00:55:29,642 --> 00:55:32,900
That we do not update
the K matrix is this
941
00:55:32,900 --> 00:55:34,250
iterative loop here.
942
00:55:34,250 --> 00:55:36,280
It remains constant.
943
00:55:36,280 --> 00:55:40,570
Remember in the initial stress
technique, initial stress
944
00:55:40,570 --> 00:55:46,020
method, we would not even update
the K matrix ever after
945
00:55:46,020 --> 00:55:49,980
the first time it has
been calculated.
946
00:55:49,980 --> 00:55:52,990
But in the modified Newton
iteration, we update it at the
947
00:55:52,990 --> 00:55:55,800
beginning of each load step.
948
00:55:55,800 --> 00:56:00,140
In the BFGS method, we
proceed as follows.
949
00:56:00,140 --> 00:56:03,420
We initial, again, our
displacements and our nodal
950
00:56:03,420 --> 00:56:05,350
point force vector corresponding
to the internal
951
00:56:05,350 --> 00:56:06,570
element stresses.
952
00:56:06,570 --> 00:56:08,470
We calculate tK.
953
00:56:08,470 --> 00:56:11,630
We calculate our incremental
displacement vector with the
954
00:56:11,630 --> 00:56:13,120
bar on top.
955
00:56:13,120 --> 00:56:18,670
We measure convergence if i
is greater or equal to.
956
00:56:18,670 --> 00:56:20,700
There's no point as I just
mentioned, measuring
957
00:56:20,700 --> 00:56:24,070
convergence when you go first
time through here.
958
00:56:24,070 --> 00:56:29,020
And we perform then the line
search having calculated our
959
00:56:29,020 --> 00:56:32,440
actual incremental in nodal
point displacement that we
960
00:56:32,440 --> 00:56:35,410
want to add to the previous
nodal point displacements to
961
00:56:35,410 --> 00:56:39,150
get to this displacement
vector.
962
00:56:39,150 --> 00:56:42,040
We, of course, also have
this force vector.
963
00:56:42,040 --> 00:56:46,960
We now can update our inverse
matrix, our inverse secret
964
00:56:46,960 --> 00:56:50,740
matrix the way we
discussed it.
965
00:56:50,740 --> 00:56:53,120
And we increment our counter.
966
00:56:53,120 --> 00:56:58,110
And we keep on cycling through
here until we find in this box
967
00:56:58,110 --> 00:57:01,300
that we have converged.
968
00:57:01,300 --> 00:57:06,870
In the full Newton iteration
with line searches, we once
969
00:57:06,870 --> 00:57:10,460
again, initialize our
displacement vector, our force
970
00:57:10,460 --> 00:57:13,490
vector, our iteration counter.
971
00:57:13,490 --> 00:57:17,710
We now calculate a new K matrix
corresponding to the
972
00:57:17,710 --> 00:57:20,720
last conditions on displacements
and internal
973
00:57:20,720 --> 00:57:22,320
forces, of course.
974
00:57:22,320 --> 00:57:25,140
We enter here to calculate our
incremental displacement
975
00:57:25,140 --> 00:57:29,460
vector with the K matrix that
was just evaluated.
976
00:57:32,110 --> 00:57:36,740
And we measure convergence
once again only when i is
977
00:57:36,740 --> 00:57:38,460
greater equal to.
978
00:57:38,460 --> 00:57:45,450
We perform a line search to
determine beta in here.
979
00:57:45,450 --> 00:57:49,010
We therefore, update our
displacement vector to the
980
00:57:49,010 --> 00:57:51,570
displacements that now
correspond to the end of
981
00:57:51,570 --> 00:57:53,230
iteration i.
982
00:57:53,230 --> 00:57:55,250
We update our iteration
counter.
983
00:57:55,250 --> 00:57:57,720
Or rather, increase our
iteration counter.
984
00:57:57,720 --> 00:58:02,580
And we are calculating a new
stiffness matrix corresponding
985
00:58:02,580 --> 00:58:05,200
to the last calculated
displacement.
986
00:58:05,200 --> 00:58:08,180
This is, of course, a
distinguishing feature of the
987
00:58:08,180 --> 00:58:10,660
full Newton-Raphson method that
we are calculating a new
988
00:58:10,660 --> 00:58:13,700
K matrix at the beginning
of each iteration.
989
00:58:16,210 --> 00:58:21,360
Well, these view graphs then,
summarize the full process of
990
00:58:21,360 --> 00:58:24,770
iteration for the modified
Newton, the BFGS, and the full
991
00:58:24,770 --> 00:58:25,440
Newton method.
992
00:58:25,440 --> 00:58:28,660
Notice that I have included the
line search as the option.
993
00:58:28,660 --> 00:58:31,480
If you don't have line
searching, of course, it would
994
00:58:31,480 --> 00:58:37,290
simply mean that you set beta
equal to 1, and skip that
995
00:58:37,290 --> 00:58:39,680
particular step of
line searching.
996
00:58:39,680 --> 00:58:43,530
But in many cases, line
searching is quite important.
997
00:58:43,530 --> 00:58:47,440
And certainly, can
be necessary in
998
00:58:47,440 --> 00:58:49,790
some types of problems.
999
00:58:49,790 --> 00:58:52,690
Let us now look at what happens
in this box here.
1000
00:58:52,690 --> 00:58:55,260
How do we measure convergence?
1001
00:58:55,260 --> 00:58:58,850
In each of these iterative
schemes, we had such a box.
1002
00:58:58,850 --> 00:59:02,730
And I like to look now more
closely at how we measure
1003
00:59:02,730 --> 00:59:05,040
convergence.
1004
00:59:05,040 --> 00:59:09,090
The measure of convergence
should, of course, tell us,
1005
00:59:09,090 --> 00:59:11,740
how well do we satisfy
equilibrium?
1006
00:59:11,740 --> 00:59:16,070
And there are basically three
items that can be used.
1007
00:59:16,070 --> 00:59:20,870
Namely energy, force, or moment,
displacement, and
1008
00:59:20,870 --> 00:59:23,110
rotation, of course,
correspondingly as well.
1009
00:59:23,110 --> 00:59:25,940
But these are the items that one
can work with to measure
1010
00:59:25,940 --> 00:59:27,280
convergence.
1011
00:59:27,280 --> 00:59:31,510
If you use an energy convergence
criteria, and we
1012
00:59:31,510 --> 00:59:33,800
are very fond of this one.
1013
00:59:33,800 --> 00:59:40,970
You might want to use this
equation here delta U bar i
1014
00:59:40,970 --> 00:59:46,800
transposed times the out of
balance load divided by delta
1015
00:59:46,800 --> 00:59:50,870
U bar 1 transposed times t
plus delta t R minus tF.
1016
00:59:50,870 --> 00:59:54,500
Now notice, this is the
increment in nodal point
1017
00:59:54,500 --> 00:59:58,020
forces, or the out of
balance load, at the
1018
00:59:58,020 --> 01:00:00,990
beginning of the load step.
1019
01:00:00,990 --> 01:00:05,850
So this here is an energy
corresponding to the beginning
1020
01:00:05,850 --> 01:00:07,090
of the load step.
1021
01:00:07,090 --> 01:00:11,460
And we are comparing here
a similar quantity, but
1022
01:00:11,460 --> 01:00:13,760
corresponding to the
current load.
1023
01:00:13,760 --> 01:00:16,140
Corresponding to the
current iteration.
1024
01:00:16,140 --> 01:00:21,190
And if this quotient here is
very small measured on that
1025
01:00:21,190 --> 01:00:25,250
convergence tolerance, then
we say we have converged.
1026
01:00:25,250 --> 01:00:28,910
An important point is that in
this convergence measure, we
1027
01:00:28,910 --> 01:00:33,400
enter with displacements and
with out of balance loads.
1028
01:00:33,400 --> 01:00:36,450
Both enter in this convergence
measure.
1029
01:00:36,450 --> 01:00:41,270
And please realize also that
if the loads don't increase
1030
01:00:41,270 --> 01:00:46,960
from time t to t plus delta t,
meaning that this here is
1031
01:00:46,960 --> 01:00:49,330
equal to 0--
1032
01:00:49,330 --> 01:00:50,410
0 vector.
1033
01:00:50,410 --> 01:00:52,860
Then, of course, we
get 0 down here.
1034
01:00:52,860 --> 01:00:55,020
And we could not apply this
convergence criteria
1035
01:00:55,020 --> 01:00:56,170
indirectly.
1036
01:00:56,170 --> 01:00:59,900
Therefore, we would do is
instead of 0 here, a
1037
01:00:59,900 --> 01:01:03,890
reasonable small number in
the computer program.
1038
01:01:03,890 --> 01:01:06,760
Also note that we apply
this test here
1039
01:01:06,760 --> 01:01:09,300
prior to line searching.
1040
01:01:09,300 --> 01:01:13,590
In line searching, after line
searching of course, this top
1041
01:01:13,590 --> 01:01:15,580
part here is small.
1042
01:01:15,580 --> 01:01:20,070
And we, therefore, do not want
to apply this scheme here or
1043
01:01:20,070 --> 01:01:22,290
this convergence measure
after line searching.
1044
01:01:22,290 --> 01:01:25,530
We should do it before
line searching.
1045
01:01:25,530 --> 01:01:31,040
To measure convergence on
forces, we can use this
1046
01:01:31,040 --> 01:01:32,420
equation here.
1047
01:01:32,420 --> 01:01:35,620
Notice here we have the
out of balance loads.
1048
01:01:35,620 --> 01:01:39,640
We take the Euclidean norm on
the out of balance loads, and
1049
01:01:39,640 --> 01:01:43,710
we compare that with
RNORM, input by the
1050
01:01:43,710 --> 01:01:45,970
user, by the analyst.
1051
01:01:45,970 --> 01:01:50,560
It is interesting, important to
note that we introduce here
1052
01:01:50,560 --> 01:01:53,860
only the components
corresponding to the
1053
01:01:53,860 --> 01:01:57,220
translational degrees
of freedom.
1054
01:01:57,220 --> 01:02:01,540
And here, of course, this
RNORM is also a force
1055
01:02:01,540 --> 01:02:03,340
corresponding, so to say, to a
1056
01:02:03,340 --> 01:02:04,840
translational degree of freedom.
1057
01:02:04,840 --> 01:02:09,430
If we have also rotational
degrees of freedom, then we
1058
01:02:09,430 --> 01:02:14,500
would take instead
of RNORM, RMNORM.
1059
01:02:14,500 --> 01:02:18,760
And we would extract from these
vectors, the components
1060
01:02:18,760 --> 01:02:20,130
corresponding to the rotational
1061
01:02:20,130 --> 01:02:21,700
degrees of freedom only.
1062
01:02:21,700 --> 01:02:24,530
The point is that we want to
have the same units up
1063
01:02:24,530 --> 01:02:27,450
here as down here.
1064
01:02:27,450 --> 01:02:31,790
Forces, forces, moments,
moments.
1065
01:02:31,790 --> 01:02:36,440
And we don't want to mix
components from these
1066
01:02:36,440 --> 01:02:38,280
quantities.
1067
01:02:38,280 --> 01:02:43,160
Of course, RMNORM and RNORM must
be chosen by the analyst.
1068
01:02:43,160 --> 01:02:48,190
The are input to the analysis as
defined in the input to the
1069
01:02:48,190 --> 01:02:49,480
computer program.
1070
01:02:49,480 --> 01:02:53,030
Typically, RTOL is 0.01.
1071
01:02:53,030 --> 01:02:58,540
RNORM is simply the maximum
of the NORM of the
1072
01:02:58,540 --> 01:03:00,920
loads that are applied.
1073
01:03:00,920 --> 01:03:06,780
Notice this Euclidean norm is
evaluated as shown down here
1074
01:03:06,780 --> 01:03:09,250
for a vector a.
1075
01:03:09,250 --> 01:03:12,580
The symbolism, two bars each
side with a 2 means nothing
1076
01:03:12,580 --> 01:03:16,440
else then taking the squares of
all the components of the
1077
01:03:16,440 --> 01:03:19,630
vector, adding those
squares up.
1078
01:03:19,630 --> 01:03:24,970
And then, taking the square
root out of that addition.
1079
01:03:24,970 --> 01:03:29,200
On displacements, we can also
measure convergence.
1080
01:03:29,200 --> 01:03:32,820
And here is the equation
that can be used.
1081
01:03:32,820 --> 01:03:34,860
The Euclidean norm on
the incremental
1082
01:03:34,860 --> 01:03:37,540
displacement vector.
1083
01:03:37,540 --> 01:03:40,970
And that is compared
with DNORM.
1084
01:03:40,970 --> 01:03:43,410
Again, DNORM.
1085
01:03:43,410 --> 01:03:47,870
Notice that this here is again,
input by the analyst to
1086
01:03:47,870 --> 01:03:49,900
the computer program.
1087
01:03:49,900 --> 01:03:53,150
It is part of the input that is
provided for the analysis.
1088
01:03:53,150 --> 01:03:55,280
And it's a reference
displacement.
1089
01:03:55,280 --> 01:03:58,300
Notice then, of course, we
should only include here the
1090
01:03:58,300 --> 01:04:01,670
displacements, translational
displacements in other words.
1091
01:04:01,670 --> 01:04:04,370
If we have also, rotational
degrees of freedom.
1092
01:04:04,370 --> 01:04:08,370
Then we would here include only
the rotational degrees of
1093
01:04:08,370 --> 01:04:13,780
freedom components and
use here DMNORM.
1094
01:04:13,780 --> 01:04:18,700
This way again, we compare one
quantity with certain units
1095
01:04:18,700 --> 01:04:21,580
with a quantity that
has the same units.
1096
01:04:21,580 --> 01:04:25,220
Once for translations and
once for rotations.
1097
01:04:25,220 --> 01:04:28,440
And so when we have rotational
degrees of freedom in the
1098
01:04:28,440 --> 01:04:32,230
analysis, we this way make
sure really, that the
1099
01:04:32,230 --> 01:04:35,710
incremental displacements and
incremental rotations are
1100
01:04:35,710 --> 01:04:38,140
small at convergence.
1101
01:04:38,140 --> 01:04:41,270
This then completes what I
wanted to share with you in
1102
01:04:41,270 --> 01:04:42,010
this lecture.
1103
01:04:42,010 --> 01:04:44,060
Thank you very much for
your attention.