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:18.825
hundreds of MIT courses, visit
MIT OpenCourseWare at
00:00:18.825 --> 00:00:21.060
ocw.mit.edu.
00:00:21.060 --> 00:00:22.760
PROFESSOR: Ladies and gentlemen,
welcome to this
00:00:22.760 --> 00:00:25.585
lecture on nonlinear finite
element analysis of solids and
00:00:25.585 --> 00:00:26.410
structures.
00:00:26.410 --> 00:00:28.840
In this lecture, I'd like to
discuss with you some solution
00:00:28.840 --> 00:00:32.759
methods that we employ to
solve the finite element
00:00:32.759 --> 00:00:37.470
equations that we encounter in
nonlinear dynamic analysis.
00:00:37.470 --> 00:00:40.830
The solution methods that are,
in general used, can be
00:00:40.830 --> 00:00:44.360
thought of to fall into these
three categories; direct
00:00:44.360 --> 00:00:48.280
integration methods, mode
superposition, and
00:00:48.280 --> 00:00:49.690
substructuring.
00:00:49.690 --> 00:00:53.500
In this lecture, I would like to
talk about the explicit and
00:00:53.500 --> 00:00:57.500
implicit techniques, direct
integration methods, that we
00:00:57.500 --> 00:01:00.720
use to solve the finite
element equations.
00:01:00.720 --> 00:01:03.020
In the next lecture, we will
talk about the mode
00:01:03.020 --> 00:01:05.410
superposition and substructuring
techniques, and
00:01:05.410 --> 00:01:07.130
we will also consider
example solutions.
00:01:09.900 --> 00:01:13.840
When we want to solve the finite
element equations in
00:01:13.840 --> 00:01:18.080
nonlinear dynamic analysis, we
recognize that basically we
00:01:18.080 --> 00:01:22.120
have to operate on this
set of equations.
00:01:22.120 --> 00:01:27.050
Here we have the inertia forces,
the damping forces,
00:01:27.050 --> 00:01:32.520
the elastic forces, and the
externally applied loads.
00:01:32.520 --> 00:01:37.820
The inertia forces and damping
forces, well, you know them
00:01:37.820 --> 00:01:38.570
quite well.
00:01:38.570 --> 00:01:41.610
They're written in terms of
mass matrices, damping
00:01:41.610 --> 00:01:45.360
matrices, accelerations,
velocities.
00:01:45.360 --> 00:01:48.140
The elastic forces, here
"elastic" in quotes, are
00:01:48.140 --> 00:01:51.660
actually the nodal point forces
that are equivalent to
00:01:51.660 --> 00:01:53.790
the internal element stresses.
00:01:53.790 --> 00:01:56.780
And we talked, in the earlier
lectures, about how we
00:01:56.780 --> 00:02:00.000
evaluate them using, in
the updated Lagrangian
00:02:00.000 --> 00:02:01.540
formulation, Cauchy stresses.
00:02:01.540 --> 00:02:04.400
In the total Lagrangian
formulation, the second
00:02:04.400 --> 00:02:06.990
Piola-Kirchhoff stresses
et cetera.
00:02:06.990 --> 00:02:09.160
Please look up all the
information that we talked
00:02:09.160 --> 00:02:13.970
about in the earlier lectures,
if you want to study how we
00:02:13.970 --> 00:02:16.900
calculate these forces in large
displacement, large
00:02:16.900 --> 00:02:21.360
deformation, large
strain analysis.
00:02:21.360 --> 00:02:25.150
This set of equations should
hold at any time t during the
00:02:25.150 --> 00:02:26.770
response solution.
00:02:26.770 --> 00:02:31.010
And in direct integration
analysis, we consider this set
00:02:31.010 --> 00:02:37.400
of equations at the discrete
time points delta t apart, as
00:02:37.400 --> 00:02:39.630
shown down here.
00:02:39.630 --> 00:02:43.380
Of course, even in static
analysis, we always worked
00:02:43.380 --> 00:02:45.940
with the data t time
increment.
00:02:45.940 --> 00:02:48.980
But in static analysis, we did
not include the inertia
00:02:48.980 --> 00:02:54.040
forces, and the damping forces
in the analysis.
00:02:54.040 --> 00:02:57.700
The issues that we like to
discuss are, what are the
00:02:57.700 --> 00:03:00.770
basic procedures for obtaining
the solutions at
00:03:00.770 --> 00:03:03.170
the discrete times?
00:03:03.170 --> 00:03:06.570
And then once we have discussed
this question, we
00:03:06.570 --> 00:03:09.640
have to address the question of
which procedure should be
00:03:09.640 --> 00:03:11.920
used for a given problem.
00:03:11.920 --> 00:03:16.090
If we use, so to say a bad
procedure, "bad" in quotes for
00:03:16.090 --> 00:03:19.160
a given problem, the costs
can be very high
00:03:19.160 --> 00:03:20.870
to solve the problem.
00:03:20.870 --> 00:03:23.230
Therefore, it can be important,
can be very
00:03:23.230 --> 00:03:26.390
important to select the
appropriate technique for the
00:03:26.390 --> 00:03:29.000
given problem.
00:03:29.000 --> 00:03:31.970
Let's talk first about the
explicit time integration
00:03:31.970 --> 00:03:33.250
techniques.
00:03:33.250 --> 00:03:36.980
And the one that I like to
focus attention on is the
00:03:36.980 --> 00:03:38.150
central difference method.
00:03:38.150 --> 00:03:42.830
The central difference method
is very widely used to solve
00:03:42.830 --> 00:03:44.410
nonlinear problems.
00:03:44.410 --> 00:03:48.240
It's an explicit technique,
as we will just now see.
00:03:48.240 --> 00:03:51.160
The basic equations that
we are working
00:03:51.160 --> 00:03:53.560
with are listed here.
00:03:53.560 --> 00:03:57.860
This is the equilibrium
equation at time t.
00:03:57.860 --> 00:04:01.920
Notice here is a tf, the force
that corresponds to the
00:04:01.920 --> 00:04:06.490
internal element stresses
at time t.
00:04:06.490 --> 00:04:10.640
Here are the damping forces, and
here the inertia forces.
00:04:10.640 --> 00:04:15.760
We expand the velocity at
time t as shown on the
00:04:15.760 --> 00:04:17.450
right hand side here.
00:04:17.450 --> 00:04:21.200
And we expand the acceleration
at time t as shown on the
00:04:21.200 --> 00:04:23.230
right hand side here.
00:04:23.230 --> 00:04:27.050
These are the central
difference formulas.
00:04:27.050 --> 00:04:31.390
This method is mainly used for
wave propagation problems.
00:04:31.390 --> 00:04:34.150
It's an explicit technique,
and this is an important
00:04:34.150 --> 00:04:37.560
point, it's an explicit
technique because we are
00:04:37.560 --> 00:04:40.930
looking at the equilibrium
equation at time t to
00:04:40.930 --> 00:04:45.210
calculate the response for
time t plus delta t.
00:04:45.210 --> 00:04:50.640
Any method that does just that,
looks at the equilibrium
00:04:50.640 --> 00:04:53.980
equation at time t to obtain the
solution for the response
00:04:53.980 --> 00:04:57.790
at time t plus delta t is termed
an explicit method, an
00:04:57.790 --> 00:04:59.220
explicit integration method.
00:04:59.220 --> 00:05:01.410
So the central difference
method is
00:05:01.410 --> 00:05:04.290
such an explicit method.
00:05:04.290 --> 00:05:07.470
The important point to recognize
right here is that
00:05:07.470 --> 00:05:12.710
we don't set up a k matrix,
a stiffness matrix.
00:05:12.710 --> 00:05:16.680
And that is one important
fact that pertains
00:05:16.680 --> 00:05:19.680
to an explicit technique.
00:05:19.680 --> 00:05:23.480
Using these three equations,
we can
00:05:23.480 --> 00:05:26.850
directly develop one equation.
00:05:26.850 --> 00:05:29.800
And that equation looks
as shown up here.
00:05:29.800 --> 00:05:33.380
Notice we have three equations
and three unknowns
00:05:33.380 --> 00:05:34.390
which can be solved.
00:05:34.390 --> 00:05:37.050
The three unknowns being the
displacement at time t plus
00:05:37.050 --> 00:05:40.340
delta t, the velocity at time
t plus delta t, and the
00:05:40.340 --> 00:05:42.180
acceleration at time
t plus delta t.
00:05:42.180 --> 00:05:43.820
Three equations and
three unknowns.
00:05:43.820 --> 00:05:46.990
We can solve these equations
and the governing equation
00:05:46.990 --> 00:05:50.840
that we obtain in that solution
is listed right here.
00:05:50.840 --> 00:05:53.680
Notice m over delta t squared.
00:05:53.680 --> 00:05:57.740
Notice c over 2 delta t here.
00:05:57.740 --> 00:06:00.500
And then the unknown
displacement corresponding to
00:06:00.500 --> 00:06:04.000
time t plus delta t, and on the
right hand side, we have
00:06:04.000 --> 00:06:06.080
what we call an effective
load vector
00:06:06.080 --> 00:06:07.780
corresponding to time t.
00:06:07.780 --> 00:06:12.700
There is a hat here, and that
hat means that there are a
00:06:12.700 --> 00:06:15.170
large number of terms to
be taken into account.
00:06:15.170 --> 00:06:18.800
Namely, those corresponding to
the inertia and the damping in
00:06:18.800 --> 00:06:19.650
the system.
00:06:19.650 --> 00:06:25.740
Here we have tr hat is equal
to tr, the actual loads
00:06:25.740 --> 00:06:27.550
applied corresponding
to a time t.
00:06:27.550 --> 00:06:29.310
These are the externally
applied loads
00:06:29.310 --> 00:06:30.860
to the nodal points.
00:06:30.860 --> 00:06:34.050
tf is the force vector
corresponding to the internal
00:06:34.050 --> 00:06:36.560
element stresses, as
I said earlier.
00:06:36.560 --> 00:06:41.260
Here we have an inertia
contribution.
00:06:41.260 --> 00:06:45.140
Here another inertia
contribution and a damping
00:06:45.140 --> 00:06:47.240
contribution.
00:06:47.240 --> 00:06:50.170
But it is important to recognize
that the right hand
00:06:50.170 --> 00:06:52.040
side is known.
00:06:52.040 --> 00:06:56.270
In other words, all the terms
here are known because these
00:06:56.270 --> 00:06:59.080
terms only involve displacements
corresponding to
00:06:59.080 --> 00:07:02.750
time t minus delta t,
and displacements
00:07:02.750 --> 00:07:04.350
corresponding to time t.
00:07:04.350 --> 00:07:10.650
And therefore we know this tr
hat, we can plug it in here on
00:07:10.650 --> 00:07:13.300
the right hand side and
calculate the unknown
00:07:13.300 --> 00:07:17.090
displacements corresponding
to time t plus delta t.
00:07:17.090 --> 00:07:22.520
The method is most useful when
m and c, the mass and damping
00:07:22.520 --> 00:07:24.340
matrices are diagonal.
00:07:24.340 --> 00:07:29.540
Because in that case, these
equations here decouple as
00:07:29.540 --> 00:07:30.960
shown right here.
00:07:30.960 --> 00:07:34.790
Notice we can calculate the
individual displacements
00:07:34.790 --> 00:07:38.180
components, or the displacement
that's individual
00:07:38.180 --> 00:07:43.870
degrees of freedom, one after
the other as shown here.
00:07:43.870 --> 00:07:52.720
Once we have evaluated tr hat,
tr hat let's remember,
00:07:52.720 --> 00:07:55.620
contains these terms here.
00:07:55.620 --> 00:08:00.880
In other words, although m and
c are diagonal, there is some
00:08:00.880 --> 00:08:06.170
coupling from the equation j
into equation j minus 1 and
00:08:06.170 --> 00:08:10.730
into equation j plus 1 right
here on the right hand side,
00:08:10.730 --> 00:08:13.790
because the element
contributions overlap the
00:08:13.790 --> 00:08:15.880
degrees of freedom.
00:08:15.880 --> 00:08:21.780
And therefore, this here is
where the element coupling at
00:08:21.780 --> 00:08:22.700
the nodal points enters.
00:08:22.700 --> 00:08:28.210
But once we have calculated
this vector here, I should
00:08:28.210 --> 00:08:30.190
point to this vector here.
00:08:30.190 --> 00:08:33.340
Once we have calculated this
vector, we know the individual
00:08:33.340 --> 00:08:37.990
components, then we enter
here with that component
00:08:37.990 --> 00:08:42.710
corresponding to a degree of
freedom i, and directly we can
00:08:42.710 --> 00:08:46.050
calculate the nodal point
displacements corresponding to
00:08:46.050 --> 00:08:48.730
a degree of freedom i.
00:08:48.730 --> 00:08:52.270
That's an important point to
recognize, that there is no
00:08:52.270 --> 00:08:55.650
need really to set
up a k matrix, to
00:08:55.650 --> 00:08:57.835
trianglize a k matrix.
00:09:00.470 --> 00:09:05.310
We don't do that here, and that
means that per time step,
00:09:05.310 --> 00:09:07.610
this solution method can
be very effective.
00:09:10.450 --> 00:09:17.120
Note that in the solution if
the damping is 0, in other
00:09:17.120 --> 00:09:20.100
words, we don't have damping in
the system, cii is equal to
00:09:20.100 --> 00:09:24.370
0, we need that mii
is greater than 0.
00:09:24.370 --> 00:09:26.860
Otherwise we can't solve
the equations.
00:09:26.860 --> 00:09:32.440
Notice that tf is calculated
by summing all the element
00:09:32.440 --> 00:09:33.670
contributions.
00:09:33.670 --> 00:09:38.620
Summing over m means summing
over all the elements.
00:09:38.620 --> 00:09:42.500
And as I mentioned earlier, this
is where coupling comes
00:09:42.500 --> 00:09:48.870
into the solution because one
nodal point, one degree of
00:09:48.870 --> 00:09:52.980
freedom at that nodal point
carries contribution from all
00:09:52.980 --> 00:09:56.320
the elements that surround
that nodal point.
00:09:56.320 --> 00:10:01.770
To start the solution, we need
a special equation because,
00:10:01.770 --> 00:10:07.120
remember in the solution for
time t plus delta t, we need
00:10:07.120 --> 00:10:09.330
the information, the
displacements corresponding to
00:10:09.330 --> 00:10:13.350
time t, and corresponding
to time t minus delta t.
00:10:13.350 --> 00:10:17.840
So, if t is equal to 0, and we
want to march ahead to delta
00:10:17.840 --> 00:10:19.825
t, in the first step--
00:10:19.825 --> 00:10:23.740
we are in the first step--
we use this equation to
00:10:23.740 --> 00:10:27.610
calculate, to evaluate
minus delta tu.
00:10:27.610 --> 00:10:32.640
And that means we can then
directly start our explicit
00:10:32.640 --> 00:10:35.680
solution process.
00:10:35.680 --> 00:10:40.020
The central difference method is
only conditionally stable.
00:10:40.020 --> 00:10:43.680
The condition that has to be
satisfied on the time step is
00:10:43.680 --> 00:10:46.030
given in this equation.
00:10:46.030 --> 00:10:51.060
Delta t critical is a critical
time step, and delta t, the
00:10:51.060 --> 00:10:54.170
actual time step that is used
in the solution, has to be
00:10:54.170 --> 00:10:55.890
smaller than that value.
00:10:55.890 --> 00:10:57.830
Smaller or equal
to that value.
00:10:57.830 --> 00:11:01.690
Notice delta t critical is given
as tn over pi, where tn
00:11:01.690 --> 00:11:06.136
is the smallest period in the
finite element assemblage.
00:11:06.136 --> 00:11:09.710
In nonlinear analysis, we have
to recognize that tn changes
00:11:09.710 --> 00:11:11.190
during the time history.
00:11:11.190 --> 00:11:15.800
Response, tn would become
smaller if the system
00:11:15.800 --> 00:11:21.820
stiffens, due for example to
large displacement effects.
00:11:21.820 --> 00:11:25.900
And tn would become larger if
the system softens due to
00:11:25.900 --> 00:11:27.340
elasto-plasticity effects.
00:11:31.440 --> 00:11:35.590
One interesting and most
important point is that tn can
00:11:35.590 --> 00:11:37.060
be estimated.
00:11:37.060 --> 00:11:43.520
And it can be estimated, or we
can find a bound on tn which
00:11:43.520 --> 00:11:47.200
we can then use to calculate
the critical time step.
00:11:47.200 --> 00:11:49.920
And that is achieved
as follows.
00:11:49.920 --> 00:11:56.120
We notice that the omega n
squared, omega n being 2 pi
00:11:56.120 --> 00:12:01.895
over tn radians per second omega
n squared is smaller or
00:12:01.895 --> 00:12:10.640
equal to the maximum of the
omega m, nms squared over all
00:12:10.640 --> 00:12:13.140
elements m.
00:12:13.140 --> 00:12:14.350
Now what does this mean?
00:12:14.350 --> 00:12:20.310
It means here that we're looking
at omega nm squared of
00:12:20.310 --> 00:12:27.140
each of the elements, and we
select the largest such value
00:12:27.140 --> 00:12:30.260
to put in here.
00:12:30.260 --> 00:12:32.520
And that then gives us an upper
00:12:32.520 --> 00:12:36.970
bound on omega n squared.
00:12:36.970 --> 00:12:40.130
The actual omega n that
we need, because
00:12:40.130 --> 00:12:43.690
omega n gives us tn.
00:12:43.690 --> 00:12:47.960
So we want to calculate the the
largest frequency of all
00:12:47.960 --> 00:12:51.190
individual elements.
00:12:51.190 --> 00:12:54.530
And with that largest frequency,
we enter in this
00:12:54.530 --> 00:13:01.930
equation here, to get the bound
on tn, noting that in
00:13:01.930 --> 00:13:07.340
nonlinear analysis omega nm
changes, just the way I said
00:13:07.340 --> 00:13:09.300
earlier tn will also change.
00:13:09.300 --> 00:13:14.190
But all the omega nm changes
and tn correspondingly
00:13:14.190 --> 00:13:16.590
changes, this equation
is always satisfied.
00:13:19.740 --> 00:13:24.410
The time integration step that
we can actually use then, is
00:13:24.410 --> 00:13:26.610
given by as this equation.
00:13:26.610 --> 00:13:30.590
We can use delta t is equal to
2 over omega nm maximum.
00:13:30.590 --> 00:13:36.580
The maximum radiant per second
frequency that we have for any
00:13:36.580 --> 00:13:39.940
one of the elements in
the complete mesh.
00:13:39.940 --> 00:13:44.760
Because that expression is
smaller or equal than the
00:13:44.760 --> 00:13:46.520
critical time step.
00:13:46.520 --> 00:13:51.310
We could call this value here 2
over omega nm , the critical
00:13:51.310 --> 00:13:53.930
time step of element m.
00:13:53.930 --> 00:13:57.240
And therefore, what we're doing
here really is, that
00:13:57.240 --> 00:14:01.430
we're taking the smallest of the
critical time steps over
00:14:01.430 --> 00:14:07.480
all elements as that time
step in the solution.
00:14:07.480 --> 00:14:10.240
And that time step we
know is smaller or
00:14:10.240 --> 00:14:11.685
equal to delta t critical.
00:14:14.200 --> 00:14:19.470
Let's look at the proof of this
relationship, because
00:14:19.470 --> 00:14:23.500
it's an interesting proof and it
gives us a bit more insight
00:14:23.500 --> 00:14:29.240
into why this relationship
really holds.
00:14:29.240 --> 00:14:32.490
So I'd like to prove to you
very briefly that this
00:14:32.490 --> 00:14:34.570
relationship here is valid.
00:14:34.570 --> 00:14:39.130
And we do so by using the
Rayleigh quotient.
00:14:39.130 --> 00:14:44.620
Please look at the textbook
where the Rayleigh quotient is
00:14:44.620 --> 00:14:48.630
discussed, and I have to now
assume that you are familiar
00:14:48.630 --> 00:14:50.310
with the Rayleigh quotient.
00:14:50.310 --> 00:14:55.750
If you write down the Rayleigh
quotient for omega n squared,
00:14:55.750 --> 00:14:57.990
it is this relationship here.
00:14:57.990 --> 00:15:04.620
Omega n being the largest
frequency of the complete
00:15:04.620 --> 00:15:05.970
finite element mesh.
00:15:05.970 --> 00:15:10.020
Omega n is equal to 2
pi divided by tn.
00:15:10.020 --> 00:15:13.520
And tn was the smallest period
in the mesh that we are
00:15:13.520 --> 00:15:15.630
interested in knowing.
00:15:15.630 --> 00:15:19.770
Let us define this relationship
here, just the
00:15:19.770 --> 00:15:23.250
definition, and this
relationship here.
00:15:23.250 --> 00:15:25.800
Once again, just
the definition.
00:15:25.800 --> 00:15:31.140
Then, by simply substituting
from here and there into this
00:15:31.140 --> 00:15:34.890
relationship here, we surely
can write this equation as
00:15:34.890 --> 00:15:35.930
shown here.
00:15:35.930 --> 00:15:39.280
We're summing over all of the
us, and we're summing
00:15:39.280 --> 00:15:40.795
over all of the is.
00:15:40.795 --> 00:15:43.620
m going over all the elements.
00:15:43.620 --> 00:15:46.580
If you have a 1,000 elements in
the mesh, then n would go
00:15:46.580 --> 00:15:50.940
from 1 to 1,000.
00:15:50.940 --> 00:15:55.640
If we now look at the Rayleigh
quotient for a single element,
00:15:55.640 --> 00:16:00.420
we can say that rho m is given
by this relationship.
00:16:00.420 --> 00:16:05.700
Which is nothing else, using our
definition of um and im,
00:16:05.700 --> 00:16:10.310
is nothing else than this
relationship here.
00:16:10.310 --> 00:16:15.300
However, we can immediately
recognize that rho m, this
00:16:15.300 --> 00:16:19.950
value here, must be smaller or
equal to omega nm squared,
00:16:19.950 --> 00:16:26.440
where omega nm is the largest
frequency radian per second of
00:16:26.440 --> 00:16:28.920
the element m.
00:16:28.920 --> 00:16:34.520
This is here a relationship
that must hold because the
00:16:34.520 --> 00:16:39.110
Rayleigh quotient statement,
or the Rayleigh quotient
00:16:39.110 --> 00:16:42.230
theorem, tells us that
this must hold.
00:16:42.230 --> 00:16:46.890
And once again, please look up
in the textbook a general
00:16:46.890 --> 00:16:51.140
proof, there's a general proof
given that in fact this is a
00:16:51.140 --> 00:16:53.760
valid statement.
00:16:53.760 --> 00:16:59.480
Where now we can use this
statement, go back to this
00:16:59.480 --> 00:17:04.800
equation and immediately
identify that this must hold.
00:17:04.800 --> 00:17:09.690
And notice once again I'm
looking here at the element m.
00:17:09.690 --> 00:17:14.319
Of course, this relationship
must hold for all elements, m
00:17:14.319 --> 00:17:17.020
going from 1 to the total
number of elements you
00:17:17.020 --> 00:17:19.619
actually have in the mesh.
00:17:19.619 --> 00:17:24.579
This is the relationship we
will be using, and we
00:17:24.579 --> 00:17:29.390
substitute that relationship in
our earlier expression, for
00:17:29.390 --> 00:17:32.730
the omega n squared and directly
arrive at this
00:17:32.730 --> 00:17:35.190
relationship here.
00:17:35.190 --> 00:17:43.730
Now we look at this part here,
and recognize that we can take
00:17:43.730 --> 00:17:50.100
this out of the summation sign
by using the largest or the
00:17:50.100 --> 00:17:53.100
effect, we can take the effect
out of the summation sign, I
00:17:53.100 --> 00:17:55.610
should've said, by taking
the largest
00:17:55.610 --> 00:17:58.530
value of omega nm squared.
00:17:58.530 --> 00:18:02.450
And that's being done here and
now we can cancel these two
00:18:02.450 --> 00:18:06.950
terms resulting into this
relationship which
00:18:06.950 --> 00:18:09.790
completes our proof.
00:18:09.790 --> 00:18:13.420
It's an interesting proof, and
once again if you're not
00:18:13.420 --> 00:18:15.740
familiar with the Rayleigh
quotient, please look up the
00:18:15.740 --> 00:18:19.180
information in the textbook,
after which I think you can
00:18:19.180 --> 00:18:22.640
quite easily follow
this proof.
00:18:22.640 --> 00:18:26.070
The important point now is that
the largest frequencies
00:18:26.070 --> 00:18:31.920
of simple elements can be
estimated analytically.
00:18:31.920 --> 00:18:35.310
Sometimes we can calculate them
exactly, and here we have
00:18:35.310 --> 00:18:39.290
a case where we can calculate
the largest frequency of the
00:18:39.290 --> 00:18:40.720
element exactly.
00:18:40.720 --> 00:18:45.430
Here we have a two nodal truss
element, m here, m here, equal
00:18:45.430 --> 00:18:47.750
to rho AL over 2.
00:18:47.750 --> 00:18:50.310
L being the length of the
truss, A being the cross
00:18:50.310 --> 00:18:53.250
sectional area, rho being
the mass density.
00:18:53.250 --> 00:18:56.440
If we write down the eigenvalue
problem for this
00:18:56.440 --> 00:19:00.030
single element, we obtain
this relationship here.
00:19:00.030 --> 00:19:05.890
Stiffness matrix times the
eigenvector phi omega squared,
00:19:05.890 --> 00:19:10.680
which will be an eigenvalue, the
mass matrix phi being the
00:19:10.680 --> 00:19:12.360
sought eigenvector.
00:19:12.360 --> 00:19:18.100
And if we solve this eigenvalue
problem, we obtain
00:19:18.100 --> 00:19:24.260
directly as the eigenvalues
these two values here.
00:19:24.260 --> 00:19:26.180
There are only two of course,
because this is
00:19:26.180 --> 00:19:27.960
a two by two matrix.
00:19:27.960 --> 00:19:28.950
And so is that.
00:19:28.950 --> 00:19:30.590
So we have only two
eigenvalues.
00:19:30.590 --> 00:19:36.460
The first eigenvalue is 0, the
rigid body mode in the system.
00:19:36.460 --> 00:19:40.330
And the second eigenvalue is
given here, which is equal to
00:19:40.330 --> 00:19:42.750
omega n squared.
00:19:42.750 --> 00:19:45.950
And there is the relationship
of that second eigenvalue
00:19:45.950 --> 00:19:51.380
obtained by some very simple
calculations, and this is a
00:19:51.380 --> 00:19:56.210
rewriting of that equation,
of that relationship here.
00:19:56.210 --> 00:20:00.160
c is a wave speed
in the system.
00:20:00.160 --> 00:20:06.230
In fact, c squared is, as you
can see here, e over rho.
00:20:06.230 --> 00:20:09.860
It's this relationship that
we can now use, very
00:20:09.860 --> 00:20:15.440
conveniently, to get a time step
in a finite element mesh
00:20:15.440 --> 00:20:21.480
that will satisfy the critical
time step criterion.
00:20:21.480 --> 00:20:25.430
In other words, if we have a
finite element mesh that is
00:20:25.430 --> 00:20:30.730
consisting of an assemblage of
such truss elements, then to
00:20:30.730 --> 00:20:34.490
get the time step for the
explicit time integration,
00:20:34.490 --> 00:20:39.040
which will be a good time step,
which will be smaller or
00:20:39.040 --> 00:20:42.730
equal to the critical time step,
what we do is we simply
00:20:42.730 --> 00:20:46.880
calculate this relationship here
for each of the elements.
00:20:46.880 --> 00:20:52.410
And this relationship will give
us then, the proper time
00:20:52.410 --> 00:20:58.710
step by evaluating the smallest
time step that we can
00:20:58.710 --> 00:21:00.240
use for any one of
the elements.
00:21:00.240 --> 00:21:05.520
In other words, we simply take
2 over omega n, which is
00:21:05.520 --> 00:21:10.840
written here, and we recognize
that L over c, where L is the
00:21:10.840 --> 00:21:17.170
length of the element, is the
time step that we can use in
00:21:17.170 --> 00:21:18.360
the analysis.
00:21:18.360 --> 00:21:23.920
Notice L over c, we want to use
the smallest value of L
00:21:23.920 --> 00:21:28.150
over c coming from all
of the elements.
00:21:28.150 --> 00:21:30.840
L would be an element length.
00:21:30.840 --> 00:21:33.880
c would be the wave speed
through the element, and we
00:21:33.880 --> 00:21:38.560
want to use the smallest value
of L over c, looking at all
00:21:38.560 --> 00:21:44.390
the elements, as the time step
for our time integration.
00:21:44.390 --> 00:21:47.710
Notice here that L over c is
the time required for wave
00:21:47.710 --> 00:21:51.130
front to travel through
the element.
00:21:51.130 --> 00:21:54.310
That is a physical
interpretation of
00:21:54.310 --> 00:21:57.980
this L over c value.
00:21:57.980 --> 00:22:00.890
Let's look at some
modelling aspects
00:22:00.890 --> 00:22:03.090
regarding this technique.
00:22:03.090 --> 00:22:09.700
Let the applied wavelengths
be Lw as shown here.
00:22:09.700 --> 00:22:12.920
Here schematically
we show a wave.
00:22:12.920 --> 00:22:18.680
And let's see that how we would
model the finite element
00:22:18.680 --> 00:22:22.930
system appropriately when
we apply to that system
00:22:22.930 --> 00:22:25.130
this kind of wave.
00:22:25.130 --> 00:22:33.330
Well, we first of all find the
value tw, which is really the
00:22:33.330 --> 00:22:37.970
time required for the total
wave to propagate across a
00:22:37.970 --> 00:22:40.320
particular point in
the mesh, with the
00:22:40.320 --> 00:22:42.870
wave speed c of course.
00:22:42.870 --> 00:22:49.380
We then choose delta t via this
relationship here, where
00:22:49.380 --> 00:22:55.290
n is the number of time steps
used to represent the wave.
00:22:55.290 --> 00:23:02.250
With delta t now known, we
directly get Le, which is the
00:23:02.250 --> 00:23:04.670
element lengths.
00:23:04.670 --> 00:23:07.360
Well, I should not quite say
it's the element lengths, it's
00:23:07.360 --> 00:23:08.950
related to the element
lengths.
00:23:08.950 --> 00:23:16.150
For very simple problems, and
particular for our problems in
00:23:16.150 --> 00:23:19.070
a model using truss elements,
two nodal truss elements,
00:23:19.070 --> 00:23:21.420
actually Le would be the
element lengths.
00:23:21.420 --> 00:23:25.370
But if we have more complex
elements to model the finite
00:23:25.370 --> 00:23:29.070
elements solution, then this
here, this Le value, would be
00:23:29.070 --> 00:23:32.250
related to the element
lengths.
00:23:32.250 --> 00:23:35.770
Note that in this formula here,
n has to be chosen by
00:23:35.770 --> 00:23:39.280
the analyst to be a reasonable
value so as to be able to
00:23:39.280 --> 00:23:43.000
represent the total wave
appropriately.
00:23:43.000 --> 00:23:45.840
Here we have some
observations.
00:23:45.840 --> 00:23:50.230
Note that in 1D analysis, c,
the wave speed, is given as
00:23:50.230 --> 00:23:52.220
the square root out
of e over rho.
00:23:52.220 --> 00:23:54.200
I mentioned that already
earlier.
00:23:54.200 --> 00:23:57.170
Notice that in nonlinear
analysis, delta t must satisfy
00:23:57.170 --> 00:24:00.250
the stability limit throughout
the solution.
00:24:00.250 --> 00:24:04.700
This is a very important point
and in the next lecture we
00:24:04.700 --> 00:24:08.650
will look at some examples that
demonstrates what happens
00:24:08.650 --> 00:24:11.980
if you don't satisfy
the stability limit
00:24:11.980 --> 00:24:14.120
throughout the solution.
00:24:14.120 --> 00:24:17.445
It may also be effective to
change the time step during
00:24:17.445 --> 00:24:19.490
the analysis.
00:24:19.490 --> 00:24:21.210
Of course, this you
would do in a
00:24:21.210 --> 00:24:22.460
nonlinear analysis, typically.
00:24:25.860 --> 00:24:27.300
For the application.
00:24:27.300 --> 00:24:31.480
or for the use of the explicit
time integration method, when
00:24:31.480 --> 00:24:34.230
you use the explicit time
integration method, one
00:24:34.230 --> 00:24:36.620
generally uses lower-order
elements.
00:24:36.620 --> 00:24:38.300
They are usually preferable.
00:24:38.300 --> 00:24:41.010
And in this particular case,
you want to have
00:24:41.010 --> 00:24:44.250
Le here and Le there.
00:24:44.250 --> 00:24:47.710
Uniform meshing, uniform
meshing.
00:24:47.710 --> 00:24:51.070
For higher-order elements, you
could also use the explicit
00:24:51.070 --> 00:24:52.600
time integration.
00:24:52.600 --> 00:24:57.580
And there again, you want to
really have that this L is
00:24:57.580 --> 00:25:00.990
equal to that L star.
00:25:00.990 --> 00:25:07.390
If you don't have that
condition, then this L here
00:25:07.390 --> 00:25:10.920
will govern the actual
time step.
00:25:10.920 --> 00:25:14.400
The smallest length through
the element will
00:25:14.400 --> 00:25:15.840
govern the time step.
00:25:15.840 --> 00:25:21.890
And if you use this value here
for Le, with 8 there, you have
00:25:21.890 --> 00:25:24.890
a conservative value for Le.
00:25:24.890 --> 00:25:29.610
This 8 here enters into the
solution because you have this
00:25:29.610 --> 00:25:34.000
midnode, and the midnode carries
more stiffness then
00:25:34.000 --> 00:25:35.250
the corner nodes.
00:25:38.750 --> 00:25:43.050
Let us look further at
some observations
00:25:43.050 --> 00:25:44.440
regarding the method.
00:25:44.440 --> 00:25:48.140
In the analysis of this problem,
where we have a beam
00:25:48.140 --> 00:25:54.480
subjected to an end step load as
shown here, we can get the
00:25:54.480 --> 00:25:58.790
exact solution by simply
choosing a finite element mesh
00:25:58.790 --> 00:26:03.070
of truss elements that satisfy
this relationship here.
00:26:03.070 --> 00:26:07.420
In other words, we represent
this beam here, using two
00:26:07.420 --> 00:26:09.280
nodal truss elements.
00:26:09.280 --> 00:26:13.600
And these two nodal truss
elements have a length equal
00:26:13.600 --> 00:26:15.440
to c times delta t.
00:26:15.440 --> 00:26:18.670
Of course, we're looking here at
a linear elastic solution,
00:26:18.670 --> 00:26:20.870
meaning that c is constant.
00:26:20.870 --> 00:26:24.040
And we will get, this is the
important point, we will get
00:26:24.040 --> 00:26:27.190
the exact solution for
any number of truss
00:26:27.190 --> 00:26:29.270
elements that are used.
00:26:29.270 --> 00:26:33.450
In fact, we will look at this
problem later on in the next
00:26:33.450 --> 00:26:38.370
lecture in more detail, and we
will use this fact to obtain
00:26:38.370 --> 00:26:42.240
the exact solution, and we'll
study also what happens if we
00:26:42.240 --> 00:26:44.900
change the time step.
00:26:44.900 --> 00:26:49.460
So I'd like to refer you to the
next lecture in which I
00:26:49.460 --> 00:26:53.910
will be talking more
about this problem.
00:26:53.910 --> 00:26:58.110
Using the explicit time
integration method, it is very
00:26:58.110 --> 00:27:01.700
important to use uniform
meshing.
00:27:01.700 --> 00:27:04.370
And the reason for it
is given right here.
00:27:04.370 --> 00:27:09.860
Namely, we want to have then a
time step that is not unduly
00:27:09.860 --> 00:27:14.620
small in any region
of the total mesh.
00:27:14.620 --> 00:27:19.170
If you don't use uniform
meshing, then the time step is
00:27:19.170 --> 00:27:24.760
governed by a particular part
in the mesh, the time step
00:27:24.760 --> 00:27:26.610
size I should say,
is governed by a
00:27:26.610 --> 00:27:28.040
particular part in the mesh.
00:27:28.040 --> 00:27:33.340
Namely, those elements which
have the smallest Le, the way
00:27:33.340 --> 00:27:34.830
I talked about it earlier.
00:27:34.830 --> 00:27:40.020
And that same time step would
be used for the whole mesh,
00:27:40.020 --> 00:27:43.370
which means then, that time
step, that small time step
00:27:43.370 --> 00:27:46.150
would actually be too
small for the other
00:27:46.150 --> 00:27:48.510
parts of the mesh.
00:27:48.510 --> 00:27:50.780
Of course, one could say, well,
why not use different
00:27:50.780 --> 00:27:53.010
time steps in the mesh?
00:27:53.010 --> 00:28:00.170
Namely, effective time steps for
the elements that are used
00:28:00.170 --> 00:28:01.940
in the various regions
of the mesh.
00:28:01.940 --> 00:28:05.410
That can be done, but then you
have to use special coupling
00:28:05.410 --> 00:28:11.040
techniques to couple the time
step integrations where you
00:28:11.040 --> 00:28:15.110
use in one mesh a small time
step, in one part of the mesh
00:28:15.110 --> 00:28:17.820
a small time step, in the other
products the mesh a
00:28:17.820 --> 00:28:18.880
larger time step.
00:28:18.880 --> 00:28:22.850
It can be done, but it needs
special considerations
00:28:22.850 --> 00:28:26.320
regarding the coupling of the
time integrations in the
00:28:26.320 --> 00:28:28.800
different parts of the mesh.
00:28:28.800 --> 00:28:32.540
We may also want to note that
a system with a very large
00:28:32.540 --> 00:28:36.740
bandwidth can be solved
frequently, quite efficiently,
00:28:36.740 --> 00:28:39.690
with the central difference
method, although the problem
00:28:39.690 --> 00:28:41.680
may not be a wave propagation
problem.
00:28:41.680 --> 00:28:47.940
The reason for this fact is
that in the explicit time
00:28:47.940 --> 00:28:53.710
integration, you do not set up a
k matrix, and therefore, the
00:28:53.710 --> 00:28:58.240
large bandwidth that you're
seeing in the k matrix is not
00:28:58.240 --> 00:29:00.510
really a hindrance.
00:29:00.510 --> 00:29:06.610
It is not being used because the
k matrix is not being set
00:29:06.610 --> 00:29:10.300
up, and you have no
factorization or LDL transpose
00:29:10.300 --> 00:29:12.550
factorization of the
coefficient matrix.
00:29:12.550 --> 00:29:16.100
And therefore, this high cost
that lies in the factorization
00:29:16.100 --> 00:29:21.280
of the coefficient matrix is not
necessary, is not used, is
00:29:21.280 --> 00:29:24.870
not expensed in the explicit
time integration.
00:29:24.870 --> 00:29:29.110
And for that reason, it can be
effective to use explicit time
00:29:29.110 --> 00:29:34.056
integration when there would
be such a large bandwidth.
00:29:34.056 --> 00:29:38.090
Explicit time integration lends
itself also very well to
00:29:38.090 --> 00:29:39.850
parallel processing.
00:29:39.850 --> 00:29:43.000
Notice that the basic equations
that we are solving
00:29:43.000 --> 00:29:47.546
for in explicit time
integration, in each time step
00:29:47.546 --> 00:29:49.740
are given as here.
00:29:49.740 --> 00:29:54.080
T plus delta tu is equal to
something that I denote here
00:29:54.080 --> 00:29:56.740
by the three dots,
times tr hat.
00:29:56.740 --> 00:30:00.390
And this tr hat vector
we defined earlier.
00:30:00.390 --> 00:30:05.120
So all we need to do is
calculate these entries in the
00:30:05.120 --> 00:30:08.570
vector here to get our unknown
displacements.
00:30:08.570 --> 00:30:11.890
Well, that can be done quite
effectively in parallel.
00:30:11.890 --> 00:30:16.560
We might want to calculate these
entries here at the same
00:30:16.560 --> 00:30:20.100
time as we calculate
say, these entries.
00:30:20.100 --> 00:30:23.130
Where this is here a particular
element group, and
00:30:23.130 --> 00:30:25.920
that are the entries
corresponding to an adjacent
00:30:25.920 --> 00:30:26.810
element group.
00:30:26.810 --> 00:30:29.260
And these are the entries
corresponding to another
00:30:29.260 --> 00:30:30.050
element group.
00:30:30.050 --> 00:30:34.460
So we can, in other words, in
parallel, process all the
00:30:34.460 --> 00:30:38.440
information that goes in here,
that goes in here, and that
00:30:38.440 --> 00:30:40.580
goes in here, and so on.
00:30:40.580 --> 00:30:44.560
Notice there's some coupling
in that information, which
00:30:44.560 --> 00:30:48.470
goes across this red
horizontal line.
00:30:48.470 --> 00:30:53.280
That coupling comes from the fm
contributions, namely the
00:30:53.280 --> 00:30:58.360
element contributions that
will be elements that lie
00:30:58.360 --> 00:31:01.220
across these equations,
so to say.
00:31:01.220 --> 00:31:04.010
And that's where a coupling
lies, and that has to be
00:31:04.010 --> 00:31:07.720
properly taken into account in
the parallel processing.
00:31:07.720 --> 00:31:12.120
But in principle, explicit time
integration because the
00:31:12.120 --> 00:31:17.110
equation that we are dealing
with, as shown here, lends
00:31:17.110 --> 00:31:22.450
itself very well to parallel
processing on computers.
00:31:22.450 --> 00:31:27.550
Let us now look at the other
technique that I wanted to
00:31:27.550 --> 00:31:30.380
discuss with you briefly
in this lecture.
00:31:30.380 --> 00:31:34.890
And this is an implicit time
integration method.
00:31:34.890 --> 00:31:38.780
It is a trapezoidal rule, which
is very widely used in
00:31:38.780 --> 00:31:41.090
implicit time integration.
00:31:41.090 --> 00:31:46.470
But let's first look at the
basic equation that we now are
00:31:46.470 --> 00:31:47.500
operating on.
00:31:47.500 --> 00:31:50.730
In implicit time integration,
in any implicit time
00:31:50.730 --> 00:31:55.340
integration scheme, we are
solving this equation here.
00:31:58.020 --> 00:32:03.690
At time t plus delta t to obtain
the solution at time t
00:32:03.690 --> 00:32:05.000
plus delta t.
00:32:05.000 --> 00:32:07.440
Notice this is the equilibrium
equation at time
00:32:07.440 --> 00:32:08.780
t plus delta t.
00:32:08.780 --> 00:32:11.470
This is quite different from
what we're doing in explicit
00:32:11.470 --> 00:32:12.300
time integration.
00:32:12.300 --> 00:32:15.050
In explicit time integration
we are looking at the
00:32:15.050 --> 00:32:18.940
equilibrium equation at time t
to obtain the solution at time
00:32:18.940 --> 00:32:20.340
t plus delta t.
00:32:20.340 --> 00:32:23.680
Here we are looking at the
equilibrium equation, we are
00:32:23.680 --> 00:32:27.630
using the equilibrium equation
at time t plus delta t to
00:32:27.630 --> 00:32:30.040
obtain the solution at
time t plus delta t.
00:32:30.040 --> 00:32:36.470
And that means that we have to
deal with a stiffness matrix.
00:32:36.470 --> 00:32:40.300
And that stiffness matrix here
is given as tk, corresponding
00:32:40.300 --> 00:32:43.230
to a modified Newton
iteration.
00:32:43.230 --> 00:32:47.590
Notice that this equation will
be solved to calculate an
00:32:47.590 --> 00:32:49.200
increment and displacement.
00:32:49.200 --> 00:32:52.280
And that increment in
displacement is added to the
00:32:52.280 --> 00:32:54.260
displacement vector that we
had from the previous
00:32:54.260 --> 00:32:59.950
iteration to obtain an updated
displacement vector.
00:32:59.950 --> 00:33:03.670
This is one equation in the
unknown displacements,
00:33:03.670 --> 00:33:07.950
velocities, and accelerations
at time t plus delta t.
00:33:07.950 --> 00:33:13.300
And we need two more
equations to solve
00:33:13.300 --> 00:33:15.200
for these three unknowns.
00:33:15.200 --> 00:33:19.750
In the trapezoidal rule, which
as I mentioned just now, is
00:33:19.750 --> 00:33:23.780
very widely used for implicit
time integration, these are
00:33:23.780 --> 00:33:31.020
the basic equations used for
displacements and for
00:33:31.020 --> 00:33:32.750
velocities.
00:33:32.750 --> 00:33:39.065
If we rewrite these equations,
we directly obtain the
00:33:39.065 --> 00:33:43.480
velocities and accelerations
in this form.
00:33:43.480 --> 00:33:46.710
Notice there is here the
increment in displacement from
00:33:46.710 --> 00:33:49.320
time t to time t plus delta t.
00:33:49.320 --> 00:33:52.370
Similarly here.
00:33:52.370 --> 00:33:55.990
In nonlinear analysis, we
are iterating for this
00:33:55.990 --> 00:34:01.420
displacement here, and therefore
we rewrite these
00:34:01.420 --> 00:34:03.240
equations a bit.
00:34:03.240 --> 00:34:05.790
And that's being shown,
that rewriting is
00:34:05.790 --> 00:34:07.870
shown on this viewgraph.
00:34:07.870 --> 00:34:11.150
Notice that we now have
the superscript k here
00:34:11.150 --> 00:34:15.679
corresponding to the iteration
k, and that we have split the
00:34:15.679 --> 00:34:19.870
displacement at time t plus
delta t into a value that we
00:34:19.870 --> 00:34:23.760
know, corresponding to iteration
k minus 1, and an
00:34:23.760 --> 00:34:28.250
increment that is unknown,
delta uk.
00:34:28.250 --> 00:34:30.739
We do the same for the
acceleration here.
00:34:33.400 --> 00:34:36.480
Of course, these terms
are as shown on
00:34:36.480 --> 00:34:38.760
the previous viewgraph.
00:34:38.760 --> 00:34:44.030
And now we have our velocity and
acceleration corresponding
00:34:44.030 --> 00:34:52.130
to iteration k, in terms of
known values and one unknown
00:34:52.130 --> 00:34:56.239
value, the incremental
displacement.
00:34:56.239 --> 00:35:00.030
We can now use these two
equations and substitute into
00:35:00.030 --> 00:35:03.950
our equilibrium equation the
equation that you saw on the
00:35:03.950 --> 00:35:07.640
first viewgraph, when we started
to talk about the
00:35:07.640 --> 00:35:09.710
implicit time integration
scheme.
00:35:09.710 --> 00:35:13.110
And the result of that
substitution is
00:35:13.110 --> 00:35:14.510
given on this viewgraph.
00:35:14.510 --> 00:35:20.010
We have a coefficient matrix
here, which we call tk hat
00:35:20.010 --> 00:35:21.590
times the incremental
displacement
00:35:21.590 --> 00:35:24.110
vector, which is unknown.
00:35:24.110 --> 00:35:27.260
And on the right hand side of
that equation, we have the
00:35:27.260 --> 00:35:31.900
externally applied loads,
the nodal point forces
00:35:31.900 --> 00:35:34.690
corresponding to the internal
element stresses at time t
00:35:34.690 --> 00:35:40.160
plus delta t, and at the end
of iteration k minus one.
00:35:40.160 --> 00:35:43.180
To evaluate these, we need the
displacements corresponding to
00:35:43.180 --> 00:35:45.570
iteration k minus 1.
00:35:45.570 --> 00:35:50.240
And here we have inertia
contributions and damping
00:35:50.240 --> 00:35:51.330
contributions.
00:35:51.330 --> 00:35:55.290
But notice all this information
here is known once
00:35:55.290 --> 00:35:58.110
you know the displacements
corresponding to t plus delta
00:35:58.110 --> 00:36:02.300
tu k minus 1.
00:36:02.300 --> 00:36:07.150
We iterate on this equation
until we satisfy certain
00:36:07.150 --> 00:36:12.440
equilibrium criteria,
convergence criteria, and much
00:36:12.440 --> 00:36:15.650
the same way as we do it in
static nonlinear analysis, and
00:36:15.650 --> 00:36:18.550
the way we discussed it
in an earlier lecture.
00:36:18.550 --> 00:36:21.490
However, let us make
a few observations.
00:36:21.490 --> 00:36:24.890
First of all, we notice that as
delta t gets smaller, the
00:36:24.890 --> 00:36:29.410
entries in this tk hat matrix,
this effective stiffness
00:36:29.410 --> 00:36:31.030
matrix, increase.
00:36:31.030 --> 00:36:32.040
Why do they increase?
00:36:32.040 --> 00:36:38.570
Because we have delta t squared
term, 1 over delta t
00:36:38.570 --> 00:36:42.490
squared term in front
of the mass matrix.
00:36:42.490 --> 00:36:45.430
And the 1 over 2 delta
t term in front
00:36:45.430 --> 00:36:46.720
of the damping matrix.
00:36:46.720 --> 00:36:54.440
As delta t becomes smaller, 1
over delta t squared becomes
00:36:54.440 --> 00:36:57.830
larger and that means that the
entries in this matrix here
00:36:57.830 --> 00:37:01.660
become larger as delta
t gets smaller.
00:37:01.660 --> 00:37:04.230
The convergence characteristics
of the
00:37:04.230 --> 00:37:06.700
equilibrium iterations
are better
00:37:06.700 --> 00:37:08.820
than in static analysis.
00:37:08.820 --> 00:37:13.200
The reason being, that the
entries in tk hat are
00:37:13.200 --> 00:37:17.060
generally larger, much larger
than what we see in static
00:37:17.060 --> 00:37:20.630
analysis because of this 1 over
delta t squared term, and
00:37:20.630 --> 00:37:24.990
1 over 2 delta t term, that
I just mentioned.
00:37:24.990 --> 00:37:27.350
We also notice that the
trapezoidal rule is
00:37:27.350 --> 00:37:31.850
unconditionally stable
in linear analysis.
00:37:31.850 --> 00:37:35.840
The analysis that shows that the
method is unconditionally
00:37:35.840 --> 00:37:42.450
stable is really performed on a
linear system, and therefore
00:37:42.450 --> 00:37:44.450
the unconditional stability
strictly
00:37:44.450 --> 00:37:46.190
holds only linear analysis.
00:37:46.190 --> 00:37:49.680
For nonlinear analysis, we want
to select delta t for
00:37:49.680 --> 00:37:53.575
accuracy, and we want to select
delta t for convergence
00:37:53.575 --> 00:37:55.630
of the iteration.
00:37:55.630 --> 00:37:59.910
If we find, for example, that
we have convergence
00:37:59.910 --> 00:38:08.260
difficulties in the iteration,
then we may want to switch to
00:38:08.260 --> 00:38:09.750
a more powerful scheme.
00:38:09.750 --> 00:38:12.260
That's one approach, of course,
just the same way as
00:38:12.260 --> 00:38:15.420
we are doing it in
static analysis.
00:38:15.420 --> 00:38:22.690
Or we may also recall, or we
may identify actually, that
00:38:22.690 --> 00:38:27.680
delta t is too large, and we
should decrease delta t.
00:38:27.680 --> 00:38:31.680
In particular, as a rule of
thumb, I'd just like to share
00:38:31.680 --> 00:38:33.010
some experience with you.
00:38:33.010 --> 00:38:36.580
We find that if you use the
BFGS method, which we
00:38:36.580 --> 00:38:40.000
discussed earlier in an earlier
lecture, if you use
00:38:40.000 --> 00:38:44.010
the BFGS method for the
equilibrium iteration, and if
00:38:44.010 --> 00:38:47.190
we have difficulties converging
there, then very
00:38:47.190 --> 00:38:50.830
frequently the time
step is too large,
00:38:50.830 --> 00:38:53.190
and should be decreased.
00:38:53.190 --> 00:38:57.870
Not just to obtain convergence
in the BFGS iterations, but to
00:38:57.870 --> 00:39:00.200
obtain an appropriate
accuracy.
00:39:00.200 --> 00:39:05.060
So we really look at the
difficulties obtaining
00:39:05.060 --> 00:39:09.440
convergence in the iterations as
also a signal to us to tell
00:39:09.440 --> 00:39:13.630
us that our time step delta
t is probably too large.
00:39:16.430 --> 00:39:21.980
The convergence criteria that we
are using are basically the
00:39:21.980 --> 00:39:25.730
same that we have seen already
in the static analysis in an
00:39:25.730 --> 00:39:27.060
earlier lecture.
00:39:27.060 --> 00:39:32.440
In static analysis we only
included this term here.
00:39:32.440 --> 00:39:37.660
In dynamic analysis we include
now the mass and damping terms
00:39:37.660 --> 00:39:40.230
in our energy criterion.
00:39:40.230 --> 00:39:45.650
We discussed this criterion
earlier, of course not
00:39:45.650 --> 00:39:48.200
including the mass and
damping effects.
00:39:48.200 --> 00:39:50.310
Now we do.
00:39:50.310 --> 00:39:57.490
Another convergence criterion is
one on forces, and earlier
00:39:57.490 --> 00:40:01.180
we had only this term here.
00:40:01.180 --> 00:40:05.460
Now we include also the inertia
and the damping terms.
00:40:07.990 --> 00:40:13.680
Once again, if we are only
having translational degrees
00:40:13.680 --> 00:40:19.900
of freedom, then here we have
only nodal point forces, nodal
00:40:19.900 --> 00:40:21.960
point forces here.
00:40:21.960 --> 00:40:26.380
And our norm would be having
the dimensions of a nodal
00:40:26.380 --> 00:40:27.990
point force.
00:40:27.990 --> 00:40:32.770
If we have rotational degrees
of freedom, then we still
00:40:32.770 --> 00:40:36.290
include here only the
nodal point forces.
00:40:36.290 --> 00:40:40.510
We extract, in other words, the
components corresponding
00:40:40.510 --> 00:40:44.230
to the translational degrees of
freedom of this vector and
00:40:44.230 --> 00:40:49.050
that vector, and include those
components only up here.
00:40:49.050 --> 00:40:53.650
We use RNORM that also has the
units of the nodal point
00:40:53.650 --> 00:40:59.430
force, and we supplement this
convergence criteria by
00:40:59.430 --> 00:41:04.580
another quotient, where we have
on top the components
00:41:04.580 --> 00:41:06.720
corresponding to the rotational
degrees of freedom
00:41:06.720 --> 00:41:10.030
only, and in the bottom,
what we called
00:41:10.030 --> 00:41:13.000
RMNORM already earlier.
00:41:13.000 --> 00:41:17.140
The point is that we want to
have up here only those
00:41:17.140 --> 00:41:23.070
components that carry the same
units as we have down here.
00:41:23.070 --> 00:41:28.080
Once for translations,
once for rotations.
00:41:28.080 --> 00:41:31.480
And of course recall we're
talking here about the
00:41:31.480 --> 00:41:32.190
Euclidean norm.
00:41:32.190 --> 00:41:34.180
which is defined
as shown here.
00:41:34.180 --> 00:41:37.440
We went through this already
earlier, when we talked about
00:41:37.440 --> 00:41:40.740
the convergence criteria
in static analysis.
00:41:40.740 --> 00:41:44.770
We have also the convergence
criterion on displacements.
00:41:44.770 --> 00:41:47.290
Notice this is what
we saw already
00:41:47.290 --> 00:41:49.050
earlier in static analysis.
00:41:49.050 --> 00:41:53.690
Once again, we use displacement
components here
00:41:53.690 --> 00:41:55.940
to compare with the
displacement.
00:41:55.940 --> 00:41:58.680
We introduced all of the
displacement components here
00:41:58.680 --> 00:42:02.370
to compare with this
displacement and if we have
00:42:02.370 --> 00:42:04.610
rotational degrees of
freedom, they would
00:42:04.610 --> 00:42:07.060
not be included here.
00:42:07.060 --> 00:42:09.670
Their components would
not be included here.
00:42:09.670 --> 00:42:13.290
When you use DNORM here, they
would, however, in a second
00:42:13.290 --> 00:42:17.960
check be included here and only
included here, and then
00:42:17.960 --> 00:42:19.710
we would use here DMNORM.
00:42:19.710 --> 00:42:24.120
We went through that
earlier already.
00:42:24.120 --> 00:42:30.180
Some modelling hints regarding
the use of the implicit time
00:42:30.180 --> 00:42:32.310
integration scheme.
00:42:32.310 --> 00:42:36.830
The integration scheme is
generally used to analyze
00:42:36.830 --> 00:42:40.120
structural vibration
problems, and we
00:42:40.120 --> 00:42:42.440
would proceed as follows.
00:42:42.440 --> 00:42:47.540
We identify the frequencies
contained in the loading by a
00:42:47.540 --> 00:42:49.770
Fourier decomposition.
00:42:49.770 --> 00:42:53.340
We chose a finite element mesh
that can accurately represent
00:42:53.340 --> 00:42:57.600
the static response and all
important frequencies.
00:42:57.600 --> 00:43:01.550
And then we perform a direct
integration using this time
00:43:01.550 --> 00:43:08.780
step here where tco is the
smallest period in seconds
00:43:08.780 --> 00:43:11.110
that we want to integrate.
00:43:11.110 --> 00:43:15.540
Now there's a lot of information
here, and I'd like
00:43:15.540 --> 00:43:20.520
to refer you to an example,
analysis of a tower that I
00:43:20.520 --> 00:43:23.330
will be presenting to you in
the next lecture, where we
00:43:23.330 --> 00:43:28.200
will talk about these items
quite a bit more.
00:43:28.200 --> 00:43:32.700
But this is basically the
approach that we use, to model
00:43:32.700 --> 00:43:37.520
the problem and then also to
perform the direct integration
00:43:37.520 --> 00:43:41.480
of the governing finite
element equations.
00:43:41.480 --> 00:43:45.120
As I said earlier already, the
method is used primarily for
00:43:45.120 --> 00:43:47.490
structural vibration problems.
00:43:47.490 --> 00:43:51.680
It is effective, quite
effective, using also
00:43:51.680 --> 00:43:53.930
higher-order elements.
00:43:53.930 --> 00:43:58.970
And if we do use higher-order
elements, we do find
00:43:58.970 --> 00:44:03.860
frequently that method, or the
total solution should be
00:44:03.860 --> 00:44:06.470
performed most effectively
using a
00:44:06.470 --> 00:44:08.680
consistent mass matrix.
00:44:08.680 --> 00:44:11.180
I'm thinking here in particular,
if you use very
00:44:11.180 --> 00:44:15.030
high-order elements such as
cubic shall elements,
00:44:15.030 --> 00:44:18.420
isoparametric cubic shell
elements, an element that we
00:44:18.420 --> 00:44:22.150
will talk about in a later
lecture, then with such
00:44:22.150 --> 00:44:25.720
element you probably want to use
a consistent mass matrix
00:44:25.720 --> 00:44:29.650
for the overall solution
of the problem.
00:44:29.650 --> 00:44:32.910
Notice that a structural
vibration problem can be
00:44:32.910 --> 00:44:36.650
thought of as a static problem
including inertia forces.
00:44:36.650 --> 00:44:42.760
And we already have noted that
in structural analysis of
00:44:42.760 --> 00:44:48.150
static problems, the
higher-order elements do quite
00:44:48.150 --> 00:44:50.490
well, and sometimes much better
00:44:50.490 --> 00:44:52.580
than the lower elements.
00:44:52.580 --> 00:44:56.730
And for this reason, they also
do quite well in dynamic
00:44:56.730 --> 00:45:00.860
analysis of structural
vibration problems.
00:45:00.860 --> 00:45:05.680
A typical problem that one might
like to solve using this
00:45:05.680 --> 00:45:07.710
implicit time integration
scheme, is
00:45:07.710 --> 00:45:09.920
schematically shown here.
00:45:09.920 --> 00:45:14.490
We have a tower that is
subjected to a blast loading
00:45:14.490 --> 00:45:18.340
as shown here schematically,
and we assume that only the
00:45:18.340 --> 00:45:20.550
structural vibrations
of this tower are
00:45:20.550 --> 00:45:21.950
required to be solved.
00:45:21.950 --> 00:45:26.210
In this case, we may very well
obtain the response of the
00:45:26.210 --> 00:45:31.250
tower using something like about
100 steps to integrate
00:45:31.250 --> 00:45:32.370
the response.
00:45:32.370 --> 00:45:36.330
In fact, we will look at a
problem of this nature in the
00:45:36.330 --> 00:45:40.080
next lecture, where we do
actually analyze a tower using
00:45:40.080 --> 00:45:41.920
the implicit time integration
scheme.
00:45:41.920 --> 00:45:45.640
And where we talk a little bit
about how the tower would be
00:45:45.640 --> 00:45:49.190
modelled, and what one would
look for to make sure that you
00:45:49.190 --> 00:45:52.720
have adequately calculated
the response.
00:45:52.720 --> 00:45:56.510
Finally, I'd like to also
point out that it can be
00:45:56.510 --> 00:46:02.330
effective to use the combined
technique of explicit and
00:46:02.330 --> 00:46:04.430
implicit integration.
00:46:04.430 --> 00:46:08.960
You might, for example, have a
problem that shows an initial
00:46:08.960 --> 00:46:12.040
wave propagation and then
a structural vibration.
00:46:12.040 --> 00:46:16.190
In this case, you might first
want to use the explicit time
00:46:16.190 --> 00:46:18.860
integration to calculate the
initial wave propagation
00:46:18.860 --> 00:46:23.380
response of the structure, and
then restart to use an
00:46:23.380 --> 00:46:26.850
implicit time integration
technique to calculate the
00:46:26.850 --> 00:46:29.960
vibrational response
of the structure.
00:46:29.960 --> 00:46:33.440
So in that case, you would
basically combine explicit and
00:46:33.440 --> 00:46:36.510
implicit time integration, but
the combination would be first
00:46:36.510 --> 00:46:39.550
explicit, and then implicit
time integration.
00:46:39.550 --> 00:46:45.570
Another kind of combination is
quite effective when you have
00:46:45.570 --> 00:46:49.370
to analyze problems where you
have very stiff regions and
00:46:49.370 --> 00:46:51.970
very flexible regions.
00:46:51.970 --> 00:46:55.550
For the flexible regions of the
problem, you might want to
00:46:55.550 --> 00:46:58.650
use explicit time integration,
and for the stiff region you
00:46:58.650 --> 00:47:01.320
might want to use implicit
time integration.
00:47:01.320 --> 00:47:05.020
A typical analysis, for example,
would be a fluid
00:47:05.020 --> 00:47:08.310
structure problem, where the
fluid is very soft and the
00:47:08.310 --> 00:47:11.160
structure is very stiff.
00:47:11.160 --> 00:47:13.660
Using the implicit time
integration for the structure
00:47:13.660 --> 00:47:16.830
and the explicit time
integration for the fluid can
00:47:16.830 --> 00:47:17.600
be effective.
00:47:17.600 --> 00:47:20.850
Of course, then you would also
have to couple these time
00:47:20.850 --> 00:47:23.940
integration schemes properly
at the boundary between the
00:47:23.940 --> 00:47:26.560
fluid and the structure.
00:47:26.560 --> 00:47:29.190
This brings me then to the end
of what I wanted to say in
00:47:29.190 --> 00:47:29.980
this lecture.
00:47:29.980 --> 00:47:33.300
In the next lecture we will
continue our discussion of
00:47:33.300 --> 00:47:37.230
solution methods for the
solution of the nonlinear
00:47:37.230 --> 00:47:39.180
dynamic equilibrium equations.
00:47:39.180 --> 00:47:42.570
And we will also show examples,
examples that
00:47:42.570 --> 00:47:45.580
demonstrate how the techniques
that I discussed in this
00:47:45.580 --> 00:47:48.750
lecture are used to
solve problems.
00:47:48.750 --> 00:47:50.160
Thank you very much for
your attention.