1
00:00:00,040 --> 00:00:02,470
The following content is
provided under a Creative
2
00:00:02,470 --> 00:00:03,880
Commons license.
3
00:00:03,880 --> 00:00:06,920
Your support will help MIT
OpenCourseWare continue to
4
00:00:06,920 --> 00:00:10,570
offer high quality educational
resources for free.
5
00:00:10,570 --> 00:00:13,470
To make a donation or view
additional materials from
6
00:00:13,470 --> 00:00:18,825
hundreds of MIT courses, visit
MIT OpenCourseWare at
7
00:00:18,825 --> 00:00:21,060
ocw.mit.edu.
8
00:00:21,060 --> 00:00:22,760
PROFESSOR: Ladies and gentlemen,
welcome to this
9
00:00:22,760 --> 00:00:25,585
lecture on nonlinear finite
element analysis of solids and
10
00:00:25,585 --> 00:00:26,410
structures.
11
00:00:26,410 --> 00:00:28,840
In this lecture, I'd like to
discuss with you some solution
12
00:00:28,840 --> 00:00:32,759
methods that we employ to
solve the finite element
13
00:00:32,759 --> 00:00:37,470
equations that we encounter in
nonlinear dynamic analysis.
14
00:00:37,470 --> 00:00:40,830
The solution methods that are,
in general used, can be
15
00:00:40,830 --> 00:00:44,360
thought of to fall into these
three categories; direct
16
00:00:44,360 --> 00:00:48,280
integration methods, mode
superposition, and
17
00:00:48,280 --> 00:00:49,690
substructuring.
18
00:00:49,690 --> 00:00:53,500
In this lecture, I would like to
talk about the explicit and
19
00:00:53,500 --> 00:00:57,500
implicit techniques, direct
integration methods, that we
20
00:00:57,500 --> 00:01:00,720
use to solve the finite
element equations.
21
00:01:00,720 --> 00:01:03,020
In the next lecture, we will
talk about the mode
22
00:01:03,020 --> 00:01:05,410
superposition and substructuring
techniques, and
23
00:01:05,410 --> 00:01:07,130
we will also consider
example solutions.
24
00:01:09,900 --> 00:01:13,840
When we want to solve the finite
element equations in
25
00:01:13,840 --> 00:01:18,080
nonlinear dynamic analysis, we
recognize that basically we
26
00:01:18,080 --> 00:01:22,120
have to operate on this
set of equations.
27
00:01:22,120 --> 00:01:27,050
Here we have the inertia forces,
the damping forces,
28
00:01:27,050 --> 00:01:32,520
the elastic forces, and the
externally applied loads.
29
00:01:32,520 --> 00:01:37,820
The inertia forces and damping
forces, well, you know them
30
00:01:37,820 --> 00:01:38,570
quite well.
31
00:01:38,570 --> 00:01:41,610
They're written in terms of
mass matrices, damping
32
00:01:41,610 --> 00:01:45,360
matrices, accelerations,
velocities.
33
00:01:45,360 --> 00:01:48,140
The elastic forces, here
"elastic" in quotes, are
34
00:01:48,140 --> 00:01:51,660
actually the nodal point forces
that are equivalent to
35
00:01:51,660 --> 00:01:53,790
the internal element stresses.
36
00:01:53,790 --> 00:01:56,780
And we talked, in the earlier
lectures, about how we
37
00:01:56,780 --> 00:02:00,000
evaluate them using, in
the updated Lagrangian
38
00:02:00,000 --> 00:02:01,540
formulation, Cauchy stresses.
39
00:02:01,540 --> 00:02:04,400
In the total Lagrangian
formulation, the second
40
00:02:04,400 --> 00:02:06,990
Piola-Kirchhoff stresses
et cetera.
41
00:02:06,990 --> 00:02:09,160
Please look up all the
information that we talked
42
00:02:09,160 --> 00:02:13,970
about in the earlier lectures,
if you want to study how we
43
00:02:13,970 --> 00:02:16,900
calculate these forces in large
displacement, large
44
00:02:16,900 --> 00:02:21,360
deformation, large
strain analysis.
45
00:02:21,360 --> 00:02:25,150
This set of equations should
hold at any time t during the
46
00:02:25,150 --> 00:02:26,770
response solution.
47
00:02:26,770 --> 00:02:31,010
And in direct integration
analysis, we consider this set
48
00:02:31,010 --> 00:02:37,400
of equations at the discrete
time points delta t apart, as
49
00:02:37,400 --> 00:02:39,630
shown down here.
50
00:02:39,630 --> 00:02:43,380
Of course, even in static
analysis, we always worked
51
00:02:43,380 --> 00:02:45,940
with the data t time
increment.
52
00:02:45,940 --> 00:02:48,980
But in static analysis, we did
not include the inertia
53
00:02:48,980 --> 00:02:54,040
forces, and the damping forces
in the analysis.
54
00:02:54,040 --> 00:02:57,700
The issues that we like to
discuss are, what are the
55
00:02:57,700 --> 00:03:00,770
basic procedures for obtaining
the solutions at
56
00:03:00,770 --> 00:03:03,170
the discrete times?
57
00:03:03,170 --> 00:03:06,570
And then once we have discussed
this question, we
58
00:03:06,570 --> 00:03:09,640
have to address the question of
which procedure should be
59
00:03:09,640 --> 00:03:11,920
used for a given problem.
60
00:03:11,920 --> 00:03:16,090
If we use, so to say a bad
procedure, "bad" in quotes for
61
00:03:16,090 --> 00:03:19,160
a given problem, the costs
can be very high
62
00:03:19,160 --> 00:03:20,870
to solve the problem.
63
00:03:20,870 --> 00:03:23,230
Therefore, it can be important,
can be very
64
00:03:23,230 --> 00:03:26,390
important to select the
appropriate technique for the
65
00:03:26,390 --> 00:03:29,000
given problem.
66
00:03:29,000 --> 00:03:31,970
Let's talk first about the
explicit time integration
67
00:03:31,970 --> 00:03:33,250
techniques.
68
00:03:33,250 --> 00:03:36,980
And the one that I like to
focus attention on is the
69
00:03:36,980 --> 00:03:38,150
central difference method.
70
00:03:38,150 --> 00:03:42,830
The central difference method
is very widely used to solve
71
00:03:42,830 --> 00:03:44,410
nonlinear problems.
72
00:03:44,410 --> 00:03:48,240
It's an explicit technique,
as we will just now see.
73
00:03:48,240 --> 00:03:51,160
The basic equations that
we are working
74
00:03:51,160 --> 00:03:53,560
with are listed here.
75
00:03:53,560 --> 00:03:57,860
This is the equilibrium
equation at time t.
76
00:03:57,860 --> 00:04:01,920
Notice here is a tf, the force
that corresponds to the
77
00:04:01,920 --> 00:04:06,490
internal element stresses
at time t.
78
00:04:06,490 --> 00:04:10,640
Here are the damping forces, and
here the inertia forces.
79
00:04:10,640 --> 00:04:15,760
We expand the velocity at
time t as shown on the
80
00:04:15,760 --> 00:04:17,450
right hand side here.
81
00:04:17,450 --> 00:04:21,200
And we expand the acceleration
at time t as shown on the
82
00:04:21,200 --> 00:04:23,230
right hand side here.
83
00:04:23,230 --> 00:04:27,050
These are the central
difference formulas.
84
00:04:27,050 --> 00:04:31,390
This method is mainly used for
wave propagation problems.
85
00:04:31,390 --> 00:04:34,150
It's an explicit technique,
and this is an important
86
00:04:34,150 --> 00:04:37,560
point, it's an explicit
technique because we are
87
00:04:37,560 --> 00:04:40,930
looking at the equilibrium
equation at time t to
88
00:04:40,930 --> 00:04:45,210
calculate the response for
time t plus delta t.
89
00:04:45,210 --> 00:04:50,640
Any method that does just that,
looks at the equilibrium
90
00:04:50,640 --> 00:04:53,980
equation at time t to obtain the
solution for the response
91
00:04:53,980 --> 00:04:57,790
at time t plus delta t is termed
an explicit method, an
92
00:04:57,790 --> 00:04:59,220
explicit integration method.
93
00:04:59,220 --> 00:05:01,410
So the central difference
method is
94
00:05:01,410 --> 00:05:04,290
such an explicit method.
95
00:05:04,290 --> 00:05:07,470
The important point to recognize
right here is that
96
00:05:07,470 --> 00:05:12,710
we don't set up a k matrix,
a stiffness matrix.
97
00:05:12,710 --> 00:05:16,680
And that is one important
fact that pertains
98
00:05:16,680 --> 00:05:19,680
to an explicit technique.
99
00:05:19,680 --> 00:05:23,480
Using these three equations,
we can
100
00:05:23,480 --> 00:05:26,850
directly develop one equation.
101
00:05:26,850 --> 00:05:29,800
And that equation looks
as shown up here.
102
00:05:29,800 --> 00:05:33,380
Notice we have three equations
and three unknowns
103
00:05:33,380 --> 00:05:34,390
which can be solved.
104
00:05:34,390 --> 00:05:37,050
The three unknowns being the
displacement at time t plus
105
00:05:37,050 --> 00:05:40,340
delta t, the velocity at time
t plus delta t, and the
106
00:05:40,340 --> 00:05:42,180
acceleration at time
t plus delta t.
107
00:05:42,180 --> 00:05:43,820
Three equations and
three unknowns.
108
00:05:43,820 --> 00:05:46,990
We can solve these equations
and the governing equation
109
00:05:46,990 --> 00:05:50,840
that we obtain in that solution
is listed right here.
110
00:05:50,840 --> 00:05:53,680
Notice m over delta t squared.
111
00:05:53,680 --> 00:05:57,740
Notice c over 2 delta t here.
112
00:05:57,740 --> 00:06:00,500
And then the unknown
displacement corresponding to
113
00:06:00,500 --> 00:06:04,000
time t plus delta t, and on the
right hand side, we have
114
00:06:04,000 --> 00:06:06,080
what we call an effective
load vector
115
00:06:06,080 --> 00:06:07,780
corresponding to time t.
116
00:06:07,780 --> 00:06:12,700
There is a hat here, and that
hat means that there are a
117
00:06:12,700 --> 00:06:15,170
large number of terms to
be taken into account.
118
00:06:15,170 --> 00:06:18,800
Namely, those corresponding to
the inertia and the damping in
119
00:06:18,800 --> 00:06:19,650
the system.
120
00:06:19,650 --> 00:06:25,740
Here we have tr hat is equal
to tr, the actual loads
121
00:06:25,740 --> 00:06:27,550
applied corresponding
to a time t.
122
00:06:27,550 --> 00:06:29,310
These are the externally
applied loads
123
00:06:29,310 --> 00:06:30,860
to the nodal points.
124
00:06:30,860 --> 00:06:34,050
tf is the force vector
corresponding to the internal
125
00:06:34,050 --> 00:06:36,560
element stresses, as
I said earlier.
126
00:06:36,560 --> 00:06:41,260
Here we have an inertia
contribution.
127
00:06:41,260 --> 00:06:45,140
Here another inertia
contribution and a damping
128
00:06:45,140 --> 00:06:47,240
contribution.
129
00:06:47,240 --> 00:06:50,170
But it is important to recognize
that the right hand
130
00:06:50,170 --> 00:06:52,040
side is known.
131
00:06:52,040 --> 00:06:56,270
In other words, all the terms
here are known because these
132
00:06:56,270 --> 00:06:59,080
terms only involve displacements
corresponding to
133
00:06:59,080 --> 00:07:02,750
time t minus delta t,
and displacements
134
00:07:02,750 --> 00:07:04,350
corresponding to time t.
135
00:07:04,350 --> 00:07:10,650
And therefore we know this tr
hat, we can plug it in here on
136
00:07:10,650 --> 00:07:13,300
the right hand side and
calculate the unknown
137
00:07:13,300 --> 00:07:17,090
displacements corresponding
to time t plus delta t.
138
00:07:17,090 --> 00:07:22,520
The method is most useful when
m and c, the mass and damping
139
00:07:22,520 --> 00:07:24,340
matrices are diagonal.
140
00:07:24,340 --> 00:07:29,540
Because in that case, these
equations here decouple as
141
00:07:29,540 --> 00:07:30,960
shown right here.
142
00:07:30,960 --> 00:07:34,790
Notice we can calculate the
individual displacements
143
00:07:34,790 --> 00:07:38,180
components, or the displacement
that's individual
144
00:07:38,180 --> 00:07:43,870
degrees of freedom, one after
the other as shown here.
145
00:07:43,870 --> 00:07:52,720
Once we have evaluated tr hat,
tr hat let's remember,
146
00:07:52,720 --> 00:07:55,620
contains these terms here.
147
00:07:55,620 --> 00:08:00,880
In other words, although m and
c are diagonal, there is some
148
00:08:00,880 --> 00:08:06,170
coupling from the equation j
into equation j minus 1 and
149
00:08:06,170 --> 00:08:10,730
into equation j plus 1 right
here on the right hand side,
150
00:08:10,730 --> 00:08:13,790
because the element
contributions overlap the
151
00:08:13,790 --> 00:08:15,880
degrees of freedom.
152
00:08:15,880 --> 00:08:21,780
And therefore, this here is
where the element coupling at
153
00:08:21,780 --> 00:08:22,700
the nodal points enters.
154
00:08:22,700 --> 00:08:28,210
But once we have calculated
this vector here, I should
155
00:08:28,210 --> 00:08:30,190
point to this vector here.
156
00:08:30,190 --> 00:08:33,340
Once we have calculated this
vector, we know the individual
157
00:08:33,340 --> 00:08:37,990
components, then we enter
here with that component
158
00:08:37,990 --> 00:08:42,710
corresponding to a degree of
freedom i, and directly we can
159
00:08:42,710 --> 00:08:46,050
calculate the nodal point
displacements corresponding to
160
00:08:46,050 --> 00:08:48,730
a degree of freedom i.
161
00:08:48,730 --> 00:08:52,270
That's an important point to
recognize, that there is no
162
00:08:52,270 --> 00:08:55,650
need really to set
up a k matrix, to
163
00:08:55,650 --> 00:08:57,835
trianglize a k matrix.
164
00:09:00,470 --> 00:09:05,310
We don't do that here, and that
means that per time step,
165
00:09:05,310 --> 00:09:07,610
this solution method can
be very effective.
166
00:09:10,450 --> 00:09:17,120
Note that in the solution if
the damping is 0, in other
167
00:09:17,120 --> 00:09:20,100
words, we don't have damping in
the system, cii is equal to
168
00:09:20,100 --> 00:09:24,370
0, we need that mii
is greater than 0.
169
00:09:24,370 --> 00:09:26,860
Otherwise we can't solve
the equations.
170
00:09:26,860 --> 00:09:32,440
Notice that tf is calculated
by summing all the element
171
00:09:32,440 --> 00:09:33,670
contributions.
172
00:09:33,670 --> 00:09:38,620
Summing over m means summing
over all the elements.
173
00:09:38,620 --> 00:09:42,500
And as I mentioned earlier, this
is where coupling comes
174
00:09:42,500 --> 00:09:48,870
into the solution because one
nodal point, one degree of
175
00:09:48,870 --> 00:09:52,980
freedom at that nodal point
carries contribution from all
176
00:09:52,980 --> 00:09:56,320
the elements that surround
that nodal point.
177
00:09:56,320 --> 00:10:01,770
To start the solution, we need
a special equation because,
178
00:10:01,770 --> 00:10:07,120
remember in the solution for
time t plus delta t, we need
179
00:10:07,120 --> 00:10:09,330
the information, the
displacements corresponding to
180
00:10:09,330 --> 00:10:13,350
time t, and corresponding
to time t minus delta t.
181
00:10:13,350 --> 00:10:17,840
So, if t is equal to 0, and we
want to march ahead to delta
182
00:10:17,840 --> 00:10:19,825
t, in the first step--
183
00:10:19,825 --> 00:10:23,740
we are in the first step--
we use this equation to
184
00:10:23,740 --> 00:10:27,610
calculate, to evaluate
minus delta tu.
185
00:10:27,610 --> 00:10:32,640
And that means we can then
directly start our explicit
186
00:10:32,640 --> 00:10:35,680
solution process.
187
00:10:35,680 --> 00:10:40,020
The central difference method is
only conditionally stable.
188
00:10:40,020 --> 00:10:43,680
The condition that has to be
satisfied on the time step is
189
00:10:43,680 --> 00:10:46,030
given in this equation.
190
00:10:46,030 --> 00:10:51,060
Delta t critical is a critical
time step, and delta t, the
191
00:10:51,060 --> 00:10:54,170
actual time step that is used
in the solution, has to be
192
00:10:54,170 --> 00:10:55,890
smaller than that value.
193
00:10:55,890 --> 00:10:57,830
Smaller or equal
to that value.
194
00:10:57,830 --> 00:11:01,690
Notice delta t critical is given
as tn over pi, where tn
195
00:11:01,690 --> 00:11:06,136
is the smallest period in the
finite element assemblage.
196
00:11:06,136 --> 00:11:09,710
In nonlinear analysis, we have
to recognize that tn changes
197
00:11:09,710 --> 00:11:11,190
during the time history.
198
00:11:11,190 --> 00:11:15,800
Response, tn would become
smaller if the system
199
00:11:15,800 --> 00:11:21,820
stiffens, due for example to
large displacement effects.
200
00:11:21,820 --> 00:11:25,900
And tn would become larger if
the system softens due to
201
00:11:25,900 --> 00:11:27,340
elasto-plasticity effects.
202
00:11:31,440 --> 00:11:35,590
One interesting and most
important point is that tn can
203
00:11:35,590 --> 00:11:37,060
be estimated.
204
00:11:37,060 --> 00:11:43,520
And it can be estimated, or we
can find a bound on tn which
205
00:11:43,520 --> 00:11:47,200
we can then use to calculate
the critical time step.
206
00:11:47,200 --> 00:11:49,920
And that is achieved
as follows.
207
00:11:49,920 --> 00:11:56,120
We notice that the omega n
squared, omega n being 2 pi
208
00:11:56,120 --> 00:12:01,895
over tn radians per second omega
n squared is smaller or
209
00:12:01,895 --> 00:12:10,640
equal to the maximum of the
omega m, nms squared over all
210
00:12:10,640 --> 00:12:13,140
elements m.
211
00:12:13,140 --> 00:12:14,350
Now what does this mean?
212
00:12:14,350 --> 00:12:20,310
It means here that we're looking
at omega nm squared of
213
00:12:20,310 --> 00:12:27,140
each of the elements, and we
select the largest such value
214
00:12:27,140 --> 00:12:30,260
to put in here.
215
00:12:30,260 --> 00:12:32,520
And that then gives us an upper
216
00:12:32,520 --> 00:12:36,970
bound on omega n squared.
217
00:12:36,970 --> 00:12:40,130
The actual omega n that
we need, because
218
00:12:40,130 --> 00:12:43,690
omega n gives us tn.
219
00:12:43,690 --> 00:12:47,960
So we want to calculate the the
largest frequency of all
220
00:12:47,960 --> 00:12:51,190
individual elements.
221
00:12:51,190 --> 00:12:54,530
And with that largest frequency,
we enter in this
222
00:12:54,530 --> 00:13:01,930
equation here, to get the bound
on tn, noting that in
223
00:13:01,930 --> 00:13:07,340
nonlinear analysis omega nm
changes, just the way I said
224
00:13:07,340 --> 00:13:09,300
earlier tn will also change.
225
00:13:09,300 --> 00:13:14,190
But all the omega nm changes
and tn correspondingly
226
00:13:14,190 --> 00:13:16,590
changes, this equation
is always satisfied.
227
00:13:19,740 --> 00:13:24,410
The time integration step that
we can actually use then, is
228
00:13:24,410 --> 00:13:26,610
given by as this equation.
229
00:13:26,610 --> 00:13:30,590
We can use delta t is equal to
2 over omega nm maximum.
230
00:13:30,590 --> 00:13:36,580
The maximum radiant per second
frequency that we have for any
231
00:13:36,580 --> 00:13:39,940
one of the elements in
the complete mesh.
232
00:13:39,940 --> 00:13:44,760
Because that expression is
smaller or equal than the
233
00:13:44,760 --> 00:13:46,520
critical time step.
234
00:13:46,520 --> 00:13:51,310
We could call this value here 2
over omega nm , the critical
235
00:13:51,310 --> 00:13:53,930
time step of element m.
236
00:13:53,930 --> 00:13:57,240
And therefore, what we're doing
here really is, that
237
00:13:57,240 --> 00:14:01,430
we're taking the smallest of the
critical time steps over
238
00:14:01,430 --> 00:14:07,480
all elements as that time
step in the solution.
239
00:14:07,480 --> 00:14:10,240
And that time step we
know is smaller or
240
00:14:10,240 --> 00:14:11,685
equal to delta t critical.
241
00:14:14,200 --> 00:14:19,470
Let's look at the proof of this
relationship, because
242
00:14:19,470 --> 00:14:23,500
it's an interesting proof and it
gives us a bit more insight
243
00:14:23,500 --> 00:14:29,240
into why this relationship
really holds.
244
00:14:29,240 --> 00:14:32,490
So I'd like to prove to you
very briefly that this
245
00:14:32,490 --> 00:14:34,570
relationship here is valid.
246
00:14:34,570 --> 00:14:39,130
And we do so by using the
Rayleigh quotient.
247
00:14:39,130 --> 00:14:44,620
Please look at the textbook
where the Rayleigh quotient is
248
00:14:44,620 --> 00:14:48,630
discussed, and I have to now
assume that you are familiar
249
00:14:48,630 --> 00:14:50,310
with the Rayleigh quotient.
250
00:14:50,310 --> 00:14:55,750
If you write down the Rayleigh
quotient for omega n squared,
251
00:14:55,750 --> 00:14:57,990
it is this relationship here.
252
00:14:57,990 --> 00:15:04,620
Omega n being the largest
frequency of the complete
253
00:15:04,620 --> 00:15:05,970
finite element mesh.
254
00:15:05,970 --> 00:15:10,020
Omega n is equal to 2
pi divided by tn.
255
00:15:10,020 --> 00:15:13,520
And tn was the smallest period
in the mesh that we are
256
00:15:13,520 --> 00:15:15,630
interested in knowing.
257
00:15:15,630 --> 00:15:19,770
Let us define this relationship
here, just the
258
00:15:19,770 --> 00:15:23,250
definition, and this
relationship here.
259
00:15:23,250 --> 00:15:25,800
Once again, just
the definition.
260
00:15:25,800 --> 00:15:31,140
Then, by simply substituting
from here and there into this
261
00:15:31,140 --> 00:15:34,890
relationship here, we surely
can write this equation as
262
00:15:34,890 --> 00:15:35,930
shown here.
263
00:15:35,930 --> 00:15:39,280
We're summing over all of the
us, and we're summing
264
00:15:39,280 --> 00:15:40,795
over all of the is.
265
00:15:40,795 --> 00:15:43,620
m going over all the elements.
266
00:15:43,620 --> 00:15:46,580
If you have a 1,000 elements in
the mesh, then n would go
267
00:15:46,580 --> 00:15:50,940
from 1 to 1,000.
268
00:15:50,940 --> 00:15:55,640
If we now look at the Rayleigh
quotient for a single element,
269
00:15:55,640 --> 00:16:00,420
we can say that rho m is given
by this relationship.
270
00:16:00,420 --> 00:16:05,700
Which is nothing else, using our
definition of um and im,
271
00:16:05,700 --> 00:16:10,310
is nothing else than this
relationship here.
272
00:16:10,310 --> 00:16:15,300
However, we can immediately
recognize that rho m, this
273
00:16:15,300 --> 00:16:19,950
value here, must be smaller or
equal to omega nm squared,
274
00:16:19,950 --> 00:16:26,440
where omega nm is the largest
frequency radian per second of
275
00:16:26,440 --> 00:16:28,920
the element m.
276
00:16:28,920 --> 00:16:34,520
This is here a relationship
that must hold because the
277
00:16:34,520 --> 00:16:39,110
Rayleigh quotient statement,
or the Rayleigh quotient
278
00:16:39,110 --> 00:16:42,230
theorem, tells us that
this must hold.
279
00:16:42,230 --> 00:16:46,890
And once again, please look up
in the textbook a general
280
00:16:46,890 --> 00:16:51,140
proof, there's a general proof
given that in fact this is a
281
00:16:51,140 --> 00:16:53,760
valid statement.
282
00:16:53,760 --> 00:16:59,480
Where now we can use this
statement, go back to this
283
00:16:59,480 --> 00:17:04,800
equation and immediately
identify that this must hold.
284
00:17:04,800 --> 00:17:09,690
And notice once again I'm
looking here at the element m.
285
00:17:09,690 --> 00:17:14,319
Of course, this relationship
must hold for all elements, m
286
00:17:14,319 --> 00:17:17,020
going from 1 to the total
number of elements you
287
00:17:17,020 --> 00:17:19,619
actually have in the mesh.
288
00:17:19,619 --> 00:17:24,579
This is the relationship we
will be using, and we
289
00:17:24,579 --> 00:17:29,390
substitute that relationship in
our earlier expression, for
290
00:17:29,390 --> 00:17:32,730
the omega n squared and directly
arrive at this
291
00:17:32,730 --> 00:17:35,190
relationship here.
292
00:17:35,190 --> 00:17:43,730
Now we look at this part here,
and recognize that we can take
293
00:17:43,730 --> 00:17:50,100
this out of the summation sign
by using the largest or the
294
00:17:50,100 --> 00:17:53,100
effect, we can take the effect
out of the summation sign, I
295
00:17:53,100 --> 00:17:55,610
should've said, by taking
the largest
296
00:17:55,610 --> 00:17:58,530
value of omega nm squared.
297
00:17:58,530 --> 00:18:02,450
And that's being done here and
now we can cancel these two
298
00:18:02,450 --> 00:18:06,950
terms resulting into this
relationship which
299
00:18:06,950 --> 00:18:09,790
completes our proof.
300
00:18:09,790 --> 00:18:13,420
It's an interesting proof, and
once again if you're not
301
00:18:13,420 --> 00:18:15,740
familiar with the Rayleigh
quotient, please look up the
302
00:18:15,740 --> 00:18:19,180
information in the textbook,
after which I think you can
303
00:18:19,180 --> 00:18:22,640
quite easily follow
this proof.
304
00:18:22,640 --> 00:18:26,070
The important point now is that
the largest frequencies
305
00:18:26,070 --> 00:18:31,920
of simple elements can be
estimated analytically.
306
00:18:31,920 --> 00:18:35,310
Sometimes we can calculate them
exactly, and here we have
307
00:18:35,310 --> 00:18:39,290
a case where we can calculate
the largest frequency of the
308
00:18:39,290 --> 00:18:40,720
element exactly.
309
00:18:40,720 --> 00:18:45,430
Here we have a two nodal truss
element, m here, m here, equal
310
00:18:45,430 --> 00:18:47,750
to rho AL over 2.
311
00:18:47,750 --> 00:18:50,310
L being the length of the
truss, A being the cross
312
00:18:50,310 --> 00:18:53,250
sectional area, rho being
the mass density.
313
00:18:53,250 --> 00:18:56,440
If we write down the eigenvalue
problem for this
314
00:18:56,440 --> 00:19:00,030
single element, we obtain
this relationship here.
315
00:19:00,030 --> 00:19:05,890
Stiffness matrix times the
eigenvector phi omega squared,
316
00:19:05,890 --> 00:19:10,680
which will be an eigenvalue, the
mass matrix phi being the
317
00:19:10,680 --> 00:19:12,360
sought eigenvector.
318
00:19:12,360 --> 00:19:18,100
And if we solve this eigenvalue
problem, we obtain
319
00:19:18,100 --> 00:19:24,260
directly as the eigenvalues
these two values here.
320
00:19:24,260 --> 00:19:26,180
There are only two of course,
because this is
321
00:19:26,180 --> 00:19:27,960
a two by two matrix.
322
00:19:27,960 --> 00:19:28,950
And so is that.
323
00:19:28,950 --> 00:19:30,590
So we have only two
eigenvalues.
324
00:19:30,590 --> 00:19:36,460
The first eigenvalue is 0, the
rigid body mode in the system.
325
00:19:36,460 --> 00:19:40,330
And the second eigenvalue is
given here, which is equal to
326
00:19:40,330 --> 00:19:42,750
omega n squared.
327
00:19:42,750 --> 00:19:45,950
And there is the relationship
of that second eigenvalue
328
00:19:45,950 --> 00:19:51,380
obtained by some very simple
calculations, and this is a
329
00:19:51,380 --> 00:19:56,210
rewriting of that equation,
of that relationship here.
330
00:19:56,210 --> 00:20:00,160
c is a wave speed
in the system.
331
00:20:00,160 --> 00:20:06,230
In fact, c squared is, as you
can see here, e over rho.
332
00:20:06,230 --> 00:20:09,860
It's this relationship that
we can now use, very
333
00:20:09,860 --> 00:20:15,440
conveniently, to get a time step
in a finite element mesh
334
00:20:15,440 --> 00:20:21,480
that will satisfy the critical
time step criterion.
335
00:20:21,480 --> 00:20:25,430
In other words, if we have a
finite element mesh that is
336
00:20:25,430 --> 00:20:30,730
consisting of an assemblage of
such truss elements, then to
337
00:20:30,730 --> 00:20:34,490
get the time step for the
explicit time integration,
338
00:20:34,490 --> 00:20:39,040
which will be a good time step,
which will be smaller or
339
00:20:39,040 --> 00:20:42,730
equal to the critical time step,
what we do is we simply
340
00:20:42,730 --> 00:20:46,880
calculate this relationship here
for each of the elements.
341
00:20:46,880 --> 00:20:52,410
And this relationship will give
us then, the proper time
342
00:20:52,410 --> 00:20:58,710
step by evaluating the smallest
time step that we can
343
00:20:58,710 --> 00:21:00,240
use for any one of
the elements.
344
00:21:00,240 --> 00:21:05,520
In other words, we simply take
2 over omega n, which is
345
00:21:05,520 --> 00:21:10,840
written here, and we recognize
that L over c, where L is the
346
00:21:10,840 --> 00:21:17,170
length of the element, is the
time step that we can use in
347
00:21:17,170 --> 00:21:18,360
the analysis.
348
00:21:18,360 --> 00:21:23,920
Notice L over c, we want to use
the smallest value of L
349
00:21:23,920 --> 00:21:28,150
over c coming from all
of the elements.
350
00:21:28,150 --> 00:21:30,840
L would be an element length.
351
00:21:30,840 --> 00:21:33,880
c would be the wave speed
through the element, and we
352
00:21:33,880 --> 00:21:38,560
want to use the smallest value
of L over c, looking at all
353
00:21:38,560 --> 00:21:44,390
the elements, as the time step
for our time integration.
354
00:21:44,390 --> 00:21:47,710
Notice here that L over c is
the time required for wave
355
00:21:47,710 --> 00:21:51,130
front to travel through
the element.
356
00:21:51,130 --> 00:21:54,310
That is a physical
interpretation of
357
00:21:54,310 --> 00:21:57,980
this L over c value.
358
00:21:57,980 --> 00:22:00,890
Let's look at some
modelling aspects
359
00:22:00,890 --> 00:22:03,090
regarding this technique.
360
00:22:03,090 --> 00:22:09,700
Let the applied wavelengths
be Lw as shown here.
361
00:22:09,700 --> 00:22:12,920
Here schematically
we show a wave.
362
00:22:12,920 --> 00:22:18,680
And let's see that how we would
model the finite element
363
00:22:18,680 --> 00:22:22,930
system appropriately when
we apply to that system
364
00:22:22,930 --> 00:22:25,130
this kind of wave.
365
00:22:25,130 --> 00:22:33,330
Well, we first of all find the
value tw, which is really the
366
00:22:33,330 --> 00:22:37,970
time required for the total
wave to propagate across a
367
00:22:37,970 --> 00:22:40,320
particular point in
the mesh, with the
368
00:22:40,320 --> 00:22:42,870
wave speed c of course.
369
00:22:42,870 --> 00:22:49,380
We then choose delta t via this
relationship here, where
370
00:22:49,380 --> 00:22:55,290
n is the number of time steps
used to represent the wave.
371
00:22:55,290 --> 00:23:02,250
With delta t now known, we
directly get Le, which is the
372
00:23:02,250 --> 00:23:04,670
element lengths.
373
00:23:04,670 --> 00:23:07,360
Well, I should not quite say
it's the element lengths, it's
374
00:23:07,360 --> 00:23:08,950
related to the element
lengths.
375
00:23:08,950 --> 00:23:16,150
For very simple problems, and
particular for our problems in
376
00:23:16,150 --> 00:23:19,070
a model using truss elements,
two nodal truss elements,
377
00:23:19,070 --> 00:23:21,420
actually Le would be the
element lengths.
378
00:23:21,420 --> 00:23:25,370
But if we have more complex
elements to model the finite
379
00:23:25,370 --> 00:23:29,070
elements solution, then this
here, this Le value, would be
380
00:23:29,070 --> 00:23:32,250
related to the element
lengths.
381
00:23:32,250 --> 00:23:35,770
Note that in this formula here,
n has to be chosen by
382
00:23:35,770 --> 00:23:39,280
the analyst to be a reasonable
value so as to be able to
383
00:23:39,280 --> 00:23:43,000
represent the total wave
appropriately.
384
00:23:43,000 --> 00:23:45,840
Here we have some
observations.
385
00:23:45,840 --> 00:23:50,230
Note that in 1D analysis, c,
the wave speed, is given as
386
00:23:50,230 --> 00:23:52,220
the square root out
of e over rho.
387
00:23:52,220 --> 00:23:54,200
I mentioned that already
earlier.
388
00:23:54,200 --> 00:23:57,170
Notice that in nonlinear
analysis, delta t must satisfy
389
00:23:57,170 --> 00:24:00,250
the stability limit throughout
the solution.
390
00:24:00,250 --> 00:24:04,700
This is a very important point
and in the next lecture we
391
00:24:04,700 --> 00:24:08,650
will look at some examples that
demonstrates what happens
392
00:24:08,650 --> 00:24:11,980
if you don't satisfy
the stability limit
393
00:24:11,980 --> 00:24:14,120
throughout the solution.
394
00:24:14,120 --> 00:24:17,445
It may also be effective to
change the time step during
395
00:24:17,445 --> 00:24:19,490
the analysis.
396
00:24:19,490 --> 00:24:21,210
Of course, this you
would do in a
397
00:24:21,210 --> 00:24:22,460
nonlinear analysis, typically.
398
00:24:25,860 --> 00:24:27,300
For the application.
399
00:24:27,300 --> 00:24:31,480
or for the use of the explicit
time integration method, when
400
00:24:31,480 --> 00:24:34,230
you use the explicit time
integration method, one
401
00:24:34,230 --> 00:24:36,620
generally uses lower-order
elements.
402
00:24:36,620 --> 00:24:38,300
They are usually preferable.
403
00:24:38,300 --> 00:24:41,010
And in this particular case,
you want to have
404
00:24:41,010 --> 00:24:44,250
Le here and Le there.
405
00:24:44,250 --> 00:24:47,710
Uniform meshing, uniform
meshing.
406
00:24:47,710 --> 00:24:51,070
For higher-order elements, you
could also use the explicit
407
00:24:51,070 --> 00:24:52,600
time integration.
408
00:24:52,600 --> 00:24:57,580
And there again, you want to
really have that this L is
409
00:24:57,580 --> 00:25:00,990
equal to that L star.
410
00:25:00,990 --> 00:25:07,390
If you don't have that
condition, then this L here
411
00:25:07,390 --> 00:25:10,920
will govern the actual
time step.
412
00:25:10,920 --> 00:25:14,400
The smallest length through
the element will
413
00:25:14,400 --> 00:25:15,840
govern the time step.
414
00:25:15,840 --> 00:25:21,890
And if you use this value here
for Le, with 8 there, you have
415
00:25:21,890 --> 00:25:24,890
a conservative value for Le.
416
00:25:24,890 --> 00:25:29,610
This 8 here enters into the
solution because you have this
417
00:25:29,610 --> 00:25:34,000
midnode, and the midnode carries
more stiffness then
418
00:25:34,000 --> 00:25:35,250
the corner nodes.
419
00:25:38,750 --> 00:25:43,050
Let us look further at
some observations
420
00:25:43,050 --> 00:25:44,440
regarding the method.
421
00:25:44,440 --> 00:25:48,140
In the analysis of this problem,
where we have a beam
422
00:25:48,140 --> 00:25:54,480
subjected to an end step load as
shown here, we can get the
423
00:25:54,480 --> 00:25:58,790
exact solution by simply
choosing a finite element mesh
424
00:25:58,790 --> 00:26:03,070
of truss elements that satisfy
this relationship here.
425
00:26:03,070 --> 00:26:07,420
In other words, we represent
this beam here, using two
426
00:26:07,420 --> 00:26:09,280
nodal truss elements.
427
00:26:09,280 --> 00:26:13,600
And these two nodal truss
elements have a length equal
428
00:26:13,600 --> 00:26:15,440
to c times delta t.
429
00:26:15,440 --> 00:26:18,670
Of course, we're looking here at
a linear elastic solution,
430
00:26:18,670 --> 00:26:20,870
meaning that c is constant.
431
00:26:20,870 --> 00:26:24,040
And we will get, this is the
important point, we will get
432
00:26:24,040 --> 00:26:27,190
the exact solution for
any number of truss
433
00:26:27,190 --> 00:26:29,270
elements that are used.
434
00:26:29,270 --> 00:26:33,450
In fact, we will look at this
problem later on in the next
435
00:26:33,450 --> 00:26:38,370
lecture in more detail, and we
will use this fact to obtain
436
00:26:38,370 --> 00:26:42,240
the exact solution, and we'll
study also what happens if we
437
00:26:42,240 --> 00:26:44,900
change the time step.
438
00:26:44,900 --> 00:26:49,460
So I'd like to refer you to the
next lecture in which I
439
00:26:49,460 --> 00:26:53,910
will be talking more
about this problem.
440
00:26:53,910 --> 00:26:58,110
Using the explicit time
integration method, it is very
441
00:26:58,110 --> 00:27:01,700
important to use uniform
meshing.
442
00:27:01,700 --> 00:27:04,370
And the reason for it
is given right here.
443
00:27:04,370 --> 00:27:09,860
Namely, we want to have then a
time step that is not unduly
444
00:27:09,860 --> 00:27:14,620
small in any region
of the total mesh.
445
00:27:14,620 --> 00:27:19,170
If you don't use uniform
meshing, then the time step is
446
00:27:19,170 --> 00:27:24,760
governed by a particular part
in the mesh, the time step
447
00:27:24,760 --> 00:27:26,610
size I should say,
is governed by a
448
00:27:26,610 --> 00:27:28,040
particular part in the mesh.
449
00:27:28,040 --> 00:27:33,340
Namely, those elements which
have the smallest Le, the way
450
00:27:33,340 --> 00:27:34,830
I talked about it earlier.
451
00:27:34,830 --> 00:27:40,020
And that same time step would
be used for the whole mesh,
452
00:27:40,020 --> 00:27:43,370
which means then, that time
step, that small time step
453
00:27:43,370 --> 00:27:46,150
would actually be too
small for the other
454
00:27:46,150 --> 00:27:48,510
parts of the mesh.
455
00:27:48,510 --> 00:27:50,780
Of course, one could say, well,
why not use different
456
00:27:50,780 --> 00:27:53,010
time steps in the mesh?
457
00:27:53,010 --> 00:28:00,170
Namely, effective time steps for
the elements that are used
458
00:28:00,170 --> 00:28:01,940
in the various regions
of the mesh.
459
00:28:01,940 --> 00:28:05,410
That can be done, but then you
have to use special coupling
460
00:28:05,410 --> 00:28:11,040
techniques to couple the time
step integrations where you
461
00:28:11,040 --> 00:28:15,110
use in one mesh a small time
step, in one part of the mesh
462
00:28:15,110 --> 00:28:17,820
a small time step, in the other
products the mesh a
463
00:28:17,820 --> 00:28:18,880
larger time step.
464
00:28:18,880 --> 00:28:22,850
It can be done, but it needs
special considerations
465
00:28:22,850 --> 00:28:26,320
regarding the coupling of the
time integrations in the
466
00:28:26,320 --> 00:28:28,800
different parts of the mesh.
467
00:28:28,800 --> 00:28:32,540
We may also want to note that
a system with a very large
468
00:28:32,540 --> 00:28:36,740
bandwidth can be solved
frequently, quite efficiently,
469
00:28:36,740 --> 00:28:39,690
with the central difference
method, although the problem
470
00:28:39,690 --> 00:28:41,680
may not be a wave propagation
problem.
471
00:28:41,680 --> 00:28:47,940
The reason for this fact is
that in the explicit time
472
00:28:47,940 --> 00:28:53,710
integration, you do not set up a
k matrix, and therefore, the
473
00:28:53,710 --> 00:28:58,240
large bandwidth that you're
seeing in the k matrix is not
474
00:28:58,240 --> 00:29:00,510
really a hindrance.
475
00:29:00,510 --> 00:29:06,610
It is not being used because the
k matrix is not being set
476
00:29:06,610 --> 00:29:10,300
up, and you have no
factorization or LDL transpose
477
00:29:10,300 --> 00:29:12,550
factorization of the
coefficient matrix.
478
00:29:12,550 --> 00:29:16,100
And therefore, this high cost
that lies in the factorization
479
00:29:16,100 --> 00:29:21,280
of the coefficient matrix is not
necessary, is not used, is
480
00:29:21,280 --> 00:29:24,870
not expensed in the explicit
time integration.
481
00:29:24,870 --> 00:29:29,110
And for that reason, it can be
effective to use explicit time
482
00:29:29,110 --> 00:29:34,056
integration when there would
be such a large bandwidth.
483
00:29:34,056 --> 00:29:38,090
Explicit time integration lends
itself also very well to
484
00:29:38,090 --> 00:29:39,850
parallel processing.
485
00:29:39,850 --> 00:29:43,000
Notice that the basic equations
that we are solving
486
00:29:43,000 --> 00:29:47,546
for in explicit time
integration, in each time step
487
00:29:47,546 --> 00:29:49,740
are given as here.
488
00:29:49,740 --> 00:29:54,080
T plus delta tu is equal to
something that I denote here
489
00:29:54,080 --> 00:29:56,740
by the three dots,
times tr hat.
490
00:29:56,740 --> 00:30:00,390
And this tr hat vector
we defined earlier.
491
00:30:00,390 --> 00:30:05,120
So all we need to do is
calculate these entries in the
492
00:30:05,120 --> 00:30:08,570
vector here to get our unknown
displacements.
493
00:30:08,570 --> 00:30:11,890
Well, that can be done quite
effectively in parallel.
494
00:30:11,890 --> 00:30:16,560
We might want to calculate these
entries here at the same
495
00:30:16,560 --> 00:30:20,100
time as we calculate
say, these entries.
496
00:30:20,100 --> 00:30:23,130
Where this is here a particular
element group, and
497
00:30:23,130 --> 00:30:25,920
that are the entries
corresponding to an adjacent
498
00:30:25,920 --> 00:30:26,810
element group.
499
00:30:26,810 --> 00:30:29,260
And these are the entries
corresponding to another
500
00:30:29,260 --> 00:30:30,050
element group.
501
00:30:30,050 --> 00:30:34,460
So we can, in other words, in
parallel, process all the
502
00:30:34,460 --> 00:30:38,440
information that goes in here,
that goes in here, and that
503
00:30:38,440 --> 00:30:40,580
goes in here, and so on.
504
00:30:40,580 --> 00:30:44,560
Notice there's some coupling
in that information, which
505
00:30:44,560 --> 00:30:48,470
goes across this red
horizontal line.
506
00:30:48,470 --> 00:30:53,280
That coupling comes from the fm
contributions, namely the
507
00:30:53,280 --> 00:30:58,360
element contributions that
will be elements that lie
508
00:30:58,360 --> 00:31:01,220
across these equations,
so to say.
509
00:31:01,220 --> 00:31:04,010
And that's where a coupling
lies, and that has to be
510
00:31:04,010 --> 00:31:07,720
properly taken into account in
the parallel processing.
511
00:31:07,720 --> 00:31:12,120
But in principle, explicit time
integration because the
512
00:31:12,120 --> 00:31:17,110
equation that we are dealing
with, as shown here, lends
513
00:31:17,110 --> 00:31:22,450
itself very well to parallel
processing on computers.
514
00:31:22,450 --> 00:31:27,550
Let us now look at the other
technique that I wanted to
515
00:31:27,550 --> 00:31:30,380
discuss with you briefly
in this lecture.
516
00:31:30,380 --> 00:31:34,890
And this is an implicit time
integration method.
517
00:31:34,890 --> 00:31:38,780
It is a trapezoidal rule, which
is very widely used in
518
00:31:38,780 --> 00:31:41,090
implicit time integration.
519
00:31:41,090 --> 00:31:46,470
But let's first look at the
basic equation that we now are
520
00:31:46,470 --> 00:31:47,500
operating on.
521
00:31:47,500 --> 00:31:50,730
In implicit time integration,
in any implicit time
522
00:31:50,730 --> 00:31:55,340
integration scheme, we are
solving this equation here.
523
00:31:58,020 --> 00:32:03,690
At time t plus delta t to obtain
the solution at time t
524
00:32:03,690 --> 00:32:05,000
plus delta t.
525
00:32:05,000 --> 00:32:07,440
Notice this is the equilibrium
equation at time
526
00:32:07,440 --> 00:32:08,780
t plus delta t.
527
00:32:08,780 --> 00:32:11,470
This is quite different from
what we're doing in explicit
528
00:32:11,470 --> 00:32:12,300
time integration.
529
00:32:12,300 --> 00:32:15,050
In explicit time integration
we are looking at the
530
00:32:15,050 --> 00:32:18,940
equilibrium equation at time t
to obtain the solution at time
531
00:32:18,940 --> 00:32:20,340
t plus delta t.
532
00:32:20,340 --> 00:32:23,680
Here we are looking at the
equilibrium equation, we are
533
00:32:23,680 --> 00:32:27,630
using the equilibrium equation
at time t plus delta t to
534
00:32:27,630 --> 00:32:30,040
obtain the solution at
time t plus delta t.
535
00:32:30,040 --> 00:32:36,470
And that means that we have to
deal with a stiffness matrix.
536
00:32:36,470 --> 00:32:40,300
And that stiffness matrix here
is given as tk, corresponding
537
00:32:40,300 --> 00:32:43,230
to a modified Newton
iteration.
538
00:32:43,230 --> 00:32:47,590
Notice that this equation will
be solved to calculate an
539
00:32:47,590 --> 00:32:49,200
increment and displacement.
540
00:32:49,200 --> 00:32:52,280
And that increment in
displacement is added to the
541
00:32:52,280 --> 00:32:54,260
displacement vector that we
had from the previous
542
00:32:54,260 --> 00:32:59,950
iteration to obtain an updated
displacement vector.
543
00:32:59,950 --> 00:33:03,670
This is one equation in the
unknown displacements,
544
00:33:03,670 --> 00:33:07,950
velocities, and accelerations
at time t plus delta t.
545
00:33:07,950 --> 00:33:13,300
And we need two more
equations to solve
546
00:33:13,300 --> 00:33:15,200
for these three unknowns.
547
00:33:15,200 --> 00:33:19,750
In the trapezoidal rule, which
as I mentioned just now, is
548
00:33:19,750 --> 00:33:23,780
very widely used for implicit
time integration, these are
549
00:33:23,780 --> 00:33:31,020
the basic equations used for
displacements and for
550
00:33:31,020 --> 00:33:32,750
velocities.
551
00:33:32,750 --> 00:33:39,065
If we rewrite these equations,
we directly obtain the
552
00:33:39,065 --> 00:33:43,480
velocities and accelerations
in this form.
553
00:33:43,480 --> 00:33:46,710
Notice there is here the
increment in displacement from
554
00:33:46,710 --> 00:33:49,320
time t to time t plus delta t.
555
00:33:49,320 --> 00:33:52,370
Similarly here.
556
00:33:52,370 --> 00:33:55,990
In nonlinear analysis, we
are iterating for this
557
00:33:55,990 --> 00:34:01,420
displacement here, and therefore
we rewrite these
558
00:34:01,420 --> 00:34:03,240
equations a bit.
559
00:34:03,240 --> 00:34:05,790
And that's being shown,
that rewriting is
560
00:34:05,790 --> 00:34:07,870
shown on this viewgraph.
561
00:34:07,870 --> 00:34:11,150
Notice that we now have
the superscript k here
562
00:34:11,150 --> 00:34:15,679
corresponding to the iteration
k, and that we have split the
563
00:34:15,679 --> 00:34:19,870
displacement at time t plus
delta t into a value that we
564
00:34:19,870 --> 00:34:23,760
know, corresponding to iteration
k minus 1, and an
565
00:34:23,760 --> 00:34:28,250
increment that is unknown,
delta uk.
566
00:34:28,250 --> 00:34:30,739
We do the same for the
acceleration here.
567
00:34:33,400 --> 00:34:36,480
Of course, these terms
are as shown on
568
00:34:36,480 --> 00:34:38,760
the previous viewgraph.
569
00:34:38,760 --> 00:34:44,030
And now we have our velocity and
acceleration corresponding
570
00:34:44,030 --> 00:34:52,130
to iteration k, in terms of
known values and one unknown
571
00:34:52,130 --> 00:34:56,239
value, the incremental
displacement.
572
00:34:56,239 --> 00:35:00,030
We can now use these two
equations and substitute into
573
00:35:00,030 --> 00:35:03,950
our equilibrium equation the
equation that you saw on the
574
00:35:03,950 --> 00:35:07,640
first viewgraph, when we started
to talk about the
575
00:35:07,640 --> 00:35:09,710
implicit time integration
scheme.
576
00:35:09,710 --> 00:35:13,110
And the result of that
substitution is
577
00:35:13,110 --> 00:35:14,510
given on this viewgraph.
578
00:35:14,510 --> 00:35:20,010
We have a coefficient matrix
here, which we call tk hat
579
00:35:20,010 --> 00:35:21,590
times the incremental
displacement
580
00:35:21,590 --> 00:35:24,110
vector, which is unknown.
581
00:35:24,110 --> 00:35:27,260
And on the right hand side of
that equation, we have the
582
00:35:27,260 --> 00:35:31,900
externally applied loads,
the nodal point forces
583
00:35:31,900 --> 00:35:34,690
corresponding to the internal
element stresses at time t
584
00:35:34,690 --> 00:35:40,160
plus delta t, and at the end
of iteration k minus one.
585
00:35:40,160 --> 00:35:43,180
To evaluate these, we need the
displacements corresponding to
586
00:35:43,180 --> 00:35:45,570
iteration k minus 1.
587
00:35:45,570 --> 00:35:50,240
And here we have inertia
contributions and damping
588
00:35:50,240 --> 00:35:51,330
contributions.
589
00:35:51,330 --> 00:35:55,290
But notice all this information
here is known once
590
00:35:55,290 --> 00:35:58,110
you know the displacements
corresponding to t plus delta
591
00:35:58,110 --> 00:36:02,300
tu k minus 1.
592
00:36:02,300 --> 00:36:07,150
We iterate on this equation
until we satisfy certain
593
00:36:07,150 --> 00:36:12,440
equilibrium criteria,
convergence criteria, and much
594
00:36:12,440 --> 00:36:15,650
the same way as we do it in
static nonlinear analysis, and
595
00:36:15,650 --> 00:36:18,550
the way we discussed it
in an earlier lecture.
596
00:36:18,550 --> 00:36:21,490
However, let us make
a few observations.
597
00:36:21,490 --> 00:36:24,890
First of all, we notice that as
delta t gets smaller, the
598
00:36:24,890 --> 00:36:29,410
entries in this tk hat matrix,
this effective stiffness
599
00:36:29,410 --> 00:36:31,030
matrix, increase.
600
00:36:31,030 --> 00:36:32,040
Why do they increase?
601
00:36:32,040 --> 00:36:38,570
Because we have delta t squared
term, 1 over delta t
602
00:36:38,570 --> 00:36:42,490
squared term in front
of the mass matrix.
603
00:36:42,490 --> 00:36:45,430
And the 1 over 2 delta
t term in front
604
00:36:45,430 --> 00:36:46,720
of the damping matrix.
605
00:36:46,720 --> 00:36:54,440
As delta t becomes smaller, 1
over delta t squared becomes
606
00:36:54,440 --> 00:36:57,830
larger and that means that the
entries in this matrix here
607
00:36:57,830 --> 00:37:01,660
become larger as delta
t gets smaller.
608
00:37:01,660 --> 00:37:04,230
The convergence characteristics
of the
609
00:37:04,230 --> 00:37:06,700
equilibrium iterations
are better
610
00:37:06,700 --> 00:37:08,820
than in static analysis.
611
00:37:08,820 --> 00:37:13,200
The reason being, that the
entries in tk hat are
612
00:37:13,200 --> 00:37:17,060
generally larger, much larger
than what we see in static
613
00:37:17,060 --> 00:37:20,630
analysis because of this 1 over
delta t squared term, and
614
00:37:20,630 --> 00:37:24,990
1 over 2 delta t term, that
I just mentioned.
615
00:37:24,990 --> 00:37:27,350
We also notice that the
trapezoidal rule is
616
00:37:27,350 --> 00:37:31,850
unconditionally stable
in linear analysis.
617
00:37:31,850 --> 00:37:35,840
The analysis that shows that the
method is unconditionally
618
00:37:35,840 --> 00:37:42,450
stable is really performed on a
linear system, and therefore
619
00:37:42,450 --> 00:37:44,450
the unconditional stability
strictly
620
00:37:44,450 --> 00:37:46,190
holds only linear analysis.
621
00:37:46,190 --> 00:37:49,680
For nonlinear analysis, we want
to select delta t for
622
00:37:49,680 --> 00:37:53,575
accuracy, and we want to select
delta t for convergence
623
00:37:53,575 --> 00:37:55,630
of the iteration.
624
00:37:55,630 --> 00:37:59,910
If we find, for example, that
we have convergence
625
00:37:59,910 --> 00:38:08,260
difficulties in the iteration,
then we may want to switch to
626
00:38:08,260 --> 00:38:09,750
a more powerful scheme.
627
00:38:09,750 --> 00:38:12,260
That's one approach, of course,
just the same way as
628
00:38:12,260 --> 00:38:15,420
we are doing it in
static analysis.
629
00:38:15,420 --> 00:38:22,690
Or we may also recall, or we
may identify actually, that
630
00:38:22,690 --> 00:38:27,680
delta t is too large, and we
should decrease delta t.
631
00:38:27,680 --> 00:38:31,680
In particular, as a rule of
thumb, I'd just like to share
632
00:38:31,680 --> 00:38:33,010
some experience with you.
633
00:38:33,010 --> 00:38:36,580
We find that if you use the
BFGS method, which we
634
00:38:36,580 --> 00:38:40,000
discussed earlier in an earlier
lecture, if you use
635
00:38:40,000 --> 00:38:44,010
the BFGS method for the
equilibrium iteration, and if
636
00:38:44,010 --> 00:38:47,190
we have difficulties converging
there, then very
637
00:38:47,190 --> 00:38:50,830
frequently the time
step is too large,
638
00:38:50,830 --> 00:38:53,190
and should be decreased.
639
00:38:53,190 --> 00:38:57,870
Not just to obtain convergence
in the BFGS iterations, but to
640
00:38:57,870 --> 00:39:00,200
obtain an appropriate
accuracy.
641
00:39:00,200 --> 00:39:05,060
So we really look at the
difficulties obtaining
642
00:39:05,060 --> 00:39:09,440
convergence in the iterations as
also a signal to us to tell
643
00:39:09,440 --> 00:39:13,630
us that our time step delta
t is probably too large.
644
00:39:16,430 --> 00:39:21,980
The convergence criteria that we
are using are basically the
645
00:39:21,980 --> 00:39:25,730
same that we have seen already
in the static analysis in an
646
00:39:25,730 --> 00:39:27,060
earlier lecture.
647
00:39:27,060 --> 00:39:32,440
In static analysis we only
included this term here.
648
00:39:32,440 --> 00:39:37,660
In dynamic analysis we include
now the mass and damping terms
649
00:39:37,660 --> 00:39:40,230
in our energy criterion.
650
00:39:40,230 --> 00:39:45,650
We discussed this criterion
earlier, of course not
651
00:39:45,650 --> 00:39:48,200
including the mass and
damping effects.
652
00:39:48,200 --> 00:39:50,310
Now we do.
653
00:39:50,310 --> 00:39:57,490
Another convergence criterion is
one on forces, and earlier
654
00:39:57,490 --> 00:40:01,180
we had only this term here.
655
00:40:01,180 --> 00:40:05,460
Now we include also the inertia
and the damping terms.
656
00:40:07,990 --> 00:40:13,680
Once again, if we are only
having translational degrees
657
00:40:13,680 --> 00:40:19,900
of freedom, then here we have
only nodal point forces, nodal
658
00:40:19,900 --> 00:40:21,960
point forces here.
659
00:40:21,960 --> 00:40:26,380
And our norm would be having
the dimensions of a nodal
660
00:40:26,380 --> 00:40:27,990
point force.
661
00:40:27,990 --> 00:40:32,770
If we have rotational degrees
of freedom, then we still
662
00:40:32,770 --> 00:40:36,290
include here only the
nodal point forces.
663
00:40:36,290 --> 00:40:40,510
We extract, in other words, the
components corresponding
664
00:40:40,510 --> 00:40:44,230
to the translational degrees of
freedom of this vector and
665
00:40:44,230 --> 00:40:49,050
that vector, and include those
components only up here.
666
00:40:49,050 --> 00:40:53,650
We use RNORM that also has the
units of the nodal point
667
00:40:53,650 --> 00:40:59,430
force, and we supplement this
convergence criteria by
668
00:40:59,430 --> 00:41:04,580
another quotient, where we have
on top the components
669
00:41:04,580 --> 00:41:06,720
corresponding to the rotational
degrees of freedom
670
00:41:06,720 --> 00:41:10,030
only, and in the bottom,
what we called
671
00:41:10,030 --> 00:41:13,000
RMNORM already earlier.
672
00:41:13,000 --> 00:41:17,140
The point is that we want to
have up here only those
673
00:41:17,140 --> 00:41:23,070
components that carry the same
units as we have down here.
674
00:41:23,070 --> 00:41:28,080
Once for translations,
once for rotations.
675
00:41:28,080 --> 00:41:31,480
And of course recall we're
talking here about the
676
00:41:31,480 --> 00:41:32,190
Euclidean norm.
677
00:41:32,190 --> 00:41:34,180
which is defined
as shown here.
678
00:41:34,180 --> 00:41:37,440
We went through this already
earlier, when we talked about
679
00:41:37,440 --> 00:41:40,740
the convergence criteria
in static analysis.
680
00:41:40,740 --> 00:41:44,770
We have also the convergence
criterion on displacements.
681
00:41:44,770 --> 00:41:47,290
Notice this is what
we saw already
682
00:41:47,290 --> 00:41:49,050
earlier in static analysis.
683
00:41:49,050 --> 00:41:53,690
Once again, we use displacement
components here
684
00:41:53,690 --> 00:41:55,940
to compare with the
displacement.
685
00:41:55,940 --> 00:41:58,680
We introduced all of the
displacement components here
686
00:41:58,680 --> 00:42:02,370
to compare with this
displacement and if we have
687
00:42:02,370 --> 00:42:04,610
rotational degrees of
freedom, they would
688
00:42:04,610 --> 00:42:07,060
not be included here.
689
00:42:07,060 --> 00:42:09,670
Their components would
not be included here.
690
00:42:09,670 --> 00:42:13,290
When you use DNORM here, they
would, however, in a second
691
00:42:13,290 --> 00:42:17,960
check be included here and only
included here, and then
692
00:42:17,960 --> 00:42:19,710
we would use here DMNORM.
693
00:42:19,710 --> 00:42:24,120
We went through that
earlier already.
694
00:42:24,120 --> 00:42:30,180
Some modelling hints regarding
the use of the implicit time
695
00:42:30,180 --> 00:42:32,310
integration scheme.
696
00:42:32,310 --> 00:42:36,830
The integration scheme is
generally used to analyze
697
00:42:36,830 --> 00:42:40,120
structural vibration
problems, and we
698
00:42:40,120 --> 00:42:42,440
would proceed as follows.
699
00:42:42,440 --> 00:42:47,540
We identify the frequencies
contained in the loading by a
700
00:42:47,540 --> 00:42:49,770
Fourier decomposition.
701
00:42:49,770 --> 00:42:53,340
We chose a finite element mesh
that can accurately represent
702
00:42:53,340 --> 00:42:57,600
the static response and all
important frequencies.
703
00:42:57,600 --> 00:43:01,550
And then we perform a direct
integration using this time
704
00:43:01,550 --> 00:43:08,780
step here where tco is the
smallest period in seconds
705
00:43:08,780 --> 00:43:11,110
that we want to integrate.
706
00:43:11,110 --> 00:43:15,540
Now there's a lot of information
here, and I'd like
707
00:43:15,540 --> 00:43:20,520
to refer you to an example,
analysis of a tower that I
708
00:43:20,520 --> 00:43:23,330
will be presenting to you in
the next lecture, where we
709
00:43:23,330 --> 00:43:28,200
will talk about these items
quite a bit more.
710
00:43:28,200 --> 00:43:32,700
But this is basically the
approach that we use, to model
711
00:43:32,700 --> 00:43:37,520
the problem and then also to
perform the direct integration
712
00:43:37,520 --> 00:43:41,480
of the governing finite
element equations.
713
00:43:41,480 --> 00:43:45,120
As I said earlier already, the
method is used primarily for
714
00:43:45,120 --> 00:43:47,490
structural vibration problems.
715
00:43:47,490 --> 00:43:51,680
It is effective, quite
effective, using also
716
00:43:51,680 --> 00:43:53,930
higher-order elements.
717
00:43:53,930 --> 00:43:58,970
And if we do use higher-order
elements, we do find
718
00:43:58,970 --> 00:44:03,860
frequently that method, or the
total solution should be
719
00:44:03,860 --> 00:44:06,470
performed most effectively
using a
720
00:44:06,470 --> 00:44:08,680
consistent mass matrix.
721
00:44:08,680 --> 00:44:11,180
I'm thinking here in particular,
if you use very
722
00:44:11,180 --> 00:44:15,030
high-order elements such as
cubic shall elements,
723
00:44:15,030 --> 00:44:18,420
isoparametric cubic shell
elements, an element that we
724
00:44:18,420 --> 00:44:22,150
will talk about in a later
lecture, then with such
725
00:44:22,150 --> 00:44:25,720
element you probably want to use
a consistent mass matrix
726
00:44:25,720 --> 00:44:29,650
for the overall solution
of the problem.
727
00:44:29,650 --> 00:44:32,910
Notice that a structural
vibration problem can be
728
00:44:32,910 --> 00:44:36,650
thought of as a static problem
including inertia forces.
729
00:44:36,650 --> 00:44:42,760
And we already have noted that
in structural analysis of
730
00:44:42,760 --> 00:44:48,150
static problems, the
higher-order elements do quite
731
00:44:48,150 --> 00:44:50,490
well, and sometimes much better
732
00:44:50,490 --> 00:44:52,580
than the lower elements.
733
00:44:52,580 --> 00:44:56,730
And for this reason, they also
do quite well in dynamic
734
00:44:56,730 --> 00:45:00,860
analysis of structural
vibration problems.
735
00:45:00,860 --> 00:45:05,680
A typical problem that one might
like to solve using this
736
00:45:05,680 --> 00:45:07,710
implicit time integration
scheme, is
737
00:45:07,710 --> 00:45:09,920
schematically shown here.
738
00:45:09,920 --> 00:45:14,490
We have a tower that is
subjected to a blast loading
739
00:45:14,490 --> 00:45:18,340
as shown here schematically,
and we assume that only the
740
00:45:18,340 --> 00:45:20,550
structural vibrations
of this tower are
741
00:45:20,550 --> 00:45:21,950
required to be solved.
742
00:45:21,950 --> 00:45:26,210
In this case, we may very well
obtain the response of the
743
00:45:26,210 --> 00:45:31,250
tower using something like about
100 steps to integrate
744
00:45:31,250 --> 00:45:32,370
the response.
745
00:45:32,370 --> 00:45:36,330
In fact, we will look at a
problem of this nature in the
746
00:45:36,330 --> 00:45:40,080
next lecture, where we do
actually analyze a tower using
747
00:45:40,080 --> 00:45:41,920
the implicit time integration
scheme.
748
00:45:41,920 --> 00:45:45,640
And where we talk a little bit
about how the tower would be
749
00:45:45,640 --> 00:45:49,190
modelled, and what one would
look for to make sure that you
750
00:45:49,190 --> 00:45:52,720
have adequately calculated
the response.
751
00:45:52,720 --> 00:45:56,510
Finally, I'd like to also
point out that it can be
752
00:45:56,510 --> 00:46:02,330
effective to use the combined
technique of explicit and
753
00:46:02,330 --> 00:46:04,430
implicit integration.
754
00:46:04,430 --> 00:46:08,960
You might, for example, have a
problem that shows an initial
755
00:46:08,960 --> 00:46:12,040
wave propagation and then
a structural vibration.
756
00:46:12,040 --> 00:46:16,190
In this case, you might first
want to use the explicit time
757
00:46:16,190 --> 00:46:18,860
integration to calculate the
initial wave propagation
758
00:46:18,860 --> 00:46:23,380
response of the structure, and
then restart to use an
759
00:46:23,380 --> 00:46:26,850
implicit time integration
technique to calculate the
760
00:46:26,850 --> 00:46:29,960
vibrational response
of the structure.
761
00:46:29,960 --> 00:46:33,440
So in that case, you would
basically combine explicit and
762
00:46:33,440 --> 00:46:36,510
implicit time integration, but
the combination would be first
763
00:46:36,510 --> 00:46:39,550
explicit, and then implicit
time integration.
764
00:46:39,550 --> 00:46:45,570
Another kind of combination is
quite effective when you have
765
00:46:45,570 --> 00:46:49,370
to analyze problems where you
have very stiff regions and
766
00:46:49,370 --> 00:46:51,970
very flexible regions.
767
00:46:51,970 --> 00:46:55,550
For the flexible regions of the
problem, you might want to
768
00:46:55,550 --> 00:46:58,650
use explicit time integration,
and for the stiff region you
769
00:46:58,650 --> 00:47:01,320
might want to use implicit
time integration.
770
00:47:01,320 --> 00:47:05,020
A typical analysis, for example,
would be a fluid
771
00:47:05,020 --> 00:47:08,310
structure problem, where the
fluid is very soft and the
772
00:47:08,310 --> 00:47:11,160
structure is very stiff.
773
00:47:11,160 --> 00:47:13,660
Using the implicit time
integration for the structure
774
00:47:13,660 --> 00:47:16,830
and the explicit time
integration for the fluid can
775
00:47:16,830 --> 00:47:17,600
be effective.
776
00:47:17,600 --> 00:47:20,850
Of course, then you would also
have to couple these time
777
00:47:20,850 --> 00:47:23,940
integration schemes properly
at the boundary between the
778
00:47:23,940 --> 00:47:26,560
fluid and the structure.
779
00:47:26,560 --> 00:47:29,190
This brings me then to the end
of what I wanted to say in
780
00:47:29,190 --> 00:47:29,980
this lecture.
781
00:47:29,980 --> 00:47:33,300
In the next lecture we will
continue our discussion of
782
00:47:33,300 --> 00:47:37,230
solution methods for the
solution of the nonlinear
783
00:47:37,230 --> 00:47:39,180
dynamic equilibrium equations.
784
00:47:39,180 --> 00:47:42,570
And we will also show examples,
examples that
785
00:47:42,570 --> 00:47:45,580
demonstrate how the techniques
that I discussed in this
786
00:47:45,580 --> 00:47:48,750
lecture are used to
solve problems.
787
00:47:48,750 --> 00:47:50,160
Thank you very much for
your attention.