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:17,400 hundreds of MIT courses, visit MIT OpenCourseWare at 7 00:00:17,400 --> 00:00:18,650 ocw.mit.edu. 8 00:00:21,980 --> 00:00:23,600 PROFESSOR: Ladies and gentlemen, welcome to this 9 00:00:23,600 --> 00:00:26,160 lecture on non-linear finite element analysis of solids and 10 00:00:26,160 --> 00:00:27,320 structures. 11 00:00:27,320 --> 00:00:30,300 In this lecture, I'd like to continue with our discussion 12 00:00:30,300 --> 00:00:34,530 of the methods of solution for the finite element equations 13 00:00:34,530 --> 00:00:38,670 that we are encountering in non-linear dynamic analysis. 14 00:00:38,670 --> 00:00:41,640 In the previous lecture, I pointed out that the 15 00:00:41,640 --> 00:00:45,100 techniques used for the solution of these finite 16 00:00:45,100 --> 00:00:47,590 element equations in non-linear dynamic analysis 17 00:00:47,590 --> 00:00:53,020 can be classified as data integration techniques, 18 00:00:53,020 --> 00:00:57,410 methods of multiple position, and substructuring. 19 00:00:57,410 --> 00:01:01,870 We talked in the previous lecture about these types of 20 00:01:01,870 --> 00:01:04,870 methods, and in particular, we discussed in detail an 21 00:01:04,870 --> 00:01:08,110 explicit integration method, namely the central difference 22 00:01:08,110 --> 00:01:11,570 method, and an implicit integration method, namely the 23 00:01:11,570 --> 00:01:13,210 trapezoidal rule. 24 00:01:13,210 --> 00:01:16,720 We talked about how these methods are applied, how time 25 00:01:16,720 --> 00:01:20,210 steps might be selected, what the restrictions of these 26 00:01:20,210 --> 00:01:21,830 methods are, and so on. 27 00:01:21,830 --> 00:01:25,080 I'd like to now in this lecture discuss with you the 28 00:01:25,080 --> 00:01:28,730 method of multiple position, and also 29 00:01:28,730 --> 00:01:31,380 substructuring technique. 30 00:01:31,380 --> 00:01:36,050 After the discussion off these two techniques, I'd like to 31 00:01:36,050 --> 00:01:40,270 talk with you about some examples. 32 00:01:40,270 --> 00:01:43,000 As a first example, we want to consider the wave 33 00:01:43,000 --> 00:01:45,320 propagation in a rod. 34 00:01:45,320 --> 00:01:48,650 As a second example, I'd like to share with you experiences 35 00:01:48,650 --> 00:01:54,000 regarding the response of a 3 degree of freedom system. 36 00:01:54,000 --> 00:01:56,830 As a third example, then, we want to consider the analysis 37 00:01:56,830 --> 00:02:00,360 of a 10-story tapered tower. 38 00:02:00,360 --> 00:02:05,300 Example 3 and example 1 are really solutions of a linear 39 00:02:05,300 --> 00:02:08,320 system, but we learn quite a bit, even regarding non-linear 40 00:02:08,320 --> 00:02:11,890 analysis, from these examples. 41 00:02:11,890 --> 00:02:14,860 As the fourth example, then, I'd like to share with you 42 00:02:14,860 --> 00:02:18,960 experiences obtained in the analysis of a pendulum that 43 00:02:18,960 --> 00:02:21,396 goes through very large displacement. 44 00:02:21,396 --> 00:02:25,660 As a fifth example, then, we consider the pipe whip 45 00:02:25,660 --> 00:02:27,500 response solution. 46 00:02:27,500 --> 00:02:30,860 This is of much practical interest. 47 00:02:30,860 --> 00:02:34,860 And after these examples, which I've prepared view 48 00:02:34,860 --> 00:02:39,220 graphs for, I'd like to consider with you some slides. 49 00:02:39,220 --> 00:02:42,840 On these slides, we will have the analysis of a control rod 50 00:02:42,840 --> 00:02:49,860 drive housing with the response of a spherical cap, a 51 00:02:49,860 --> 00:02:52,410 highly non-linear problem, large displacement, 52 00:02:52,410 --> 00:02:56,940 elastoplastic response, and then the analysis of a fluid 53 00:02:56,940 --> 00:03:01,890 structure interaction problem, a pipe test in which we are 54 00:03:01,890 --> 00:03:04,030 comparing our computer's results with 55 00:03:04,030 --> 00:03:06,610 experimental results. 56 00:03:06,610 --> 00:03:09,970 I should point out right now that the details regarding 57 00:03:09,970 --> 00:03:14,650 these example solutions are given in some papers that I 58 00:03:14,650 --> 00:03:16,220 had also earlier. 59 00:03:16,220 --> 00:03:20,530 And please look at the study guide, where the references to 60 00:03:20,530 --> 00:03:22,560 these papers are given. 61 00:03:22,560 --> 00:03:25,580 Let me now walk over here to the view graph that we have 62 00:03:25,580 --> 00:03:28,310 prepared regarding this lecture. 63 00:03:28,310 --> 00:03:30,940 And the first method of solution that I wanted to 64 00:03:30,940 --> 00:03:35,080 discuss with you is the method of multiple position. 65 00:03:35,080 --> 00:03:39,410 Quite a number of you probably have used this technique 66 00:03:39,410 --> 00:03:41,240 already in linear analysis. 67 00:03:41,240 --> 00:03:43,910 In fact, the word "superposition" somehow 68 00:03:43,910 --> 00:03:47,440 implies that you are dealing with a linear analysis. 69 00:03:47,440 --> 00:03:50,670 However, we can extend that method directly to non-linear 70 00:03:50,670 --> 00:03:55,660 analysis if we recognize that the mode shape vectors of the 71 00:03:55,660 --> 00:04:00,160 system at a particular time can be used as generalized 72 00:04:00,160 --> 00:04:04,670 displacements, or basis vectors, to express the 73 00:04:04,670 --> 00:04:08,100 response to any other time, as well. 74 00:04:08,100 --> 00:04:11,370 The method is effective when, in non-linear analysis, the 75 00:04:11,370 --> 00:04:16,410 response lies in only a few vibration modes, just the same 76 00:04:16,410 --> 00:04:19,670 way as we are applying multiple position also in 77 00:04:19,670 --> 00:04:24,840 linear analysis, namely, it's only effective really when a 78 00:04:24,840 --> 00:04:28,840 few vibration modes capture the total dynamic response. 79 00:04:28,840 --> 00:04:30,080 The same holds here, of course, 80 00:04:30,080 --> 00:04:32,220 in non-linear analysis. 81 00:04:32,220 --> 00:04:34,420 And when the system has possibly only local 82 00:04:34,420 --> 00:04:36,350 non-linearities. 83 00:04:36,350 --> 00:04:40,550 Let's look at is equations that we're using in this 84 00:04:40,550 --> 00:04:42,250 multiple position technique. 85 00:04:42,250 --> 00:04:45,810 The basic starting point are the equations that we have 86 00:04:45,810 --> 00:04:49,310 seen before regarding dynamic analysis. 87 00:04:49,310 --> 00:04:53,060 But notice in this particular equation, I have assumed now 88 00:04:53,060 --> 00:04:57,280 that the damping matrix is not present, so we have no damping 89 00:04:57,280 --> 00:05:00,240 in the system coming from the damping matrix. 90 00:05:00,240 --> 00:05:03,110 There could, of course, be damping due to the material 91 00:05:03,110 --> 00:05:05,200 non-linearities, which are captured in 92 00:05:05,200 --> 00:05:07,450 the K and the F matrix. 93 00:05:07,450 --> 00:05:09,680 Well, this is the set of equations 94 00:05:09,680 --> 00:05:11,200 that we want to solve. 95 00:05:11,200 --> 00:05:16,230 And usually, we just iterate on these equations until the 96 00:05:16,230 --> 00:05:19,940 right-hand side is 0, where of course we also would bring 97 00:05:19,940 --> 00:05:22,330 this one here onto the right-hand side. 98 00:05:22,330 --> 00:05:27,210 So think of this one here as a minus on the right-hand side, 99 00:05:27,210 --> 00:05:30,980 and if, then, this right-hand side, with M on that 100 00:05:30,980 --> 00:05:38,230 right-hand side is 0, meaning the data UK is 0, no more 101 00:05:38,230 --> 00:05:41,040 updating to the incremental displacements, then we have 102 00:05:41,040 --> 00:05:44,020 established equilibrium configuration corresponding to 103 00:05:44,020 --> 00:05:49,040 time t, of course, including inertia effects. 104 00:05:49,040 --> 00:05:54,000 In multiple position, we don't want to directly to deal with 105 00:05:54,000 --> 00:05:56,210 these equations in this form. 106 00:05:56,210 --> 00:06:00,750 Instead, we want to transform them into a new form, and into 107 00:06:00,750 --> 00:06:08,670 a form such that we obtain less unknowns to deal with. 108 00:06:08,670 --> 00:06:13,730 Let's assume for the moment that tau is equal to 0 so that 109 00:06:13,730 --> 00:06:17,610 we are dealing with the solution of the equations 110 00:06:17,610 --> 00:06:19,610 using the initial stress method. 111 00:06:19,610 --> 00:06:23,270 We talked about initial stress method in an earlier lecture. 112 00:06:23,270 --> 00:06:29,990 Then using this transformation here, where phi R are the 113 00:06:29,990 --> 00:06:37,150 eigenvectors of this problem down here, And the xi are the 114 00:06:37,150 --> 00:06:39,870 generalized displacements corresponding to these 115 00:06:39,870 --> 00:06:42,000 eigenvectors. 116 00:06:42,000 --> 00:06:45,350 Notice this equation is exactly the same that we're 117 00:06:45,350 --> 00:06:48,440 used to seeing in linear analysis. 118 00:06:48,440 --> 00:06:52,080 Of course, in linear analysis, we are usually not putting 119 00:06:52,080 --> 00:06:53,500 this t plus delta t up there. 120 00:06:53,500 --> 00:06:58,710 It's understood that the time dependency is in xi, but the 121 00:06:58,710 --> 00:07:02,040 space dependency is in phi i. 122 00:07:02,040 --> 00:07:05,784 And notice we are summing here from r to s. 123 00:07:05,784 --> 00:07:09,060 r would have to be selected by the analyst, 124 00:07:09,060 --> 00:07:12,580 and s the same way. 125 00:07:12,580 --> 00:07:17,010 Having introduced this transformation with this 126 00:07:17,010 --> 00:07:23,710 eigenvalue problem here, and substituting for the U 127 00:07:23,710 --> 00:07:27,880 displacements in the general equilibrium equation, and 128 00:07:27,880 --> 00:07:31,400 using the same approach as we are used to in linear 129 00:07:31,400 --> 00:07:35,990 analysis, we directly obtain this set of equations. 130 00:07:35,990 --> 00:07:42,980 Now, notice here that in this vector x double dot, and in 131 00:07:42,980 --> 00:07:47,340 this vector x, delta meaning increment, we have only the 132 00:07:47,340 --> 00:07:53,760 components corresponding to the r mode up to the s mode. 133 00:07:53,760 --> 00:07:56,530 These are the generalized displacements corresponding to 134 00:07:56,530 --> 00:07:59,540 phi r up to phi s. 135 00:07:59,540 --> 00:08:03,470 Here we have phi r to phi s listed in this 136 00:08:03,470 --> 00:08:05,800 matrix, capital phi. 137 00:08:05,800 --> 00:08:09,070 Lower-case phi for these here. 138 00:08:09,070 --> 00:08:13,980 These are the columns in this matrix that is denoted as a 139 00:08:13,980 --> 00:08:15,650 capital phi. 140 00:08:15,650 --> 00:08:21,890 Notice here we define an omega squared matrix, 0 elements in 141 00:08:21,890 --> 00:08:23,370 the off-diagonals. 142 00:08:23,370 --> 00:08:28,160 And on the diagonal, we have the frequencies of the system, 143 00:08:28,160 --> 00:08:32,880 of the actually linear system, because we use 0K in the 144 00:08:32,880 --> 00:08:35,590 eigenvalue solution. 145 00:08:35,590 --> 00:08:40,750 These are the entries for this equation up here. 146 00:08:40,750 --> 00:08:45,400 And notice that we have now here much less equations than 147 00:08:45,400 --> 00:08:49,650 the total number of finite element degrees of freedom. 148 00:08:49,650 --> 00:08:55,330 Typically, r might be 1, s might be 20, and then we would 149 00:08:55,330 --> 00:08:59,330 have 20 equations here, 20 equations only. 150 00:08:59,330 --> 00:09:06,280 Of course, r could also be 25, and s 30, in which case we 151 00:09:06,280 --> 00:09:10,470 have here a multiple position corresponding to the 25th up 152 00:09:10,470 --> 00:09:13,240 to the 30th mode. 153 00:09:13,240 --> 00:09:18,050 In any case, in general, we will have much less equations 154 00:09:18,050 --> 00:09:22,830 then finite element degrees of freedom that are listed, for 155 00:09:22,830 --> 00:09:26,970 example, in this vector here, R, I should say. 156 00:09:26,970 --> 00:09:30,920 The externally applied load corresponding to these degrees 157 00:09:30,920 --> 00:09:33,280 of freedom are listed in here. 158 00:09:33,280 --> 00:09:37,940 Of course, the F vector is the vector of nodal point forces 159 00:09:37,940 --> 00:09:41,480 corresponding to the element stresses at time t plus delta 160 00:09:41,480 --> 00:09:47,180 t, and at the end of iteration, k minus 1. 161 00:09:47,180 --> 00:09:54,160 This matrix times that vector gives us a vector which has 162 00:09:54,160 --> 00:09:59,930 the same length as this vector has and that vector has. 163 00:09:59,930 --> 00:10:03,930 So we are not dealing anymore with all the degrees of 164 00:10:03,930 --> 00:10:10,010 freedom that are listed for which we have the forces here 165 00:10:10,010 --> 00:10:13,440 in the solution of the equation, but with much less 166 00:10:13,440 --> 00:10:16,620 degrees of freedom, the degrees of freedom being equal 167 00:10:16,620 --> 00:10:20,660 to the number of mode shapes that are being employed. 168 00:10:20,660 --> 00:10:23,760 An interesting and most important point is that in 169 00:10:23,760 --> 00:10:27,660 non-linear analysis, we have a coupling, due to this vector 170 00:10:27,660 --> 00:10:31,530 appearing in here, between the different mode shape 171 00:10:31,530 --> 00:10:35,390 equations, or modal equations. 172 00:10:35,390 --> 00:10:40,450 We cannot solve, as we usually do in linear analysis, one 173 00:10:40,450 --> 00:10:44,510 equation of the modes of a position after the other. 174 00:10:44,510 --> 00:10:48,340 In linear analysis, you would typically solve the dynamic 175 00:10:48,340 --> 00:10:52,440 response equation in mode 1, then the one in mode 2, then 176 00:10:52,440 --> 00:10:57,700 the one in mode 3, and like that, you'd sweep over all the 177 00:10:57,700 --> 00:11:02,100 generalized degrees of freedom, one after the other, 178 00:11:02,100 --> 00:11:07,260 for the total response from time 0 to the end of the 179 00:11:07,260 --> 00:11:10,030 response time that you're interested in. 180 00:11:10,030 --> 00:11:13,260 In non-linear analysis, we have to go step by step still 181 00:11:13,260 --> 00:11:22,100 forward in time, solving all equations, all modal equations 182 00:11:22,100 --> 00:11:26,360 at time t, then all modal equations at time t plus delta 183 00:11:26,360 --> 00:11:31,670 t, then at time t plus 2 delta t, and so on, because there is 184 00:11:31,670 --> 00:11:36,620 this coupling, and because we don't know f unless we have 185 00:11:36,620 --> 00:11:39,015 all the information coming from all the mode shapes. 186 00:11:41,900 --> 00:11:46,670 Well, a typical problem that you might want to solve with 187 00:11:46,670 --> 00:11:49,440 this approach is the pipe whip problem. 188 00:11:49,440 --> 00:11:55,160 Here we have a very long, slender pipe. 189 00:11:55,160 --> 00:11:57,980 There is here a restraint that has a gap. 190 00:11:57,980 --> 00:12:01,580 And suddenly, this end of the pipe is subjected to a very 191 00:12:01,580 --> 00:12:02,780 large force. 192 00:12:02,780 --> 00:12:05,750 The pipe undergoes dynamic motion, a wave travels along 193 00:12:05,750 --> 00:12:09,520 here, goes plastic, hits this stop, et cetera. 194 00:12:09,520 --> 00:12:12,470 It's a very complicated problem, and this approach of 195 00:12:12,470 --> 00:12:16,760 multiple position is it cheap way of obtaining a solution to 196 00:12:16,760 --> 00:12:18,380 this problem. 197 00:12:18,380 --> 00:12:21,210 Of course, the accuracy of the solution depends on how many 198 00:12:21,210 --> 00:12:23,570 mode shapes you include. 199 00:12:23,570 --> 00:12:26,340 If you only include a few mode shapes, the solution is quite 200 00:12:26,340 --> 00:12:29,610 inexpensive, but the accuracy is not very high. 201 00:12:29,610 --> 00:12:34,500 If you include as many mode shapes as there are finite 202 00:12:34,500 --> 00:12:37,190 element degrees of freedom, then of course you get the 203 00:12:37,190 --> 00:12:38,130 exact response. 204 00:12:38,130 --> 00:12:41,250 All you have done is transformed from one basis, 205 00:12:41,250 --> 00:12:43,560 the finite element degrees of freedom basis, to the mode 206 00:12:43,560 --> 00:12:44,640 shape basis. 207 00:12:44,640 --> 00:12:47,030 But then, of course, the solution is very expensive, 208 00:12:47,030 --> 00:12:49,160 and you don't want to ever do is that. 209 00:12:49,160 --> 00:12:52,950 The point is that you want to perform such analysis if you 210 00:12:52,950 --> 00:12:55,750 only need to include a few modes in 211 00:12:55,750 --> 00:12:57,980 the response solution. 212 00:12:57,980 --> 00:13:02,120 We would actually look at this problem, or such a problem, a 213 00:13:02,120 --> 00:13:04,930 little later, when we discuss example analyses. 214 00:13:07,540 --> 00:13:11,660 As a second technique that can also be very useful in 215 00:13:11,660 --> 00:13:14,750 non-linear dynamic analysis, I'd like to discuss with you 216 00:13:14,750 --> 00:13:17,130 the substructuring technique. 217 00:13:17,130 --> 00:13:20,810 This procedure is used is implicit time integration, and 218 00:13:20,810 --> 00:13:26,140 the important aspect is that prior to the analysis, we are 219 00:13:26,140 --> 00:13:30,310 setting u this effective stiffness matrix corresponding 220 00:13:30,310 --> 00:13:32,510 to the linear degrees of freedom only, and we are 221 00:13:32,510 --> 00:13:36,640 statically condensing out all the degrees of freedom that 222 00:13:36,640 --> 00:13:41,130 will remain linear throughout the response solution. 223 00:13:41,130 --> 00:13:44,700 It is quite effective for local non-linearities. 224 00:13:44,700 --> 00:13:48,140 Here, for example, we have a tower structure, and the only 225 00:13:48,140 --> 00:13:51,210 non-linearity that this structure sees, say, is the 226 00:13:51,210 --> 00:13:52,460 slip down here. 227 00:13:55,270 --> 00:13:58,360 This typically would mean that we have a large number of 228 00:13:58,360 --> 00:14:01,620 linear degrees of freedom that we can deal with very 229 00:14:01,620 --> 00:14:04,640 effectively, using the substructuring approach. 230 00:14:04,640 --> 00:14:07,480 The non-linear degrees of freedom here, of course, 231 00:14:07,480 --> 00:14:11,430 require iteration and updating of matrices, and those we deal 232 00:14:11,430 --> 00:14:14,000 with in the usual way. 233 00:14:14,000 --> 00:14:19,100 Let me show you in more detail what I mean here. 234 00:14:19,100 --> 00:14:21,260 Let's look schematically at an example. 235 00:14:21,260 --> 00:14:25,990 Here we have a 10-story building, and that 10-story 236 00:14:25,990 --> 00:14:29,740 building is discretized, idealized, as shown here. 237 00:14:29,740 --> 00:14:34,340 We also want to include the elasticity of the foundation, 238 00:14:34,340 --> 00:14:37,100 and that is modeled very simply here as an assemblage 239 00:14:37,100 --> 00:14:39,320 of springs. 240 00:14:39,320 --> 00:14:45,700 We deal here with two substructures, shown in red 241 00:14:45,700 --> 00:14:49,390 for the substructure degrees of freedom and in black for 242 00:14:49,390 --> 00:14:52,660 the master degrees of freedom. 243 00:14:52,660 --> 00:14:57,710 This is also part of the master degrees of freedom. 244 00:14:57,710 --> 00:15:03,080 Notice that here we have one substructure-- 245 00:15:03,080 --> 00:15:07,650 in other words, this part here is nothing else than that part 246 00:15:07,650 --> 00:15:08,900 right here. 247 00:15:10,840 --> 00:15:12,905 That's what is one substructure. 248 00:15:15,450 --> 00:15:18,750 And of course, you see here another substructure. 249 00:15:18,750 --> 00:15:22,730 But we have a repetition here-- in other words, this 250 00:15:22,730 --> 00:15:25,250 substructure here is identically equal to that 251 00:15:25,250 --> 00:15:27,370 substructure. 252 00:15:27,370 --> 00:15:31,570 This, of course, is the fact because of the particular 253 00:15:31,570 --> 00:15:33,440 design of this building. 254 00:15:33,440 --> 00:15:36,610 The upper stories look the same way as the lower stories. 255 00:15:36,610 --> 00:15:39,870 That's the reason why we can say, this is one substructure 256 00:15:39,870 --> 00:15:44,980 model which can model this part of the structure and that 257 00:15:44,980 --> 00:15:47,420 part of the structure, as well. 258 00:15:47,420 --> 00:15:51,850 If we look at the stiffness matrix, or the effective 259 00:15:51,850 --> 00:15:56,420 stiffness matrix of this structure that he would 260 00:15:56,420 --> 00:15:59,750 encounter in a direct integration solution using 261 00:15:59,750 --> 00:16:02,050 implicit time integration, 262 00:16:02,050 --> 00:16:04,930 schematically, we see the following. 263 00:16:04,930 --> 00:16:08,490 Notice here we have the entries coming from the 264 00:16:08,490 --> 00:16:13,630 substructures, and here we have entries 265 00:16:13,630 --> 00:16:16,720 from the master structure. 266 00:16:16,720 --> 00:16:20,380 This coupling here corresponds also to master degrees of 267 00:16:20,380 --> 00:16:23,540 freedom, and similarly, this coupling corresponds to master 268 00:16:23,540 --> 00:16:25,700 degrees of freedom. 269 00:16:25,700 --> 00:16:28,440 Here we have the displacement vector, the incremental 270 00:16:28,440 --> 00:16:29,770 displacement vector. 271 00:16:29,770 --> 00:16:34,720 And let's see here that these master degrees of freedom of 272 00:16:34,720 --> 00:16:38,710 course correspond to these element entries here. 273 00:16:38,710 --> 00:16:41,570 These master degrees of freedom correspond to these 274 00:16:41,570 --> 00:16:44,650 entries, and these master degrees of freedom correspond 275 00:16:44,650 --> 00:16:46,140 to these entries. 276 00:16:46,140 --> 00:16:51,230 Now what we can do is condense out all the linear degrees of 277 00:16:51,230 --> 00:16:53,360 freedom effects. 278 00:16:53,360 --> 00:16:57,220 And that is achieved as follows. 279 00:16:57,220 --> 00:17:01,390 We recognize that the total effective stiffness matrix 280 00:17:01,390 --> 00:17:06,740 corresponding to time t is made up of this part here, 281 00:17:06,740 --> 00:17:13,720 which is the linear stiffness matrix plus 282 00:17:13,720 --> 00:17:17,410 the total mass matrix. 283 00:17:17,410 --> 00:17:21,190 I should point out here that this K matrix contains all the 284 00:17:21,190 --> 00:17:25,790 linear contributions, where that M matrix here corresponds 285 00:17:25,790 --> 00:17:29,990 to all the contributions from the masses at the master 286 00:17:29,990 --> 00:17:33,440 degrees of freedom and the 287 00:17:33,440 --> 00:17:35,360 substructure degrees of freedom. 288 00:17:35,360 --> 00:17:38,790 Because we're assuming here that the total mass matrix is 289 00:17:38,790 --> 00:17:42,170 linear, we can put all the mass effects right in there. 290 00:17:42,170 --> 00:17:46,540 And then since this is not the total effect of stiffness, we 291 00:17:46,540 --> 00:17:49,060 add the effect of the non-linear element 292 00:17:49,060 --> 00:17:52,110 stiffnesses, right here. 293 00:17:52,110 --> 00:17:55,480 And in this particular analysis, this might be, for 294 00:17:55,480 --> 00:17:58,140 example, the non-linear springs down there. 295 00:17:58,140 --> 00:18:00,510 The non-linear springs modelling the foundation on 296 00:18:00,510 --> 00:18:01,760 the structure. 297 00:18:04,290 --> 00:18:09,560 Well, this can be written in this form. 298 00:18:09,560 --> 00:18:15,260 And if we perform static condensation on this combined 299 00:18:15,260 --> 00:18:21,020 matrix, we arrive at this matrix here. 300 00:18:21,020 --> 00:18:24,810 Notice here we have a condensed degree of freedom 301 00:18:24,810 --> 00:18:29,360 vector carrying only the master degrees of freedom, and 302 00:18:29,360 --> 00:18:33,280 here we have the entries corresponding to the master 303 00:18:33,280 --> 00:18:36,720 degrees of freedom in the tK hat matrix. 304 00:18:36,720 --> 00:18:39,750 We put a c on there, meaning condensed. 305 00:18:39,750 --> 00:18:42,290 After static condensation, this is a matrix 306 00:18:42,290 --> 00:18:45,660 that has been reached. 307 00:18:45,660 --> 00:18:50,530 These black squares here correspond to as the black 308 00:18:50,530 --> 00:18:52,570 squares that you saw earlier. 309 00:18:52,570 --> 00:18:57,110 Of course, they have now been changed a bit because of the 310 00:18:57,110 --> 00:19:00,480 static condensation, because the effect of the linear 311 00:19:00,480 --> 00:19:03,900 degrees of freedom have been carried into 312 00:19:03,900 --> 00:19:05,150 these entries here. 313 00:19:08,130 --> 00:19:12,760 The major steps in the solution are as follows. 314 00:19:12,760 --> 00:19:16,220 Prior to the step-by-step solution, we establish the K 315 00:19:16,220 --> 00:19:22,140 hat, the linear effective stiffness matrix containing K 316 00:19:22,140 --> 00:19:26,510 plus 4 over delta t squared m. 317 00:19:26,510 --> 00:19:31,400 And this K hat matrix is then used to statically condense 318 00:19:31,400 --> 00:19:36,310 out all internal substructure degrees of freedom to obtain 319 00:19:36,310 --> 00:19:40,810 this K hat c, c standing for condensed matrix. 320 00:19:40,810 --> 00:19:45,390 We note that the actual effective stiffness matrix 321 00:19:45,390 --> 00:19:49,490 including the non-linear effects is then obtained by 322 00:19:49,490 --> 00:19:53,510 taking the linear effective stiffness matrix and adding 323 00:19:53,510 --> 00:19:56,940 the non-linear stiffness effects to it. 324 00:19:56,940 --> 00:19:59,960 Notice, once again, this matrix is obtained from this 325 00:19:59,960 --> 00:20:02,770 matrix here, as I mentioned earlier. 326 00:20:02,770 --> 00:20:06,590 We use this matrix to statically condense out all 327 00:20:06,590 --> 00:20:09,930 the linear degrees of freedom. 328 00:20:09,930 --> 00:20:15,950 The solution is then obtained in each time step as follows. 329 00:20:15,950 --> 00:20:20,830 We update the condensed matrix K hat c-- 330 00:20:20,830 --> 00:20:23,070 there should be no t here-- 331 00:20:23,070 --> 00:20:26,090 K hat c for non-linearities. 332 00:20:26,090 --> 00:20:30,440 That means we're adding the tK matrix to it, corresponding to 333 00:20:30,440 --> 00:20:34,520 the non-linearities that have occurred in the system. 334 00:20:34,520 --> 00:20:37,670 We establish the complete load vector for all degrees of 335 00:20:37,670 --> 00:20:41,620 freedom, and condense out the substructure internal degrees 336 00:20:41,620 --> 00:20:43,190 of freedom. 337 00:20:43,190 --> 00:20:47,000 This is important to recognize. 338 00:20:47,000 --> 00:20:51,920 We have to calculate the complete load vector because 339 00:20:51,920 --> 00:20:57,970 of the intertia effects on the internal substructure degrees 340 00:20:57,970 --> 00:20:58,710 of freedom. 341 00:20:58,710 --> 00:21:01,240 The intertia forces on the internal substructure degrees 342 00:21:01,240 --> 00:21:05,300 of freedom carry forces to the master degrees of freedom that 343 00:21:05,300 --> 00:21:10,660 we have to take into account in each time step and in each 344 00:21:10,660 --> 00:21:13,020 equilibrium iteration. 345 00:21:13,020 --> 00:21:16,810 We then solve for the master degrees of freedom, using the 346 00:21:16,810 --> 00:21:23,650 small system size, the tK hat c, the velocities, 347 00:21:23,650 --> 00:21:27,680 accelerations, and calculate all the substructural degrees 348 00:21:27,680 --> 00:21:28,650 of freedom-- 349 00:21:28,650 --> 00:21:32,450 velocities, displacements, and accelerations. 350 00:21:32,450 --> 00:21:38,080 These are required here to calculate the load vector 351 00:21:38,080 --> 00:21:41,460 corresponding to all degrees of freedom. 352 00:21:41,460 --> 00:21:44,460 Schematically, the solution procedure, therefore, proceeds 353 00:21:44,460 --> 00:21:45,620 as follows. 354 00:21:45,620 --> 00:21:49,690 For each time step and each iteration, we start with the 355 00:21:49,690 --> 00:21:53,320 displacement velocities and accelerations at time t. 356 00:21:53,320 --> 00:21:56,090 We calculate the effective load vector, corresponding to 357 00:21:56,090 --> 00:21:59,170 time t plus delta t, the way I've been 358 00:21:59,170 --> 00:22:02,060 talking about it earlier. 359 00:22:02,060 --> 00:22:06,160 And then we condense that load vector, as shown here, to the 360 00:22:06,160 --> 00:22:11,800 condensed load vector, c meaning condensation. 361 00:22:11,800 --> 00:22:17,080 We then, from this load vector, calculate with the 362 00:22:17,080 --> 00:22:21,270 effective tension stiffness matrix, the t plus delta t Uc 363 00:22:21,270 --> 00:22:25,810 vector, and this vector is used to calculate the 364 00:22:25,810 --> 00:22:28,750 displacement, velocity, and accelerations of all the 365 00:22:28,750 --> 00:22:31,860 degrees of freedom at time t plus delta t. 366 00:22:31,860 --> 00:22:34,750 Of course here I'm showing only the process for a 367 00:22:34,750 --> 00:22:36,210 typical time step. 368 00:22:36,210 --> 00:22:40,990 If you iterate, as we usually do in non-linear analysis, of 369 00:22:40,990 --> 00:22:43,820 course, then you would have to also introduce the iteration 370 00:22:43,820 --> 00:22:47,150 counter here in these vectors. 371 00:22:47,150 --> 00:22:51,570 Notice that using this solution procedure, we get the 372 00:22:51,570 --> 00:22:57,190 same results as if we had not performed the substructuring. 373 00:22:57,190 --> 00:23:00,210 In other words, there is no additional assumption that 374 00:23:00,210 --> 00:23:02,940 we're introducing into the solution. 375 00:23:02,940 --> 00:23:04,280 You get the same results. 376 00:23:04,280 --> 00:23:08,890 However, the solution cost is reduced, because we are 377 00:23:08,890 --> 00:23:13,120 dealing with one substructure, a typical substructure, only 378 00:23:13,120 --> 00:23:15,630 once regarding the factorization-- 379 00:23:15,630 --> 00:23:18,410 the LDL transpose factorization and static 380 00:23:18,410 --> 00:23:22,290 condensation that has to be performed. 381 00:23:22,290 --> 00:23:25,780 Let us now look at some example solutions. 382 00:23:25,780 --> 00:23:29,570 And the first example that I'd like to consider with you is 383 00:23:29,570 --> 00:23:33,360 the solution of the wave propagation in a rod. 384 00:23:33,360 --> 00:23:38,560 Here we have a uniform freely floating rod, and we apply to 385 00:23:38,560 --> 00:23:43,100 the left end of that rod a force that varies as shown 386 00:23:43,100 --> 00:23:47,140 here as step load of 1000 newtons. 387 00:23:47,140 --> 00:23:50,840 The material data for the rod are given here. 388 00:23:50,840 --> 00:23:54,440 We want to perform a linear elastic analysis, but although 389 00:23:54,440 --> 00:23:58,050 this is a linear analysis, we still will learn quite a bit 390 00:23:58,050 --> 00:23:59,850 from it regarding non-linear analysis. 391 00:24:02,490 --> 00:24:05,820 We want to consider the compressive force at a point 392 00:24:05,820 --> 00:24:10,980 at the center of the rod, right there. 393 00:24:10,980 --> 00:24:16,120 And the exact solution for this force at that point is 394 00:24:16,120 --> 00:24:18,600 plotted here. 395 00:24:18,600 --> 00:24:24,030 You have, once again, apply 1000 newtons, and here we have 396 00:24:24,030 --> 00:24:25,920 the time scale. 397 00:24:25,920 --> 00:24:31,320 Notice up to 1/2 t star, there's no force at point a, 398 00:24:31,320 --> 00:24:34,180 because the applied force in fact has to 399 00:24:34,180 --> 00:24:36,230 propagate through the rod. 400 00:24:36,230 --> 00:24:41,380 Then we have the full force up to 3/2 t star, and then no 401 00:24:41,380 --> 00:24:45,300 force again, where t star is the time for the stress waves 402 00:24:45,300 --> 00:24:48,460 to travel through the complete rod. 403 00:24:48,460 --> 00:24:51,500 This is the analytical solution, well-known, 404 00:24:51,500 --> 00:24:55,290 well-established, and we want to compare with this solution 405 00:24:55,290 --> 00:24:58,330 when we perform our finite element analysis. 406 00:24:58,330 --> 00:25:03,180 For the finite element analysis, we use a mesh off 407 00:25:03,180 --> 00:25:10,670 2-noded truss elements, as shown here, and we measure the 408 00:25:10,670 --> 00:25:14,560 force at the midpoint of the total structure by the 409 00:25:14,560 --> 00:25:19,230 compressive force in this element. 410 00:25:19,230 --> 00:25:23,450 We use in this analysis first essential difference method, 411 00:25:23,450 --> 00:25:27,740 and recognize that the time step that we use has to be 412 00:25:27,740 --> 00:25:30,850 small or equal to the critical time step, as we discussed in 413 00:25:30,850 --> 00:25:32,450 the previous lecture. 414 00:25:32,450 --> 00:25:35,820 And that critical time step is given by Le/c. 415 00:25:35,820 --> 00:25:38,270 Please refer to as the previous lecture notes, and 416 00:25:38,270 --> 00:25:41,860 you would see this formula there, which is nothing else 417 00:25:41,860 --> 00:25:45,460 than t star divided by the number of elements. 418 00:25:45,460 --> 00:25:49,520 We recall that if delta t is greater than delta t critical, 419 00:25:49,520 --> 00:25:52,730 we would produce an unstable solution. 420 00:25:52,730 --> 00:25:54,190 We also need to deal with the initial 421 00:25:54,190 --> 00:25:56,190 conditions for this problem. 422 00:25:56,190 --> 00:26:02,140 And we recognize that if the initial displacements are 0 423 00:26:02,140 --> 00:26:06,360 and R is given, then we can calculate the initial 424 00:26:06,360 --> 00:26:08,280 acceleration, as shown here. 425 00:26:11,540 --> 00:26:16,170 Using now a time step equal to the critical time step, we 426 00:26:16,170 --> 00:26:19,160 obtain the correct result, and I pointed that out in the 427 00:26:19,160 --> 00:26:22,360 earlier lecture already, that we in fact obtain the exact 428 00:26:22,360 --> 00:26:25,530 solution for any number of finite elements 429 00:26:25,530 --> 00:26:28,270 to model the rod. 430 00:26:28,270 --> 00:26:32,650 And here is a finite element solution, compared with the 431 00:26:32,650 --> 00:26:33,780 exact solution. 432 00:26:33,780 --> 00:26:38,860 A wonderful solution, right on the exact solution. 433 00:26:38,860 --> 00:26:41,920 Using, however, now a smaller time step-- 434 00:26:41,920 --> 00:26:43,480 and this is frequently of surprise-- 435 00:26:46,300 --> 00:26:51,070 using a smaller time step, we get a much worse result. 436 00:26:51,070 --> 00:26:55,230 Here we use a time step of 1/2 delta t critical. 437 00:26:55,230 --> 00:26:57,730 Half the time step that we used before. 438 00:26:57,730 --> 00:27:01,150 And just solution looks actually very bad right here. 439 00:27:04,280 --> 00:27:08,820 Frequently one sort of intuitively thinks, well, if I 440 00:27:08,820 --> 00:27:12,020 use a smaller time step, I should get better results. 441 00:27:12,020 --> 00:27:16,880 My costs go up, and surely I should get better results. 442 00:27:16,880 --> 00:27:20,250 Well, it turns out that using the explicit time integration 443 00:27:20,250 --> 00:27:23,930 scheme for this type of problem, this wave propagation 444 00:27:23,930 --> 00:27:28,490 problem, with a smaller time step, more cost, higher cost, 445 00:27:28,490 --> 00:27:30,410 you get worse results. 446 00:27:30,410 --> 00:27:33,580 That is really something to keep in mind. 447 00:27:33,580 --> 00:27:35,482 Now let us consider the trapezoidal rule. 448 00:27:37,990 --> 00:27:42,130 The trapezoidal rule-- we discussed this method in the 449 00:27:42,130 --> 00:27:43,020 earlier lecture. 450 00:27:43,020 --> 00:27:46,320 And I'd like to refer you once again to the notes 451 00:27:46,320 --> 00:27:49,410 corresponding to that lecture, where we have identified that 452 00:27:49,410 --> 00:27:51,810 a stable solution would be obtained with any 453 00:27:51,810 --> 00:27:53,440 choice of delta t. 454 00:27:53,440 --> 00:27:56,500 Remember, this is a linear analysis, and the trapezoidal 455 00:27:56,500 --> 00:27:59,170 rule is unconditionally stable. 456 00:27:59,170 --> 00:28:02,790 We could use a consistent or lumped mass matrix. 457 00:28:02,790 --> 00:28:05,560 And in this particular case, we use a lumped mass matrix, 458 00:28:05,560 --> 00:28:08,470 because we also like to compare our results with those 459 00:28:08,470 --> 00:28:12,610 that we obtained using the central difference method. 460 00:28:12,610 --> 00:28:17,350 If we use the trapezoidal rule with the time step that we 461 00:28:17,350 --> 00:28:20,170 employ in the central difference method-- 462 00:28:20,170 --> 00:28:23,420 in fact, the critical time step that we used there-- 463 00:28:23,420 --> 00:28:26,880 and we also use the initial conditions computed as shown 464 00:28:26,880 --> 00:28:30,930 here, the way we did it for the central difference method 465 00:28:30,930 --> 00:28:35,820 solution, we find that the solution is very inaccurate 466 00:28:35,820 --> 00:28:38,660 when compared to the exact solution. 467 00:28:38,660 --> 00:28:43,720 Here, of course, we also use the 10-element mesh. 468 00:28:43,720 --> 00:28:50,590 If we use 0 initial conditions for the trapezoidal rule, we 469 00:28:50,590 --> 00:28:52,430 get this solution. 470 00:28:52,430 --> 00:28:55,390 Not much different, really, from the solution that I just 471 00:28:55,390 --> 00:28:56,620 showed you. 472 00:28:56,620 --> 00:29:00,120 So whether we use 0 initial conditions or the initial 473 00:29:00,120 --> 00:29:04,210 conditions corresponding to the equation MU double dot is 474 00:29:04,210 --> 00:29:07,420 equal to R, which we have seen on the previous view graph, 475 00:29:07,420 --> 00:29:09,390 makes very little difference when using 476 00:29:09,390 --> 00:29:10,570 the trapezoidal rule. 477 00:29:10,570 --> 00:29:13,480 It makes a lot of difference using the explicit central 478 00:29:13,480 --> 00:29:14,750 difference method. 479 00:29:14,750 --> 00:29:17,940 And for that reason, in fact, we have to use the proper 480 00:29:17,940 --> 00:29:22,540 initial conditions obtained from MU double dot equals R at 481 00:29:22,540 --> 00:29:26,640 time 0 when we use the central difference method. 482 00:29:26,640 --> 00:29:30,370 Using now twice the time step of what we just used, namely, 483 00:29:30,370 --> 00:29:34,220 delta t equal to 2 times delta t, critical of the central 484 00:29:34,220 --> 00:29:37,360 difference method, this solution is stable, but very 485 00:29:37,360 --> 00:29:38,210 inaccurate. 486 00:29:38,210 --> 00:29:41,120 It is stable because we're using an unconditionally 487 00:29:41,120 --> 00:29:42,540 stable operator. 488 00:29:42,540 --> 00:29:44,450 We discussed it in the previous lectures. 489 00:29:44,450 --> 00:29:48,790 Note that this is the solution we're obtaining when compared 490 00:29:48,790 --> 00:29:51,600 to the exact solution. 491 00:29:51,600 --> 00:29:55,130 Now we look at the solution where we have half the time 492 00:29:55,130 --> 00:30:00,620 step off the critical size, corresponding to the central 493 00:30:00,620 --> 00:30:01,920 difference method. 494 00:30:01,920 --> 00:30:05,630 And you see here, the solution, again, very 495 00:30:05,630 --> 00:30:09,790 inaccurate when compared to the exact solution, but of 496 00:30:09,790 --> 00:30:13,300 course stable, because the method is always stable for 497 00:30:13,300 --> 00:30:15,980 any time step delta t. 498 00:30:15,980 --> 00:30:18,730 The question now could, of course, be, and should be, 499 00:30:18,730 --> 00:30:22,440 what happens if we take more elements to discretize the 500 00:30:22,440 --> 00:30:23,160 whole structure? 501 00:30:23,160 --> 00:30:25,480 Here we use the 10-element model. 502 00:30:25,480 --> 00:30:27,480 Why not use more elements? 503 00:30:27,480 --> 00:30:32,060 And here we now look at the solution obtained with the 504 00:30:32,060 --> 00:30:35,870 essential difference method when we use 100 elements to 505 00:30:35,870 --> 00:30:38,170 model the total bar. 506 00:30:38,170 --> 00:30:43,670 Here the exact solution in red, and that exact solution 507 00:30:43,670 --> 00:30:48,040 is reached if you integrate the equilibrium equations with 508 00:30:48,040 --> 00:30:49,400 this time step. 509 00:30:49,400 --> 00:30:50,790 We discussed that earlier. 510 00:30:50,790 --> 00:30:54,460 That exact solution is always obtained for any number of 511 00:30:54,460 --> 00:30:58,030 elements that you use to discretize a structure, 512 00:30:58,030 --> 00:31:01,950 provided you use delta t equals delta t critical. 513 00:31:01,950 --> 00:31:06,180 If you then use a smaller time step, in fact, delta t equals 514 00:31:06,180 --> 00:31:11,120 1/2 delta t critical was used here, we obtain a solution 515 00:31:11,120 --> 00:31:15,700 that is not the exact solution, but as you can see 516 00:31:15,700 --> 00:31:19,490 here, it oscillates for this time step this amount about 517 00:31:19,490 --> 00:31:22,700 the exact solution. 518 00:31:22,700 --> 00:31:27,060 So we have much the same observations that we have seen 519 00:31:27,060 --> 00:31:30,780 already earlier, namely, when the time step is decreased in 520 00:31:30,780 --> 00:31:32,515 the central difference method solution-- 521 00:31:35,340 --> 00:31:38,850 from the critical time step, here we see a decrease-- 522 00:31:38,850 --> 00:31:43,220 the solution that one might expect to be better actually 523 00:31:43,220 --> 00:31:44,010 deteriorates. 524 00:31:44,010 --> 00:31:45,060 Of course, it couldn't be better 525 00:31:45,060 --> 00:31:46,240 than the exact solution. 526 00:31:46,240 --> 00:31:49,140 But we see here a deterioration, and a quite 527 00:31:49,140 --> 00:31:51,620 significant one. 528 00:31:51,620 --> 00:31:55,170 If you use the trapezoidal rule, this is the solution 529 00:31:55,170 --> 00:32:00,190 that you obtain with the time step equal to the critical 530 00:32:00,190 --> 00:32:02,550 time step of the central difference method. 531 00:32:02,550 --> 00:32:06,980 Again, here, we're using the 100-element mesh. 532 00:32:06,980 --> 00:32:10,250 The solution is certainly better than what we have seen 533 00:32:10,250 --> 00:32:14,300 for the 10-element mesh, but it does not give 534 00:32:14,300 --> 00:32:16,080 us the exact solution-- 535 00:32:16,080 --> 00:32:19,060 in other words, a [UNINTELLIGIBLE] wavefront as 536 00:32:19,060 --> 00:32:21,850 we would like to obtain. 537 00:32:21,850 --> 00:32:25,630 Of course, this is a very difficult problem to solve, 538 00:32:25,630 --> 00:32:29,960 because we have infinitely steep wave fronts here and 539 00:32:29,960 --> 00:32:32,230 there, basically, to model. 540 00:32:32,230 --> 00:32:35,940 And an overshoot would obviously obtain it, but we 541 00:32:35,940 --> 00:32:38,570 would like to have said overshoot and undershoot, with 542 00:32:38,570 --> 00:32:41,590 this oscillation here as small as possible, of course, with 543 00:32:41,590 --> 00:32:45,060 the solution scheme that we employ. 544 00:32:45,060 --> 00:32:48,400 There is one more important consideration that I'd like to 545 00:32:48,400 --> 00:32:50,470 bring to your attention. 546 00:32:50,470 --> 00:32:56,030 Here we use a 2-dimensional element mesh to model the bar. 547 00:32:56,030 --> 00:33:01,440 Notice these are 4-node elements, plane stress 548 00:33:01,440 --> 00:33:03,640 elements here to model the bar. 549 00:33:03,640 --> 00:33:06,180 And in this particular case, we use the 550 00:33:06,180 --> 00:33:08,500 central difference methods. 551 00:33:08,500 --> 00:33:11,440 The thickness is given here, the material 552 00:33:11,440 --> 00:33:14,690 data are given here. 553 00:33:14,690 --> 00:33:21,350 For this solution, notice that this distance here is half the 554 00:33:21,350 --> 00:33:25,480 distance that we have here. 555 00:33:25,480 --> 00:33:28,490 And that is something that leads to a phenomenon that one 556 00:33:28,490 --> 00:33:29,810 has to be aware of. 557 00:33:29,810 --> 00:33:33,660 Notice our the critical time step that we identified for 558 00:33:33,660 --> 00:33:37,720 the one-dimensional wave propagation here would, of 559 00:33:37,720 --> 00:33:41,650 course, be calculated depending on this length. 560 00:33:41,650 --> 00:33:45,900 However, we have to allow for the fact that there is also a 561 00:33:45,900 --> 00:33:50,170 possibility of wave propagation this way, and the 562 00:33:50,170 --> 00:33:54,970 actual critical time step that we have to use in this problem 563 00:33:54,970 --> 00:33:58,560 is given by this length, by the shortest element distance. 564 00:33:58,560 --> 00:34:02,590 Please refer to the previous lecture, in which I discussed 565 00:34:02,590 --> 00:34:03,440 this already. 566 00:34:03,440 --> 00:34:08,940 You have to use the shortest distance between nodal points 567 00:34:08,940 --> 00:34:15,880 within an element to determine the critical time step. 568 00:34:15,880 --> 00:34:21,639 So this relationship here does not hold anymore, because this 569 00:34:21,639 --> 00:34:23,699 distance is less than that. 570 00:34:23,699 --> 00:34:28,560 And if you were to use delta t equal to t star over 10 571 00:34:28,560 --> 00:34:34,139 elements, which was a critical time step that we used when we 572 00:34:34,139 --> 00:34:38,080 solved the truss element model, then we would find that 573 00:34:38,080 --> 00:34:39,560 the solution diverges. 574 00:34:39,560 --> 00:34:42,510 In fact, in element 5, which was the element that we always 575 00:34:42,510 --> 00:34:46,420 used to measure our stress, we find that tau zz, the absolute 576 00:34:46,420 --> 00:34:51,280 value of tau zz, is greater than 1000 newtons divided by 577 00:34:51,280 --> 00:34:55,440 the area at this time. 578 00:34:55,440 --> 00:34:58,550 Notice that it was 1000 newton that we applied. 579 00:34:58,550 --> 00:35:04,490 So a stress that should be 0, really, because tau zz should 580 00:35:04,490 --> 00:35:07,210 really be 0, it's the stress in the transverse direction, 581 00:35:07,210 --> 00:35:12,010 the wave travels this way, and tau zz is acting that way. 582 00:35:12,010 --> 00:35:13,260 This stress should be 0. 583 00:35:13,260 --> 00:35:17,790 We find, actually, that it grows, and it grows because we 584 00:35:17,790 --> 00:35:21,720 are using a time step that is larger than the critical time 585 00:35:21,720 --> 00:35:28,180 step that should be used for this problem. 586 00:35:28,180 --> 00:35:32,460 Let us now look at another example, a second example from 587 00:35:32,460 --> 00:35:34,390 which we can also learn quite a bit. 588 00:35:34,390 --> 00:35:38,500 Here we have a very simple 3 degree of freedom system. 589 00:35:38,500 --> 00:35:42,540 These are the degrees of freedom, x1, x2, x3. 590 00:35:42,540 --> 00:35:46,730 Notice we have masses m here, and we have a linear spring 591 00:35:46,730 --> 00:35:49,660 connecting these masses, a non-linear spring connecting 592 00:35:49,660 --> 00:35:53,460 these masses, and a linear spring connecting this mass to 593 00:35:53,460 --> 00:35:56,400 the rigid support. 594 00:35:56,400 --> 00:36:02,640 For this problem, the linear spring has is distance here, 595 00:36:02,640 --> 00:36:06,050 and the non-linear spring is characterized by this force 596 00:36:06,050 --> 00:36:08,420 displacement relationship. 597 00:36:08,420 --> 00:36:11,660 Notice there's a kink going from 1 to 100. 598 00:36:11,660 --> 00:36:18,000 So the spring gets suddenly, at a particular distance, or 599 00:36:18,000 --> 00:36:20,870 rather deformations, I should say-- 600 00:36:20,870 --> 00:36:23,110 displacement x2 minus x3-- 601 00:36:23,110 --> 00:36:27,170 at a particular deformation, the spring becomes 100 times 602 00:36:27,170 --> 00:36:31,940 as stiff as prior to that point. 603 00:36:31,940 --> 00:36:37,050 Well, we can calculate a critical time step based on 604 00:36:37,050 --> 00:36:42,770 this spring stiffness here, the linear part, or rather, 605 00:36:42,770 --> 00:36:46,960 this spring stiffness for small displacements. 606 00:36:46,960 --> 00:36:50,710 And that critical time step for this problem is given as 607 00:36:50,710 --> 00:36:53,740 1.11 seconds. 608 00:36:53,740 --> 00:36:59,830 As long as the response of this system is such that we 609 00:36:59,830 --> 00:37:03,250 would, in this non-linear spring, never go around this 610 00:37:03,250 --> 00:37:06,990 kink, we would deal with a linear system, and this would 611 00:37:06,990 --> 00:37:10,250 be the proper critical time step to use. 612 00:37:10,250 --> 00:37:14,860 However, as soon as we are going around this kink here, 613 00:37:14,860 --> 00:37:18,640 due to the fact that the forces that I applied, or the 614 00:37:18,640 --> 00:37:21,310 initial conditions that we are applying to these degrees of 615 00:37:21,310 --> 00:37:24,370 freedom, are large enough, as soon as we go around this 616 00:37:24,370 --> 00:37:28,700 kink, our new critical time step is based on this spring 617 00:37:28,700 --> 00:37:33,600 stiffness, and if this given right down here. 618 00:37:33,600 --> 00:37:35,560 This is a new critical time step. 619 00:37:35,560 --> 00:37:40,750 So if we now perform essential difference method solutions of 620 00:37:40,750 --> 00:37:45,590 this problem, and if the solution brings us into this 621 00:37:45,590 --> 00:37:50,510 regime, it is very important to use, as a critical time 622 00:37:50,510 --> 00:37:54,080 step, 0.14 seconds. 623 00:37:54,080 --> 00:37:58,480 It is very important, even if we're only in this regime for 624 00:37:58,480 --> 00:38:00,800 just a few steps. 625 00:38:00,800 --> 00:38:04,120 You see, what might happen in the response is that most of 626 00:38:04,120 --> 00:38:06,310 the time, you're in this regime. 627 00:38:06,310 --> 00:38:10,350 But then all of a sudden, for just a few steps, you're 628 00:38:10,350 --> 00:38:12,870 getting into this regime, and because of this high 629 00:38:12,870 --> 00:38:16,260 stiffness, you're thrown back into this regime. 630 00:38:16,260 --> 00:38:20,470 If your time step that you use in the time integration is 631 00:38:20,470 --> 00:38:28,960 then larger than that value, you are having an instability 632 00:38:28,960 --> 00:38:32,090 for just these few time steps, and that introduces an error 633 00:38:32,090 --> 00:38:36,505 in the solution which, if you don't have a comparison with 634 00:38:36,505 --> 00:38:40,430 an exact analytical solution, will of course be an error 635 00:38:40,430 --> 00:38:43,770 that you not necessarily will see, directly. 636 00:38:43,770 --> 00:38:45,620 And that's one of the important points. 637 00:38:45,620 --> 00:38:49,410 You see, in linear analysis, when you use the essential 638 00:38:49,410 --> 00:38:53,630 difference method in linear analysis with a time step that 639 00:38:53,630 --> 00:38:57,680 is larger than the critical time step, you will, after a 640 00:38:57,680 --> 00:39:01,320 few time step integrations, see that the response becomes 641 00:39:01,320 --> 00:39:05,930 very large, and the instability is obvious. 642 00:39:05,930 --> 00:39:08,230 The instability is really obvious when you apply the 643 00:39:08,230 --> 00:39:10,850 central difference method with a time step that is larger 644 00:39:10,850 --> 00:39:14,730 than the critical time step to the solution of the system. 645 00:39:14,730 --> 00:39:17,400 However, in non-linear analysis, as I mentioned in 646 00:39:17,400 --> 00:39:21,990 the earlier lectures, and you are now well aware of it, in 647 00:39:21,990 --> 00:39:26,260 non-linear analysis, the frequencies change. 648 00:39:26,260 --> 00:39:30,230 The critical frequency being the largest frequency in the 649 00:39:30,230 --> 00:39:32,220 system changes. 650 00:39:32,220 --> 00:39:35,580 The smallest period of the system changes, as we are 651 00:39:35,580 --> 00:39:38,160 seeing here in this view graph, also. 652 00:39:38,160 --> 00:39:44,770 And it is important that your time step is always smaller 653 00:39:44,770 --> 00:39:48,540 than the critical time step based on the highest frequency 654 00:39:48,540 --> 00:39:51,560 that you can ever encounter in the solution. 655 00:39:51,560 --> 00:39:54,150 If that is not the case, then you might have only an 656 00:39:54,150 --> 00:39:57,820 instability for a few time steps, and the solution may 657 00:39:57,820 --> 00:40:02,390 not have time to blow up, or give you very large 658 00:40:02,390 --> 00:40:03,500 displacements. 659 00:40:03,500 --> 00:40:06,590 So the instability will not be obvious. 660 00:40:06,590 --> 00:40:09,680 Instead, the instability over these few time steps only 661 00:40:09,680 --> 00:40:15,220 introduces an error in the solution, an error of which 662 00:40:15,220 --> 00:40:17,140 you really don't know the magnitude. 663 00:40:17,140 --> 00:40:21,390 And after these few steps have been completed, you're using 664 00:40:21,390 --> 00:40:23,890 now a time step that is again stable, 665 00:40:23,890 --> 00:40:25,410 giving a stable solution. 666 00:40:25,410 --> 00:40:27,400 That error which we have accumulated, of course, is 667 00:40:27,400 --> 00:40:30,110 carried forward, but you have not seen a 668 00:40:30,110 --> 00:40:31,670 blow-up of the solution. 669 00:40:31,670 --> 00:40:35,560 And that's one important consideration when you apply 670 00:40:35,560 --> 00:40:39,210 the central difference method in non-linear analysis. 671 00:40:39,210 --> 00:40:45,420 Of course, the important point, the sort of ground rule 672 00:40:45,420 --> 00:40:49,780 that you have to observe, is that delta t must always be 673 00:40:49,780 --> 00:40:54,520 smaller than the critical time step of the mesh where that 674 00:40:54,520 --> 00:40:58,020 critical time step changes throughout the solution. 675 00:40:58,020 --> 00:40:59,690 Let's look at this example here. 676 00:40:59,690 --> 00:41:05,080 What happens if we violate that condition? 677 00:41:05,080 --> 00:41:08,850 I show you here the solution with delta t equals 0.15 678 00:41:08,850 --> 00:41:12,740 seconds, compared to the solution with 679 00:41:12,740 --> 00:41:16,230 delta t 0.05 seconds. 680 00:41:16,230 --> 00:41:22,330 Notice that the response of the right math here for both 681 00:41:22,330 --> 00:41:27,860 of these time steps is slightly different, but we can 682 00:41:27,860 --> 00:41:30,720 say that for this time step, we don't get displacements 683 00:41:30,720 --> 00:41:33,760 that are very much different than we obtain 684 00:41:33,760 --> 00:41:35,650 with this time step. 685 00:41:35,650 --> 00:41:41,380 If we look at the response of the center mass, there are 686 00:41:41,380 --> 00:41:44,920 more differences in displacements. 687 00:41:44,920 --> 00:41:50,650 And if we look at the response of the left mass, in fact, we 688 00:41:50,650 --> 00:41:52,550 have even larger differences. 689 00:41:52,550 --> 00:41:58,000 But the important point is that if you only have this 690 00:41:58,000 --> 00:42:05,510 solution here corresponding to the time step 0.15 seconds, 691 00:42:05,510 --> 00:42:08,750 which is larger than the critical time step, you would 692 00:42:08,750 --> 00:42:12,730 certainly not expect that this solution was unstable. 693 00:42:12,730 --> 00:42:16,585 It in fact was only unstable for a very few time steps. 694 00:42:16,585 --> 00:42:19,960 There is an error accumulation, 695 00:42:19,960 --> 00:42:23,190 fairly large here. 696 00:42:23,190 --> 00:42:24,130 In fact, very large. 697 00:42:24,130 --> 00:42:26,880 You can see that these displacements are only about 698 00:42:26,880 --> 00:42:30,170 50% of those displacements. 699 00:42:30,170 --> 00:42:34,860 But you cannot see directly an instability here if you don't 700 00:42:34,860 --> 00:42:36,110 know this solution. 701 00:42:39,770 --> 00:42:43,460 Let's look at the force in the center truss. 702 00:42:43,460 --> 00:42:48,130 Here we list the time, here we list the solution for delta t, 703 00:42:48,130 --> 00:42:51,230 0.05 seconds. 704 00:42:51,230 --> 00:42:55,940 And here we list the solution for delta t 0.15 seconds. 705 00:42:55,940 --> 00:42:59,380 And you can see here that at the beginning, the forces 706 00:42:59,380 --> 00:43:02,480 compare quite well, and then the error 707 00:43:02,480 --> 00:43:04,790 becomes larger and larger. 708 00:43:04,790 --> 00:43:08,480 In fact, here we have an opposite sign in the force. 709 00:43:08,480 --> 00:43:11,230 This being here, say, the accurate solution with a small 710 00:43:11,230 --> 00:43:14,280 enough time step, and this being, of course, a grossly 711 00:43:14,280 --> 00:43:15,310 inaccurate solution. 712 00:43:15,310 --> 00:43:20,350 But it is not evident that this is a grossly inaccurate 713 00:43:20,350 --> 00:43:23,670 solution if you don't know this solution. 714 00:43:23,670 --> 00:43:26,170 And that's the point I wanted to make with this example. 715 00:43:28,870 --> 00:43:31,790 The next example that I like to discuss with you Is the 716 00:43:31,790 --> 00:43:34,260 analysis of a tower. 717 00:43:34,260 --> 00:43:38,770 The tower is shown here, subjected to a pressure that 718 00:43:38,770 --> 00:43:40,970 is induced by a blast. 719 00:43:40,970 --> 00:43:43,550 Here we have some geometric and material data 720 00:43:43,550 --> 00:43:46,510 regarding the tower. 721 00:43:46,510 --> 00:43:51,670 The applied load on the tower is shown here. 722 00:43:51,670 --> 00:43:56,970 Notice it's a load that we assume to be instantaneously 723 00:43:56,970 --> 00:44:01,570 there, and then to decay, as shown here, over 200 724 00:44:01,570 --> 00:44:03,920 milliseconds. 725 00:44:03,920 --> 00:44:08,280 The purpose of the analysis is to determine the displacement 726 00:44:08,280 --> 00:44:12,150 and velocities at the top off the tower, to determine the 727 00:44:12,150 --> 00:44:14,960 moments at the base of the tower. 728 00:44:14,960 --> 00:44:17,180 And in this particular analysis, we use the 729 00:44:17,180 --> 00:44:22,000 trapezoidal rule and the lumped mass matrix to solve 730 00:44:22,000 --> 00:44:24,010 the problem. 731 00:44:24,010 --> 00:44:27,030 We must make a number of decisions here, and two 732 00:44:27,030 --> 00:44:29,610 decisions are most important. 733 00:44:29,610 --> 00:44:34,760 We must choose a mesh that is, in other words, we have to 734 00:44:34,760 --> 00:44:37,560 establish an appropriate finite element discretization 735 00:44:37,560 --> 00:44:41,370 for the tower, and we have to choose a time step for the 736 00:44:41,370 --> 00:44:43,570 time integration. 737 00:44:43,570 --> 00:44:48,330 The choice of mesh and the choice of the time step are 738 00:44:48,330 --> 00:44:53,670 two very closely related items, and that's what I would 739 00:44:53,670 --> 00:44:58,050 like to discuss with you in this example in more depth. 740 00:44:58,050 --> 00:45:01,610 Of course, the time step and the mesh must both depend on 741 00:45:01,610 --> 00:45:02,860 the loading that is supplied. 742 00:45:05,520 --> 00:45:08,720 Let's look at some observations. 743 00:45:08,720 --> 00:45:12,330 The point is that the choice of mesh will determines the 744 00:45:12,330 --> 00:45:16,370 highest natural frequency, and of course, the corresponding 745 00:45:16,370 --> 00:45:21,520 mode shape, that is properly, accurately represented in the 746 00:45:21,520 --> 00:45:23,960 finite element analysis. 747 00:45:23,960 --> 00:45:30,430 The more elements you use, just as a rule of thumb, the 748 00:45:30,430 --> 00:45:34,590 higher the frequency that can be properly represented, that 749 00:45:34,590 --> 00:45:37,530 can be accurately picked up, with the choice 750 00:45:37,530 --> 00:45:39,440 of mesh we're using. 751 00:45:39,440 --> 00:45:44,200 The time step that you're using will of course determine 752 00:45:44,200 --> 00:45:49,030 how high a frequency you're integrating accurately. 753 00:45:49,030 --> 00:45:54,760 In other words, you want to use really enough elements to 754 00:45:54,760 --> 00:46:00,440 have a mesh that is accurate, and you want to use is small 755 00:46:00,440 --> 00:46:05,430 enough time step to integrate the frequency of the mesh that 756 00:46:05,430 --> 00:46:08,530 you need to integrate in order to have an accurate response 757 00:46:08,530 --> 00:46:09,510 prediction. 758 00:46:09,510 --> 00:46:12,810 I used the word "accurate," and we have to now of course 759 00:46:12,810 --> 00:46:16,160 look closely at what I mean by that. 760 00:46:16,160 --> 00:46:20,430 It is most effective, really, to choose the mesh and the 761 00:46:20,430 --> 00:46:24,320 time step such that the highest frequency, accurately 762 00:46:24,320 --> 00:46:29,180 integrated, is equal to the highest frequency accurately 763 00:46:29,180 --> 00:46:32,790 represented by the mesh. 764 00:46:32,790 --> 00:46:37,500 To choose the mesh and the time step, we look at the 765 00:46:37,500 --> 00:46:40,700 applied loading, in other words, the physics of the 766 00:46:40,700 --> 00:46:45,340 problem, and we represent this applied loading 767 00:46:45,340 --> 00:46:47,570 as a Fourier series. 768 00:46:47,570 --> 00:46:51,360 If we do so, that Fourier series will display the 769 00:46:51,360 --> 00:46:55,430 important frequencies that are to be accurately represented 770 00:46:55,430 --> 00:46:58,210 by the mesh. 771 00:46:58,210 --> 00:47:01,670 Here for this particular loading that we're looking at, 772 00:47:01,670 --> 00:47:07,530 we did a Fourier analysis, and you can see that with this 773 00:47:07,530 --> 00:47:12,930 Fourier analysis, we include terms now in our first 774 00:47:12,930 --> 00:47:19,320 analysis, in case one, up to fn 17 hertz, and in case two, 775 00:47:19,320 --> 00:47:24,700 in our second analysis, up fn equal to 30 hertz. 776 00:47:24,700 --> 00:47:30,860 If we truncate our Fourier representation of the loading 777 00:47:30,860 --> 00:47:38,960 at 17 hertz, we obtain this approximation here to the 778 00:47:38,960 --> 00:47:41,560 actual applied noting. 779 00:47:41,560 --> 00:47:47,520 If we truncate our Fourier representation at 30 hertz, we 780 00:47:47,520 --> 00:47:51,390 obtain this Fourier approximation 781 00:47:51,390 --> 00:47:53,280 to the actual loading. 782 00:47:53,280 --> 00:47:58,650 As you can see, surely, as expected, if you include more 783 00:47:58,650 --> 00:48:01,790 terms in the Fourier approximation, your 784 00:48:01,790 --> 00:48:07,060 approximation to the actual loading becomes better. 785 00:48:07,060 --> 00:48:10,940 Now having identified that we want to perform the analysis 786 00:48:10,940 --> 00:48:14,370 with two meshes, first there's a mesh that gives us 787 00:48:14,370 --> 00:48:18,280 frequencies accurate up to 17 hertz, and then with a mesh 788 00:48:18,280 --> 00:48:21,640 that is accurate up to 30 hertz, we have to look for 789 00:48:21,640 --> 00:48:22,810 those meshes. 790 00:48:22,810 --> 00:48:26,900 And we do so by establishing three meshes 791 00:48:26,900 --> 00:48:28,410 that are shown here. 792 00:48:28,410 --> 00:48:30,480 Here we have a 30-element mesh. 793 00:48:30,480 --> 00:48:34,590 Notice that one element spans from the bottom of a floor to 794 00:48:34,590 --> 00:48:37,020 the top off the floor, and across from one 795 00:48:37,020 --> 00:48:39,210 column to the other. 796 00:48:39,210 --> 00:48:44,000 A 60-element mesh, in which we simply have subdivided each 797 00:48:44,000 --> 00:48:48,510 element here in this more coarse mesh into two elements, 798 00:48:48,510 --> 00:48:52,700 and a 120-element mesh, in which we again have taken one 799 00:48:52,700 --> 00:48:56,970 element here and represent it as two elements. 800 00:48:56,970 --> 00:49:02,100 We use these meshes here to calculate frequencies. 801 00:49:02,100 --> 00:49:06,750 And on this view graph here, you see the frequencies 802 00:49:06,750 --> 00:49:10,750 calculated for the 30-element mesh and the 60-element mesh. 803 00:49:10,750 --> 00:49:13,840 Notice mode numbers in this column. 804 00:49:13,840 --> 00:49:17,970 30-element mesh results here, 60-element mesh results here, 805 00:49:17,970 --> 00:49:19,770 in each case, hertz. 806 00:49:19,770 --> 00:49:26,550 We notice that up to this red line, the difference between 807 00:49:26,550 --> 00:49:32,210 the frequencies is less than 1 hertz. 808 00:49:32,210 --> 00:49:34,250 In other words, the difference between the frequencies 809 00:49:34,250 --> 00:49:38,050 obtained from the 30-element mesh and the 60-element mesh 810 00:49:38,050 --> 00:49:41,820 is less than 1 hertz, and that is also right at 811 00:49:41,820 --> 00:49:45,830 the 17 hertz level. 812 00:49:45,830 --> 00:49:49,520 So we can say that the 30-element mesh really 813 00:49:49,520 --> 00:49:55,090 represents the frequencies up to 17 hertz quite accurately. 814 00:49:55,090 --> 00:49:57,980 You should also remember, please, that since we are 815 00:49:57,980 --> 00:50:02,280 using a lumped mass matrix in the analysis, there's no 816 00:50:02,280 --> 00:50:07,150 bounding theorem saying that the frequencies of the 817 00:50:07,150 --> 00:50:13,030 60-element mesh must always be below the frequencies of the 818 00:50:13,030 --> 00:50:15,540 30-element mesh. 819 00:50:15,540 --> 00:50:20,170 For a larger mesh, we can actually have frequencies that 820 00:50:20,170 --> 00:50:27,370 are smaller or larger than for the more coarse mesh. 821 00:50:27,370 --> 00:50:30,530 That of course also applies for when we compare the 822 00:50:30,530 --> 00:50:33,280 60-element mesh result with the 120-element mesh 823 00:50:33,280 --> 00:50:34,980 result later on. 824 00:50:34,980 --> 00:50:37,520 So we say that for the 30-element mesh, we have 825 00:50:37,520 --> 00:50:42,240 accurately represented now the frequencies up to 17 hertz, 826 00:50:42,240 --> 00:50:45,830 and therefore, there is the mesh to use. 827 00:50:45,830 --> 00:50:51,480 We calculate the time step via is this relationship here as 828 00:50:51,480 --> 00:50:57,990 being 0.033 seconds, being the cutoff frequency, and delta t 829 00:50:57,990 --> 00:51:04,130 then is given here as 0.0017 seconds. 830 00:51:04,130 --> 00:51:09,420 This mesh now should be integrated with this time 831 00:51:09,420 --> 00:51:17,440 step, and then we can say is that the mesh up to 17 hertz 832 00:51:17,440 --> 00:51:22,110 represents the frequencies quite accurately, and we have 833 00:51:22,110 --> 00:51:26,620 integrated the frequencies up to 17 hertz quite accurately, 834 00:51:26,620 --> 00:51:30,490 and therefore, we have an accurate solution of the 835 00:51:30,490 --> 00:51:33,210 problem up to 17 hertz-- 836 00:51:33,210 --> 00:51:36,750 and what I mean by this "up to 17 hertz" is that also in the 837 00:51:36,750 --> 00:51:41,310 Fourier representation, we in effect only have represented 838 00:51:41,310 --> 00:51:44,600 the load by the Fourier representation up to 839 00:51:44,600 --> 00:51:47,050 17 hertz, in effect. 840 00:51:47,050 --> 00:51:50,540 Of course, in the actual step-by-step solution, we also 841 00:51:50,540 --> 00:51:54,540 will carry terms above 17 hertz. 842 00:51:54,540 --> 00:52:00,070 But these terms are assumed not to contribute to very much 843 00:52:00,070 --> 00:52:04,180 to the response, because we are not picking them up very 844 00:52:04,180 --> 00:52:07,100 well in the mesh and in the time step selection. 845 00:52:09,780 --> 00:52:16,140 If we proceed now with the same ideas to also look at the 846 00:52:16,140 --> 00:52:20,960 60-element mesh, then we perform there also the 847 00:52:20,960 --> 00:52:25,610 frequency analysis for the 60-element mesh and the 848 00:52:25,610 --> 00:52:34,190 120-element mesh, and we find that up to 30 hertz in the 849 00:52:34,190 --> 00:52:39,680 60-element mesh, we gain represent all the frequencies 850 00:52:39,680 --> 00:52:40,710 accurately. 851 00:52:40,710 --> 00:52:49,300 Notice that above the red line, below these frequencies 852 00:52:49,300 --> 00:52:52,820 here, the difference between these frequencies 853 00:52:52,820 --> 00:52:56,560 is less than 1 hertz. 854 00:52:56,560 --> 00:52:59,350 Above this line, the difference between the 855 00:52:59,350 --> 00:53:03,320 frequencies is larger than 1 hertz, and that 856 00:53:03,320 --> 00:53:05,500 is our cutoff line. 857 00:53:05,500 --> 00:53:09,090 Well, if we want to perform now the direct time 858 00:53:09,090 --> 00:53:12,630 integration of the 60-element mesh, we proceed as for the 859 00:53:12,630 --> 00:53:13,790 30-element mesh. 860 00:53:13,790 --> 00:53:17,770 We calculate our cut-off frequency, our delta t, the 861 00:53:17,770 --> 00:53:23,680 proper time step, 20 steps for the highest frequency or the 862 00:53:23,680 --> 00:53:28,100 smallest period that we want to integrate accurately. 863 00:53:28,100 --> 00:53:30,600 That's the rule of thumb we're using. 864 00:53:30,600 --> 00:53:33,760 And with this time step, we perform the analysis. 865 00:53:33,760 --> 00:53:37,410 Notice a smaller time step would be accurately integrate 866 00:53:37,410 --> 00:53:39,740 frequencies which are not accurately 867 00:53:39,740 --> 00:53:41,130 represented by the mesh. 868 00:53:41,130 --> 00:53:45,120 Therefore, there is no point choosing a smaller time step. 869 00:53:45,120 --> 00:53:48,840 A larger time step would not accurately integrate all 870 00:53:48,840 --> 00:53:50,750 frequencies which are accurately 871 00:53:50,750 --> 00:53:53,810 represented by the mesh. 872 00:53:53,810 --> 00:53:57,540 It is not good to choose a larger time step, because your 873 00:53:57,540 --> 00:54:03,930 mesh is chosen to represent the frequencies up to 30 hertz 874 00:54:03,930 --> 00:54:08,480 quite accurately, and to get all of what you can get out of 875 00:54:08,480 --> 00:54:12,890 the mesh, so to say, you want to take a time step that is in 876 00:54:12,890 --> 00:54:16,060 accordance selected with the mesh. 877 00:54:16,060 --> 00:54:19,790 And that times step selection gives you this result here, 878 00:54:19,790 --> 00:54:22,240 0.003 seconds. 879 00:54:22,240 --> 00:54:26,210 Let's look now at the actual response that is calculated. 880 00:54:26,210 --> 00:54:29,350 Considering the horizontal displacement at the top of the 881 00:54:29,350 --> 00:54:35,780 tower, denoted as u, here, we find that in fact the 882 00:54:35,780 --> 00:54:39,160 30-element mesh and the 60-element mesh give virtually 883 00:54:39,160 --> 00:54:40,410 the same response. 884 00:54:43,060 --> 00:54:47,470 If we look at the horizontal velocity at the top off the 885 00:54:47,470 --> 00:54:53,620 tower, we find that there are some small differences between 886 00:54:53,620 --> 00:54:57,330 the 30-element and 60-element mesh. 887 00:54:57,330 --> 00:54:59,640 If we were to look at the accelerations, these 888 00:54:59,640 --> 00:55:03,480 differences would of course be even larger. 889 00:55:03,480 --> 00:55:07,570 But the purpose of this analysis was to calculate the 890 00:55:07,570 --> 00:55:12,890 displacement and velocities at the top of the tower. 891 00:55:12,890 --> 00:55:16,270 Let's look at the moment, which was also one requirement 892 00:55:16,270 --> 00:55:19,990 for our analysis, the moment right down here at the bottom 893 00:55:19,990 --> 00:55:21,810 off the tower. 894 00:55:21,810 --> 00:55:26,280 And here we see that the 30-element mesh gives the 895 00:55:26,280 --> 00:55:27,380 smooth result. 896 00:55:27,380 --> 00:55:32,950 The 60-element mesh gives the zig-zagged result. 897 00:55:32,950 --> 00:55:37,220 But overall, really, basically a good response prediction 898 00:55:37,220 --> 00:55:40,370 already with the 30-element mesh. 899 00:55:40,370 --> 00:55:46,160 If we look at the displacements of the tower, we 900 00:55:46,160 --> 00:55:48,510 have amplified these displacements here. 901 00:55:48,510 --> 00:55:52,740 The 30-element mesh at 200 millisecond looks like that. 902 00:55:52,740 --> 00:55:54,100 The 60-element mesh. 903 00:55:54,100 --> 00:55:57,380 looks like that. 904 00:55:57,380 --> 00:56:00,770 At 400 milliseconds, the 30-element mesh 905 00:56:00,770 --> 00:56:02,390 looks as shown here. 906 00:56:02,390 --> 00:56:05,061 The 60-element mesh looks as shown there. 907 00:56:07,710 --> 00:56:10,840 Regarding these analyses, we now can 908 00:56:10,840 --> 00:56:12,950 summarize some comments. 909 00:56:12,950 --> 00:56:16,820 We saw some high-frequency oscillations when we looked at 910 00:56:16,820 --> 00:56:21,360 the moment reaction response of the 60-element mesh. 911 00:56:21,360 --> 00:56:25,570 However, this high-frequency response, which actually can 912 00:56:25,570 --> 00:56:30,310 be identified to be 110 hertz, is probably inaccurate, 913 00:56:30,310 --> 00:56:34,030 because our mesh did not represent frequencies 914 00:56:34,030 --> 00:56:38,100 accurately up to that level. 915 00:56:38,100 --> 00:56:42,160 The obtained solutions for the horizontal displacement and 916 00:56:42,160 --> 00:56:46,380 the velocity at the top off the tower are really virtually 917 00:56:46,380 --> 00:56:49,210 identical for the 30- and 60- element mesh. 918 00:56:51,760 --> 00:56:54,970 Let us now look at another example. 919 00:56:54,970 --> 00:56:59,090 This is the example of a pendulum. 920 00:56:59,090 --> 00:57:03,660 A large displacement analysis of this pendulum is performed, 921 00:57:03,660 --> 00:57:08,280 and basically what we're looking at here is, if I can 922 00:57:08,280 --> 00:57:10,990 just have your attention here, is at the pendulum that 923 00:57:10,990 --> 00:57:16,460 originally is horizontal, it's pinned right here. 924 00:57:16,460 --> 00:57:21,000 And we let this pendulum go through very large motions, 925 00:57:21,000 --> 00:57:24,170 and pendle around as shown here. 926 00:57:24,170 --> 00:57:27,630 If always should come back to the horizontal position, 927 00:57:27,630 --> 00:57:30,720 unless we have damping in the system. 928 00:57:30,720 --> 00:57:34,190 Since we don't have damping in the system, the pendulum 929 00:57:34,190 --> 00:57:36,470 should just keep on going forever. 930 00:57:36,470 --> 00:57:39,690 We have a frictionless hinge right here, and this is 931 00:57:39,690 --> 00:57:42,240 pictorially represented here on the view graph, a 932 00:57:42,240 --> 00:57:44,960 frictionless hinge right there. 933 00:57:44,960 --> 00:57:52,050 The pendulum is modelled by a truss element of stiffness EA. 934 00:57:52,050 --> 00:57:55,350 The length is given here. 935 00:57:55,350 --> 00:58:00,000 The mass is concentrated right there. 936 00:58:00,000 --> 00:58:03,100 So there's no distributed mass here. 937 00:58:03,100 --> 00:58:06,460 All of the mass of the pendulum is concentrated here, 938 00:58:06,460 --> 00:58:10,400 half of it here, and half of it there. 939 00:58:10,400 --> 00:58:13,240 The initial conditions of the pendulum are listed here. 940 00:58:13,240 --> 00:58:17,760 The initial angle is 90 degrees theta, so the pendulum 941 00:58:17,760 --> 00:58:20,110 is up here, as I explained already. 942 00:58:20,110 --> 00:58:23,480 And initially, it has 0 velocity. 943 00:58:23,480 --> 00:58:26,530 And now we let the pendulum go, and as I said, it should 944 00:58:26,530 --> 00:58:31,260 spring up and down here, of course, 945 00:58:31,260 --> 00:58:32,510 under the gravity loading. 946 00:58:35,600 --> 00:58:41,510 The dynamic response that we want to calculate is obtained 947 00:58:41,510 --> 00:58:44,280 using the trapezoidal rule. 948 00:58:44,280 --> 00:58:47,770 And we are using full Newton iterations to establish the 949 00:58:47,770 --> 00:58:51,270 equilibrium during every time step. 950 00:58:51,270 --> 00:58:55,340 The convergence tolerance that we use, and that I have to 951 00:58:55,340 --> 00:59:01,480 refer you now to our earlier view graph regarding 952 00:59:01,480 --> 00:59:02,970 convergence tolerances. 953 00:59:02,970 --> 00:59:05,720 The tolerance that we're using here is ETOL equal to 10 to 954 00:59:05,720 --> 00:59:06,660 the minus seventh. 955 00:59:06,660 --> 00:59:11,686 We talked about what ETOL is in an earlier lecture. 956 00:59:11,686 --> 00:59:15,110 In our first analysis, we chose the time step to be 957 00:59:15,110 --> 00:59:17,220 equal to 0.1 seconds. 958 00:59:17,220 --> 00:59:19,200 And in that case, the following 959 00:59:19,200 --> 00:59:20,850 response is obtained here. 960 00:59:23,870 --> 00:59:28,120 We expect that solution is shown as a solid line, and 961 00:59:28,120 --> 00:59:31,890 these circles here show the calculated 962 00:59:31,890 --> 00:59:35,510 finite element solution. 963 00:59:35,510 --> 00:59:40,010 Notice 90 degrees theta here minus 90 degrees here, and 964 00:59:40,010 --> 00:59:41,505 then we return to 90 degrees. 965 00:59:44,280 --> 00:59:46,670 The solution looks quite fine, looks stable. 966 00:59:46,670 --> 00:59:49,900 And in fact, if you look at the textbook, you would find 967 00:59:49,900 --> 00:59:53,280 that there I'm giving the solution up to here. 968 00:59:53,280 --> 00:59:54,760 Everything looks good. 969 00:59:54,760 --> 00:59:59,580 However, if we plot the strain in the truss, we find that at 970 00:59:59,580 --> 01:00:06,130 larger times, the strain oscillates very much, and in 971 01:00:06,130 --> 01:00:09,890 fact, this oscillation grows. 972 01:00:09,890 --> 01:00:13,760 This clearly signifies an instability in 973 01:00:13,760 --> 01:00:16,270 the solution process. 974 01:00:16,270 --> 01:00:19,640 The instability is unchanged when we tighten our 975 01:00:19,640 --> 01:00:25,060 convergence tolerance or when we use the BFGS method instead 976 01:00:25,060 --> 01:00:27,880 of the full Newton method. 977 01:00:27,880 --> 01:00:31,390 We should recall here that the trapezoidal rule is only 978 01:00:31,390 --> 01:00:34,920 unconditionally stable in a linear analysis. 979 01:00:34,920 --> 01:00:38,490 Here we are talking about a non-linear analysis, and you 980 01:00:38,490 --> 01:00:44,320 cannot select any time step magnitude delta t and expect a 981 01:00:44,320 --> 01:00:47,420 stable solution. 982 01:00:47,420 --> 01:00:52,070 So what we have to do now is to use a smaller time step, 983 01:00:52,070 --> 01:00:57,210 and we choose delta t equal to 0.025 seconds, 984 01:00:57,210 --> 01:00:59,850 and repeat the analysis. 985 01:00:59,850 --> 01:01:03,060 Theta degrees here as a function of time. 986 01:01:03,060 --> 01:01:09,010 And once again, a very nice, smooth response solution, 987 01:01:09,010 --> 01:01:12,860 starting with 90 degrees, going to minus 90 degrees, 988 01:01:12,860 --> 01:01:16,930 coming back to 90 degrees, and so on. 989 01:01:16,930 --> 01:01:20,330 Of course, even with a larger time step, you got a smooth 990 01:01:20,330 --> 01:01:24,100 solution for theta. 991 01:01:24,100 --> 01:01:28,090 What really told us about the instability in the solution 992 01:01:28,090 --> 01:01:33,030 was our look at the strain in the truss. 993 01:01:33,030 --> 01:01:38,750 So we look now also at this strain, and in this particular 994 01:01:38,750 --> 01:01:42,070 case, with this time step, we see a smooth 995 01:01:42,070 --> 01:01:44,460 solution for the strain. 996 01:01:44,460 --> 01:01:49,990 You should please observe that the scale on the strain is 997 01:01:49,990 --> 01:01:54,020 here different from the strain scale that we had in the 998 01:01:54,020 --> 01:01:55,100 earlier view graph. 999 01:01:55,100 --> 01:01:58,230 If you refer back to the earlier view graph, you see 1000 01:01:58,230 --> 01:02:00,711 that we had a different scale there. 1001 01:02:00,711 --> 01:02:04,450 But we have a smooth solution here, in fact, a solution the 1002 01:02:04,450 --> 01:02:07,780 way we would expect it. 1003 01:02:07,780 --> 01:02:16,870 If we don't iterate, then we obtain a solution that looks 1004 01:02:16,870 --> 01:02:20,430 as shown here. 1005 01:02:20,430 --> 01:02:23,150 In other words, instead of getting this oscillatory 1006 01:02:23,150 --> 01:02:28,070 response that we would like to obtain, we obtain really a 1007 01:02:28,070 --> 01:02:31,440 highly damped solution response. 1008 01:02:31,440 --> 01:02:35,050 This, of course, shows the importance of iteration. 1009 01:02:35,050 --> 01:02:40,860 We have to iterate in the solution of a dynamic problem. 1010 01:02:40,860 --> 01:02:44,070 It's very important to integrate in order to 1011 01:02:44,070 --> 01:02:48,640 establish equilibrium for every time step, for every 1012 01:02:48,640 --> 01:02:52,050 discrete time that we're considering. 1013 01:02:52,050 --> 01:02:55,020 In this particular case, we see that we obtained a highly 1014 01:02:55,020 --> 01:02:56,160 damped solution. 1015 01:02:56,160 --> 01:02:58,700 In other cases, the solution might blow up. 1016 01:03:01,460 --> 01:03:04,670 So it is important to iterate in the nonlinear dynamic 1017 01:03:04,670 --> 01:03:06,520 solutions when you use an implicit 1018 01:03:06,520 --> 01:03:09,950 time integration scheme. 1019 01:03:09,950 --> 01:03:13,370 Although the solution obtained with the equilibrium iteration 1020 01:03:13,370 --> 01:03:15,510 is highly inaccurate, the solution 1021 01:03:15,510 --> 01:03:17,010 is, of course, stable. 1022 01:03:17,010 --> 01:03:22,020 And that is also seen when one looks at the strain plot. 1023 01:03:22,020 --> 01:03:25,050 Notice here we're plotting strains in the truss element 1024 01:03:25,050 --> 01:03:27,050 as a function of time. 1025 01:03:27,050 --> 01:03:31,740 No blow-up of that strain, no violent oscillation. 1026 01:03:31,740 --> 01:03:35,800 It is, in fact, going up here, and then becoming this rather 1027 01:03:35,800 --> 01:03:40,260 small oscillating, as shown here. 1028 01:03:40,260 --> 01:03:43,040 So we don't have an unstable solution, but a highly 1029 01:03:43,040 --> 01:03:46,140 inaccurate solution. 1030 01:03:46,140 --> 01:03:49,500 Let's look at the next example now, the analysis of a pipe 1031 01:03:49,500 --> 01:03:51,190 whip problem. 1032 01:03:51,190 --> 01:03:59,390 Here we have a long pipe, dimension given here and here, 1033 01:03:59,390 --> 01:04:01,700 that is subjected suddenly to a 1034 01:04:01,700 --> 01:04:04,570 concentrated load at its end. 1035 01:04:04,570 --> 01:04:07,990 There is also, below the pipe, a restraint. 1036 01:04:07,990 --> 01:04:12,750 That restraint is only active once this gap is closed. 1037 01:04:12,750 --> 01:04:16,310 Notice there is a 3-inch gap here between the 1038 01:04:16,310 --> 01:04:19,040 pipe and the restraint. 1039 01:04:19,040 --> 01:04:22,050 Of course, this restraint, in an engineering design, is 1040 01:04:22,050 --> 01:04:25,850 installed to prevent violent emotions of the 1041 01:04:25,850 --> 01:04:29,190 pipe at this point. 1042 01:04:29,190 --> 01:04:34,070 What our objective here is, to determine the transient 1043 01:04:34,070 --> 01:04:37,580 response of this pipe. 1044 01:04:37,580 --> 01:04:41,050 Notice that under this very large load, the pipe very 1045 01:04:41,050 --> 01:04:44,590 rapidly becomes plastic. 1046 01:04:44,590 --> 01:04:47,820 And the difficulty of this problem really lies in that 1047 01:04:47,820 --> 01:04:52,520 the pipe going plastic becomes very flexible here, but the 1048 01:04:52,520 --> 01:04:55,560 restraint suddenly hitting the pipe, or the pipe hitting the 1049 01:04:55,560 --> 01:04:59,020 restraint, rather, means that a very large difference is 1050 01:04:59,020 --> 01:05:02,730 suddenly introduced at this end. 1051 01:05:02,730 --> 01:05:06,160 The finite element model we used is shown here. 1052 01:05:06,160 --> 01:05:12,020 6 summation beam elements that of course can go plastic. 1053 01:05:12,020 --> 01:05:13,560 There's more of the elastoplasticity 1054 01:05:13,560 --> 01:05:15,130 in the system properly. 1055 01:05:15,130 --> 01:05:18,640 And a truss element from this node to that node. 1056 01:05:18,640 --> 01:05:21,370 This truss element contains a 3-inch gap. 1057 01:05:24,600 --> 01:05:30,290 The material properties for the pipes are listed here. 1058 01:05:30,290 --> 01:05:33,110 Notice the material property corresponding to the 1059 01:05:33,110 --> 01:05:34,430 elastoplasticity. 1060 01:05:34,430 --> 01:05:38,950 We assume no strain hardening and for the restraint. 1061 01:05:38,950 --> 01:05:42,090 Again, the restraint also is modeled as an 1062 01:05:42,090 --> 01:05:45,480 elastoplastic material. 1063 01:05:45,480 --> 01:05:49,280 The analysis performed, in this particular case, using 1064 01:05:49,280 --> 01:05:53,580 multiple position with two modes and direct time 1065 01:05:53,580 --> 01:05:55,180 integration. 1066 01:05:55,180 --> 01:05:58,770 We want to compare the response that we predict using 1067 01:05:58,770 --> 01:06:01,960 the multiple position procedure, which I discussed 1068 01:06:01,960 --> 01:06:05,150 earlier in this lecture, with the response that we obtain 1069 01:06:05,150 --> 01:06:08,350 using a time integration analysis. 1070 01:06:08,350 --> 01:06:10,990 For the time integration analysis, we use the 1071 01:06:10,990 --> 01:06:12,810 trapezoidal rule. 1072 01:06:12,810 --> 01:06:16,410 And in fact, even in the multiple position, we also use 1073 01:06:16,410 --> 01:06:18,870 the trapezoidal rule. 1074 01:06:18,870 --> 01:06:22,120 And for the mass representation, we have chosen 1075 01:06:22,120 --> 01:06:24,440 a consistent mass matrix. 1076 01:06:24,440 --> 01:06:27,690 In this particular case of convergence tolerance, ETOL 1077 01:06:27,690 --> 01:06:30,140 was used with this value down here. 1078 01:06:33,210 --> 01:06:37,230 If we perform the eigenvalue solution of the initial 1079 01:06:37,230 --> 01:06:43,490 system, in other words, of the elastic pipe in the initial 1080 01:06:43,490 --> 01:06:49,240 state, we obtain these natural frequencies, 8.5 1081 01:06:49,240 --> 01:06:52,960 hertz and 53 hertz. 1082 01:06:52,960 --> 01:06:54,440 These are also the mode shapes. 1083 01:06:57,230 --> 01:07:02,350 And of course, these frequencies will change as 1084 01:07:02,350 --> 01:07:05,920 time progresses in the elastoplastic response, and 1085 01:07:05,920 --> 01:07:09,390 also the mode shapes will change, particularly once 1086 01:07:09,390 --> 01:07:11,090 their restraint is hit. 1087 01:07:11,090 --> 01:07:18,510 However, what we are saying is that we believe we can capture 1088 01:07:18,510 --> 01:07:23,510 the elastoplastic response of the pipe, approximately, by 1089 01:07:23,510 --> 01:07:28,380 using this mode shape and that mode shape as generalized 1090 01:07:28,380 --> 01:07:32,300 displacement patterns to represent 1091 01:07:32,300 --> 01:07:34,150 the total pipe response. 1092 01:07:34,150 --> 01:07:39,480 Of course, there are the heavy approximations that we reduce 1093 01:07:39,480 --> 01:07:43,090 all of the pipe degrees of freedom to adjust 2 1094 01:07:43,090 --> 01:07:44,520 generalized degrees of freedom. 1095 01:07:47,470 --> 01:07:51,670 To integrate a response, we have to choose a time step. 1096 01:07:51,670 --> 01:07:58,280 And here we have delta t 1/20 of the cutoff period that we 1097 01:07:58,280 --> 01:08:00,000 want to use here. 1098 01:08:00,000 --> 01:08:04,200 Notice that we have only two modes in this problem that we 1099 01:08:04,200 --> 01:08:08,120 want to consider, with the frequencies that I 1100 01:08:08,120 --> 01:08:09,440 showed you just now. 1101 01:08:09,440 --> 01:08:12,850 And this equation gives us, as a time 1102 01:08:12,850 --> 01:08:16,140 step, 1/1000 of a second. 1103 01:08:16,140 --> 01:08:20,529 Of course, once again, this estimate, or this choice of 1104 01:08:20,529 --> 01:08:25,510 time step, is based solely on looking at the linear response 1105 01:08:25,510 --> 01:08:26,800 of the pipe. 1106 01:08:26,800 --> 01:08:32,170 In other words, assuming that the pipe remains elastic and 1107 01:08:32,170 --> 01:08:36,390 that the restraint is not there either. 1108 01:08:36,390 --> 01:08:39,330 If we look at the results obtained from these two 1109 01:08:39,330 --> 01:08:43,970 analyses, we find that the multiple position response for 1110 01:08:43,970 --> 01:08:49,160 the displacements looks as shown here, the direct time 1111 01:08:49,160 --> 01:08:53,569 integration response looks as shown here. 1112 01:08:53,569 --> 01:08:58,670 Notice that this is here the distance of the gap, so at 1113 01:08:58,670 --> 01:09:02,760 this particular time, the gap is closed and the restraint 1114 01:09:02,760 --> 01:09:05,430 comes into action. 1115 01:09:05,430 --> 01:09:08,359 We note from this graph that the multiple position solution 1116 01:09:08,359 --> 01:09:12,740 is really very close to the direct integration solution 1117 01:09:12,740 --> 01:09:16,550 for the displacement at the tip off the pipe. 1118 01:09:16,550 --> 01:09:19,850 Now let us look at the moment at the 1119 01:09:19,850 --> 01:09:21,899 built-in end of the pipe. 1120 01:09:21,899 --> 01:09:25,990 And here we have the moment response plotted as a function 1121 01:09:25,990 --> 01:09:30,770 of time for the direct integration solution, as shown 1122 01:09:30,770 --> 01:09:34,890 here, and for the multiple position 1123 01:09:34,890 --> 01:09:38,170 solution, as shown here. 1124 01:09:38,170 --> 01:09:43,990 We notice that the maximum moment occurring at the 1125 01:09:43,990 --> 01:09:49,410 built-in end in the pipe is almost the same when predicted 1126 01:09:49,410 --> 01:09:52,260 by the multiple position solution and the direct 1127 01:09:52,260 --> 01:09:54,650 integration solution. 1128 01:09:54,650 --> 01:09:57,150 There's some small difference, but the difference is 1129 01:09:57,150 --> 01:09:59,370 certainly not very large. 1130 01:09:59,370 --> 01:10:03,790 However, if we look at what happened at earlier times, we 1131 01:10:03,790 --> 01:10:07,930 see that there is a large difference between the results 1132 01:10:07,930 --> 01:10:10,820 obtained using the multiple position solution and the 1133 01:10:10,820 --> 01:10:13,300 direct integration solution. 1134 01:10:13,300 --> 01:10:16,990 The reason is that the response here at earlier times 1135 01:10:16,990 --> 01:10:21,100 is very complex, and the multiple position solution 1136 01:10:21,100 --> 01:10:26,680 cannot pick that complex response up very well. 1137 01:10:26,680 --> 01:10:31,610 However, it's one important point you should keep in mind. 1138 01:10:31,610 --> 01:10:35,610 And that is, if we had used more modes in the multiple 1139 01:10:35,610 --> 01:10:39,920 position solution, we would come closer to the response 1140 01:10:39,920 --> 01:10:44,760 obtained, predicted with the direct integration solution. 1141 01:10:44,760 --> 01:10:49,900 In fact, if the number of modes used in the multiple 1142 01:10:49,900 --> 01:10:53,250 position solution is equal to the number of degrees of 1143 01:10:53,250 --> 01:10:56,490 freedom in the finite element mesh, then the multiple 1144 01:10:56,490 --> 01:11:00,310 position solution gives you the same response as you 1145 01:11:00,310 --> 01:11:03,220 calculate with the direct integration solution. 1146 01:11:03,220 --> 01:11:06,120 I mentioned that earlier in the lecture already, and the 1147 01:11:06,120 --> 01:11:10,300 reason, of course, is that in that case, all you would have 1148 01:11:10,300 --> 01:11:13,620 done is a change of bases from the finite element 1149 01:11:13,620 --> 01:11:17,580 displacements to the mode shape bases. 1150 01:11:17,580 --> 01:11:24,200 So if you use, even with this initial stiffness matrix that 1151 01:11:24,200 --> 01:11:27,600 you have used to calculate the mode shapes, even with that 1152 01:11:27,600 --> 01:11:32,890 initial stiffness matrix, if use all the modes, as many 1153 01:11:32,890 --> 01:11:35,650 modes as there are degrees of freedom, in this multiple 1154 01:11:35,650 --> 01:11:38,620 position solution, you would get the same response 1155 01:11:38,620 --> 01:11:41,340 prediction with the multiple position solution as we're 1156 01:11:41,340 --> 01:11:46,350 seeing here, with the direct integration solution. 1157 01:11:46,350 --> 01:11:50,840 And that means we can be quite assured that as we increase 1158 01:11:50,840 --> 01:11:54,500 the number of modes in the multiple position solution, we 1159 01:11:54,500 --> 01:11:59,420 will get closer and closer to this complex response that is 1160 01:11:59,420 --> 01:12:02,590 predicted using the direct integration solution. 1161 01:12:02,590 --> 01:12:07,010 In general, one, of course, does not know the error when 1162 01:12:07,010 --> 01:12:11,340 you use the multiple position solution, and there must 1163 01:12:11,340 --> 01:12:15,490 always be a question of how large is the error, for 1164 01:12:15,490 --> 01:12:18,170 example, that we're seeing here. 1165 01:12:18,170 --> 01:12:22,510 And that can only be assessed by a more accurate solution. 1166 01:12:22,510 --> 01:12:25,360 However, I believe that's the multiple position solution is 1167 01:12:25,360 --> 01:12:29,450 of much value in the initial design analysis stages when 1168 01:12:29,450 --> 01:12:33,920 you want to obtain fairly inexpensive solutions for the 1169 01:12:33,920 --> 01:12:37,850 response, and you want to just understand how the system acts 1170 01:12:37,850 --> 01:12:41,120 that you are analyzing and designing. 1171 01:12:41,120 --> 01:12:44,700 The multiple position solution really yields relatively 1172 01:12:44,700 --> 01:12:49,400 inexpensive response solutions of the system which may be 1173 01:12:49,400 --> 01:12:52,920 very helpful in the initial design stages of that 1174 01:12:52,920 --> 01:12:54,650 structural system. 1175 01:12:54,650 --> 01:12:58,420 This completes, then, the discussion of our examples for 1176 01:12:58,420 --> 01:13:00,360 which I have prepared view graphs. 1177 01:13:00,360 --> 01:13:03,360 I'd like to now share with you some further experiences that 1178 01:13:03,360 --> 01:13:05,570 are summarized on this slide. 1179 01:13:05,570 --> 01:13:11,110 And let me go over here, where I have 1180 01:13:11,110 --> 01:13:16,560 the first slide projected. 1181 01:13:16,560 --> 01:13:19,250 I should point out, before we discuss the information on 1182 01:13:19,250 --> 01:13:22,800 these slides once more, that we will go fairly rapidly 1183 01:13:22,800 --> 01:13:26,540 through these slides, and that not all details pertaining to 1184 01:13:26,540 --> 01:13:30,110 the solutions that we are looking at here are given on 1185 01:13:30,110 --> 01:13:32,180 these slides. 1186 01:13:32,180 --> 01:13:35,230 Please refer to the papers-- 1187 01:13:35,230 --> 01:13:37,530 the references are given in the study guide-- in order to 1188 01:13:37,530 --> 01:13:40,630 find all of the details pertaining to the solutions 1189 01:13:40,630 --> 01:13:43,850 that I'm now discussing with the slides. 1190 01:13:43,850 --> 01:13:46,100 Here we consider now the analysis of a control rod 1191 01:13:46,100 --> 01:13:49,310 drive housing, which is modelled using four beam 1192 01:13:49,310 --> 01:13:52,290 elements and point masses. 1193 01:13:52,290 --> 01:13:56,740 Notice that the structure is subjected to motion up here. 1194 01:13:56,740 --> 01:13:59,790 The non-linearity in the system lies in the fact that 1195 01:13:59,790 --> 01:14:03,770 there's a gap between the tip of that beam here and the 1196 01:14:03,770 --> 01:14:08,180 spring on the right-hand side and on the left-hand side. 1197 01:14:08,180 --> 01:14:12,760 On the next slide, we see now the solutions obtained. 1198 01:14:12,760 --> 01:14:17,900 We obtain a solution using direct integration, and using 1199 01:14:17,900 --> 01:14:21,720 a multiple position analysis, including two modes. 1200 01:14:21,720 --> 01:14:24,640 These two solutions are compared with a solution that 1201 01:14:24,640 --> 01:14:29,000 was reported at some earlier time by Peterson and myself. 1202 01:14:29,000 --> 01:14:34,190 Here you see the comparison of the results obtained. 1203 01:14:34,190 --> 01:14:37,860 Notice there is a difference here between the solution that 1204 01:14:37,860 --> 01:14:40,440 we reported earlier, that Peterson and myself reported 1205 01:14:40,440 --> 01:14:44,610 earlier, and the solution that we generated now. 1206 01:14:44,610 --> 01:14:48,520 But that difference is due to the fact that Peterson and 1207 01:14:48,520 --> 01:14:53,220 myself used, at that time, a different solution scheme. 1208 01:14:53,220 --> 01:14:55,560 The response solution that we obtained with the multiple 1209 01:14:55,560 --> 01:14:59,330 position analysis is very close to the response solution 1210 01:14:59,330 --> 01:15:02,700 that we calculated using direct integration. 1211 01:15:02,700 --> 01:15:06,660 Now the next analysis that I'd like to discuss this you 1212 01:15:06,660 --> 01:15:12,230 involves the analysis of a spherical shell subjected to a 1213 01:15:12,230 --> 01:15:13,710 uniform pressure. 1214 01:15:13,710 --> 01:15:17,630 And that pressure is suddenly applied as shown down here. 1215 01:15:17,630 --> 01:15:22,260 Notice this is a step pressure constant in time. 1216 01:15:22,260 --> 01:15:26,270 The material data of the shell are given here. 1217 01:15:26,270 --> 01:15:30,230 We want to calculate the elastoplastic response. 1218 01:15:30,230 --> 01:15:34,130 And some geometric data are given here. 1219 01:15:34,130 --> 01:15:37,150 Notice the shell was idealized using 10 8-node 1220 01:15:37,150 --> 01:15:38,920 axis-symmetric elements. 1221 01:15:38,920 --> 01:15:42,290 It's a spherical shell, so we use an axis-symmetric model to 1222 01:15:42,290 --> 01:15:43,870 model the shell. 1223 01:15:43,870 --> 01:15:47,860 We use Newmark integration with delta and 1224 01:15:47,860 --> 01:15:49,560 alpha as given here. 1225 01:15:52,230 --> 01:15:57,380 We chose these values because Nagarajan and Popov had chosen 1226 01:15:57,380 --> 01:15:58,920 these same values, and we wanted to 1227 01:15:58,920 --> 01:16:01,460 compare with their solution. 1228 01:16:01,460 --> 01:16:05,050 We used 2-by-2 Gauss integration and a consistent 1229 01:16:05,050 --> 01:16:07,340 mass matrix for the problem. 1230 01:16:07,340 --> 01:16:10,840 The time step was 10 microseconds, and we used the 1231 01:16:10,840 --> 01:16:13,810 total Lagrangian formulation in the analysis. 1232 01:16:13,810 --> 01:16:17,670 The next slide now shows the response that we predicted. 1233 01:16:17,670 --> 01:16:22,190 Here we see the ADINA solution, a dashed line, 1234 01:16:22,190 --> 01:16:27,990 compared with the solution by Nagarajan and Popov given by 1235 01:16:27,990 --> 01:16:29,870 the solid line. 1236 01:16:29,870 --> 01:16:36,780 Notice we are plotting here the apex displacement of the 1237 01:16:36,780 --> 01:16:41,280 shell as a function of time. 1238 01:16:41,280 --> 01:16:44,060 Now in this particular case, we use the consistent mass 1239 01:16:44,060 --> 01:16:47,750 matrix, 2-by-2 integration, and specific solution 1240 01:16:47,750 --> 01:16:51,860 parameter, and we want to do now change these parameters in 1241 01:16:51,860 --> 01:16:55,730 order to see the effect of these parameters on the 1242 01:16:55,730 --> 01:16:57,940 solution that is predicted. 1243 01:16:57,940 --> 01:17:02,400 The next slide now shows again the apex displacement of the 1244 01:17:02,400 --> 01:17:07,420 shell as a function of time, but using once a consistent 1245 01:17:07,420 --> 01:17:14,790 mass matrix, and once using a lumped mass matrix. 1246 01:17:14,790 --> 01:17:18,190 And in this particular case, we use the trapezoidal rule, 1247 01:17:18,190 --> 01:17:21,020 which corresponds to the Newmark integration with delta 1248 01:17:21,020 --> 01:17:25,030 equals 0.5 and alpha equals 1/4. 1249 01:17:25,030 --> 01:17:28,150 Notice that there are only slide differences between the 1250 01:17:28,150 --> 01:17:31,020 lumped mass solution results and the consistent mass 1251 01:17:31,020 --> 01:17:33,940 solution results. 1252 01:17:33,940 --> 01:17:37,690 Another parameter that one might want to vary is the 1253 01:17:37,690 --> 01:17:39,090 integration order. 1254 01:17:39,090 --> 01:17:43,160 Here we have the response, using 2-by-2 integration, 1255 01:17:43,160 --> 01:17:45,480 which we saw already earlier. 1256 01:17:45,480 --> 01:17:48,150 Actually, we didn't quite see this response earlier, because 1257 01:17:48,150 --> 01:17:51,490 earlier we used different delta and alpha parameters. 1258 01:17:51,490 --> 01:17:54,920 Very close to these parameters, however. 1259 01:17:54,920 --> 01:17:59,070 Anyway, here we have the 2-by-2 integration results, 1260 01:17:59,070 --> 01:18:01,870 and then we have the 3-by-3 integration results and the 1261 01:18:01,870 --> 01:18:04,010 4-by-4 integration results here. 1262 01:18:04,010 --> 01:18:08,650 Very small differences between these solutions, but it's 1263 01:18:08,650 --> 01:18:11,940 interesting to perform these solutions, and really assure 1264 01:18:11,940 --> 01:18:14,220 oneself that there are small differences. 1265 01:18:14,220 --> 01:18:18,150 We know that shells are rather sensitive structures, and we 1266 01:18:18,150 --> 01:18:21,730 want to make sure about the modeling that we're using for 1267 01:18:21,730 --> 01:18:26,090 the representation of these types of structures. 1268 01:18:26,090 --> 01:18:32,670 The next slide now shows the analysis of a problem for 1269 01:18:32,670 --> 01:18:36,310 which we had experimental results to compare with. 1270 01:18:36,310 --> 01:18:42,130 Here we have a pipe, a very rigid pipe, and a flexible 1271 01:18:42,130 --> 01:18:45,740 pipe, a nickel pipe here. 1272 01:18:45,740 --> 01:18:52,710 And these pipes, both being straight, are totally filled 1273 01:18:52,710 --> 01:18:53,740 with water. 1274 01:18:53,740 --> 01:18:58,850 At the end off this pipe here, a pulse gun sets off a pulse, 1275 01:18:58,850 --> 01:19:03,570 and that pulse travels through the water, all the way along. 1276 01:19:03,570 --> 01:19:07,160 And of course, there's an interaction between the pipe 1277 01:19:07,160 --> 01:19:08,430 and the water. 1278 01:19:08,430 --> 01:19:12,670 In other words, it's a fluid structure interaction problem. 1279 01:19:12,670 --> 01:19:15,940 This pipe here is very rigid, so that we model it as 1280 01:19:15,940 --> 01:19:17,620 infinitely rigid. 1281 01:19:17,620 --> 01:19:20,810 You will see that just now on the next slide. 1282 01:19:20,810 --> 01:19:24,790 Notice here we have pressure gauges indicated. 1283 01:19:24,790 --> 01:19:28,200 The pressures and the strains in the pipe, in this nickel 1284 01:19:28,200 --> 01:19:31,720 pipe, the flexible pipe, were measured by the Stanford 1285 01:19:31,720 --> 01:19:34,460 Research Institute. 1286 01:19:34,460 --> 01:19:38,480 Here we have the nickel properties, we assume an 1287 01:19:38,480 --> 01:19:45,750 elastoplastic material, and the water properties. 1288 01:19:45,750 --> 01:19:50,720 On the next slide now, we see the pressure pulse as a 1289 01:19:50,720 --> 01:19:53,510 function of time. 1290 01:19:53,510 --> 01:19:58,690 Of course, this was given to us for our analysis. 1291 01:19:58,690 --> 01:20:00,650 And this slide now shows the finite element 1292 01:20:00,650 --> 01:20:01,810 model that we use. 1293 01:20:01,810 --> 01:20:04,830 We used an axis-symmetric model. 1294 01:20:04,830 --> 01:20:08,050 The center line is right here. 1295 01:20:08,050 --> 01:20:11,980 And notice that the water here in the rigid pipe was 1296 01:20:11,980 --> 01:20:14,840 represented by 81 elements. 1297 01:20:14,840 --> 01:20:18,790 Notice the rigid pipe was modelled as infinitely rigid, 1298 01:20:18,790 --> 01:20:22,440 because these are rollers on an infinitely rigid 1299 01:20:22,440 --> 01:20:25,200 surface, so to say. 1300 01:20:25,200 --> 01:20:29,630 Notice here is a nickel pipe representation, and here we 1301 01:20:29,630 --> 01:20:32,825 have water elements again, and a thin layer of 1302 01:20:32,825 --> 01:20:35,590 water elements here. 1303 01:20:35,590 --> 01:20:38,800 On the next slide now, we see the response predicted by 1304 01:20:38,800 --> 01:20:43,280 ADINA and the experimental results, the ADINA response 1305 01:20:43,280 --> 01:20:47,970 shown as a solid line and the experimental result shown as a 1306 01:20:47,970 --> 01:20:48,610 dashed line. 1307 01:20:48,610 --> 01:20:50,740 Of course, the experimental results reported by the 1308 01:20:50,740 --> 01:20:52,700 Stanford Research Institute. 1309 01:20:52,700 --> 01:20:56,470 Notice we are plotting here the pressure as a function of 1310 01:20:56,470 --> 01:21:00,560 time for one of the pressure stations. 1311 01:21:00,560 --> 01:21:03,100 And here you'll see the same information for another 1312 01:21:03,100 --> 01:21:08,280 pressure station, and for another pressure station, and 1313 01:21:08,280 --> 01:21:09,970 one more pressure station-- 1314 01:21:09,970 --> 01:21:13,970 the pressure stations that I pointed out to you earlier. 1315 01:21:13,970 --> 01:21:18,010 Notice that here finally we compare the hoop stress 1316 01:21:18,010 --> 01:21:22,390 predicted by ADINA with the experimentally observed hoop 1317 01:21:22,390 --> 01:21:25,720 stress as a function of time at one of the 1318 01:21:25,720 --> 01:21:29,450 stations within the pipe. 1319 01:21:29,450 --> 01:21:32,920 Notice that for the hoop stress here, we have four 1320 01:21:32,920 --> 01:21:35,680 experimental results corresponding to the 1321 01:21:35,680 --> 01:21:38,680 measurements that were taken at the north pole, the south 1322 01:21:38,680 --> 01:21:43,000 pole, east pole, and west pole on the pipe. 1323 01:21:43,000 --> 01:21:45,620 In the ADINA solution, of course, we use an 1324 01:21:45,620 --> 01:21:48,980 axis-symmetric model, and that axis-symmetric model only 1325 01:21:48,980 --> 01:21:53,180 gives one hoop strain. 1326 01:21:53,180 --> 01:21:57,190 And that hoop strain, as you can see here, lies quite 1327 01:21:57,190 --> 01:22:01,850 nicely within the experimental results obtained by the 1328 01:22:01,850 --> 01:22:04,020 Stanford Research Institute. 1329 01:22:04,020 --> 01:22:08,340 I should also point out that the comparison of the ADINA 1330 01:22:08,340 --> 01:22:12,740 results with the experimental results is really very good, 1331 01:22:12,740 --> 01:22:15,610 except that we don't pick up these high-frequency 1332 01:22:15,610 --> 01:22:20,410 oscillations that have been observed in the experiment. 1333 01:22:20,410 --> 01:22:24,040 Well, the ADINA finite element mesh that we used in this 1334 01:22:24,040 --> 01:22:27,370 analysis is not able, really, to represent these 1335 01:22:27,370 --> 01:22:30,330 high-frequency oscillations, and of course, that is 1336 01:22:30,330 --> 01:22:33,660 certainly one of the reasons why we could not hope to pick 1337 01:22:33,660 --> 01:22:36,690 these frequency oscillations up. 1338 01:22:36,690 --> 01:22:39,330 This, then, brings me to the last slide that I wanted to 1339 01:22:39,330 --> 01:22:41,870 share with you in this lecture, and also to the end 1340 01:22:41,870 --> 01:22:42,740 of this lecture. 1341 01:22:42,740 --> 01:22:44,090 Thank you very much for your attention.