WEBVTT
00:00:00.040 --> 00:00:02.470
The following content is
provided under a Creative
00:00:02.470 --> 00:00:03.880
Commons license.
00:00:03.880 --> 00:00:06.920
Your support will help MIT
OpenCourseWare continue to
00:00:06.920 --> 00:00:10.570
offer high-quality educational
resources for free.
00:00:10.570 --> 00:00:13.470
To make a donation, or view
additional materials from
00:00:13.470 --> 00:00:17.872
hundreds of MIT courses, visit
MIT OpenCourseWare at
00:00:17.872 --> 00:00:19.122
ocw.mit.edu.
00:00:21.230 --> 00:00:22.980
PROFESSOR: Ladies and gentlemen,
welcome to this
00:00:22.980 --> 00:00:25.550
lecture on nonlinear finite
element analysis of solids and
00:00:25.550 --> 00:00:26.820
structures.
00:00:26.820 --> 00:00:29.390
In this lecture, I like to
discuss with you solution
00:00:29.390 --> 00:00:32.580
methods that we use to solve the
finite element equations
00:00:32.580 --> 00:00:34.790
in nonlinear static analysis.
00:00:34.790 --> 00:00:38.520
In the previous lectures, we
derived this set of equations,
00:00:38.520 --> 00:00:42.020
where tK is a tangent
stiffness matrix.
00:00:42.020 --> 00:00:46.360
Delta Ui is a nodal point
vector of incremental
00:00:46.360 --> 00:00:49.760
displacements corresponding
to iteration i.
00:00:49.760 --> 00:00:54.210
t plus delta t R is the
externally applied load vector
00:00:54.210 --> 00:00:58.260
corresponding to time
t plus delta t.
00:00:58.260 --> 00:01:02.920
And t plus delta t F i minus 1
is equal to the nodal point
00:01:02.920 --> 00:01:06.530
force vector corresponding to
the internal elements stresses
00:01:06.530 --> 00:01:09.120
at time t plus delta t,
and at the end of
00:01:09.120 --> 00:01:12.060
iteration i minus 1.
00:01:12.060 --> 00:01:15.770
The displacements are updated
as shown here.
00:01:15.770 --> 00:01:17.620
Delta Ui, of course,
is calculated
00:01:17.620 --> 00:01:19.450
from this set of equations.
00:01:19.450 --> 00:01:22.850
We add delta Ui to the previous
displacements that
00:01:22.850 --> 00:01:25.210
corresponded to iteration
t plus delta --
00:01:25.210 --> 00:01:28.720
to time t plus delta t and
iteration i minus 1.
00:01:28.720 --> 00:01:31.880
It's the end of iteration
i minus 1.
00:01:31.880 --> 00:01:37.110
And this right-hand side gives
us the new displacement vector
00:01:37.110 --> 00:01:42.340
corresponding to iteration
i, end of iteration i.
00:01:42.340 --> 00:01:45.680
Now this is one scheme that it's
used to solve the finite
00:01:45.680 --> 00:01:47.770
element equations.
00:01:47.770 --> 00:01:50.220
We refer to this as
the mortified
00:01:50.220 --> 00:01:51.940
Newton-Raphson method.
00:01:51.940 --> 00:01:53.790
However, there are other
schemes as well.
00:01:53.790 --> 00:01:55.600
And schemes that are in certain
00:01:55.600 --> 00:01:58.550
problems much more effective.
00:01:58.550 --> 00:02:01.770
Since it is so important to use
an effective method for
00:02:01.770 --> 00:02:04.430
the solution of the finite
element equations because the
00:02:04.430 --> 00:02:07.470
costs can otherwise be very
high, we should be familiar
00:02:07.470 --> 00:02:09.360
with those other techniques.
00:02:09.360 --> 00:02:13.980
And I'd like to share
those now with you.
00:02:13.980 --> 00:02:16.400
Therefore, we will look in
this lecture at other
00:02:16.400 --> 00:02:18.810
techniques to solve the finite
element equations.
00:02:18.810 --> 00:02:22.240
And we need, of course, to
discuss convergence criteria.
00:02:22.240 --> 00:02:23.980
It's very important to use
00:02:23.980 --> 00:02:25.880
appropriate convergence criteria.
00:02:25.880 --> 00:02:30.850
Otherwise, you iterate too long
and/or on the other side,
00:02:30.850 --> 00:02:33.930
you may actually stop iterating
before you have
00:02:33.930 --> 00:02:36.250
achieved proper convergence.
00:02:36.250 --> 00:02:41.310
So let me then go over to my
view graphs and let us start
00:02:41.310 --> 00:02:42.580
the discussion here.
00:02:42.580 --> 00:02:47.316
The basic set of equations that
we would like to solve
00:02:47.316 --> 00:02:48.530
are given here.
00:02:48.530 --> 00:02:52.510
Notice R is the vector of
externally applied nodal point
00:02:52.510 --> 00:02:55.390
forces at time t plus delta t.
00:02:55.390 --> 00:03:00.580
And F is the vector of nodal
point forces corresponding to
00:03:00.580 --> 00:03:04.560
the internal element stresses
at time t plus delta t.
00:03:04.560 --> 00:03:07.870
Of course, this vector is
unknown, and we want to
00:03:07.870 --> 00:03:13.560
iterate somehow to find it and
make sure, of course, that at
00:03:13.560 --> 00:03:17.150
the end of that iteration, if we
have the proper vector F, R
00:03:17.150 --> 00:03:21.670
is equal to F, or R minus
F is equal to 0.
00:03:21.670 --> 00:03:24.570
We assume that the loading is
deformation-independent.
00:03:24.570 --> 00:03:29.820
If say, loading is
deformation-dependent, we can
00:03:29.820 --> 00:03:31.250
also deal with that situation.
00:03:31.250 --> 00:03:34.890
But we will have to add
additional terms to our
00:03:34.890 --> 00:03:35.980
iterative scheme.
00:03:35.980 --> 00:03:38.650
And I will get back to
that a little later.
00:03:38.650 --> 00:03:41.940
Notice that F in the total
Lagrangian formulation the way
00:03:41.940 --> 00:03:47.030
we have discussed it earlier is
evaluated by this product
00:03:47.030 --> 00:03:50.290
here integrated over the
original volume for a single
00:03:50.290 --> 00:03:54.270
element and for the UL
formulation, F is calculated
00:03:54.270 --> 00:03:55.310
as shown here.
00:03:55.310 --> 00:03:58.500
Here, of course, we integrate
over the current volume, or
00:03:58.500 --> 00:04:02.580
the volume that we actually
want to calculate.
00:04:02.580 --> 00:04:05.360
Notice in the iteration, of
course, this volume here would
00:04:05.360 --> 00:04:12.020
be t plus delta t V i minus 1
as we will discuss just now.
00:04:12.020 --> 00:04:15.980
Because in the iteration, we
always update this integral
00:04:15.980 --> 00:04:20.089
and that integral with the
iteration counter that we are
00:04:20.089 --> 00:04:24.950
having on the right-hand side
of the system of equations.
00:04:24.950 --> 00:04:31.170
The methods that we use are
based on the Newton-Raphson
00:04:31.170 --> 00:04:38.020
method, which is really used
very abundantly to find roots
00:04:38.020 --> 00:04:40.260
of an equation.
00:04:40.260 --> 00:04:41.960
A small historical note:
00:04:41.960 --> 00:04:46.360
Newton gave a first version
of the method in 1669.
00:04:46.360 --> 00:04:49.560
Raphson then generalized and
simplified the method
00:04:49.560 --> 00:04:52.860
actually, in 1690.
00:04:52.860 --> 00:04:55.940
Both mathematicians used
the same concept.
00:04:55.940 --> 00:04:59.290
And both algorithms gave really,
the same results.
00:04:59.290 --> 00:05:03.090
But it is very appropriate to
refer to the methods that
00:05:03.090 --> 00:05:06.490
we're using as the
Newton-Raphson method because
00:05:06.490 --> 00:05:11.010
Raphson really contributed quite
a bit to that method.
00:05:11.010 --> 00:05:14.110
So let us now consider a single
00:05:14.110 --> 00:05:16.220
Newton-Raphson iteration.
00:05:16.220 --> 00:05:19.030
The way you may have actually
encountered it already earlier
00:05:19.030 --> 00:05:23.220
in your studies of solution
methods for the roots of an
00:05:23.220 --> 00:05:24.730
algebraic equation.
00:05:24.730 --> 00:05:27.070
Basically, what you're
doing is this.
00:05:27.070 --> 00:05:32.510
You're saying if you have xi
minus 1, you calculate
00:05:32.510 --> 00:05:34.660
f at xi minus 1.
00:05:34.660 --> 00:05:39.900
You divide this value by the
f prime at xi minus 1.
00:05:39.900 --> 00:05:43.910
And then this right-hand side
gives you better value, a
00:05:43.910 --> 00:05:47.340
better approximation to the
root of the equation.
00:05:47.340 --> 00:05:50.540
Once you have xi, you repeat
the process and
00:05:50.540 --> 00:05:53.090
you get xi plus 1.
00:05:53.090 --> 00:05:58.630
You keep on repeating the
process until, basically f xi
00:05:58.630 --> 00:06:01.070
here is close to 0.
00:06:01.070 --> 00:06:03.980
Because then you have a root.
00:06:03.980 --> 00:06:17.470
Well, if we use that
Newton-Raphson formula, it is
00:06:17.470 --> 00:06:21.010
quite interesting to see how
it has been derived.
00:06:21.010 --> 00:06:25.130
We can write for any point xi,
a neighboring point xi minus
00:06:25.130 --> 00:06:29.970
1, directly this equation here
by a Taylor series expansion.
00:06:29.970 --> 00:06:33.450
And if we neglect the higher
order terms, we get that f of
00:06:33.450 --> 00:06:37.180
xi is approximately equal to
this relationship on the
00:06:37.180 --> 00:06:39.220
right-hand side.
00:06:39.220 --> 00:06:44.335
Well, since f of xi is supposed
to be 0, because we
00:06:44.335 --> 00:06:47.060
are looking for a 0 of
the equation, we
00:06:47.060 --> 00:06:48.430
set it equal to 0.
00:06:48.430 --> 00:06:50.860
And here you see directly
the formula that
00:06:50.860 --> 00:06:52.120
I showed you earlier.
00:06:52.120 --> 00:06:55.360
So this is a very quick
derivation of the
00:06:55.360 --> 00:06:59.020
Newton-Raphson procedure.
00:06:59.020 --> 00:07:03.040
Let us look at a mathematical
example to see how the method
00:07:03.040 --> 00:07:07.370
works, and to get some insight
into the technique.
00:07:07.370 --> 00:07:10.960
Let us choose a very simple
example, not much to do with
00:07:10.960 --> 00:07:12.680
finite element analysis.
00:07:12.680 --> 00:07:16.280
Let us say that f of x equals
sine x and that our starting
00:07:16.280 --> 00:07:19.080
value is equal to 2
for the iteration.
00:07:19.080 --> 00:07:21.770
You always have to choose, of
course, a starting value for
00:07:21.770 --> 00:07:23.030
the iteration.
00:07:23.030 --> 00:07:28.750
Then, in this column you see the
values that are calculated
00:07:28.750 --> 00:07:30.830
in the successive iterations.
00:07:30.830 --> 00:07:34.020
And we also show
here the error.
00:07:34.020 --> 00:07:37.700
It is interesting to observe
that when you are close to the
00:07:37.700 --> 00:07:40.210
root, you have quadratic
convergence.
00:07:40.210 --> 00:07:45.410
Meaning that the error here,
epsilon becomes an error.
00:07:45.410 --> 00:07:48.570
Epsilon squared here.
00:07:48.570 --> 00:07:51.880
This is summarized on
this view graph.
00:07:51.880 --> 00:07:58.790
Mathematically, if we have that
the error, Ei minus 1, is
00:07:58.790 --> 00:08:00.220
given by this equation here.
00:08:00.220 --> 00:08:02.200
Of course, there's an
approximate sign here.
00:08:02.200 --> 00:08:04.340
Then the error in the
next iteration
00:08:04.340 --> 00:08:07.906
will be of this order.
00:08:07.906 --> 00:08:12.590
So the convergence is really
very rapid once we are close
00:08:12.590 --> 00:08:15.660
to the root.
00:08:15.660 --> 00:08:19.170
Here we can see, however,
that if we are not
00:08:19.170 --> 00:08:21.300
close to the root--
00:08:21.300 --> 00:08:26.690
in other words, if we are too
far from the root, and we pick
00:08:26.690 --> 00:08:31.390
a bad value, then we don't
get to the desired
00:08:31.390 --> 00:08:34.080
root, which is pi.
00:08:34.080 --> 00:08:38.919
But we get a value that
is far away from pi.
00:08:38.919 --> 00:08:42.159
It is also root, but certainly
not the root that we are
00:08:42.159 --> 00:08:43.635
interested in.
00:08:43.635 --> 00:08:47.550
So Newton-Raphson iteration
is not always convergent
00:08:47.550 --> 00:08:51.420
directly, converging directly
to the result that we would
00:08:51.420 --> 00:08:52.960
like to obtain.
00:08:52.960 --> 00:08:57.090
There has to be some care when
we use that iterative scheme.
00:08:57.090 --> 00:09:02.720
And in fact, we have to be close
enough just to basically
00:09:02.720 --> 00:09:04.190
say what has to happen.
00:09:04.190 --> 00:09:07.110
We have to be close enough to
the root in order to have
00:09:07.110 --> 00:09:09.930
convergence to the root that
we're looking for.
00:09:09.930 --> 00:09:12.750
And ultimately, if we do
converge to that root, we will
00:09:12.750 --> 00:09:15.060
get quadratic convergence
to the root.
00:09:15.060 --> 00:09:17.420
And that is, of course,
very desirable.
00:09:17.420 --> 00:09:21.040
But the starting value is
critical, as you can see here
00:09:21.040 --> 00:09:22.400
from this view graph.
00:09:22.400 --> 00:09:27.050
Let us look pictorially at
the solution process.
00:09:27.050 --> 00:09:30.020
Here we have our
sine function.
00:09:30.020 --> 00:09:31.770
And we are looking
for this root.
00:09:34.280 --> 00:09:41.170
In iteration 1, we basically set
up a tangent to that sine
00:09:41.170 --> 00:09:44.380
function at the starting
value, x0.
00:09:44.380 --> 00:09:52.140
And we calculate x1 using this
tangent as another estimate to
00:09:52.140 --> 00:09:55.820
the value that we're
interested in.
00:09:55.820 --> 00:10:03.540
Now in the iteration 2, we lay
a tangent to the f function.
00:10:03.540 --> 00:10:06.520
Which, of course, is sine
x, at the new value
00:10:06.520 --> 00:10:09.050
corresponding to x1.
00:10:09.050 --> 00:10:12.520
And with this tangent,
we get to this value.
00:10:12.520 --> 00:10:15.150
And you can see already
that we got closer
00:10:15.150 --> 00:10:16.220
to the desired value.
00:10:16.220 --> 00:10:21.070
We were this far away
originally.
00:10:21.070 --> 00:10:23.460
Then that far away.
00:10:23.460 --> 00:10:26.910
And now we are only
that far away.
00:10:26.910 --> 00:10:29.950
Notice our sine function
has this value here.
00:10:29.950 --> 00:10:32.680
And certainly, it's not 0 yet.
00:10:32.680 --> 00:10:37.540
So we lay another tangent
at that point.
00:10:37.540 --> 00:10:42.120
And that's being done on
this view graph here.
00:10:42.120 --> 00:10:45.890
Another tangent, and now
we get much closer
00:10:45.890 --> 00:10:46.990
to the desired value.
00:10:46.990 --> 00:10:50.000
Right here, x3 is already
quite close.
00:10:50.000 --> 00:10:54.340
But not close enough as
shown by this value.
00:10:54.340 --> 00:10:59.380
It's not yet close enough
to 0, this length here.
00:10:59.380 --> 00:11:03.830
And we lay another tangent
to that point now.
00:11:03.830 --> 00:11:08.500
And that's shown on this
view graph now.
00:11:08.500 --> 00:11:12.110
And we can see with the blue
line being the tangent at this
00:11:12.110 --> 00:11:18.470
point, we get x4 very close
to the 0 value that we are
00:11:18.470 --> 00:11:20.770
actually interested in.
00:11:20.770 --> 00:11:23.840
Certainly on this view graph,
you can't see any difference
00:11:23.840 --> 00:11:28.290
anymore between x4 and
the desired value.
00:11:28.290 --> 00:11:31.220
So this is basically the process
that we are following
00:11:31.220 --> 00:11:34.510
through when we go from
iteration 1 to iteration 2 to
00:11:34.510 --> 00:11:38.540
iteration 3 and to
iteration 4.
00:11:38.540 --> 00:11:41.580
Beautiful convergence in this
particular example.
00:11:41.580 --> 00:11:47.290
And the reason being that x0 is
close enough to the desired
00:11:47.290 --> 00:11:50.730
value, the desired root.
00:11:50.730 --> 00:11:52.650
Which, itself, of course,
is close to x4.
00:11:52.650 --> 00:11:53.840
Four
00:11:53.840 --> 00:12:01.420
Well, if however, we were to
take a value that is not close
00:12:01.420 --> 00:12:07.420
enough and then we, as I showed
earlier, we do not
00:12:07.420 --> 00:12:11.050
converge to the actual root.
00:12:11.050 --> 00:12:13.020
And this is shown here.
00:12:13.020 --> 00:12:18.500
For example, if we take a
tangent at a point where f of
00:12:18.500 --> 00:12:24.920
prime is almost 0, then of
course, this tangent throws us
00:12:24.920 --> 00:12:27.080
far away from this root.
00:12:27.080 --> 00:12:31.220
And we will converge
to some other root.
00:12:31.220 --> 00:12:34.120
Of course, we don't want to use
an exact 0 here because we
00:12:34.120 --> 00:12:37.310
have f of x divided
by f prime of x.
00:12:37.310 --> 00:12:40.080
So with an exact 0, we could
not [UNINTELLIGIBLE].
00:12:40.080 --> 00:12:45.330
But a tangent that is close to
0, an f prime value that is
00:12:45.330 --> 00:12:49.220
close to 0, would throw us far
away from the desired root.
00:12:49.220 --> 00:12:53.180
And that shows some of the
difficulties using--
00:12:53.180 --> 00:12:59.040
in fact, almost all iterative
schemes that we may not get
00:12:59.040 --> 00:13:00.750
fast enough to our root.
00:13:00.750 --> 00:13:03.520
And in fact, we may never get
to the root that we are
00:13:03.520 --> 00:13:05.010
looking for.
00:13:05.010 --> 00:13:09.380
Well, let us then look at the
Newton-Raphson iteration for
00:13:09.380 --> 00:13:12.130
multiple degrees of freedom.
00:13:12.130 --> 00:13:15.350
Here we have our R minus
F that we want to
00:13:15.350 --> 00:13:17.700
be setting to 0.
00:13:17.700 --> 00:13:21.130
Or we want to rather, solve for
the displacements t plus
00:13:21.130 --> 00:13:27.400
delta t U, which are the
solution because this t plus
00:13:27.400 --> 00:13:32.060
delta t u displacements would
give us this force vector.
00:13:32.060 --> 00:13:36.610
And this force vector
equilibrates the R vector.
00:13:36.610 --> 00:13:41.580
Meaning f of U, this little
f of U, is equal to 0.
00:13:41.580 --> 00:13:44.990
And that, of course, is what
we want to achieve.
00:13:44.990 --> 00:13:47.130
Notice that we have here,
of course, now
00:13:47.130 --> 00:13:48.720
n degrees of freedom.
00:13:48.720 --> 00:13:51.660
Therefore, n equations
to solve.
00:13:51.660 --> 00:13:55.290
Previously, we only looked
at one equation.
00:13:55.290 --> 00:14:01.070
If we apply the same principle
as before to this little f,
00:14:01.070 --> 00:14:06.840
namely we expand f as a Taylor
series about t plus delta t U
00:14:06.840 --> 00:14:15.000
i minus 1, then we get
this equation on
00:14:15.000 --> 00:14:17.770
the right-hand side.
00:14:17.770 --> 00:14:20.300
And of course, on the left-hand
side, we have f of t
00:14:20.300 --> 00:14:22.120
plus delta t Ui.
00:14:22.120 --> 00:14:26.660
Notice they are higher terms,
which we will neglect just as
00:14:26.660 --> 00:14:28.630
we have done earlier for
the single equation.
00:14:28.630 --> 00:14:31.650
Because we are looking for
Taylor series expansion about
00:14:31.650 --> 00:14:34.790
t plus delta t U i minus 1.
00:14:34.790 --> 00:14:41.560
If we neglect this part here,
we obtain directly this
00:14:41.560 --> 00:14:44.010
equation here.
00:14:44.010 --> 00:14:47.020
0 on the left-hand side, because
once again, we are
00:14:47.020 --> 00:14:54.440
looking for the displacements
values for which f is 0.
00:14:54.440 --> 00:14:56.770
So therefore, we set this
deliberately equal to 0.
00:14:56.770 --> 00:15:00.750
And on the right-hand side, we
have f t plus delta t U i
00:15:00.750 --> 00:15:04.540
minus 1 plus partial f with
respect to U times the delta
00:15:04.540 --> 00:15:07.190
U. This delta U is unknown,
of course,
00:15:07.190 --> 00:15:08.420
that we want to calculate.
00:15:08.420 --> 00:15:12.910
And this delta U is nothing
else then the difference
00:15:12.910 --> 00:15:16.610
between the U value in iteration
i and the U value in
00:15:16.610 --> 00:15:19.140
iteration i minus 1.
00:15:19.140 --> 00:15:22.670
Notice that this f is
evaluated at t plus
00:15:22.670 --> 00:15:24.350
delta t U i minus 1.
00:15:27.510 --> 00:15:29.530
Let us look now at
this equation.
00:15:29.530 --> 00:15:36.620
And if we write it a little bit
more out, we find that on
00:15:36.620 --> 00:15:40.300
the left-hand side, we have,
of course, a vector of 0's.
00:15:40.300 --> 00:15:43.990
On the right-hand side, we
have a vector of f i
00:15:43.990 --> 00:15:47.070
components, f1 to fn.
00:15:47.070 --> 00:15:52.600
All of which are evaluated at
t plus delta t U i minus 1.
00:15:52.600 --> 00:15:55.940
And then we have here a matrix,
which is a square
00:15:55.940 --> 00:16:02.560
matrix in which the individual
elements are partials of f i
00:16:02.560 --> 00:16:04.240
with respect to Uj.
00:16:04.240 --> 00:16:09.020
For example, f1, partial f1
with respect to U1 here.
00:16:09.020 --> 00:16:13.990
Partial fn with respect to
U n here, et cetera.
00:16:13.990 --> 00:16:18.130
This matrix here will give us
a tangent stiffness matrix.
00:16:18.130 --> 00:16:21.120
Notice this matrix is multiplied
by the vector of
00:16:21.120 --> 00:16:24.090
incremental nodal point
displacements corresponding to
00:16:24.090 --> 00:16:25.340
iteration i.
00:16:28.840 --> 00:16:37.040
If we remember that f, the
little f, at t plus delta t U
00:16:37.040 --> 00:16:41.450
i minus 1 is equal to R minus
2 plus delta t F i minus 1.
00:16:41.450 --> 00:16:44.980
Capital F here, little
f there, remember?
00:16:44.980 --> 00:16:50.230
Then, the partial of little f
with respect to U at t plus
00:16:50.230 --> 00:16:53.733
delta t U i minus 1 is given
via this equation.
00:16:53.733 --> 00:16:58.720
Now, please recognize that
partial of R with respect to U
00:16:58.720 --> 00:17:03.870
is equal to 0 if the loads are
deformation-independent.
00:17:03.870 --> 00:17:06.770
Of course, if the loads depend
on the deformations, in other
00:17:06.770 --> 00:17:10.670
words, if the loads depend on
the U displacements, then this
00:17:10.670 --> 00:17:15.010
would not be 0, and we would
have to actually put a term in
00:17:15.010 --> 00:17:17.979
here, carry that term along
in the solution of
00:17:17.979 --> 00:17:20.069
the nonlinear equations.
00:17:20.069 --> 00:17:24.050
But for the moment,
we neglect this.
00:17:24.050 --> 00:17:27.210
We assume that the loads are
deformation-independent.
00:17:27.210 --> 00:17:30.230
And then, this term
is equal to 0.
00:17:30.230 --> 00:17:33.770
This here is now giving us the
tangent stiffness matrix.
00:17:33.770 --> 00:17:37.430
And it is written down here.
00:17:37.430 --> 00:17:39.870
The tangent stiffness matrix
is nothing else than the
00:17:39.870 --> 00:17:45.440
partials of capital F with
respect to the U's, to the
00:17:45.440 --> 00:17:46.690
displacements.
00:17:48.520 --> 00:17:51.990
The final result then is given
on this view graph.
00:17:51.990 --> 00:17:55.540
We substitute all the
information that I shared with
00:17:55.540 --> 00:18:00.790
you into the Taylor series
expansion around t plus delta
00:18:00.790 --> 00:18:02.470
t U i minus 1.
00:18:02.470 --> 00:18:05.110
We get directly this
equation here.
00:18:05.110 --> 00:18:10.300
Notice tangent stiffness matrix
corresponding to time t
00:18:10.300 --> 00:18:14.360
plus delta t and iteration
i minus 1.
00:18:14.360 --> 00:18:20.580
Delta U i corresponding to
iteration i we have to be very
00:18:20.580 --> 00:18:23.330
clear about this, that this is
an increment in displacement
00:18:23.330 --> 00:18:28.460
from time t plus delta t
iteration i minus 1 to
00:18:28.460 --> 00:18:30.220
iteration i.
00:18:30.220 --> 00:18:33.380
And of course, the nodal
point force vector.
00:18:33.380 --> 00:18:36.020
Externally applied nodal point
forces go in here.
00:18:36.020 --> 00:18:38.470
And this is the nodal point
force vector corresponding to
00:18:38.470 --> 00:18:43.090
the internal element stresses at
time t plus delta t and at
00:18:43.090 --> 00:18:45.090
the end of iteration
i minus 1.
00:18:45.090 --> 00:18:49.100
Evaluated differently in the
total Lagrangian formulation
00:18:49.100 --> 00:18:52.480
from the way we're evaluating
it in the updated Lagrangian
00:18:52.480 --> 00:18:53.070
formulation.
00:18:53.070 --> 00:18:55.690
We talked about this vector
abundantly in
00:18:55.690 --> 00:18:57.330
the previous lectures.
00:18:57.330 --> 00:19:00.510
We talked about this vector,
of course, also in the
00:19:00.510 --> 00:19:01.630
previous lectures.
00:19:01.630 --> 00:19:05.970
We had in the previous lectures
here at tK, a
00:19:05.970 --> 00:19:10.260
constant tangent stiffness
matrix, which was set up at
00:19:10.260 --> 00:19:14.210
the beginning of this whole
interactive process, and was
00:19:14.210 --> 00:19:16.030
never updated.
00:19:16.030 --> 00:19:19.850
Well, here we now update
that stiffness matrix.
00:19:19.850 --> 00:19:23.650
And I think if you look through
the information that I
00:19:23.650 --> 00:19:26.080
just discussed with you,
you recognize why
00:19:26.080 --> 00:19:27.360
we're updating it.
00:19:27.360 --> 00:19:29.840
We are doing so because we are
always starting with the new
00:19:29.840 --> 00:19:37.120
Taylor series expansion about
the point i minus 1.
00:19:37.120 --> 00:19:41.990
This set of equations if of
course solved for delta Ui.
00:19:41.990 --> 00:19:44.070
We add delta Ui to what we
had already in terms of
00:19:44.070 --> 00:19:48.020
displacements to get our new
estimate for the nodal point
00:19:48.020 --> 00:19:50.160
displacements.
00:19:50.160 --> 00:19:54.450
It is important to realize that
the K matrix, which we
00:19:54.450 --> 00:19:57.440
are using in the solution
process, is symmetric.
00:19:57.440 --> 00:20:01.190
Because first of all, we use
symmetric stress and strain
00:20:01.190 --> 00:20:03.900
measures in our governing
equation.
00:20:03.900 --> 00:20:06.990
Remember that when we apply the
principle of virtual work,
00:20:06.990 --> 00:20:10.760
we use Cauchy stresses
and infinitesimally
00:20:10.760 --> 00:20:12.250
small virtual strains.
00:20:12.250 --> 00:20:16.850
Now notice that both of these
tenses are symmetric tenses.
00:20:16.850 --> 00:20:19.980
Of course, these are the tenses
we're using in the UL,
00:20:19.980 --> 00:20:22.280
in the updated Lagrangian
formulation.
00:20:22.280 --> 00:20:25.180
In the total Lagrangian
formulation, we use the second
00:20:25.180 --> 00:20:28.340
Piola-Kirchhoff stress in
Green-Lagrange strain.
00:20:28.340 --> 00:20:32.350
Both of these measures are
again, symmetric measures.
00:20:32.350 --> 00:20:34.650
Both tenses are symmetric
tenses.
00:20:34.650 --> 00:20:38.880
If we had introduced for a
formulation non-symmetric
00:20:38.880 --> 00:20:42.210
tenses, for example a
non-symmetric stress tenser,
00:20:42.210 --> 00:20:44.960
and of course, an energy
conjugate non-symmetric strain
00:20:44.960 --> 00:20:48.780
tenser, we would have obtained
a non-symmetric K matrix,
00:20:48.780 --> 00:20:53.500
which would be much more
difficult to deal with.
00:20:53.500 --> 00:20:57.680
Much less cost effective in
the solution process.
00:20:57.680 --> 00:21:01.470
Also, please recognize that
we interpolated the real
00:21:01.470 --> 00:21:05.810
displacements and the virtual
displacements with exactly the
00:21:05.810 --> 00:21:07.820
same functions.
00:21:07.820 --> 00:21:11.020
Whereas this point here is a
continuum mechanics point,
00:21:11.020 --> 00:21:13.680
this is really a finite
element point.
00:21:13.680 --> 00:21:17.590
If we had used different
interpolation functions for
00:21:17.590 --> 00:21:21.930
the real displacements, then for
the virtual displacements,
00:21:21.930 --> 00:21:26.300
or vice-versa, then we would
have not necessarily obtained
00:21:26.300 --> 00:21:27.860
a symmetric K matrix.
00:21:27.860 --> 00:21:32.110
Of course, the most natural
procedure is to use the same
00:21:32.110 --> 00:21:35.460
kind of interpolation functions
for the virtual
00:21:35.460 --> 00:21:37.940
displacements as we used for
the real displacements, and
00:21:37.940 --> 00:21:39.340
that's what we did.
00:21:39.340 --> 00:21:42.460
Finally, we also assumed
that the loading was
00:21:42.460 --> 00:21:43.710
deformation-independent.
00:21:45.880 --> 00:21:49.000
If we have deformation-dependent
loading,
00:21:49.000 --> 00:21:52.680
then if you go more back to the
earlier view graph, then
00:21:52.680 --> 00:21:57.020
of course, this right-hand
side vector R here would
00:21:57.020 --> 00:21:59.750
depend on the displacements.
00:21:59.750 --> 00:22:02.630
And there are basically two
different approaches
00:22:02.630 --> 00:22:04.310
that one can take.
00:22:04.310 --> 00:22:08.480
In the first approach, one
simply updates this vector
00:22:08.480 --> 00:22:11.770
with the iteration, or with
taking into account the
00:22:11.770 --> 00:22:15.910
current or last calculated
volume and surface areas for
00:22:15.910 --> 00:22:16.920
the elements.
00:22:16.920 --> 00:22:21.430
So we would simply put here
an i minus 1 up there.
00:22:21.430 --> 00:22:24.660
The left-hand side matrix,
we leave unchanged.
00:22:24.660 --> 00:22:27.630
And if that converges fast,
certainly it's a very
00:22:27.630 --> 00:22:29.540
effective approach to use.
00:22:29.540 --> 00:22:32.950
We use that approach
abundantly.
00:22:32.950 --> 00:22:36.490
However, another approach
is to actually take the
00:22:36.490 --> 00:22:39.980
differentiation of R with
respect to the U's, the way we
00:22:39.980 --> 00:22:41.890
have been looking at it earlier
on an earlier view
00:22:41.890 --> 00:22:45.260
graph, and then you get some
components that you
00:22:45.260 --> 00:22:47.110
add to the K matrix.
00:22:47.110 --> 00:22:51.200
And that K matrix may then
out to be non-symmetric.
00:22:51.200 --> 00:22:54.400
However, if the loading is
deformation-independent, then
00:22:54.400 --> 00:22:57.910
the differentiation of R with
respect to the displacements U
00:22:57.910 --> 00:23:01.770
is equal to 0, and there's
no component from that
00:23:01.770 --> 00:23:04.250
differentiation coming
into the K matrix.
00:23:04.250 --> 00:23:07.550
And our K, of course, provided
these are also
00:23:07.550 --> 00:23:10.850
satisfied, is symmetric.
00:23:10.850 --> 00:23:13.840
The iterative scheme that we
just discussed is referred to
00:23:13.840 --> 00:23:16.780
as the full Newton-Raphson
method.
00:23:16.780 --> 00:23:19.630
Full because we are setting
up a new K matrix at the
00:23:19.630 --> 00:23:22.340
beginning of each iteration.
00:23:22.340 --> 00:23:26.120
The full Newton-Raphson method
shows mathematically quadratic
00:23:26.120 --> 00:23:28.650
convergence the way we discussed
it a bit earlier in
00:23:28.650 --> 00:23:31.640
the lecture.
00:23:31.640 --> 00:23:36.400
And that, of course, is always
the case provided you are
00:23:36.400 --> 00:23:38.050
close enough to the root.
00:23:38.050 --> 00:23:41.550
This quadratic convergence only
holds provided you are
00:23:41.550 --> 00:23:45.780
close enough to the root when
you solve your equations.
00:23:45.780 --> 00:23:49.420
In finite element analysis, it
is also important to recognize
00:23:49.420 --> 00:23:52.410
that a number of requirements
must be fulfilled.
00:23:52.410 --> 00:23:57.410
For example, in elasto-plastic
analysis, the stresses and
00:23:57.410 --> 00:23:59.020
strains must be properly--
00:23:59.020 --> 00:24:03.520
plastic strains must be
properly updated.
00:24:03.520 --> 00:24:07.050
And similarly, the rotations
in a shell analysis must be
00:24:07.050 --> 00:24:07.870
properly updated.
00:24:07.870 --> 00:24:11.380
So it is not necessarily the
case that you automatically
00:24:11.380 --> 00:24:14.780
get quadratic convergence when
you do finite element
00:24:14.780 --> 00:24:18.250
announces with a full
Newton-Raphson method.
00:24:18.250 --> 00:24:22.120
It is very important to also,
on the level of updating the
00:24:22.120 --> 00:24:27.800
stresses and the rotations,
really do things properly in
00:24:27.800 --> 00:24:33.040
quotes in order to obtain the
full quadratic convergence of
00:24:33.040 --> 00:24:36.900
the Newton-Raphson method.
00:24:36.900 --> 00:24:39.580
We can depict the iteration
process in
00:24:39.580 --> 00:24:41.350
two equivalent ways.
00:24:41.350 --> 00:24:45.300
The first way is shown
here, the left.
00:24:45.300 --> 00:24:48.520
And it's really the way we've
been discussing just now the
00:24:48.520 --> 00:24:54.280
solution of f equals R minus
capital F. Where we want to
00:24:54.280 --> 00:24:58.150
solve for the root, the
0 of this equation.
00:24:58.150 --> 00:25:02.980
Notice here we have in red
the little f depicted.
00:25:02.980 --> 00:25:07.280
Notice that at this point here,
we have t plus delta t
00:25:07.280 --> 00:25:10.650
capital F i minus 1.
00:25:10.650 --> 00:25:14.180
And t plus delta t
R is this value.
00:25:14.180 --> 00:25:18.130
Now, as we get closer to the
root, which is of course, the
00:25:18.130 --> 00:25:23.800
point where the red curve
crosses this U-axis, as we get
00:25:23.800 --> 00:25:28.620
closer to the root, this capital
F will get closer and
00:25:28.620 --> 00:25:33.140
closer to the R. And that is,
of course, when the little f
00:25:33.140 --> 00:25:34.310
is equal to 0.
00:25:34.310 --> 00:25:36.080
That's all we are
showing here.
00:25:36.080 --> 00:25:41.830
Notice that we are setting up
a slope f prime at the point
00:25:41.830 --> 00:25:44.330
corresponding to the i minus
first iteration.
00:25:44.330 --> 00:25:49.180
And that slope brings us into
this point, which will be the
00:25:49.180 --> 00:25:51.240
next point for--
00:25:51.240 --> 00:25:55.330
or the next starting point
for the next iteration.
00:25:55.330 --> 00:25:57.930
I should say, the point
of starting
00:25:57.930 --> 00:25:59.600
with the next iteration.
00:25:59.600 --> 00:26:02.890
Now, this is one way
of looking at
00:26:02.890 --> 00:26:04.420
the iterative process.
00:26:04.420 --> 00:26:06.960
We can also look at
it as shown here.
00:26:06.960 --> 00:26:10.830
Notice we have here the R
plotted now horizontally.
00:26:10.830 --> 00:26:13.240
The displacement vertically.
00:26:13.240 --> 00:26:19.680
Notice the R at a particular
time t plus delta t is shown
00:26:19.680 --> 00:26:22.160
by this dashed line.
00:26:22.160 --> 00:26:27.110
And we really want to find this
particular point here and
00:26:27.110 --> 00:26:30.970
the corresponding U displacement
of course.
00:26:30.970 --> 00:26:34.380
At this point, the little f is
0 and the corresponding U
00:26:34.380 --> 00:26:36.250
displacement is down here.
00:26:36.250 --> 00:26:40.240
Notice that at the iteration
i minus 1, it's the end of
00:26:40.240 --> 00:26:44.570
iteration i minus 1, we have
obtained this point here.
00:26:44.570 --> 00:26:47.850
The U displacement corresponding
to that point is
00:26:47.850 --> 00:26:52.180
obtained by projecting down
on the displacement axes.
00:26:52.180 --> 00:26:57.800
And this slope here, the blue
line, gives us a tangent
00:26:57.800 --> 00:27:00.350
stiffness matrix slope.
00:27:00.350 --> 00:27:02.700
These are two equivalent ways.
00:27:02.700 --> 00:27:06.070
This is here more like a force
deflection curve, and we will
00:27:06.070 --> 00:27:10.270
use this representation
now abundantly.
00:27:10.270 --> 00:27:13.500
But keep in mind, it's really
the same thing if you look at
00:27:13.500 --> 00:27:15.300
it closely.
00:27:15.300 --> 00:27:17.890
Well, modifications
to the iterative
00:27:17.890 --> 00:27:21.410
scheme are the following.
00:27:21.410 --> 00:27:25.100
If we leave the stiffness matrix
constant throughout a
00:27:25.100 --> 00:27:26.670
complete solution process--
00:27:26.670 --> 00:27:30.070
in other words, tau is
equal to 0 here.
00:27:30.070 --> 00:27:31.390
It's never updated.
00:27:31.390 --> 00:27:34.570
We talk about the initial
stress method.
00:27:34.570 --> 00:27:39.610
If tau is equal to t, where t,
of course, corresponds to a
00:27:39.610 --> 00:27:44.670
particular solution step and
meaning that the K matrix is
00:27:44.670 --> 00:27:48.620
constant, but it is updated
as the beginning
00:27:48.620 --> 00:27:50.010
of a solution step.
00:27:50.010 --> 00:27:53.160
Then we talk about the modified
Newton method.
00:27:53.160 --> 00:27:55.860
Modified Newton-Raphson
method.
00:27:55.860 --> 00:27:59.670
Or, we may also find it
effective to update the K
00:27:59.670 --> 00:28:04.040
matrix only at particular
solution steps, at certain
00:28:04.040 --> 00:28:06.570
times only.
00:28:06.570 --> 00:28:09.590
We note that the initial stress
method and the modified
00:28:09.590 --> 00:28:13.810
Newton methods are, of course,
much less expensive than the
00:28:13.810 --> 00:28:16.646
full Newton method per
iteration, however.
00:28:16.646 --> 00:28:21.070
We should add that many more
iterations are necessary to
00:28:21.070 --> 00:28:27.150
achieve the same accuracy if we
don't set up a new K matrix
00:28:27.150 --> 00:28:29.210
in every iteration.
00:28:29.210 --> 00:28:32.030
The initial stress method and
modified Newton iteration
00:28:32.030 --> 00:28:36.240
technique do not exhibit
quadratic convergence because
00:28:36.240 --> 00:28:38.770
to obtain quadratic convergence,
you need to set
00:28:38.770 --> 00:28:41.750
up a new K matrix in
each iteration.
00:28:41.750 --> 00:28:44.920
And of course, you have to be
close enough to the root, to
00:28:44.920 --> 00:28:46.860
the solution that you're
looking for.
00:28:46.860 --> 00:28:51.200
Let us now look at
a simple example.
00:28:51.200 --> 00:28:55.250
An example where we have just
one degree of freedom.
00:28:55.250 --> 00:28:57.790
And where we deal with
two load steps.
00:28:57.790 --> 00:29:01.580
Here we show the force
applied to the
00:29:01.580 --> 00:29:04.210
problem, or to the structure.
00:29:04.210 --> 00:29:07.460
In the first load step,
we apply 1R.
00:29:07.460 --> 00:29:10.530
And then in the second
load step, 2R.
00:29:10.530 --> 00:29:14.570
Notice horizontally here, I'm
plotting displacements.
00:29:14.570 --> 00:29:17.460
And the force displacement
curve is
00:29:17.460 --> 00:29:20.130
shown by the red line.
00:29:20.130 --> 00:29:26.020
Now, to obtain the solution for
this displacement, U1, and
00:29:26.020 --> 00:29:31.340
that displacement U2, we, as
we discussed, set up a K
00:29:31.340 --> 00:29:34.600
matrix, which corresponds
to the slope, a
00:29:34.600 --> 00:29:36.940
slope on the red curve.
00:29:36.940 --> 00:29:40.490
And we iterate towards
the sought solutions.
00:29:40.490 --> 00:29:44.680
Here, we have the iterative
process using the initial
00:29:44.680 --> 00:29:45.950
stress method.
00:29:45.950 --> 00:29:49.540
In other words, where tau is
equal to 0 in our governing
00:29:49.540 --> 00:29:52.680
finite element equations.
00:29:52.680 --> 00:29:58.380
This means that we are setting
up a K matrix initially, and
00:29:58.380 --> 00:30:04.960
that K matrix is depicted here
by the slope, the blue slope
00:30:04.960 --> 00:30:06.280
that you're seeing.
00:30:06.280 --> 00:30:12.380
Now with this slope and this
load applied, we get the
00:30:12.380 --> 00:30:14.390
displacement shown down here.
00:30:14.390 --> 00:30:19.125
Notice the left superscript
means load step 1.
00:30:19.125 --> 00:30:24.200
The right superscript means
solution after iteration 1.
00:30:24.200 --> 00:30:28.980
And this is the solution
obtained from this triangle.
00:30:28.980 --> 00:30:33.990
Notice there's a blue line up
there, going up there, and
00:30:33.990 --> 00:30:39.650
where that blue line crosses the
dashed line, we pick up a
00:30:39.650 --> 00:30:41.730
vertical fine, black line.
00:30:41.730 --> 00:30:45.890
And that vertical black line
cuts through the displacement
00:30:45.890 --> 00:30:48.570
axes at this particular value.
00:30:48.570 --> 00:30:51.600
So this is the solution of
our first iteration.
00:30:51.600 --> 00:30:55.740
Now, having calculated this
displacement value, we go up
00:30:55.740 --> 00:31:00.145
vertically and get to the
red line, the red curve.
00:31:00.145 --> 00:31:04.710
And that red curve corresponds,
of course, to the
00:31:04.710 --> 00:31:06.560
internal forces.
00:31:06.560 --> 00:31:09.340
To the internal force I should
say, corresponding to this
00:31:09.340 --> 00:31:11.320
displacement.
00:31:11.320 --> 00:31:16.200
There's a red dot right
on that red curve.
00:31:16.200 --> 00:31:20.690
And at that red dot now, we've
set up a K matrix again.
00:31:20.690 --> 00:31:24.440
Now notice, in the initial
stress method, we keep the K
00:31:24.440 --> 00:31:26.000
matrix constant.
00:31:26.000 --> 00:31:32.410
This means that the slope here
of this blue line that goes
00:31:32.410 --> 00:31:37.830
through this red point is the
same as the slope that we had
00:31:37.830 --> 00:31:40.060
earlier right here.
00:31:40.060 --> 00:31:44.670
With that slope and the out of
balance load-- and let's look
00:31:44.670 --> 00:31:46.330
now very closely here--
00:31:46.330 --> 00:31:50.190
that corresponds to the distance
between this red
00:31:50.190 --> 00:31:56.330
point and that dashed line with
this out of balance load
00:31:56.330 --> 00:32:00.200
and this blue slope, we get an
increment in displacement
00:32:00.200 --> 00:32:02.730
shown right here.
00:32:02.730 --> 00:32:05.260
And that increment in
displacement added to the
00:32:05.260 --> 00:32:11.780
previous displacement gives
us this value down here.
00:32:11.780 --> 00:32:16.240
And in fact, this value is
very close to the correct
00:32:16.240 --> 00:32:20.380
solution, the exact solution,
which is the solution
00:32:20.380 --> 00:32:24.750
corresponding to the dashed
vertical line here.
00:32:24.750 --> 00:32:27.180
Notice that the dashed vertical
line is the exact
00:32:27.180 --> 00:32:30.160
solution to that dashed
horizontal line.
00:32:30.160 --> 00:32:34.920
And our vertical black line,
which we are calculating, is
00:32:34.920 --> 00:32:37.970
virtually on that dashed
black line.
00:32:37.970 --> 00:32:40.900
So we can accept now,
convergence.
00:32:40.900 --> 00:32:42.960
Let's now go into the
next load step.
00:32:42.960 --> 00:32:46.330
In the next load step,
we want to, again,
00:32:46.330 --> 00:32:47.620
deal with a K matrix.
00:32:47.620 --> 00:32:50.640
And in the initial stress
method, we use the same K
00:32:50.640 --> 00:32:52.970
matrix throughout
the solution.
00:32:52.970 --> 00:32:56.980
Meaning that the slope now,
which we are laying onto the
00:32:56.980 --> 00:33:02.860
red curve, that slope being this
blue line, is the same
00:33:02.860 --> 00:33:08.470
slope that we used before in
the first two iterations.
00:33:08.470 --> 00:33:13.750
With this slope and the out of
balance load corresponding to
00:33:13.750 --> 00:33:17.720
the distance between this
horizontal line, the dashed
00:33:17.720 --> 00:33:21.780
horizontal line, and that dashed
horizontal line because
00:33:21.780 --> 00:33:25.040
we conversed in the
first load step.
00:33:25.040 --> 00:33:27.450
With that out of balance load,
we get an incremental
00:33:27.450 --> 00:33:34.610
displacement shown from
here to there.
00:33:34.610 --> 00:33:37.220
And that incremental
displacement then, added to
00:33:37.220 --> 00:33:39.990
the previous displacement,
gives us this value in
00:33:39.990 --> 00:33:40.720
displacement.
00:33:40.720 --> 00:33:46.390
Notice time t plus delta t or
load step t plus delta t,
00:33:46.390 --> 00:33:47.600
which is here.
00:33:47.600 --> 00:33:52.620
Load step 2 and iteration
1, end of iteration 1.
00:33:52.620 --> 00:33:54.960
This is now our current
displacement.
00:33:54.960 --> 00:33:58.670
With this current displacement,
once again,
00:33:58.670 --> 00:34:05.760
which comes from this triangle
here, with this current
00:34:05.760 --> 00:34:10.610
displacement, we go vertically
up and get to that red point,
00:34:10.610 --> 00:34:13.940
which lies on the force
displacement curve--
00:34:13.940 --> 00:34:17.500
internal force displacement
curve more accurately.
00:34:17.500 --> 00:34:24.510
And now we have this red point,
and we set up, or we
00:34:24.510 --> 00:34:27.170
use now, again, a K matrix.
00:34:27.170 --> 00:34:29.110
Of course, in the initial
stress method we use the
00:34:29.110 --> 00:34:30.110
constant K matrix.
00:34:30.110 --> 00:34:32.949
So this slope is the
same as that slope.
00:34:32.949 --> 00:34:35.320
And now watch closely.
00:34:35.320 --> 00:34:42.960
From this triangle here, which
corresponds to this out of
00:34:42.960 --> 00:34:45.130
balance load, which is
the same as that
00:34:45.130 --> 00:34:47.820
out of balance load.
00:34:47.820 --> 00:34:50.580
Out of balance load obtained
by putting here
00:34:50.580 --> 00:34:53.469
a horizontal line.
00:34:53.469 --> 00:34:57.100
Notice and the difference
between the dashed line up
00:34:57.100 --> 00:35:00.965
there and the horizontal line
given by my pointer is the out
00:35:00.965 --> 00:35:02.270
of balance load.
00:35:02.270 --> 00:35:05.820
Which of course, is the same as
the distance between this
00:35:05.820 --> 00:35:08.430
red point and that
dashed line.
00:35:08.430 --> 00:35:12.580
This out of balance load with
this slope gives us another
00:35:12.580 --> 00:35:15.190
increment in displacement,
which brings us
00:35:15.190 --> 00:35:17.930
to this value here.
00:35:17.930 --> 00:35:22.660
We repeat that process and
iteratively, we get closer and
00:35:22.660 --> 00:35:26.340
closer to the correct solution,
the exact solution,
00:35:26.340 --> 00:35:30.050
which is the dashed vertical
line here for the
00:35:30.050 --> 00:35:31.700
displacement.
00:35:31.700 --> 00:35:36.040
Because at that dashed vertical
line, we are getting
00:35:36.040 --> 00:35:38.660
to that red curve.
00:35:38.660 --> 00:35:41.070
Now, this is the initial
stress method.
00:35:41.070 --> 00:35:46.030
And once again, the point here
is that we kept the same K
00:35:46.030 --> 00:35:48.900
matrix throughout the
solution process.
00:35:48.900 --> 00:35:55.710
Surely, if you now were to look
at this iterative scheme
00:35:55.710 --> 00:35:59.120
to obtain a faster solution,
you would want to
00:35:59.120 --> 00:36:00.770
set up a new K matrix.
00:36:00.770 --> 00:36:04.490
In other words, once you have
this displacement value
00:36:04.490 --> 00:36:12.670
calculated, you want to update
the slope corresponding to the
00:36:12.670 --> 00:36:16.180
slope of the red curve.
00:36:16.180 --> 00:36:18.980
And if you do so--
00:36:18.980 --> 00:36:22.470
in other words, you update this
blue slope corresponding
00:36:22.470 --> 00:36:26.980
to the red dots that you
have obtained in
00:36:26.980 --> 00:36:28.640
your iterative scheme.
00:36:28.640 --> 00:36:31.300
These red dots, of course,
being different from what
00:36:31.300 --> 00:36:32.740
you're seeing here now.
00:36:32.740 --> 00:36:37.440
Then you would have the full
Newton-Raphson method.
00:36:37.440 --> 00:36:42.750
If you only update the K matrix
or the blue tangent
00:36:42.750 --> 00:36:47.230
that we are seeing here whenever
you have converged to
00:36:47.230 --> 00:36:50.680
an equilibrium point.
00:36:50.680 --> 00:36:54.870
Meaning when you have converged
corresponding for
00:36:54.870 --> 00:36:58.710
the first load step in
this particular case.
00:36:58.710 --> 00:37:01.120
If you only update the stiffness
matrix once you have
00:37:01.120 --> 00:37:03.560
converged to this load--
00:37:03.560 --> 00:37:06.690
for this load step, then you
would have the modified
00:37:06.690 --> 00:37:09.060
Newton-Raphson method.
00:37:09.060 --> 00:37:13.420
I think if you look at this
picture with these
00:37:13.420 --> 00:37:16.700
possibilities, you will realize
that surely, the full
00:37:16.700 --> 00:37:21.080
Newton-Raphson method with
updating the slope after each
00:37:21.080 --> 00:37:23.970
iteration will converge
fastest.
00:37:23.970 --> 00:37:26.690
And the modified Newton-Raphson
method will
00:37:26.690 --> 00:37:29.150
converge a little slower than
the full Newton-Raphson
00:37:29.150 --> 00:37:34.260
method, but still faster than
the initial stress method.
00:37:34.260 --> 00:37:38.210
Well, this then brings us to the
end of the discussion of
00:37:38.210 --> 00:37:39.080
this example.
00:37:39.080 --> 00:37:42.800
I like to now introduce you to
another scheme that can be
00:37:42.800 --> 00:37:46.430
very important when we solve
nonlinear equations, finite
00:37:46.430 --> 00:37:47.770
element equations.
00:37:47.770 --> 00:37:52.290
And that scheme is the scheme
of line searches.
00:37:52.290 --> 00:37:55.480
The basic equations that we are
looking at that we want to
00:37:55.480 --> 00:37:57.380
solve are once more
depicted here.
00:37:57.380 --> 00:37:59.990
Tangent stiffness matrix,
incremental displacement
00:37:59.990 --> 00:38:03.350
vector, externally applied
loads, nodal point forces
00:38:03.350 --> 00:38:05.790
corresponding to the internal
element stresses at time t
00:38:05.790 --> 00:38:10.140
plus delta t, and at the end
of iteration i minus 1,
00:38:10.140 --> 00:38:12.640
beginning of iteration
i, of course.
00:38:12.640 --> 00:38:15.070
Notice I don't have an i on
here, I just want to denote
00:38:15.070 --> 00:38:18.220
this for the moment as an
incremental nodal point
00:38:18.220 --> 00:38:22.200
displacement vector, which
carries, however, a bar.
00:38:22.200 --> 00:38:29.560
Let us consider forming Fi using
this equation where beta
00:38:29.560 --> 00:38:31.800
is an unknown.
00:38:31.800 --> 00:38:35.460
We want to choose beta
now judiciously.
00:38:35.460 --> 00:38:40.380
We want to choose beta such
that, in some sense, R minus
00:38:40.380 --> 00:38:44.580
Fi is small.
00:38:44.580 --> 00:38:48.515
And we have to, of course, look
at the mechanism that we
00:38:48.515 --> 00:38:53.560
use to make R minus F
small in some sense.
00:38:53.560 --> 00:38:59.040
Well, as a side note, we
recognize that when this
00:38:59.040 --> 00:39:06.180
equation here is 0 for all
possible U's, then clearly
00:39:06.180 --> 00:39:08.640
this must be 0.
00:39:08.640 --> 00:39:14.380
In other words, if we were to
say to ourselves, let's make
00:39:14.380 --> 00:39:17.410
this equation 0 for
all possible U's.
00:39:17.410 --> 00:39:22.770
Then really, we would have
solved this equal to 0.
00:39:22.770 --> 00:39:25.640
The reason, of course, for it
is that we could choose here
00:39:25.640 --> 00:39:31.310
for U a vector, which carries
0's everywhere except a 1 in
00:39:31.310 --> 00:39:32.610
one location.
00:39:32.610 --> 00:39:36.880
And if that is a particular
row, then corresponding to
00:39:36.880 --> 00:39:44.170
that row, surely this vector
here, this R minus F
00:39:44.170 --> 00:39:46.900
corresponding to that row
must be equal to 0.
00:39:46.900 --> 00:39:50.390
So if we were simply to use here
the identity matrix, or
00:39:50.390 --> 00:39:54.660
if we were to construct the
identity matrix by allowing n
00:39:54.660 --> 00:39:57.730
U's that are linearly
independent, of course,
00:39:57.730 --> 00:40:03.480
because they all make up such
vector but with a 1 in a
00:40:03.480 --> 00:40:04.500
different location.
00:40:04.500 --> 00:40:08.950
As a matter of fact, we want to
use U1 with a 1 here, 0's
00:40:08.950 --> 00:40:09.850
everywhere.
00:40:09.850 --> 00:40:14.650
U2, 1 here, 0 here, everywhere
else 0, and so on.
00:40:14.650 --> 00:40:16.870
To construct an identity
matrix here.
00:40:16.870 --> 00:40:20.500
Then with that identity matrix,
clearly R minus F
00:40:20.500 --> 00:40:23.090
would be 0 everywhere.
00:40:23.090 --> 00:40:27.250
Of course, that is not a
viable scheme really.
00:40:27.250 --> 00:40:34.170
But what we want to do is to
choose a particular vector
00:40:34.170 --> 00:40:40.430
here so as to make R minus
F still 0 in some sense.
00:40:40.430 --> 00:40:44.760
And the vector that is quite
effectively chosen is the
00:40:44.760 --> 00:40:49.460
vector that we obtain from the
solution of K delta U bar
00:40:49.460 --> 00:40:51.430
equals R minus F.
00:40:51.430 --> 00:40:58.070
If we use that vector here,
basically projecting the R,
00:40:58.070 --> 00:41:01.650
capital R, minus F
onto that vector.
00:41:01.650 --> 00:41:06.770
And we look for this to be 0.
00:41:06.770 --> 00:41:12.350
Then we have so to say, searched
in the direction of
00:41:12.350 --> 00:41:20.210
delta U bar and made delta U
bar transposed times this
00:41:20.210 --> 00:41:22.360
vector here equal to 0.
00:41:22.360 --> 00:41:23.970
How does this scheme work?
00:41:23.970 --> 00:41:29.240
Well, we add beta times delta U
bar to the value that we had
00:41:29.240 --> 00:41:31.600
already from the previous
iteration.
00:41:31.600 --> 00:41:39.340
And we search along beta such
that this top here is much
00:41:39.340 --> 00:41:41.360
smaller than the bottom.
00:41:41.360 --> 00:41:45.650
In other words, STOL shall be
a convergence tolerance.
00:41:45.650 --> 00:41:47.940
Notice the way it works.
00:41:47.940 --> 00:41:49.630
We choose a value of beta.
00:41:49.630 --> 00:41:53.550
With beta known, we add this
quantity to this one.
00:41:53.550 --> 00:41:55.470
We get a value here.
00:41:55.470 --> 00:42:00.470
We take this value, calculate
Fi corresponding to that
00:42:00.470 --> 00:42:02.510
value, see whether
the convergence
00:42:02.510 --> 00:42:03.520
tolerance is satisfied.
00:42:03.520 --> 00:42:06.090
If not, we have to choose
a new beta.
00:42:06.090 --> 00:42:13.590
And like that, we choose new
betas until we have satisfied
00:42:13.590 --> 00:42:15.360
this convergence tolerance.
00:42:15.360 --> 00:42:17.320
That's what we call
line search.
00:42:17.320 --> 00:42:21.700
We are searching along the
line given by delta U bar
00:42:21.700 --> 00:42:25.110
until this is satisfied.
00:42:25.110 --> 00:42:29.210
It's a very important scheme,
and this scheme can be
00:42:29.210 --> 00:42:32.590
combined with modified Newton
iteration, with the initial
00:42:32.590 --> 00:42:35.590
stress method, and with
full Newton iteration.
00:42:35.590 --> 00:42:40.680
And it adds a lot to the
stability of the solution, of
00:42:40.680 --> 00:42:43.660
the overall solution or the
nonlinear equations.
00:42:43.660 --> 00:42:47.470
We will find, or discuss
applications just now.
00:42:47.470 --> 00:42:52.890
But at this point, we need to
stop for a minute because we
00:42:52.890 --> 00:42:54.390
need to go onto a new tape.
00:42:58.160 --> 00:43:00.970
Finally, I would like to discuss
with you a method that
00:43:00.970 --> 00:43:03.540
we also find to be very
effective in nonlinear finite
00:43:03.540 --> 00:43:05.050
element computations.
00:43:05.050 --> 00:43:07.590
Namely, the BFGS method.
00:43:07.590 --> 00:43:08.835
BFGS stands for Broyden-Fletcher
00:43:08.835 --> 00:43:10.085
-Goldfarb-Shanno.
00:43:12.200 --> 00:43:15.930
And in this method, we
proceed as follows.
00:43:15.930 --> 00:43:20.850
We calculate or define
a delta i as shown
00:43:20.850 --> 00:43:23.070
here in this equation.
00:43:23.070 --> 00:43:28.560
And we define a gamma i as
shown in this equation.
00:43:28.560 --> 00:43:32.200
We want to calculate the
coefficient matrix such that
00:43:32.200 --> 00:43:34.110
this equation is
here satisfied.
00:43:34.110 --> 00:43:37.720
In other words, given delta i
and gamma i, we want to use
00:43:37.720 --> 00:43:41.040
the coefficient matrix that
satisfies this equation.
00:43:41.040 --> 00:43:45.680
And that's basically what's
being done in the BFGS scheme.
00:43:45.680 --> 00:43:50.610
Pictorially then, we can show
for a single degree of freedom
00:43:50.610 --> 00:43:55.260
system, say the F plotted
along this axis.
00:43:55.260 --> 00:43:58.730
The U displacements plotted
along this axis.
00:43:58.730 --> 00:44:01.770
Notice here in green
we show delta i.
00:44:01.770 --> 00:44:04.520
And we show here in
green gamma i.
00:44:04.520 --> 00:44:08.890
And notice that the matrix
that we want to use is
00:44:08.890 --> 00:44:11.760
indicated by this blue slope.
00:44:11.760 --> 00:44:16.560
Which is, in other words, given
by delta i and gamma i.
00:44:16.560 --> 00:44:21.300
Well, the BFGS method is an
iterative algorithm, which
00:44:21.300 --> 00:44:26.090
produces successive
approximations to efficient
00:44:26.090 --> 00:44:27.980
stiffness matrix.
00:44:27.980 --> 00:44:30.410
Actually, we actually work
with the inverse of that
00:44:30.410 --> 00:44:34.450
stiffness matrix as I will
elaborate upon just now.
00:44:34.450 --> 00:44:36.770
It's really a compromise
between the full Newton
00:44:36.770 --> 00:44:38.600
iteration method and
the modified
00:44:38.600 --> 00:44:40.590
Newton iteration method.
00:44:40.590 --> 00:44:43.350
It shows very excellent
convergence characteristics in
00:44:43.350 --> 00:44:45.950
the solution of many
types of problems.
00:44:45.950 --> 00:44:48.680
Let's look at the individual
steps that we use when we
00:44:48.680 --> 00:44:50.900
apply the BFGS method.
00:44:50.900 --> 00:44:54.530
We, first of all, calculate
the direction of the
00:44:54.530 --> 00:44:56.160
displacement increment.
00:44:56.160 --> 00:45:03.090
Here we have delta U bar i is
equal to K inverse, the
00:45:03.090 --> 00:45:06.120
inverse of the K matrix,
corresponding to
00:45:06.120 --> 00:45:08.050
iteration i minus 1.
00:45:08.050 --> 00:45:12.150
And here we have the out
of balance load vector.
00:45:12.150 --> 00:45:15.170
Notice once again we do not
calculate actually an inverse
00:45:15.170 --> 00:45:16.080
of a matrix.
00:45:16.080 --> 00:45:20.730
We calculate the LDL transpose
factorization of a matrix, and
00:45:20.730 --> 00:45:25.690
then we, as I will just now
show, update that inverse in
00:45:25.690 --> 00:45:31.445
the BFGS iteration by specific
matrices of rank 2.
00:45:31.445 --> 00:45:33.610
We get to that just now.
00:45:33.610 --> 00:45:39.360
Well, having calculated this
delta U bar i, we use that
00:45:39.360 --> 00:45:44.850
delta U bar i in the line search
to obtain a better
00:45:44.850 --> 00:45:49.650
displacement vector, a
displacement vector
00:45:49.650 --> 00:45:51.550
corresponding to iteration i.
00:45:51.550 --> 00:45:56.450
A better value then just adding
the delta U bar i to it
00:45:56.450 --> 00:45:57.910
by adjusting beta.
00:45:57.910 --> 00:46:02.960
In other words, we choose beta
in such way as to satisfy this
00:46:02.960 --> 00:46:04.250
equation here.
00:46:04.250 --> 00:46:06.060
This is, of course, the equation
that we are now
00:46:06.060 --> 00:46:11.680
familiar with as being equation
that shows that
00:46:11.680 --> 00:46:16.300
convergence has been reached
in a line search.
00:46:16.300 --> 00:46:22.650
Having now calculated the t
plus delta t Ui and the
00:46:22.650 --> 00:46:28.390
corresponding t plus delta t Fi,
we can calculate the delta
00:46:28.390 --> 00:46:30.660
i and gamma i.
00:46:30.660 --> 00:46:33.930
And therefore, apply
the BFGS method.
00:46:33.930 --> 00:46:36.630
And that's now step three,
calculation of
00:46:36.630 --> 00:46:38.760
the new secant matrix.
00:46:38.760 --> 00:46:46.420
K inverse of iteration i is
equal to Ai transposed k
00:46:46.420 --> 00:46:50.100
inverse of iteration
i minus 1 times Ai.
00:46:50.100 --> 00:46:55.393
This A matrix is obtained
as shown here.
00:46:55.393 --> 00:46:57.910
v wi transposed.
00:46:57.910 --> 00:47:04.390
The vi and wi's are given as a
function of the delta i, gamma
00:47:04.390 --> 00:47:10.450
i, and K i minus 1 value.
00:47:10.450 --> 00:47:15.120
You should look at the textbook
to see the details of
00:47:15.120 --> 00:47:18.590
actually evaluating
the vi and wi.
00:47:18.590 --> 00:47:23.430
But once you have vi, wi, you
can calculate A. And clearly,
00:47:23.430 --> 00:47:29.100
by this equation then, you get
a new stiffness matrix.
00:47:29.100 --> 00:47:33.670
It is important to note that in
this iterative scheme, only
00:47:33.670 --> 00:47:38.550
vector products are needed
to obtain vi and wi.
00:47:38.550 --> 00:47:41.680
And it's also then important
to note that only vector
00:47:41.680 --> 00:47:45.520
products are used to calculate
delta U bar i.
00:47:45.520 --> 00:47:50.040
If this would not be the case,
particularly here.
00:47:50.040 --> 00:47:51.540
Then of course, the
iterative scheme
00:47:51.540 --> 00:47:53.470
would be very expensive.
00:47:53.470 --> 00:47:59.960
Let me show you a bit of the
details why we only need to
00:47:59.960 --> 00:48:03.230
use vector products to
get the solution.
00:48:03.230 --> 00:48:09.230
Here, we have one typical
step of the iteration.
00:48:09.230 --> 00:48:13.730
Notice on the left-hand side,
we have delta U bar i.
00:48:13.730 --> 00:48:15.890
The value of displacement
increment
00:48:15.890 --> 00:48:17.880
that we want to calculate.
00:48:17.880 --> 00:48:22.460
On the right-hand side, we
have i plus wi minus 1 vi
00:48:22.460 --> 00:48:23.990
minus 1 transposed.
00:48:23.990 --> 00:48:28.890
In other words, the Ai
minus 1 transposed.
00:48:28.890 --> 00:48:33.860
And we would have such matrices
going on here until
00:48:33.860 --> 00:48:37.670
we are coming to the A1
transposed matrix.
00:48:37.670 --> 00:48:40.830
Then we have the inverse
of a stiffness matrix.
00:48:40.830 --> 00:48:43.860
Here, corresponding
to time tau.
00:48:43.860 --> 00:48:50.590
And then we have the Ai matrices
as shown here.
00:48:50.590 --> 00:48:57.590
Now, notice that this total
amount then is multiplied by
00:48:57.590 --> 00:49:02.690
the R minus F i minus 1, the
out of balance load vector.
00:49:02.690 --> 00:49:06.370
Let's now go through the
computations that are required
00:49:06.370 --> 00:49:10.500
to get delta U bar i.
00:49:10.500 --> 00:49:14.850
Here we have delta R. This delta
R, the out of balance
00:49:14.850 --> 00:49:17.030
load vector, which I
just call delta R,
00:49:17.030 --> 00:49:19.110
multiplies this matrix.
00:49:19.110 --> 00:49:22.010
But if you look at this
multiplication, we find that
00:49:22.010 --> 00:49:29.430
this delta R times w i minus 1
transposed is simply a number.
00:49:29.430 --> 00:49:31.460
So it's one vector
multiplication
00:49:31.460 --> 00:49:33.020
that gives us a number.
00:49:33.020 --> 00:49:39.110
This number multiplied by
v i minus 1 is just the
00:49:39.110 --> 00:49:41.100
multiplication of a number
times a vector.
00:49:41.100 --> 00:49:42.930
It's very simple.
00:49:42.930 --> 00:49:48.230
And notice, of course, this
delta R i minus 1 is also
00:49:48.230 --> 00:49:49.850
multiplying the identity
matrix.
00:49:49.850 --> 00:49:53.910
So in this step here, we have
really nothing else then a
00:49:53.910 --> 00:49:55.080
vector multiplication.
00:49:55.080 --> 00:49:58.040
The only vector multiplication
that, in fact, is there, is
00:49:58.040 --> 00:50:02.030
this value, this vector here,
times that vector there,
00:50:02.030 --> 00:50:03.510
transposed.
00:50:03.510 --> 00:50:04.610
That gives us a number.
00:50:04.610 --> 00:50:07.120
This number times a vector
[UNINTELLIGIBLE]
00:50:07.120 --> 00:50:08.040
operation.
00:50:08.040 --> 00:50:11.780
And of course, this one here
times the i matrix is just
00:50:11.780 --> 00:50:13.500
this vector by itself.
00:50:13.500 --> 00:50:16.690
So in this step, what
we obtain is
00:50:16.690 --> 00:50:19.950
simply another vector.
00:50:19.950 --> 00:50:24.240
We take this vector and
repeat the process.
00:50:24.240 --> 00:50:28.500
And since the repetition is
just what we discussed
00:50:28.500 --> 00:50:34.470
already, we find that we only
need vector multiplications to
00:50:34.470 --> 00:50:38.790
roll up all of these
multiplications up to here.
00:50:38.790 --> 00:50:43.045
Then we have a vector multiplied
by K inverse.
00:50:43.045 --> 00:50:47.380
Well, that involves only vector
multiplications again,
00:50:47.380 --> 00:50:49.810
because we have already
the LDL transposed
00:50:49.810 --> 00:50:51.460
factorization of tau k.
00:50:54.140 --> 00:50:57.080
These vector multiplications
result again, in a vector.
00:50:57.080 --> 00:51:00.880
And that vector multiplied by
these matrices here in the
00:51:00.880 --> 00:51:04.330
same way as we have proceeded
here, means again, only vector
00:51:04.330 --> 00:51:05.930
multiplications.
00:51:05.930 --> 00:51:10.440
So this very important step
here only requires vector
00:51:10.440 --> 00:51:11.140
multiplications.
00:51:11.140 --> 00:51:16.520
And that makes this whole
process really efficient.
00:51:16.520 --> 00:51:20.740
In summary then, we have the
following solution procedures
00:51:20.740 --> 00:51:23.130
that we feel are
very effective.
00:51:23.130 --> 00:51:27.020
The mortified Newton-Raphson
iteration with line searches.
00:51:27.020 --> 00:51:29.540
Here we have the basic modified
00:51:29.540 --> 00:51:31.190
Newton iteration equation.
00:51:31.190 --> 00:51:34.110
And here we have the line
search equation.
00:51:34.110 --> 00:51:38.930
We, of course, after each such
solution, evaluate an
00:51:38.930 --> 00:51:39.810
efficient beta.
00:51:39.810 --> 00:51:43.600
A beta that makes U--
00:51:43.600 --> 00:51:48.660
t plus delta U i a good vector
when measured on the line
00:51:48.660 --> 00:51:52.010
search criteria the way
I discussed it.
00:51:52.010 --> 00:51:56.210
And these two equations then,
summarize really the modified
00:51:56.210 --> 00:51:58.475
Newton iteration with
the line search.
00:51:58.475 --> 00:52:01.750
The BFGS method is another
effective scheme.
00:52:01.750 --> 00:52:05.210
We use the BFGS method always
with line searches.
00:52:05.210 --> 00:52:07.950
And then, as a third option,
the full Newton-Raphson
00:52:07.950 --> 00:52:10.740
iteration with or without
line searches.
00:52:10.740 --> 00:52:13.220
The full Newton-Raphson
iteration with line searches
00:52:13.220 --> 00:52:14.890
is, of course, most powerful.
00:52:14.890 --> 00:52:18.210
But it is also most expensive
per iteration.
00:52:18.210 --> 00:52:21.020
We should point out that these
techniques that I discussed so
00:52:21.020 --> 00:52:24.760
far cannot be applied for
post-buckling analysis.
00:52:24.760 --> 00:52:29.160
We will consider the solution of
the response after buckling
00:52:29.160 --> 00:52:32.520
has been occurring, after the
item load level has been
00:52:32.520 --> 00:52:36.020
reached in a forthcoming
lecture.
00:52:36.020 --> 00:52:41.610
The next view graphs summarize
these schemes: the modified
00:52:41.610 --> 00:52:44.450
Newton-Raphson, the BFGS
method, and the full
00:52:44.450 --> 00:52:46.840
Newton-Raphson method.
00:52:46.840 --> 00:52:50.030
And I'd like to just very
quickly go through the whole
00:52:50.030 --> 00:52:54.070
solution process for each
of these methods.
00:52:54.070 --> 00:52:57.700
In the modified Newton iteration
technique, we
00:52:57.700 --> 00:53:03.580
initialize our displacements
corresponding to the beginning
00:53:03.580 --> 00:53:09.610
of the first iteration, the end
of the 0 iteration as tU.
00:53:09.610 --> 00:53:12.920
And similarly, we initialize our
force vector corresponding
00:53:12.920 --> 00:53:14.470
to the internal element
stresses.
00:53:14.470 --> 00:53:19.010
We set i is equal to 1, the
iteration counter equal to 1.
00:53:19.010 --> 00:53:22.610
We calculate a K matrix at the
beginning of the load step,
00:53:22.610 --> 00:53:25.510
and we keep that K
matrix constant.
00:53:25.510 --> 00:53:29.450
We then go into this step here
where we calculate the out of
00:53:29.450 --> 00:53:30.700
balance loads.
00:53:33.210 --> 00:53:36.790
Calculate a new displacement
vector.
00:53:36.790 --> 00:53:41.240
And in this box here, we
measure whether we have
00:53:41.240 --> 00:53:42.160
converged already.
00:53:42.160 --> 00:53:45.390
Of course, this convergence
would only be measured after
00:53:45.390 --> 00:53:49.560
we have gone say, at
least twice through
00:53:49.560 --> 00:53:51.070
the iteration process.
00:53:51.070 --> 00:53:54.650
Or rather, we have calculated
this increment in displacement
00:53:54.650 --> 00:53:56.020
vector twice.
00:53:56.020 --> 00:53:59.560
Because at the beginning
corresponding to i equal to 1,
00:53:59.560 --> 00:54:03.400
we have now this fairly large
out of balance load.
00:54:03.400 --> 00:54:05.870
We have just incremented
the load vector.
00:54:05.870 --> 00:54:10.880
Therefore, this right inside
should be non-zero.
00:54:10.880 --> 00:54:13.580
And we would get an increment
in displacements.
00:54:13.580 --> 00:54:17.500
And to measure how large that
increment in displacement is,
00:54:17.500 --> 00:54:20.350
we really want to go through
this whole cycle at
00:54:20.350 --> 00:54:22.170
least one more time.
00:54:22.170 --> 00:54:25.610
We then, having calculated this
incremental displacement
00:54:25.610 --> 00:54:30.620
vector with the bar on there,
we perform a line search the
00:54:30.620 --> 00:54:34.230
way we discussed it to update
our displacements
00:54:34.230 --> 00:54:36.720
corresponding to the
iteration i.
00:54:39.960 --> 00:54:42.300
Of course, we only have come so
far once through it, so i
00:54:42.300 --> 00:54:43.520
is equal to 1 still.
00:54:43.520 --> 00:54:48.360
And this gives us the accepted
value of displacements
00:54:48.360 --> 00:54:50.936
corresponding to the
first iteration.
00:54:50.936 --> 00:54:55.560
We now increment our iteration
counter, and go in to the
00:54:55.560 --> 00:54:57.370
second iteration.
00:54:57.370 --> 00:55:01.060
Notice now in the second
iteration, we have R minus t
00:55:01.060 --> 00:55:03.670
plus delta t F1 here.
00:55:03.670 --> 00:55:06.250
We calculate delta U bar 2.
00:55:06.250 --> 00:55:10.040
And now we would certainly
measure convergence.
00:55:10.040 --> 00:55:13.360
And I will just now talk about
how we measure convergence.
00:55:13.360 --> 00:55:17.240
If we have not converged yet,
we keep on cycling through
00:55:17.240 --> 00:55:22.080
here until we do actually
get convergence.
00:55:22.080 --> 00:55:25.840
The distinguishing feature of
the modified Newton iteration
00:55:25.840 --> 00:55:29.642
is that we are having
a constant K matrix.
00:55:29.642 --> 00:55:32.900
That we do not update
the K matrix is this
00:55:32.900 --> 00:55:34.250
iterative loop here.
00:55:34.250 --> 00:55:36.280
It remains constant.
00:55:36.280 --> 00:55:40.570
Remember in the initial stress
technique, initial stress
00:55:40.570 --> 00:55:46.020
method, we would not even update
the K matrix ever after
00:55:46.020 --> 00:55:49.980
the first time it has
been calculated.
00:55:49.980 --> 00:55:52.990
But in the modified Newton
iteration, we update it at the
00:55:52.990 --> 00:55:55.800
beginning of each load step.
00:55:55.800 --> 00:56:00.140
In the BFGS method, we
proceed as follows.
00:56:00.140 --> 00:56:03.420
We initial, again, our
displacements and our nodal
00:56:03.420 --> 00:56:05.350
point force vector corresponding
to the internal
00:56:05.350 --> 00:56:06.570
element stresses.
00:56:06.570 --> 00:56:08.470
We calculate tK.
00:56:08.470 --> 00:56:11.630
We calculate our incremental
displacement vector with the
00:56:11.630 --> 00:56:13.120
bar on top.
00:56:13.120 --> 00:56:18.670
We measure convergence if i
is greater or equal to.
00:56:18.670 --> 00:56:20.700
There's no point as I just
mentioned, measuring
00:56:20.700 --> 00:56:24.070
convergence when you go first
time through here.
00:56:24.070 --> 00:56:29.020
And we perform then the line
search having calculated our
00:56:29.020 --> 00:56:32.440
actual incremental in nodal
point displacement that we
00:56:32.440 --> 00:56:35.410
want to add to the previous
nodal point displacements to
00:56:35.410 --> 00:56:39.150
get to this displacement
vector.
00:56:39.150 --> 00:56:42.040
We, of course, also have
this force vector.
00:56:42.040 --> 00:56:46.960
We now can update our inverse
matrix, our inverse secret
00:56:46.960 --> 00:56:50.740
matrix the way we
discussed it.
00:56:50.740 --> 00:56:53.120
And we increment our counter.
00:56:53.120 --> 00:56:58.110
And we keep on cycling through
here until we find in this box
00:56:58.110 --> 00:57:01.300
that we have converged.
00:57:01.300 --> 00:57:06.870
In the full Newton iteration
with line searches, we once
00:57:06.870 --> 00:57:10.460
again, initialize our
displacement vector, our force
00:57:10.460 --> 00:57:13.490
vector, our iteration counter.
00:57:13.490 --> 00:57:17.710
We now calculate a new K matrix
corresponding to the
00:57:17.710 --> 00:57:20.720
last conditions on displacements
and internal
00:57:20.720 --> 00:57:22.320
forces, of course.
00:57:22.320 --> 00:57:25.140
We enter here to calculate our
incremental displacement
00:57:25.140 --> 00:57:29.460
vector with the K matrix that
was just evaluated.
00:57:32.110 --> 00:57:36.740
And we measure convergence
once again only when i is
00:57:36.740 --> 00:57:38.460
greater equal to.
00:57:38.460 --> 00:57:45.450
We perform a line search to
determine beta in here.
00:57:45.450 --> 00:57:49.010
We therefore, update our
displacement vector to the
00:57:49.010 --> 00:57:51.570
displacements that now
correspond to the end of
00:57:51.570 --> 00:57:53.230
iteration i.
00:57:53.230 --> 00:57:55.250
We update our iteration
counter.
00:57:55.250 --> 00:57:57.720
Or rather, increase our
iteration counter.
00:57:57.720 --> 00:58:02.580
And we are calculating a new
stiffness matrix corresponding
00:58:02.580 --> 00:58:05.200
to the last calculated
displacement.
00:58:05.200 --> 00:58:08.180
This is, of course, a
distinguishing feature of the
00:58:08.180 --> 00:58:10.660
full Newton-Raphson method that
we are calculating a new
00:58:10.660 --> 00:58:13.700
K matrix at the beginning
of each iteration.
00:58:16.210 --> 00:58:21.360
Well, these view graphs then,
summarize the full process of
00:58:21.360 --> 00:58:24.770
iteration for the modified
Newton, the BFGS, and the full
00:58:24.770 --> 00:58:25.440
Newton method.
00:58:25.440 --> 00:58:28.660
Notice that I have included the
line search as the option.
00:58:28.660 --> 00:58:31.480
If you don't have line
searching, of course, it would
00:58:31.480 --> 00:58:37.290
simply mean that you set beta
equal to 1, and skip that
00:58:37.290 --> 00:58:39.680
particular step of
line searching.
00:58:39.680 --> 00:58:43.530
But in many cases, line
searching is quite important.
00:58:43.530 --> 00:58:47.440
And certainly, can
be necessary in
00:58:47.440 --> 00:58:49.790
some types of problems.
00:58:49.790 --> 00:58:52.690
Let us now look at what happens
in this box here.
00:58:52.690 --> 00:58:55.260
How do we measure convergence?
00:58:55.260 --> 00:58:58.850
In each of these iterative
schemes, we had such a box.
00:58:58.850 --> 00:59:02.730
And I like to look now more
closely at how we measure
00:59:02.730 --> 00:59:05.040
convergence.
00:59:05.040 --> 00:59:09.090
The measure of convergence
should, of course, tell us,
00:59:09.090 --> 00:59:11.740
how well do we satisfy
equilibrium?
00:59:11.740 --> 00:59:16.070
And there are basically three
items that can be used.
00:59:16.070 --> 00:59:20.870
Namely energy, force, or moment,
displacement, and
00:59:20.870 --> 00:59:23.110
rotation, of course,
correspondingly as well.
00:59:23.110 --> 00:59:25.940
But these are the items that one
can work with to measure
00:59:25.940 --> 00:59:27.280
convergence.
00:59:27.280 --> 00:59:31.510
If you use an energy convergence
criteria, and we
00:59:31.510 --> 00:59:33.800
are very fond of this one.
00:59:33.800 --> 00:59:40.970
You might want to use this
equation here delta U bar i
00:59:40.970 --> 00:59:46.800
transposed times the out of
balance load divided by delta
00:59:46.800 --> 00:59:50.870
U bar 1 transposed times t
plus delta t R minus tF.
00:59:50.870 --> 00:59:54.500
Now notice, this is the
increment in nodal point
00:59:54.500 --> 00:59:58.020
forces, or the out of
balance load, at the
00:59:58.020 --> 01:00:00.990
beginning of the load step.
01:00:00.990 --> 01:00:05.850
So this here is an energy
corresponding to the beginning
01:00:05.850 --> 01:00:07.090
of the load step.
01:00:07.090 --> 01:00:11.460
And we are comparing here
a similar quantity, but
01:00:11.460 --> 01:00:13.760
corresponding to the
current load.
01:00:13.760 --> 01:00:16.140
Corresponding to the
current iteration.
01:00:16.140 --> 01:00:21.190
And if this quotient here is
very small measured on that
01:00:21.190 --> 01:00:25.250
convergence tolerance, then
we say we have converged.
01:00:25.250 --> 01:00:28.910
An important point is that in
this convergence measure, we
01:00:28.910 --> 01:00:33.400
enter with displacements and
with out of balance loads.
01:00:33.400 --> 01:00:36.450
Both enter in this convergence
measure.
01:00:36.450 --> 01:00:41.270
And please realize also that
if the loads don't increase
01:00:41.270 --> 01:00:46.960
from time t to t plus delta t,
meaning that this here is
01:00:46.960 --> 01:00:49.330
equal to 0--
01:00:49.330 --> 01:00:50.410
0 vector.
01:00:50.410 --> 01:00:52.860
Then, of course, we
get 0 down here.
01:00:52.860 --> 01:00:55.020
And we could not apply this
convergence criteria
01:00:55.020 --> 01:00:56.170
indirectly.
01:00:56.170 --> 01:00:59.900
Therefore, we would do is
instead of 0 here, a
01:00:59.900 --> 01:01:03.890
reasonable small number in
the computer program.
01:01:03.890 --> 01:01:06.760
Also note that we apply
this test here
01:01:06.760 --> 01:01:09.300
prior to line searching.
01:01:09.300 --> 01:01:13.590
In line searching, after line
searching of course, this top
01:01:13.590 --> 01:01:15.580
part here is small.
01:01:15.580 --> 01:01:20.070
And we, therefore, do not want
to apply this scheme here or
01:01:20.070 --> 01:01:22.290
this convergence measure
after line searching.
01:01:22.290 --> 01:01:25.530
We should do it before
line searching.
01:01:25.530 --> 01:01:31.040
To measure convergence on
forces, we can use this
01:01:31.040 --> 01:01:32.420
equation here.
01:01:32.420 --> 01:01:35.620
Notice here we have the
out of balance loads.
01:01:35.620 --> 01:01:39.640
We take the Euclidean norm on
the out of balance loads, and
01:01:39.640 --> 01:01:43.710
we compare that with
RNORM, input by the
01:01:43.710 --> 01:01:45.970
user, by the analyst.
01:01:45.970 --> 01:01:50.560
It is interesting, important to
note that we introduce here
01:01:50.560 --> 01:01:53.860
only the components
corresponding to the
01:01:53.860 --> 01:01:57.220
translational degrees
of freedom.
01:01:57.220 --> 01:02:01.540
And here, of course, this
RNORM is also a force
01:02:01.540 --> 01:02:03.340
corresponding, so to say, to a
01:02:03.340 --> 01:02:04.840
translational degree of freedom.
01:02:04.840 --> 01:02:09.430
If we have also rotational
degrees of freedom, then we
01:02:09.430 --> 01:02:14.500
would take instead
of RNORM, RMNORM.
01:02:14.500 --> 01:02:18.760
And we would extract from these
vectors, the components
01:02:18.760 --> 01:02:20.130
corresponding to the rotational
01:02:20.130 --> 01:02:21.700
degrees of freedom only.
01:02:21.700 --> 01:02:24.530
The point is that we want to
have the same units up
01:02:24.530 --> 01:02:27.450
here as down here.
01:02:27.450 --> 01:02:31.790
Forces, forces, moments,
moments.
01:02:31.790 --> 01:02:36.440
And we don't want to mix
components from these
01:02:36.440 --> 01:02:38.280
quantities.
01:02:38.280 --> 01:02:43.160
Of course, RMNORM and RNORM must
be chosen by the analyst.
01:02:43.160 --> 01:02:48.190
The are input to the analysis as
defined in the input to the
01:02:48.190 --> 01:02:49.480
computer program.
01:02:49.480 --> 01:02:53.030
Typically, RTOL is 0.01.
01:02:53.030 --> 01:02:58.540
RNORM is simply the maximum
of the NORM of the
01:02:58.540 --> 01:03:00.920
loads that are applied.
01:03:00.920 --> 01:03:06.780
Notice this Euclidean norm is
evaluated as shown down here
01:03:06.780 --> 01:03:09.250
for a vector a.
01:03:09.250 --> 01:03:12.580
The symbolism, two bars each
side with a 2 means nothing
01:03:12.580 --> 01:03:16.440
else then taking the squares of
all the components of the
01:03:16.440 --> 01:03:19.630
vector, adding those
squares up.
01:03:19.630 --> 01:03:24.970
And then, taking the square
root out of that addition.
01:03:24.970 --> 01:03:29.200
On displacements, we can also
measure convergence.
01:03:29.200 --> 01:03:32.820
And here is the equation
that can be used.
01:03:32.820 --> 01:03:34.860
The Euclidean norm on
the incremental
01:03:34.860 --> 01:03:37.540
displacement vector.
01:03:37.540 --> 01:03:40.970
And that is compared
with DNORM.
01:03:40.970 --> 01:03:43.410
Again, DNORM.
01:03:43.410 --> 01:03:47.870
Notice that this here is again,
input by the analyst to
01:03:47.870 --> 01:03:49.900
the computer program.
01:03:49.900 --> 01:03:53.150
It is part of the input that is
provided for the analysis.
01:03:53.150 --> 01:03:55.280
And it's a reference
displacement.
01:03:55.280 --> 01:03:58.300
Notice then, of course, we
should only include here the
01:03:58.300 --> 01:04:01.670
displacements, translational
displacements in other words.
01:04:01.670 --> 01:04:04.370
If we have also, rotational
degrees of freedom.
01:04:04.370 --> 01:04:08.370
Then we would here include only
the rotational degrees of
01:04:08.370 --> 01:04:13.780
freedom components and
use here DMNORM.
01:04:13.780 --> 01:04:18.700
This way again, we compare one
quantity with certain units
01:04:18.700 --> 01:04:21.580
with a quantity that
has the same units.
01:04:21.580 --> 01:04:25.220
Once for translations and
once for rotations.
01:04:25.220 --> 01:04:28.440
And so when we have rotational
degrees of freedom in the
01:04:28.440 --> 01:04:32.230
analysis, we this way make
sure really, that the
01:04:32.230 --> 01:04:35.710
incremental displacements and
incremental rotations are
01:04:35.710 --> 01:04:38.140
small at convergence.
01:04:38.140 --> 01:04:41.270
This then completes what I
wanted to share with you in
01:04:41.270 --> 01:04:42.010
this lecture.
01:04:42.010 --> 01:04:44.060
Thank you very much for
your attention.