1
00:00:14,450 --> 00:00:20,350
We'll continue with our discussion on methods
for integrating equations of equilibrium in
2
00:00:20,350 --> 00:00:28,080
dynamics. In today's lecture we will talk
about the concept of energy conservation and
3
00:00:28,080 --> 00:00:34,290
also address few issues related to nonlinear
systems, so what we have been discussing we
4
00:00:34,290 --> 00:00:38,300
have
developed, discussed this Newmark's method
5
00:00:38,300 --> 00:00:46,250
which is a second-order implicit method, the
three stage formulary for this is displayed
6
00:00:46,250 --> 00:00:52,610
here, we have shown that the method is conditionally
stable if we select this parameter delta and
7
00:00:52,610 --> 00:00:58,810
alpha to satisfy this requirement, and these
requirements are independent of system natural
8
00:00:58,810 --> 00:01:05,059
frequency and damping. So Newmark's method
is implicit, self-starting, and it is a single
9
00:01:05,059 --> 00:01:11,149
step method.
Now we can also obtain the other methods like
10
00:01:11,149 --> 00:01:15,530
average acceleration method, linear acceleration
method, as special cases of Newmark's method
11
00:01:15,530 --> 00:01:22,810
by assigning the parameters alpha and delta
values as shown here, this also we have seen.
12
00:01:22,810 --> 00:01:29,969
Now for delta greater than half the method
displays high frequency dissipation characteristic
13
00:01:29,969 --> 00:01:39,170
which are desirable, but the global error
would be of order delta T which is not desirable.
14
00:01:39,170 --> 00:01:47,350
We have also seen questions about error estimates
and convergence of forward Euler method, we
15
00:01:47,350 --> 00:01:55,689
have shown that the error is of the, is bounded
by this quantity, so we saw that local, although
16
00:01:55,689 --> 00:02:02,149
the local error is order delta T square the
global error is order delta T, this is a common
17
00:02:02,149 --> 00:02:06,789
feature in these algorithms intuitively the
explanation for this can be this, global error
18
00:02:06,789 --> 00:02:14,120
at nth step, if we write it as, if EN is a
local error the global error at the end of
19
00:02:14,120 --> 00:02:20,590
nth step will be N into EN, so N can be written
as TN/delta T, so if EN is of order delta
20
00:02:20,590 --> 00:02:26,220
T square as it is here when divided by delta
T it becomes order delta T, that means you
21
00:02:26,220 --> 00:02:31,159
have to sum up to all the errors up to that
time instant that reduces the order of accuracy,
22
00:02:31,159 --> 00:02:35,930
so global error typically would be an order
less than the local error.
23
00:02:35,930 --> 00:02:47,959
The HHT alpha method and generalizations overcome
one weakness of Newmark's method that to be,
24
00:02:47,959 --> 00:02:55,370
to possess high frequency numerical dissipation
the method would become only order delta T
25
00:02:55,370 --> 00:03:05,430
accurate, so what this HHT alpha method does
is, it introduces additional algorithmic parameters,
26
00:03:05,430 --> 00:03:13,430
for example in a generalized version we have
alpha M, alpha F, which are additional parameters,
27
00:03:13,430 --> 00:03:18,200
in addition to delta and alpha that are used
in Newmark's method, so we have now 4 parameters,
28
00:03:18,200 --> 00:03:25,319
so we have greater flexibility to you know
design the algorithm to display desired properties,
29
00:03:25,319 --> 00:03:30,300
so it can be, we have seen that the method
possesses desirable high frequency dissipation
30
00:03:30,300 --> 00:03:36,569
characteristics it is unconditionally stable
and it has second-order accuracy.
31
00:03:36,569 --> 00:03:41,370
So we can summarize what are the desirable
features of a numerical integration scheme,
32
00:03:41,370 --> 00:03:47,860
so it should at least have second-order accuracy
and it should be unconditionally stable when
33
00:03:47,860 --> 00:03:52,730
applied to linear time-invariant systems,
and we want to have controllable algorithmic
34
00:03:52,730 --> 00:04:00,340
damping in higher modes that means by changing
parameter other than the step size, that means
35
00:04:00,340 --> 00:04:05,120
these parameters are algorithmic parameter
not the system parameters we should be able
36
00:04:05,120 --> 00:04:10,799
to ensure that the numerical damping that
will be present in the system you know does
37
00:04:10,799 --> 00:04:15,930
not distort the lower order contributing modes,
and it suppresses the spurious higher-order
38
00:04:15,930 --> 00:04:23,250
modes, so again this requires investigation
into spectral radius of the amplification
39
00:04:23,250 --> 00:04:28,180
matrix as frequency tends to infinity, for
large frequency the spectral radius should
40
00:04:28,180 --> 00:04:33,640
go to zero so that higher order, high-frequency
components are dissipated and the spectral
41
00:04:33,640 --> 00:04:39,480
radius should be as close to 1 for low frequency
region where there will be participating modes
42
00:04:39,480 --> 00:04:44,820
which are contributing significantly to the
response, and also there should be no overshoot
43
00:04:44,820 --> 00:04:49,480
that means excessive oscillations during the
first few steps should not be present ideally
44
00:04:49,480 --> 00:04:54,320
the scheme should be self-starting and when
we deal with nonlinear systems as we will
45
00:04:54,320 --> 00:05:01,160
do today later, no more than one set of implicit
equations need to be solved at each step,
46
00:05:01,160 --> 00:05:07,370
so the concerns here are let me summarize
if you are carrying out a dynamic analysis
47
00:05:07,370 --> 00:05:13,250
of a structure you make a finite element model
while selecting the size of the model you
48
00:05:13,250 --> 00:05:17,570
have to keep in mind the highest frequency
that is participating, highest frequency that
49
00:05:17,570 --> 00:05:24,140
is present in the excitation then suppose
that highest frequency is say 20 hertz as
50
00:05:24,140 --> 00:05:29,800
perhaps in the case of an earthquake engineering
problem then we should ensure that the finite
51
00:05:29,800 --> 00:05:37,170
element model that we make should have acceptable
accuracies on Eigen characteristics maybe
52
00:05:37,170 --> 00:05:43,730
up to say 50 hertz or 80 hertz 3 to 4 times
the highest frequency present in the excitation,
53
00:05:43,730 --> 00:05:49,120
over that frequency range we should capture
all the modes that are present in the structure
54
00:05:49,120 --> 00:05:57,180
with good accuracy, suppose between say 0
to 50 hertz there are say 12 modes in a structure
55
00:05:57,180 --> 00:06:03,031
which contributes, whose natural frequency
lie in 0 to 50 hertz say for example, then
56
00:06:03,031 --> 00:06:07,660
our model should have at least, the degrees
of freedom in our model should be at least
57
00:06:07,660 --> 00:06:12,630
10 times the number of modes that you would
like to retain in the model expansion.
58
00:06:12,630 --> 00:06:21,570
So in order that we get acceptable accuracy
on the lower order model characteristics,
59
00:06:21,570 --> 00:06:27,100
our model size need to be large, consequently
there will be higher order spurious modes
60
00:06:27,100 --> 00:06:33,410
which are numerically not accurate and in
a physical sense they don't really contribute
61
00:06:33,410 --> 00:06:38,350
to the response, so the algorithm that we
use should ensure that those spurious modes
62
00:06:38,350 --> 00:06:42,920
are dissipated through numerical dissipation,
and the contributing lower order modes are
63
00:06:42,920 --> 00:06:51,450
not distorted, so that is the main issue when
it comes to discussion on algorithmic dissipation.
64
00:06:51,450 --> 00:06:59,300
Now, we will now consider what is the relationship
between stability and energy conservation?
65
00:06:59,300 --> 00:07:05,350
We have talked about the so-called spectral
stability, you know the integration schemes
66
00:07:05,350 --> 00:07:10,860
unapplied to linear time-invariant systems,
let's ask the question now, those concerns
67
00:07:10,860 --> 00:07:16,250
are related to conservation of energy.
So let's focus our discussion on Newmark's
68
00:07:16,250 --> 00:07:22,220
method, so in the Newmark method we have this
set of 3 formulae, this is for velocity, this
69
00:07:22,220 --> 00:07:27,720
is for displacement, and this is the equilibrium
equation at time TN+1, now for purpose of
70
00:07:27,720 --> 00:07:33,800
discussion let us take delta to be 0.5, and
alpha to be 0.25, we can do a general discussion
71
00:07:33,800 --> 00:07:40,150
but let us simplify our study a bit, we will
use this representation, so in which case
72
00:07:40,150 --> 00:07:48,280
UN, velocity is given by this, and displacement
would be given by this. Now we can reorganize
73
00:07:48,280 --> 00:07:54,070
some of these terms, suppose I write UN+ 1
dot - UN dot, this
74
00:07:54,070 --> 00:08:06,550
is given by this. Similarly UN + 1 - UN is
given by this. Now, this we can simplify and
75
00:08:06,550 --> 00:08:18,150
we can show that UN+1 - UN is given by this.
Now we have the equilibrium equation at time
76
00:08:18,150 --> 00:08:26,940
T, now what I will do is I will multiply by
U dot transpose, pre multiply by U dot transpose
77
00:08:26,940 --> 00:08:34,550
so this equation is what we get. Now if U
dot transpose into F(t) is actually the work
78
00:08:34,550 --> 00:08:39,320
done by external force per unit time, for
example in a single degree freedom system
79
00:08:39,320 --> 00:08:44,330
velocity into force will be work done per
unit time, force into displacement is a work
80
00:08:44,330 --> 00:08:49,020
done, therefore force into velocity will be
work done per unit time which is the power
81
00:08:49,020 --> 00:08:54,440
input, so this is basically the power balance
equation where we are equating the power input
82
00:08:54,440 --> 00:09:02,370
to the system to the contributions from inertial
forces, stiffness and damping forces. Now
83
00:09:02,370 --> 00:09:08,700
that sum of kinetic energy and potential energy
will now consider that it will be T+V, and
84
00:09:08,700 --> 00:09:15,950
T is 1/2 U dot transpose MU dot, V is 1/2
U transpose KU this we have seen. Now the
85
00:09:15,950 --> 00:09:21,260
time variation of the total energy that D/DT
(T+V) if you are interested we can differentiate
86
00:09:21,260 --> 00:09:28,710
this, so we will get U double dot transpose
MU dot + 1/2 U dot transpose MU double dot,
87
00:09:28,710 --> 00:09:35,140
similarly the other two terms.
Now this is a scalar quantity U dot transpose
88
00:09:35,140 --> 00:09:41,380
KU is a scalar quantity, so is U transpose
KU dot, U is N cross 1, K is N cross N, this
89
00:09:41,380 --> 00:09:48,530
is N cross 1,. so U transpose will be 1 cross
N, so this will be a scalar quantity, so for
90
00:09:48,530 --> 00:09:53,810
this we can as well write U dot transpose
KU, you can't replace it by this transpose,
91
00:09:53,810 --> 00:09:59,430
so we can or in other words these two can
be added and these two can be added and we
92
00:09:59,430 --> 00:10:09,740
can write in this form. So the time variation
of total energy is therefore given by this.
93
00:10:09,740 --> 00:10:25,850
Now we will now see if I now integrate this
you know T+V from TN to TN+1 we will get this
94
00:10:25,850 --> 00:10:37,640
integral. Now at N+1 the kinetic energy will
be given by this, and at TN I have this expression
95
00:10:37,640 --> 00:10:42,840
with you replace N+1/N, I get this. So if
I now consider the increment in kinetic energy
96
00:10:42,840 --> 00:10:51,410
this is the increment in kinetic, so now I
will simplify this, I am adding and subtracting
97
00:10:51,410 --> 00:10:57,700
some terms and rearranging these terms to
show that this incrementing kinetic energy
98
00:10:57,700 --> 00:11:07,980
can be given by this quantity.
Now similarly we get VN+1 - VN as this quantity,
99
00:11:07,980 --> 00:11:15,660
so we have TN+1 -TN as this, VN+1 - VN as
this, and increments in velocity and displacements
100
00:11:15,660 --> 00:11:20,560
are given by this, so these terms are here
now this is according to Newmark's method,
101
00:11:20,560 --> 00:11:25,340
so now I will substitute these terms into
the expression for increment in kinetic energy
102
00:11:25,340 --> 00:11:31,480
and potential energy, there is a Newmark beta,
Newmark's approximation, so upon doing that
103
00:11:31,480 --> 00:11:39,940
I get this as the increment in kinetic energy,
this is the increment in strain energy, so
104
00:11:39,940 --> 00:11:47,550
consequently the D/DT of total energy is given
by this.
105
00:11:47,550 --> 00:11:56,779
Now therefore if we now recall MU double dot
+ U dot + K = F from this we can write this
106
00:11:56,779 --> 00:12:09,971
expression, okay that means I am looking for
T+V at TN+1 and TN, I get this, so this is
107
00:12:09,971 --> 00:12:13,320
the
expression we have got, and now this is again
108
00:12:13,320 --> 00:12:18,240
Newmark beta approximation for displacement
if I now substitute that I get the change
109
00:12:18,240 --> 00:12:26,120
in total energy from time step TN to TN+1
is given by this.
110
00:12:26,120 --> 00:12:33,480
Now if we talk about free vibration F will
be 0, and C will be 0, undamped free vibration,
111
00:12:33,480 --> 00:12:41,130
then I see that the change in total energy
is 0, which is what we expect from a conservative
112
00:12:41,130 --> 00:12:47,240
system, so for free vibration of undamped
system that is C = 0 and F = 0 we get T +
113
00:12:47,240 --> 00:12:53,370
V, the change in T+V from TN to TN+1 is 0,
that is the Newmark algorithm conserves the
114
00:12:53,370 --> 00:12:59,430
total energy, this is true for any time step
size.
115
00:12:59,430 --> 00:13:04,880
Now we have seen that the method is unconditionally
stable if this condition is satisfied, these
116
00:13:04,880 --> 00:13:09,420
requirements are independent of system natural
frequency and damping that also we have seen,
117
00:13:09,420 --> 00:13:15,870
we can show that for this case the energy
balance or one time steps is negative, okay
118
00:13:15,870 --> 00:13:24,670
in general when there is a damping as well
as forcing, in regions in the alpha delta
119
00:13:24,670 --> 00:13:29,560
space where the method is not stable the energy
balance or a time step turns out to be positive,
120
00:13:29,560 --> 00:13:32,980
please bear in mind that the formulary I have
derived is for specific choice of delta and
121
00:13:32,980 --> 00:13:40,430
alpha, so you could re-derive the whole story
by keeping delta and alpha as it is and you
122
00:13:40,430 --> 00:13:47,589
can see how the question about change in energy
can be tackled.
123
00:13:47,589 --> 00:13:52,770
So some quick illustrations a single degree
freedom system with natural frequency 2 pi,
124
00:13:52,770 --> 00:13:58,770
this is damping is 0 free vibration problems
were considering, and I am considering step
125
00:13:58,770 --> 00:14:04,650
size of one tenth of critical step size and
this is the algorithmic parameters, so blue
126
00:14:04,650 --> 00:14:10,000
line here is the exact solution, and the red
one is a Newmark solution, obviously it's
127
00:14:10,000 --> 00:14:15,232
not accurate but it is stable. So if I now
plot the energy balance we see that this is
128
00:14:15,232 --> 00:14:23,920
10 to the power of -14 so you can see that
the energy balance criteria is satisfied.
129
00:14:23,920 --> 00:14:31,730
Similarly if I now take same data with higher
step size now T critical is delta T is T critical/100
130
00:14:31,730 --> 00:14:37,610
so we see that the exact and the Newmark solution
match very well and the error is again quite
131
00:14:37,610 --> 00:14:45,350
small, the energy balance. Now how about damped
free vibration? So I introduce 3% damping
132
00:14:45,350 --> 00:14:51,470
and again start with T critical/10 although
the solution is not accurate, the error is
133
00:14:51,470 --> 00:14:58,589
still that the energy balance is satisfied,
so energy balance is satisfied it does not
134
00:14:58,589 --> 00:15:02,460
mean that you have got accurate answers, okay
it's like telling if you are getting stable
135
00:15:02,460 --> 00:15:10,339
solution you need not ensure that the solution
has acceptable accuracy. Now if I increase
136
00:15:10,339 --> 00:15:17,260
step size to T critical/100 I get you
know the two exact solution and Newmark solution
137
00:15:17,260 --> 00:15:27,560
math, and this is the behavior of the error.
Now I deliberately select the algorithmic
138
00:15:27,560 --> 00:15:35,970
parameter in regions where the method is not
stable, so we see that the Newmark solution
139
00:15:35,970 --> 00:15:43,770
diverges as time increases, and you can see
that the energy balance condition is also
140
00:15:43,770 --> 00:15:48,459
violated, now this is, numbers are very large
is 0.15, if you look at the previous graph
141
00:15:48,459 --> 00:15:53,100
these numbers were 10 to the power of - 12
this is very large number and it is also blowing
142
00:15:53,100 --> 00:16:01,900
up.
The same issue about algorithmic energy conservation
143
00:16:01,900 --> 00:16:06,550
can also be studied for other integration
scheme, so I have shown some illustrations
144
00:16:06,550 --> 00:16:14,740
here this is the HHT alpha method with these
parameters a damped free vibration with step
145
00:16:14,740 --> 00:16:20,500
size of 0.1 second and we see that the energy
balance condition is met and it so happens
146
00:16:20,500 --> 00:16:25,730
that the step size is in the ensure stability
of the solution.
147
00:16:25,730 --> 00:16:34,230
Now with improved step size here for this
step size there are 10 points in a period
148
00:16:34,230 --> 00:16:39,200
because omega is 2 pi that means period is
1 second, I have taken 10 points in a cycle
149
00:16:39,200 --> 00:16:43,300
that does not lead to acceptable solutions.
On the other hand if you take 100 points the
150
00:16:43,300 --> 00:16:51,200
two, the exact solution and HHT alpha solution
matches and the error is acceptable.
151
00:16:51,200 --> 00:16:56,240
Now the discussion on energy conservation
in Newmark based time integration algorithms
152
00:16:56,240 --> 00:17:01,709
has been presented in this work by Krenk,
I have given the reference so in this paper
153
00:17:01,709 --> 00:17:07,851
he has derived that discussed this issue for
both Newmark's method HHT alpha and generalized
154
00:17:07,851 --> 00:17:20,169
alpha methods, okay.
Now, we'll now start discussing some issues
155
00:17:20,169 --> 00:17:28,119
about integrating equations of motion if we
are dealing with nonlinear systems, we will
156
00:17:28,119 --> 00:17:32,620
start by discussing some preliminary ideas
perhaps later in the course we will address
157
00:17:32,620 --> 00:17:39,739
these questions in a more elaborate way. Now
why do we get nonlinear systems? The sources
158
00:17:39,739 --> 00:17:45,679
of non-linearity in structural dynamics could
be due to nonlinear strain displacement relations,
159
00:17:45,679 --> 00:17:50,720
or due two non-linear stress-strain relations,
or non-linear energy dissipation mechanisms
160
00:17:50,720 --> 00:17:56,940
for example friction at joins, aerodynamic
damping, hysteresis and so on and so forth.
161
00:17:56,940 --> 00:18:03,299
Now for the purpose of illustration what we
will do is we will consider a nonlinear equation
162
00:18:03,299 --> 00:18:09,619
of motion of this form and at this stage let
us not inquire into the question on how to
163
00:18:09,619 --> 00:18:14,119
arrive at this finite element model from a
given continuum problem, that question will
164
00:18:14,119 --> 00:18:21,070
consider later. Now the question that we wish
to discuss now is if you were to get this
165
00:18:21,070 --> 00:18:29,830
type of equation of motion how will you integrate
the equation of motion, what are the issues?
166
00:18:29,830 --> 00:18:36,721
So let's see what happens to Newmark's method,
so the representation for velocity and displacement
167
00:18:36,721 --> 00:18:45,529
remains the same, and the equilibrium equation
at T = N+1 will have now this new term, so
168
00:18:45,529 --> 00:18:51,580
upon simplifying these equations from equation
2 I get if I solve for U double dot N+1 I
169
00:18:51,580 --> 00:18:59,049
get this, and if I substitute now this and
these quantities into the equation of motion
170
00:18:59,049 --> 00:19:07,119
or in the first equation I get U dot of N+1
is UN dot plus these terms, and this is what
171
00:19:07,119 --> 00:19:14,549
we got earlier also.
Now if we now go back to, we have now acceleration,
172
00:19:14,549 --> 00:19:19,510
displacement and I mean velocity and displacement
all of that if we now substitute into the
173
00:19:19,510 --> 00:19:26,989
governing equation of motion we get this equation
for UN+1.
174
00:19:26,989 --> 00:19:35,950
Now this equation also has this term GN+1,
UN+1 so at time TN+1 to obtain UN+1 I have
175
00:19:35,950 --> 00:19:42,270
to solve this non-linear equation, this is
a non-linear algebraic equation. Earlier in
176
00:19:42,270 --> 00:19:46,289
the absence of this term we would have inverted
this matrix and got the solution, but now
177
00:19:46,289 --> 00:19:52,419
these are a set of nonlinear algebraic equations
which at every time step we need to solve,
178
00:19:52,419 --> 00:19:59,019
okay this is the additional complexity, so
how do we tackle this? One of the popular
179
00:19:59,019 --> 00:20:03,929
algorithms for this is the
Newton-Raphson method so we will quickly recall
180
00:20:03,929 --> 00:20:09,150
what is the, how do we implement that method?
Suppose we consider a scalar nonlinear equation
181
00:20:09,150 --> 00:20:16,549
F(x) = 0, now our objective is to find roots
of this equation, we follow an iterative procedure,
182
00:20:16,549 --> 00:20:22,419
so what we do is XK, we consider XK to be
approximation to the root at K-th iteration
183
00:20:22,419 --> 00:20:30,549
step and XK+1 that is the improved solution
is XK plus a correction. Now the objective
184
00:20:30,549 --> 00:20:37,139
is to find what this correction should be,
so I will consider therefore F(xk+hk) which
185
00:20:37,139 --> 00:20:43,869
is since we are truncating this, this will
be F(xk) + HK DF/DX at X = XK. Now if XK +
186
00:20:43,869 --> 00:20:50,100
HK is a true solution the left hand side would
be 0 from which I get HK to be this, so that
187
00:20:50,100 --> 00:20:56,389
XK+1 is given by this, so the one quick observation
that we need to make is to obtain the rules
188
00:20:56,389 --> 00:21:09,740
now I need gradients of function F.
Now we are talking about Newmark's beta method,
189
00:21:09,740 --> 00:21:15,730
now if we are going to use Newton-Raphson
scheme within that, the scheme would be implicit
190
00:21:15,730 --> 00:21:22,389
now it also needs derivatives of the nonlinear
functions, and if I take X naught to be initial
191
00:21:22,389 --> 00:21:30,740
guess on the root, then we have to make this
guess to start the solution. Now in integrating
192
00:21:30,740 --> 00:21:34,960
equations of motion the initial guess could
be the solution from the linear system response
193
00:21:34,960 --> 00:21:40,690
or the nonlinear response at the preceding
time instant or some other reasonable approximation.
194
00:21:40,690 --> 00:21:44,970
Step size for implicit schemes it typically
is larger than what is needed for explicit
195
00:21:44,970 --> 00:21:49,690
schemes, the trade-off is large step size
and need to solve nonlinear algebraic equations
196
00:21:49,690 --> 00:21:53,929
at every time step, when we are talking about
implicit schemes.
197
00:21:53,929 --> 00:22:00,989
Now we will now consider to complete the discussion
on Newton-Raphson method a set of nonlinear
198
00:22:00,989 --> 00:22:07,250
equations in SN variables, suppose I have
FS (X1, X2, XN) = 0 for S = 1, 2, to N, there
199
00:22:07,250 --> 00:22:13,029
are N equations in N unknowns. Now how do
we solve these equations? So what I do, I
200
00:22:13,029 --> 00:22:19,720
take XKS from S = 1, 2 to N to the approximation
to the root at the K-th iteration step, and
201
00:22:19,720 --> 00:22:27,260
K+1, S is what I need to find. This is XKS
+ a correction, so I need to find this correction,
202
00:22:27,260 --> 00:22:37,970
again we use Taylor's expansion so this is
the function FS at the improved at XK+1S and
203
00:22:37,970 --> 00:22:44,429
this is given by this. So now by using a first-order
Taylor's expansion I get this equation and
204
00:22:44,429 --> 00:22:50,210
if this is indeed the true solution the left
hand side would be 0, so I get a set of algebraic
205
00:22:50,210 --> 00:22:59,510
equations for the unknown increments in X
case, so if you solve that I get HK as inverse
206
00:22:59,510 --> 00:23:06,330
of this matrix of gradients, this has to be
evaluated to implement the method.
207
00:23:06,330 --> 00:23:12,960
Now we have discussed the implicit scheme,
now let us see if what happens if I now tackle
208
00:23:12,960 --> 00:23:16,889
the non-linear equation using explicit scheme,
so let us consider the central difference
209
00:23:16,889 --> 00:23:21,320
method. So I have this equilibrium equation
at time T and these are the initial condition,
210
00:23:21,320 --> 00:23:26,539
and this is a time duration for which I want
to solve the equation of motion, so as before
211
00:23:26,539 --> 00:23:31,970
we approximate U dot(t) by a central difference
approximation, and U double dot(t) we get
212
00:23:31,970 --> 00:23:37,419
this, this we have done before, now I substitute
for these 2 in the governing equation and
213
00:23:37,419 --> 00:23:45,539
I get this equation and here the nonlinear
term I am writing this equation at T and the
214
00:23:45,539 --> 00:23:52,409
nonlinear term is here, the unknown is U(t)
+ delta T and that is here so on the right
215
00:23:52,409 --> 00:24:01,359
hand side I would have had G(U(t)) which would
be known to me, so here as you can see we
216
00:24:01,359 --> 00:24:07,669
need not solve any nonlinear equations to
advance from one step size to the next step
217
00:24:07,669 --> 00:24:14,100
size, so this is explicit.
So we need not have to solve nonlinear equations
218
00:24:14,100 --> 00:24:18,200
at every time step, we need to invert this
matrix all this is come, these discussions,
219
00:24:18,200 --> 00:24:23,239
these issues are common to what we already
discussed this is an explicit scheme it requires
220
00:24:23,239 --> 00:24:29,359
a special starting scheme, and also for linear
systems we already seen that this method is
221
00:24:29,359 --> 00:24:34,050
conditionally stable and there is a critical
step size which depends on the natural frequency
222
00:24:34,050 --> 00:24:39,269
of the, highest natural frequency of the system
and we expect that similar requirement also
223
00:24:39,269 --> 00:24:46,629
remains valid even when there is a non-linearity
may be more stringent requirement, so in any
224
00:24:46,629 --> 00:24:53,690
sense, in any case we can expect that to get
acceptable solutions from this we need to
225
00:24:53,690 --> 00:24:58,580
use a very small step size.
Now the idea of spectral stability does not
226
00:24:58,580 --> 00:25:04,869
apply to nonlinear systems, because of the
presence of non-linear terms we cannot use
227
00:25:04,869 --> 00:25:09,669
the Eigenvalues of the, spectral radius of
the amplification matrix of the linear system
228
00:25:09,669 --> 00:25:14,960
to infer stability of schemes when applied
to nonlinear system, so that issue has to
229
00:25:14,960 --> 00:25:20,470
be addressed separately. Now however the concept
of energy balance provides a means to get
230
00:25:20,470 --> 00:25:24,950
an idea about algorithmic damping and stability,
so when you are integrating this we could
231
00:25:24,950 --> 00:25:31,359
as well check for energy balance and see whether,
how it is behaving, so that may give you some
232
00:25:31,359 --> 00:25:36,570
idea about your choice of step size, if energy
balance conditions are violated then you are
233
00:25:36,570 --> 00:25:46,950
using wrong step size.
Now we have talked about explicit scheme and
234
00:25:46,950 --> 00:25:52,820
implicit scheme. Implicit schemes permits
you to use larger step size, but at every
235
00:25:52,820 --> 00:25:58,980
time step you have to solve non-linear algebraic
equation, explicit schemes there is no need
236
00:25:58,980 --> 00:26:05,740
to solve nonlinear algebraic equations but
they require small step size. Now the question
237
00:26:05,740 --> 00:26:11,850
is can we combine explicit and implicit approaches
simultaneously in one numerical algorithm,
238
00:26:11,850 --> 00:26:18,360
okay? For example if we have a structural
system in which there is local non-linearity,
239
00:26:18,360 --> 00:26:22,720
so where there is non-linearity you could
use explicit schemes and where there is linearity
240
00:26:22,720 --> 00:26:29,899
you can parts of the structures are linear
you could use implicit scheme, so that large
241
00:26:29,899 --> 00:26:40,200
step sizes could be used for certain parts
of the system and explicit schemes the smaller
242
00:26:40,200 --> 00:26:43,789
step size can be used for certain other parts,
so these ideas have been discussed in papers
243
00:26:43,789 --> 00:26:48,749
by Hughes, Pister, and Taylor and there is
another paper I have given the reference,
244
00:26:48,749 --> 00:26:56,210
so what I will do now is I will quickly outline
what the main issue is in this discussion.
245
00:26:56,210 --> 00:27:00,989
Now we want to discuss now implicit, explicit
and implicit explicit methods, actually we
246
00:27:00,989 --> 00:27:05,570
need to discuss how to combine implicit and
explicit approaches in a single, in a simultaneous
247
00:27:05,570 --> 00:27:10,929
linear same algorithm, so we will review that
we'll start with an implicit formulation and
248
00:27:10,929 --> 00:27:15,299
then we will look at explicit formulation
in some detail and then see what are the issues
249
00:27:15,299 --> 00:27:21,320
when we combine the two. So let's consider
the equation of motion MU double dot + N(U,U
250
00:27:21,320 --> 00:27:28,440
dot) = F(t), this N is for nonlinear terms,
so with certain initial condition and these
251
00:27:28,440 --> 00:27:34,450
are time duration of interest, so this U,
N and F are S cross 1 matrices, M is S cross
252
00:27:34,450 --> 00:27:42,999
S mass matrix, it could be diagonal and N
is S cross 1 means actually what we have,
253
00:27:42,999 --> 00:27:47,229
we have to understand that there are S number
of functions, each one of which is function
254
00:27:47,229 --> 00:27:53,369
of U1, U2, US and U1 dot, U2 dot, and US dot,
so I runs from 1, 2, S.
255
00:27:53,369 --> 00:27:58,919
Now we define what is known as tangent stiffness
matrix where we differentiate the nonlinear
256
00:27:58,919 --> 00:28:04,099
term with respect to UJ, so i-th function
d function differentiated with respect to
257
00:28:04,099 --> 00:28:09,309
J will give me the IJ-th element of the tangent
stiffness matrix. Similarly we talk about
258
00:28:09,309 --> 00:28:17,080
a tangent damping matrix where I differentiate
NI with respect to UJ dot and I get CT IJ,
259
00:28:17,080 --> 00:28:23,159
so this is tangent damping matrix. Now we
can see that mass matrix is symmetric, and
260
00:28:23,159 --> 00:28:28,950
tangent stiffness and tangent damping matrices
are also symmetric, now we denote by UN double
261
00:28:28,950 --> 00:28:36,229
dot ,UN dot and UN, the approximations to
U double dot(tn), and U dot(tn) and U(tn)
262
00:28:36,229 --> 00:28:42,630
respectively, so these are exact which are
not known these are the approximations.
263
00:28:42,630 --> 00:28:50,610
Now let's start by discussing the solution
for an implicit scheme, let's use in Newmark's
264
00:28:50,610 --> 00:28:55,009
method for illustration, it could be other
method as well, it could be HHT alpha method
265
00:28:55,009 --> 00:28:59,379
or generalized alpha method, it could be anything
else, but for purpose of discussion we will
266
00:28:59,379 --> 00:29:04,919
stick to Newmark's method. So the equilibrium
equation at time T is given by this, and at
267
00:29:04,919 --> 00:29:16,659
time N+1 the approximate equation will be
this. Now the Newmark method we assume that
268
00:29:16,659 --> 00:29:22,499
displacement is given by this where alpha
is an algorithmic parameter, and similarly
269
00:29:22,499 --> 00:29:30,340
what I will do now is I will collect the term
U double dot N+1 and delta T square alpha
270
00:29:30,340 --> 00:29:38,349
this is one term, this term makes the algorithm
implicit so all other terms which are at N
271
00:29:38,349 --> 00:29:45,690
are lumped in this U tilde N+1. Similarly
for velocity U dot of N+1 is given by this,
272
00:29:45,690 --> 00:29:55,619
delta is algorithmic parameter, again here
I will collect all the terms involving velocities
273
00:29:55,619 --> 00:30:03,119
and accelerations at TN and call this as U
tilde dot, N+1 and the other term which makes
274
00:30:03,119 --> 00:30:13,429
the algorithm implicit I will collect it here,
so this tilde are U tilde N+1 is this, and
275
00:30:13,429 --> 00:30:19,899
U tilde dot N+1 is given by this, we call
these values as predictor values, as you enter
276
00:30:19,899 --> 00:30:25,210
any time you would be knowing this at any
given time step you will be knowing this from
277
00:30:25,210 --> 00:30:31,529
your calculations at the preceding steps.
Now this is what we should find, these are
278
00:30:31,529 --> 00:30:38,999
the corrector values.
Now how do we implement this method? We set
279
00:30:38,999 --> 00:30:45,279
T = 0, and N = 0, so T = 0 the initial displacement
and velocity are known, and using the equation
280
00:30:45,279 --> 00:30:53,229
of motion I derive the initial acceleration.
Now predictor value is, now N is 0 therefore
281
00:30:53,229 --> 00:30:58,369
UN tilde will be U naught + U naught dot delta
T + delta T square U naught double dot, this
282
00:30:58,369 --> 00:31:04,899
into this term, so all these terms are known.
Similarly U tilde dot 1 because N is 0 is
283
00:31:04,899 --> 00:31:11,059
U dot naught + delta T 1 - delta U double
dot naught that is known, so that's what I
284
00:31:11,059 --> 00:31:15,289
meant predictor values will be known when
you enter this step.
285
00:31:15,289 --> 00:31:22,090
Now the corrector steps would be given by
this, now for U double dot N+1 I can use the
286
00:31:22,090 --> 00:31:31,099
equation of motion and write this terms there,
okay, this term, this M inverse of this. Similarly
287
00:31:31,099 --> 00:31:38,549
U dot N+1 is the predictor term + this term,
now you see here the unknowns are UN+1 and
288
00:31:38,549 --> 00:31:44,130
U dot N+1, and on the right hand side the
unknowns are contained in this nonlinear terms,
289
00:31:44,130 --> 00:31:48,239
so actually these two equation therefore form
a set of nonlinear equations for the unknowns
290
00:31:48,239 --> 00:31:55,529
UN+1 and U dot N+1, we use Newton-Raphson's
method with predictors as starting solution,
291
00:31:55,529 --> 00:32:03,700
we use this tilde as quantities with tilde
relations as the initial guess, and use Newton
292
00:32:03,700 --> 00:32:10,139
Raphson, so how does that work out? So I have
the first equation is given by this, I have
293
00:32:10,139 --> 00:32:13,739
taken the terms on the
right hand side to left hand side and I equate
294
00:32:13,739 --> 00:32:20,700
it to 0, let us call it as some capital G(UN+1,
U dot N+1) is 0. The second equation for velocity
295
00:32:20,700 --> 00:32:29,070
leads to this equation. Now these two are
pairs of nonlinear equations in UN+1 and U
296
00:32:29,070 --> 00:32:35,570
dot N+1, so there are 2S number of unknowns,
and there are 2S number of equations.
297
00:32:35,570 --> 00:32:40,619
Now Newton-Raphson solution now we have to
set up an iteration, and that iteration count
298
00:32:40,619 --> 00:32:48,059
is I, so I will first give a few mathematical
details and then we'll return to the algorithmic
299
00:32:48,059 --> 00:32:53,220
implementation, so at I-th step if this is
the guess and delta I is the correction, then
300
00:32:53,220 --> 00:32:57,299
this is the improved solution. So similarly
on velocity this is the improved solution,
301
00:32:57,299 --> 00:33:02,340
now you substitute this into these two equations
whatever we discussed till now use first-order
302
00:33:02,340 --> 00:33:08,049
Taylor's expansion I get this expression,
so what are unknowns here? This delta IJ and
303
00:33:08,049 --> 00:33:13,179
delta dot IJ are the unknowns which we need
to find. Similarly using Taylor's expansion
304
00:33:13,179 --> 00:33:20,389
on H, I get this equation, so I have now,
there are 2S number of unknowns and 2S number
305
00:33:20,389 --> 00:33:26,830
of linear algebraic equations they can be
cast in this form so the work involves evaluation
306
00:33:26,830 --> 00:33:32,899
of this matrix and inversion, so you have
this.
307
00:33:32,899 --> 00:33:43,479
Now G function is consequently this, and H
function, where G and H are this, I am just
308
00:33:43,479 --> 00:33:49,679
repeating for sake of completion here. Now
I need these gradients dou G/dou U, dou G/dou
309
00:33:49,679 --> 00:33:58,619
dot, dou H/dou U and dou H/U dot so they can
be evaluated from the Newton-Raphson, sorry
310
00:33:58,619 --> 00:34:03,960
Newmark beta, Newmark's as a model, and we
get these quantities in terms of tangent stiffness
311
00:34:03,960 --> 00:34:12,210
matrix and tangent damping matrix as shown
here. So in terms of, we can write now this
312
00:34:12,210 --> 00:34:15,750
matrix in terms of KT and CT which we will
do now.
313
00:34:15,750 --> 00:34:23,150
So let us now return to you know we started
with the 3 steps, we completed the predictor
314
00:34:23,150 --> 00:34:28,849
step and we are now formulating the corrector
steps, so corrector steps now involves Newton-Raphson
315
00:34:28,849 --> 00:34:36,220
steps, so a 4.1 we initialize, I = 0 that's
a iteration and these are the initial guesses
316
00:34:36,220 --> 00:34:41,740
which are the predictors which we have already
derived. Now iterations begin so I use this
317
00:34:41,740 --> 00:34:49,059
relation which is Newton-Raphson relation,
and increment I and then I will see whether
318
00:34:49,059 --> 00:34:56,950
this algorithm has converged by using epsilon
1 and epsilon 2 as measures of you know convergence,
319
00:34:56,950 --> 00:35:03,210
if this is satisfied we exit out of this iteration
loop, otherwise we go to 4.2 with, continue
320
00:35:03,210 --> 00:35:13,069
with these iterations. Now if the convergence
requirements are met we will be here and I
321
00:35:13,069 --> 00:35:20,520
will assign now to the displacement at N+1
and velocity at N+1 these converge quantities,
322
00:35:20,520 --> 00:35:27,550
then I have to compute U double dot N+1 for
which again I use the equation of motion.
323
00:35:27,550 --> 00:35:35,220
Now I increment N the time marching and if
I cross the final step size I exit, otherwise
324
00:35:35,220 --> 00:35:44,490
I go to step 2 which is that? We again start
with predictor for the next step, okay? So
325
00:35:44,490 --> 00:35:54,070
this is implementation for, implementation
of Newmark's method for a nonlinear system,
326
00:35:54,070 --> 00:35:59,549
this is an implicit scheme because as we have
seen average time step we are solving a set
327
00:35:59,549 --> 00:36:05,240
of non-linear equations.
Now we look at the same problem formulation
328
00:36:05,240 --> 00:36:14,910
using explicit scheme, now what we do is we
introduce this additional equation, see we
329
00:36:14,910 --> 00:36:20,859
got nonlinear equations in the preceding step
because for U double dot N+1 I replaced it
330
00:36:20,859 --> 00:36:29,560
by, if you see here carefully for U double
dot I have N+1, I am replacing this by F naught
331
00:36:29,560 --> 00:36:38,410
- N(UN+1, U dot N+1) and we are solving these
nonlinear equations, if instead of that for
332
00:36:38,410 --> 00:36:48,359
these terms inside this bracket if I use the
predictors, then this becomes explicit, okay
333
00:36:48,359 --> 00:36:55,140
so that leads us to the next scheme, the idea
is we will use this as, we introduce this
334
00:36:55,140 --> 00:37:02,380
equation MU double dot N+1 as N U tilde N+1,
U dot this also tilde, now use this equation
335
00:37:02,380 --> 00:37:06,750
implementing the corrector step, that means
at the stage of implementing the corrector
336
00:37:06,750 --> 00:37:12,680
step we will not use, we need not have to
solve a set of nonlinear algebraic equation,
337
00:37:12,680 --> 00:37:17,260
but we still need to iterate, there will be
multiple passes that will be needed and iteration
338
00:37:17,260 --> 00:37:21,839
will be needed, this avoids the solution to
the set of simultaneous nonlinear equations
339
00:37:21,839 --> 00:37:25,480
at every time step and the method becomes
explicit.
340
00:37:25,480 --> 00:37:31,079
So how does it work? You start with T = 0,
you have initial displacement and velocity
341
00:37:31,079 --> 00:37:39,990
use equation of motion get initial acceleration,
then you use predictor values you set I = 0,
342
00:37:39,990 --> 00:37:46,720
and your UN(i) is U naught, U1 dot(i) is U
naught dot and U double dot(i) is U double
343
00:37:46,720 --> 00:37:56,619
dot naught when you are at N = 0. Next I need
U tilde of I N+1, this is the predictor, this
344
00:37:56,619 --> 00:38:01,799
we have, this is for displacement, this is
for velocity. Now after having found out I
345
00:38:01,799 --> 00:38:08,600
want now an improved estimate of U double
dot N+1, to do that I will use the equation
346
00:38:08,600 --> 00:38:14,890
of motion where U and U dot are replaced by
the predictor values, so I get consequently
347
00:38:14,890 --> 00:38:19,920
U double dot I+1, N+1 is given by this, this
is an explicit equation there is no nonlinear
348
00:38:19,920 --> 00:38:27,650
equation solving here, so I get acceleration
from which I will again improve upon the displacement
349
00:38:27,650 --> 00:38:33,079
and velocity and check whether the convergence
has occurred, if it has occurred we will exit
350
00:38:33,079 --> 00:38:37,440
the iteration loop otherwise we iterate once
again with an improved, the predictor values
351
00:38:37,440 --> 00:38:46,289
are now replaced by these improved values.
So we if we exit after satisfactory convergence
352
00:38:46,289 --> 00:38:53,181
I get UN+1 which is the, displacement will
be what has converged here UI N+1, similarly
353
00:38:53,181 --> 00:38:59,039
velocity and acceleration are obtained, and
here I will use the, not the predictors but
354
00:38:59,039 --> 00:39:03,930
the converge solutions to find acceleration.
Then I will increment time and I will exit
355
00:39:03,930 --> 00:39:09,299
if we have cross the final time instant, otherwise
we go back and continue with the algorithm,
356
00:39:09,299 --> 00:39:16,720
so this is an explicit scheme, it involves
iterations, multiple passes but there is no
357
00:39:16,720 --> 00:39:25,050
nonlinear equation solving okay. Now the idea
is the question that we are trying to ask
358
00:39:25,050 --> 00:39:35,579
is now can we combine these 2 approaches?
So the question is can we use explicit and
359
00:39:35,579 --> 00:39:40,890
implicit concepts simultaneously in one algorithm,
the basic idea is we divide the finite elements
360
00:39:40,890 --> 00:39:46,099
into implicit and explicit sets, that means
the mesh, finite element mesh part of the
361
00:39:46,099 --> 00:39:53,620
mesh is designated as implicit, part of the
mesh is designated as explicit, so where you
362
00:39:53,620 --> 00:40:00,589
expect for instance local nonlinearities you
could or stress concentrations or things like
363
00:40:00,589 --> 00:40:08,400
that you use finer mesh or you use explicit
mesh, where you expect your understand, I
364
00:40:08,400 --> 00:40:15,990
mean linear behavior or things like that you
could use implicit mesh.
365
00:40:15,990 --> 00:40:21,490
Now let us yeah the superscripts I and E denote
the two sets respectively, that is implicit
366
00:40:21,490 --> 00:40:28,790
and explicit, so associated with this partitioning
I will have a MI, NI, FI as assembled mass,
367
00:40:28,790 --> 00:40:33,039
internal force, and an external force for
the implicit elements, and similarly ME, NE,
368
00:40:33,039 --> 00:40:38,099
FE as assembled mass, internal force and external
force for the explicit elements. So the equation
369
00:40:38,099 --> 00:40:46,630
of motion now can be written as (MI+ME) U
double dot + NI (U,U dot) + NE (U,U dot) is
370
00:40:46,630 --> 00:40:51,150
FI(t) FE(t), there is no, this is simply rewriting
the equation of motion without making any
371
00:40:51,150 --> 00:40:59,029
approximation. The approximation is now is
displayed here, for implicit part of the mesh
372
00:40:59,029 --> 00:41:04,160
I will use the current states UN+ 1 and U
dot N+1, for the explicit parts I will use
373
00:41:04,160 --> 00:41:13,690
the predictors, okay. So the idea is when
you solve for the nonlinear equations you
374
00:41:13,690 --> 00:41:25,390
will be using NI naught N okay, and this will
be, the structure of the matrices that you
375
00:41:25,390 --> 00:41:30,589
encounter and there banded ness features etc
is controlled by NI not by N, and NE makes
376
00:41:30,589 --> 00:41:36,049
no contribution to that in your solution to
nonlinear equations, so the effort in solving
377
00:41:36,049 --> 00:41:40,549
the nonlinear equations can come down substantially,
so these are the
378
00:41:40,549 --> 00:41:48,029
predictor values as before, and I get UN+1
is the predictor plus this, which now I will
379
00:41:48,029 --> 00:41:54,740
write for U double dot N+1, in implicit scheme
I wrote simply the governing equation here,
380
00:41:54,740 --> 00:42:00,380
in the explicit scheme I used the predictor
equation, so here what I am doing I am splitting
381
00:42:00,380 --> 00:42:06,780
for implicit elements I am writing NI with
current UN+1 and U dot N+1 retained as it
382
00:42:06,780 --> 00:42:13,890
is, for explicit elements I am using the predictors.
So the unknown nonlinear equation, effort
383
00:42:13,890 --> 00:42:20,069
to solve nonlinear equation is now related
to this NI, so similarly I get an equation
384
00:42:20,069 --> 00:42:26,480
for U dot N+1 which is again written here.
The tangent stiffness and tangent damping
385
00:42:26,480 --> 00:42:33,820
matrices are now evaluated only with respect
to NI, you don't need the tangent KT and CT
386
00:42:33,820 --> 00:42:41,760
for complete N vector, so you need only for
NI, and only for NI, so consequently the band
387
00:42:41,760 --> 00:42:47,319
profile of these matrices would correspond
only to the connectivity of the implicit elements,
388
00:42:47,319 --> 00:42:51,930
so they can be quite sparse and well-structured
and maybe easier to solve.
389
00:42:51,930 --> 00:42:57,650
Now I will not get into the details of this
implementation but we can make some observations,
390
00:42:57,650 --> 00:43:02,380
now what we have done is we have discussed
the implicit, explicit, and implicit-explicit
391
00:43:02,380 --> 00:43:06,809
methods in the context of Newmark's method.
Now same discussion or similar discussion
392
00:43:06,809 --> 00:43:12,140
can also be developed for other schemes such
as HHT alpha and generalized alpha methods,
393
00:43:12,140 --> 00:43:19,050
so the predictor equations that you derive
for HHT alpha should originate from HHT alpha
394
00:43:19,050 --> 00:43:26,029
representation, here the predictor equation
you know was from the Newmark beta method,
395
00:43:26,029 --> 00:43:31,290
Newmark's method, so that's one thing that
we have to bear in mind. Now the partitioning
396
00:43:31,290 --> 00:43:36,779
of the mesh into explicit and implicit elements
can also be interpreted as splitting the operator
397
00:43:36,779 --> 00:43:42,500
into explicit and implicit parts, so this
is so-called operator splitting implicit-explicit
398
00:43:42,500 --> 00:43:45,869
method.
Now what I wish to do in the remaining few
399
00:43:45,869 --> 00:43:51,039
minutes of this lecture is to give a brief
introduction to what are the issues that we
400
00:43:51,039 --> 00:43:57,690
would like to address in the next module.
The next module of this course is on model
401
00:43:57,690 --> 00:44:05,079
reduction and sub-structuring techniques.
Why do we need that? The need arises because
402
00:44:05,079 --> 00:44:13,069
of, for example treatment of large scale problems
I may like to reduce the size of the problem.
403
00:44:13,069 --> 00:44:19,119
Now the need also arises whenever you have
to deal with situations when results from
404
00:44:19,119 --> 00:44:24,640
experiments need to be discussed in conjunction
with predictions from mathematical models.
405
00:44:24,640 --> 00:44:30,650
For example if you have an instrumented structure
there, if you say that if you have used say
406
00:44:30,650 --> 00:44:37,339
S number of sensors, the experimental model
will have S degrees of freedom if we call
407
00:44:37,339 --> 00:44:42,520
the number of sensors as degrees of freedom,
it has S degrees of freedom, but the computational
408
00:44:42,520 --> 00:44:48,589
model for the same structure can have different
degrees of freedom, typically the size of
409
00:44:48,589 --> 00:44:54,470
a mathematical model could be fairly large
compared with number of sensors that we use
410
00:44:54,470 --> 00:44:59,500
in experimental part, so there is a mismatch
of degrees of freedom in computational model
411
00:44:59,500 --> 00:45:08,180
and experimental model, so how do we deal
with that? Or the other situation is different
412
00:45:08,180 --> 00:45:14,500
parts of a structure could be developed by
different teams, possibly working independently,
413
00:45:14,500 --> 00:45:18,099
working independently and possibly by using
both experimental and computational tools,
414
00:45:18,099 --> 00:45:22,650
for example this type of issues can arise
in dealing with automotive systems or space
415
00:45:22,650 --> 00:45:27,460
vehicles and things like that where different
parts of the structure can be designed and
416
00:45:27,460 --> 00:45:31,869
developed by different teams, so each team
will make a valid finite element model for
417
00:45:31,869 --> 00:45:37,180
the particular sub-structure that they are
dealing with and they may also conduct experiments
418
00:45:37,180 --> 00:45:39,990
on each of the sub-structures that they are
studying.
419
00:45:39,990 --> 00:45:44,190
And finally when the structure is assembled
then computational model need to be assembled,
420
00:45:44,190 --> 00:45:48,470
the experimental model needs to be assembled,
so what are the issues? So these issues take
421
00:45:48,470 --> 00:45:53,810
us to a discussion on sub-structuring methods.
Finally another context in which these issues
422
00:45:53,810 --> 00:45:59,619
will become important is in the context of
what are known as hybrid simulations, here
423
00:45:59,619 --> 00:46:06,461
the word hybrid means in testing of structures
under dynamic loads we combine both experimental
424
00:46:06,461 --> 00:46:11,180
and computational methods, that is why they
are called hybrid simulation. Hybrid simulations
425
00:46:11,180 --> 00:46:18,710
are essentially experimental testing techniques
to qualify structures for specified dynamic
426
00:46:18,710 --> 00:46:25,940
loads, for example earthquake qualification
testing of a structure, here we combine both
427
00:46:25,940 --> 00:46:29,690
experimental and computational models for
the same structure and the same team works
428
00:46:29,690 --> 00:46:34,680
on both of them, so a part of the structure
is studied experimentally, and a part of the
429
00:46:34,680 --> 00:46:41,390
structure studied computationally. So how
do we combine in the mathematical model that
430
00:46:41,390 --> 00:46:46,300
we use in our computational model with the
part of the structure which is actually studied
431
00:46:46,300 --> 00:46:52,010
experimentally and for which no mathematical
model is done, so this takes us to different
432
00:46:52,010 --> 00:46:57,750
types of sub-structuring problems.
Now I'll briefly introduce the problem of
433
00:46:57,750 --> 00:47:04,299
model reduction, then we'll take up discussions
in subsequent lectures. Now suppose we consider
434
00:47:04,299 --> 00:47:08,720
a N degree of freedom finite element model
for a linear system, say the linear system
435
00:47:08,720 --> 00:47:15,890
is governed by MX double dot + CX + KX is
F(t) with some initial conditions. What is
436
00:47:15,890 --> 00:47:21,880
the problem of model reduction here? The objective
of model reduction is to replace this N degree
437
00:47:21,880 --> 00:47:27,589
of freedom system by an equivalent lower order
system, a lower case N degree of freedom system,
438
00:47:27,589 --> 00:47:32,099
where the reduced system will have much lesser
number of degrees of freedom than the larger
439
00:47:32,099 --> 00:47:37,260
system, so the objective is to derive the
equation of motion for another system which
440
00:47:37,260 --> 00:47:42,760
has less number of degrees of freedom as shown
here, and these two should be related this
441
00:47:42,760 --> 00:47:51,540
MR, CR, KR are the reduced N cross N structural
matrix. XM is N cross 1 vector of degrees
442
00:47:51,540 --> 00:47:59,730
of freedom which have been retained in the
reduced model, in all the model reduction
443
00:47:59,730 --> 00:48:05,950
techniques the displacement vector X(t) is
related to, I mean is partitioned as XM(t)
444
00:48:05,950 --> 00:48:11,369
and XS(t), XM(t) are the degrees of freedom
that we want to retain in the reduced model
445
00:48:11,369 --> 00:48:17,930
XS(t) are the degrees of freedom which we
would like to eliminate, so this M and S refer
446
00:48:17,930 --> 00:48:25,210
to the subscripts M refers to master SS slave,
so those degrees of freedom which are retained
447
00:48:25,210 --> 00:48:28,960
are called master degrees of freedom, those
degrees of freedom which are eliminated are
448
00:48:28,960 --> 00:48:39,500
called slave degrees of freedom. So XM(t)
is N cross 1, master degrees of freedom and
449
00:48:39,500 --> 00:48:54,660
XS(t) is N - N slave degrees of freedom.
Now we represent X(t) as some sai into XM(t),
450
00:48:54,660 --> 00:49:00,660
okay, this is how we eliminate the slave degrees
of freedom, where sai is an N cross M transformation
451
00:49:00,660 --> 00:49:12,630
matrix, so we substitute this into the governing
equation I have M sai XM double dot + C sai
452
00:49:12,630 --> 00:49:23,740
XM dot + and
for X I am writing sai into XM, so this is
453
00:49:23,740 --> 00:49:30,020
M sai XM double dot + C sai XM dot + K sai
XM is F(t), I'll pre multiply by sai transpose,
454
00:49:30,020 --> 00:49:36,289
so I get this equation. So now I will call
sai transpose M sai as MR which is a reduced
455
00:49:36,289 --> 00:49:43,079
mass matrix, and CR which is sai transpose
C sai is the reduced damping matrix, KR is
456
00:49:43,079 --> 00:49:47,450
sai transpose K sai which is a reduced stiffness
matrix, you can quickly verify that these
457
00:49:47,450 --> 00:49:54,500
matrices are all symmetric, and FR(t) is sai
transpose F(t) which is a reduced force vector.
458
00:49:54,500 --> 00:50:02,000
Now the question is how do we select sai?
Okay, that is the main question. When we want
459
00:50:02,000 --> 00:50:08,109
to, when we reduce the size of the model for
example the original model would have N pairs
460
00:50:08,109 --> 00:50:15,299
of natural frequencies and Eigen vectors,
the reduced model will have lowercase N Eigen
461
00:50:15,299 --> 00:50:20,673
pairs, should these be equivalent? For example
100 degrees of freedom is represented in terms
462
00:50:20,673 --> 00:50:25,019
of 10 degrees of freedom model, so the 100
degrees of freedom will have 100 natural frequencies
463
00:50:25,019 --> 00:50:29,690
possible, and a 10 degree of freedom will
have 10 possible natural frequency, should
464
00:50:29,690 --> 00:50:34,410
these 10 natural frequencies be related to
the natural frequencies of the larger model,
465
00:50:34,410 --> 00:50:39,070
okay, that's the question. Similarly mode
shapes, and we can also ask this question
466
00:50:39,070 --> 00:50:43,730
on the frequency response function between
2 coordinates, for example should the frequency
467
00:50:43,730 --> 00:50:49,259
response function over a given frequency range
of the reduced system serve as acceptable
468
00:50:49,259 --> 00:50:52,770
approximations to the corresponding FRF's
of the original system, what is the objective
469
00:50:52,770 --> 00:50:57,799
of model reduction that we have to specify.
Similarly if you are interested in transient
470
00:50:57,799 --> 00:51:02,609
behavior of the system, should transient response
to dynamic excitations for the reduced system
471
00:51:02,609 --> 00:51:09,089
serve as acceptable approximation to the response
of the original system, so the question on
472
00:51:09,089 --> 00:51:15,339
how to select sai, therefore gets answered
in a broad way that it depends on the situation
473
00:51:15,339 --> 00:51:20,230
that context in which you want to reduce the
model, so we have to discuss different options
474
00:51:20,230 --> 00:51:28,170
that are possible.
Now actually we had talked about model reduction
475
00:51:28,170 --> 00:51:33,529
but there can be a problem of model expansion
also, how does it come about? Suppose if you
476
00:51:33,529 --> 00:51:38,340
consider a structural system that is being
studied both experimentally and computationally,
477
00:51:38,340 --> 00:51:42,420
if lowercase n is the number of measured degrees
of freedom, and capital N is the degrees of
478
00:51:42,420 --> 00:51:46,431
freedom in computational model, typically
the degrees of freedom in computational model
479
00:51:46,431 --> 00:51:50,990
is much greater than the degrees of freedom
in an experiment. Now when we are reconciling
480
00:51:50,990 --> 00:51:56,059
the prediction from the computational model
with what we have measured, then there are
481
00:51:56,059 --> 00:52:01,690
2 options available to us, one is reduce the
size of computational model so that only the
482
00:52:01,690 --> 00:52:06,990
degrees of freedom which are common to both
experimental and computational models are
483
00:52:06,990 --> 00:52:11,160
retained, or alternately expand the size of
the experimental model so that the degrees
484
00:52:11,160 --> 00:52:16,430
of freedom not present in the experimental
model are approximately represented, okay
485
00:52:16,430 --> 00:52:20,200
so that there is a matching of degrees of
freedom between experimental and computational
486
00:52:20,200 --> 00:52:25,279
models.
Now some of these issues will be discussing
487
00:52:25,279 --> 00:52:30,160
what we will do is, I'd propose to discuss
three methods, one is what is known as static
488
00:52:30,160 --> 00:52:34,579
condensation, that is also known as Guyan's
reduction, and there is a method known as
489
00:52:34,579 --> 00:52:40,039
dynamic condensation, and third method is
what is known as system equivalent reduction
490
00:52:40,039 --> 00:52:51,309
expansion process and it is abbreviated as
SEREP, we will be discussing this in the lectures
491
00:52:51,309 --> 00:52:55,230
to follow, we will conclude this lecture at
this stage.