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.