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,872 hundreds of MIT courses, visit MIT OpenCourseWare at 7 00:00:17,872 --> 00:00:19,122 ocw.mit.edu. 8 00:00:21,230 --> 00:00:22,980 PROFESSOR: Ladies and gentlemen, welcome to this 9 00:00:22,980 --> 00:00:25,550 lecture on nonlinear finite element analysis of solids and 10 00:00:25,550 --> 00:00:26,820 structures. 11 00:00:26,820 --> 00:00:29,390 In this lecture, I like to discuss with you solution 12 00:00:29,390 --> 00:00:32,580 methods that we use to solve the finite element equations 13 00:00:32,580 --> 00:00:34,790 in nonlinear static analysis. 14 00:00:34,790 --> 00:00:38,520 In the previous lectures, we derived this set of equations, 15 00:00:38,520 --> 00:00:42,020 where tK is a tangent stiffness matrix. 16 00:00:42,020 --> 00:00:46,360 Delta Ui is a nodal point vector of incremental 17 00:00:46,360 --> 00:00:49,760 displacements corresponding to iteration i. 18 00:00:49,760 --> 00:00:54,210 t plus delta t R is the externally applied load vector 19 00:00:54,210 --> 00:00:58,260 corresponding to time t plus delta t. 20 00:00:58,260 --> 00:01:02,920 And t plus delta t F i minus 1 is equal to the nodal point 21 00:01:02,920 --> 00:01:06,530 force vector corresponding to the internal elements stresses 22 00:01:06,530 --> 00:01:09,120 at time t plus delta t, and at the end of 23 00:01:09,120 --> 00:01:12,060 iteration i minus 1. 24 00:01:12,060 --> 00:01:15,770 The displacements are updated as shown here. 25 00:01:15,770 --> 00:01:17,620 Delta Ui, of course, is calculated 26 00:01:17,620 --> 00:01:19,450 from this set of equations. 27 00:01:19,450 --> 00:01:22,850 We add delta Ui to the previous displacements that 28 00:01:22,850 --> 00:01:25,210 corresponded to iteration t plus delta -- 29 00:01:25,210 --> 00:01:28,720 to time t plus delta t and iteration i minus 1. 30 00:01:28,720 --> 00:01:31,880 It's the end of iteration i minus 1. 31 00:01:31,880 --> 00:01:37,110 And this right-hand side gives us the new displacement vector 32 00:01:37,110 --> 00:01:42,340 corresponding to iteration i, end of iteration i. 33 00:01:42,340 --> 00:01:45,680 Now this is one scheme that it's used to solve the finite 34 00:01:45,680 --> 00:01:47,770 element equations. 35 00:01:47,770 --> 00:01:50,220 We refer to this as the mortified 36 00:01:50,220 --> 00:01:51,940 Newton-Raphson method. 37 00:01:51,940 --> 00:01:53,790 However, there are other schemes as well. 38 00:01:53,790 --> 00:01:55,600 And schemes that are in certain 39 00:01:55,600 --> 00:01:58,550 problems much more effective. 40 00:01:58,550 --> 00:02:01,770 Since it is so important to use an effective method for 41 00:02:01,770 --> 00:02:04,430 the solution of the finite element equations because the 42 00:02:04,430 --> 00:02:07,470 costs can otherwise be very high, we should be familiar 43 00:02:07,470 --> 00:02:09,360 with those other techniques. 44 00:02:09,360 --> 00:02:13,980 And I'd like to share those now with you. 45 00:02:13,980 --> 00:02:16,400 Therefore, we will look in this lecture at other 46 00:02:16,400 --> 00:02:18,810 techniques to solve the finite element equations. 47 00:02:18,810 --> 00:02:22,240 And we need, of course, to discuss convergence criteria. 48 00:02:22,240 --> 00:02:23,980 It's very important to use 49 00:02:23,980 --> 00:02:25,880 appropriate convergence criteria. 50 00:02:25,880 --> 00:02:30,850 Otherwise, you iterate too long and/or on the other side, 51 00:02:30,850 --> 00:02:33,930 you may actually stop iterating before you have 52 00:02:33,930 --> 00:02:36,250 achieved proper convergence. 53 00:02:36,250 --> 00:02:41,310 So let me then go over to my view graphs and let us start 54 00:02:41,310 --> 00:02:42,580 the discussion here. 55 00:02:42,580 --> 00:02:47,316 The basic set of equations that we would like to solve 56 00:02:47,316 --> 00:02:48,530 are given here. 57 00:02:48,530 --> 00:02:52,510 Notice R is the vector of externally applied nodal point 58 00:02:52,510 --> 00:02:55,390 forces at time t plus delta t. 59 00:02:55,390 --> 00:03:00,580 And F is the vector of nodal point forces corresponding to 60 00:03:00,580 --> 00:03:04,560 the internal element stresses at time t plus delta t. 61 00:03:04,560 --> 00:03:07,870 Of course, this vector is unknown, and we want to 62 00:03:07,870 --> 00:03:13,560 iterate somehow to find it and make sure, of course, that at 63 00:03:13,560 --> 00:03:17,150 the end of that iteration, if we have the proper vector F, R 64 00:03:17,150 --> 00:03:21,670 is equal to F, or R minus F is equal to 0. 65 00:03:21,670 --> 00:03:24,570 We assume that the loading is deformation-independent. 66 00:03:24,570 --> 00:03:29,820 If say, loading is deformation-dependent, we can 67 00:03:29,820 --> 00:03:31,250 also deal with that situation. 68 00:03:31,250 --> 00:03:34,890 But we will have to add additional terms to our 69 00:03:34,890 --> 00:03:35,980 iterative scheme. 70 00:03:35,980 --> 00:03:38,650 And I will get back to that a little later. 71 00:03:38,650 --> 00:03:41,940 Notice that F in the total Lagrangian formulation the way 72 00:03:41,940 --> 00:03:47,030 we have discussed it earlier is evaluated by this product 73 00:03:47,030 --> 00:03:50,290 here integrated over the original volume for a single 74 00:03:50,290 --> 00:03:54,270 element and for the UL formulation, F is calculated 75 00:03:54,270 --> 00:03:55,310 as shown here. 76 00:03:55,310 --> 00:03:58,500 Here, of course, we integrate over the current volume, or 77 00:03:58,500 --> 00:04:02,580 the volume that we actually want to calculate. 78 00:04:02,580 --> 00:04:05,360 Notice in the iteration, of course, this volume here would 79 00:04:05,360 --> 00:04:12,020 be t plus delta t V i minus 1 as we will discuss just now. 80 00:04:12,020 --> 00:04:15,980 Because in the iteration, we always update this integral 81 00:04:15,980 --> 00:04:20,089 and that integral with the iteration counter that we are 82 00:04:20,089 --> 00:04:24,950 having on the right-hand side of the system of equations. 83 00:04:24,950 --> 00:04:31,170 The methods that we use are based on the Newton-Raphson 84 00:04:31,170 --> 00:04:38,020 method, which is really used very abundantly to find roots 85 00:04:38,020 --> 00:04:40,260 of an equation. 86 00:04:40,260 --> 00:04:41,960 A small historical note: 87 00:04:41,960 --> 00:04:46,360 Newton gave a first version of the method in 1669. 88 00:04:46,360 --> 00:04:49,560 Raphson then generalized and simplified the method 89 00:04:49,560 --> 00:04:52,860 actually, in 1690. 90 00:04:52,860 --> 00:04:55,940 Both mathematicians used the same concept. 91 00:04:55,940 --> 00:04:59,290 And both algorithms gave really, the same results. 92 00:04:59,290 --> 00:05:03,090 But it is very appropriate to refer to the methods that 93 00:05:03,090 --> 00:05:06,490 we're using as the Newton-Raphson method because 94 00:05:06,490 --> 00:05:11,010 Raphson really contributed quite a bit to that method. 95 00:05:11,010 --> 00:05:14,110 So let us now consider a single 96 00:05:14,110 --> 00:05:16,220 Newton-Raphson iteration. 97 00:05:16,220 --> 00:05:19,030 The way you may have actually encountered it already earlier 98 00:05:19,030 --> 00:05:23,220 in your studies of solution methods for the roots of an 99 00:05:23,220 --> 00:05:24,730 algebraic equation. 100 00:05:24,730 --> 00:05:27,070 Basically, what you're doing is this. 101 00:05:27,070 --> 00:05:32,510 You're saying if you have xi minus 1, you calculate 102 00:05:32,510 --> 00:05:34,660 f at xi minus 1. 103 00:05:34,660 --> 00:05:39,900 You divide this value by the f prime at xi minus 1. 104 00:05:39,900 --> 00:05:43,910 And then this right-hand side gives you better value, a 105 00:05:43,910 --> 00:05:47,340 better approximation to the root of the equation. 106 00:05:47,340 --> 00:05:50,540 Once you have xi, you repeat the process and 107 00:05:50,540 --> 00:05:53,090 you get xi plus 1. 108 00:05:53,090 --> 00:05:58,630 You keep on repeating the process until, basically f xi 109 00:05:58,630 --> 00:06:01,070 here is close to 0. 110 00:06:01,070 --> 00:06:03,980 Because then you have a root. 111 00:06:03,980 --> 00:06:17,470 Well, if we use that Newton-Raphson formula, it is 112 00:06:17,470 --> 00:06:21,010 quite interesting to see how it has been derived. 113 00:06:21,010 --> 00:06:25,130 We can write for any point xi, a neighboring point xi minus 114 00:06:25,130 --> 00:06:29,970 1, directly this equation here by a Taylor series expansion. 115 00:06:29,970 --> 00:06:33,450 And if we neglect the higher order terms, we get that f of 116 00:06:33,450 --> 00:06:37,180 xi is approximately equal to this relationship on the 117 00:06:37,180 --> 00:06:39,220 right-hand side. 118 00:06:39,220 --> 00:06:44,335 Well, since f of xi is supposed to be 0, because we 119 00:06:44,335 --> 00:06:47,060 are looking for a 0 of the equation, we 120 00:06:47,060 --> 00:06:48,430 set it equal to 0. 121 00:06:48,430 --> 00:06:50,860 And here you see directly the formula that 122 00:06:50,860 --> 00:06:52,120 I showed you earlier. 123 00:06:52,120 --> 00:06:55,360 So this is a very quick derivation of the 124 00:06:55,360 --> 00:06:59,020 Newton-Raphson procedure. 125 00:06:59,020 --> 00:07:03,040 Let us look at a mathematical example to see how the method 126 00:07:03,040 --> 00:07:07,370 works, and to get some insight into the technique. 127 00:07:07,370 --> 00:07:10,960 Let us choose a very simple example, not much to do with 128 00:07:10,960 --> 00:07:12,680 finite element analysis. 129 00:07:12,680 --> 00:07:16,280 Let us say that f of x equals sine x and that our starting 130 00:07:16,280 --> 00:07:19,080 value is equal to 2 for the iteration. 131 00:07:19,080 --> 00:07:21,770 You always have to choose, of course, a starting value for 132 00:07:21,770 --> 00:07:23,030 the iteration. 133 00:07:23,030 --> 00:07:28,750 Then, in this column you see the values that are calculated 134 00:07:28,750 --> 00:07:30,830 in the successive iterations. 135 00:07:30,830 --> 00:07:34,020 And we also show here the error. 136 00:07:34,020 --> 00:07:37,700 It is interesting to observe that when you are close to the 137 00:07:37,700 --> 00:07:40,210 root, you have quadratic convergence. 138 00:07:40,210 --> 00:07:45,410 Meaning that the error here, epsilon becomes an error. 139 00:07:45,410 --> 00:07:48,570 Epsilon squared here. 140 00:07:48,570 --> 00:07:51,880 This is summarized on this view graph. 141 00:07:51,880 --> 00:07:58,790 Mathematically, if we have that the error, Ei minus 1, is 142 00:07:58,790 --> 00:08:00,220 given by this equation here. 143 00:08:00,220 --> 00:08:02,200 Of course, there's an approximate sign here. 144 00:08:02,200 --> 00:08:04,340 Then the error in the next iteration 145 00:08:04,340 --> 00:08:07,906 will be of this order. 146 00:08:07,906 --> 00:08:12,590 So the convergence is really very rapid once we are close 147 00:08:12,590 --> 00:08:15,660 to the root. 148 00:08:15,660 --> 00:08:19,170 Here we can see, however, that if we are not 149 00:08:19,170 --> 00:08:21,300 close to the root-- 150 00:08:21,300 --> 00:08:26,690 in other words, if we are too far from the root, and we pick 151 00:08:26,690 --> 00:08:31,390 a bad value, then we don't get to the desired 152 00:08:31,390 --> 00:08:34,080 root, which is pi. 153 00:08:34,080 --> 00:08:38,919 But we get a value that is far away from pi. 154 00:08:38,919 --> 00:08:42,159 It is also root, but certainly not the root that we are 155 00:08:42,159 --> 00:08:43,635 interested in. 156 00:08:43,635 --> 00:08:47,550 So Newton-Raphson iteration is not always convergent 157 00:08:47,550 --> 00:08:51,420 directly, converging directly to the result that we would 158 00:08:51,420 --> 00:08:52,960 like to obtain. 159 00:08:52,960 --> 00:08:57,090 There has to be some care when we use that iterative scheme. 160 00:08:57,090 --> 00:09:02,720 And in fact, we have to be close enough just to basically 161 00:09:02,720 --> 00:09:04,190 say what has to happen. 162 00:09:04,190 --> 00:09:07,110 We have to be close enough to the root in order to have 163 00:09:07,110 --> 00:09:09,930 convergence to the root that we're looking for. 164 00:09:09,930 --> 00:09:12,750 And ultimately, if we do converge to that root, we will 165 00:09:12,750 --> 00:09:15,060 get quadratic convergence to the root. 166 00:09:15,060 --> 00:09:17,420 And that is, of course, very desirable. 167 00:09:17,420 --> 00:09:21,040 But the starting value is critical, as you can see here 168 00:09:21,040 --> 00:09:22,400 from this view graph. 169 00:09:22,400 --> 00:09:27,050 Let us look pictorially at the solution process. 170 00:09:27,050 --> 00:09:30,020 Here we have our sine function. 171 00:09:30,020 --> 00:09:31,770 And we are looking for this root. 172 00:09:34,280 --> 00:09:41,170 In iteration 1, we basically set up a tangent to that sine 173 00:09:41,170 --> 00:09:44,380 function at the starting value, x0. 174 00:09:44,380 --> 00:09:52,140 And we calculate x1 using this tangent as another estimate to 175 00:09:52,140 --> 00:09:55,820 the value that we're interested in. 176 00:09:55,820 --> 00:10:03,540 Now in the iteration 2, we lay a tangent to the f function. 177 00:10:03,540 --> 00:10:06,520 Which, of course, is sine x, at the new value 178 00:10:06,520 --> 00:10:09,050 corresponding to x1. 179 00:10:09,050 --> 00:10:12,520 And with this tangent, we get to this value. 180 00:10:12,520 --> 00:10:15,150 And you can see already that we got closer 181 00:10:15,150 --> 00:10:16,220 to the desired value. 182 00:10:16,220 --> 00:10:21,070 We were this far away originally. 183 00:10:21,070 --> 00:10:23,460 Then that far away. 184 00:10:23,460 --> 00:10:26,910 And now we are only that far away. 185 00:10:26,910 --> 00:10:29,950 Notice our sine function has this value here. 186 00:10:29,950 --> 00:10:32,680 And certainly, it's not 0 yet. 187 00:10:32,680 --> 00:10:37,540 So we lay another tangent at that point. 188 00:10:37,540 --> 00:10:42,120 And that's being done on this view graph here. 189 00:10:42,120 --> 00:10:45,890 Another tangent, and now we get much closer 190 00:10:45,890 --> 00:10:46,990 to the desired value. 191 00:10:46,990 --> 00:10:50,000 Right here, x3 is already quite close. 192 00:10:50,000 --> 00:10:54,340 But not close enough as shown by this value. 193 00:10:54,340 --> 00:10:59,380 It's not yet close enough to 0, this length here. 194 00:10:59,380 --> 00:11:03,830 And we lay another tangent to that point now. 195 00:11:03,830 --> 00:11:08,500 And that's shown on this view graph now. 196 00:11:08,500 --> 00:11:12,110 And we can see with the blue line being the tangent at this 197 00:11:12,110 --> 00:11:18,470 point, we get x4 very close to the 0 value that we are 198 00:11:18,470 --> 00:11:20,770 actually interested in. 199 00:11:20,770 --> 00:11:23,840 Certainly on this view graph, you can't see any difference 200 00:11:23,840 --> 00:11:28,290 anymore between x4 and the desired value. 201 00:11:28,290 --> 00:11:31,220 So this is basically the process that we are following 202 00:11:31,220 --> 00:11:34,510 through when we go from iteration 1 to iteration 2 to 203 00:11:34,510 --> 00:11:38,540 iteration 3 and to iteration 4. 204 00:11:38,540 --> 00:11:41,580 Beautiful convergence in this particular example. 205 00:11:41,580 --> 00:11:47,290 And the reason being that x0 is close enough to the desired 206 00:11:47,290 --> 00:11:50,730 value, the desired root. 207 00:11:50,730 --> 00:11:52,650 Which, itself, of course, is close to x4. 208 00:11:52,650 --> 00:11:53,840 Four 209 00:11:53,840 --> 00:12:01,420 Well, if however, we were to take a value that is not close 210 00:12:01,420 --> 00:12:07,420 enough and then we, as I showed earlier, we do not 211 00:12:07,420 --> 00:12:11,050 converge to the actual root. 212 00:12:11,050 --> 00:12:13,020 And this is shown here. 213 00:12:13,020 --> 00:12:18,500 For example, if we take a tangent at a point where f of 214 00:12:18,500 --> 00:12:24,920 prime is almost 0, then of course, this tangent throws us 215 00:12:24,920 --> 00:12:27,080 far away from this root. 216 00:12:27,080 --> 00:12:31,220 And we will converge to some other root. 217 00:12:31,220 --> 00:12:34,120 Of course, we don't want to use an exact 0 here because we 218 00:12:34,120 --> 00:12:37,310 have f of x divided by f prime of x. 219 00:12:37,310 --> 00:12:40,080 So with an exact 0, we could not [UNINTELLIGIBLE]. 220 00:12:40,080 --> 00:12:45,330 But a tangent that is close to 0, an f prime value that is 221 00:12:45,330 --> 00:12:49,220 close to 0, would throw us far away from the desired root. 222 00:12:49,220 --> 00:12:53,180 And that shows some of the difficulties using-- 223 00:12:53,180 --> 00:12:59,040 in fact, almost all iterative schemes that we may not get 224 00:12:59,040 --> 00:13:00,750 fast enough to our root. 225 00:13:00,750 --> 00:13:03,520 And in fact, we may never get to the root that we are 226 00:13:03,520 --> 00:13:05,010 looking for. 227 00:13:05,010 --> 00:13:09,380 Well, let us then look at the Newton-Raphson iteration for 228 00:13:09,380 --> 00:13:12,130 multiple degrees of freedom. 229 00:13:12,130 --> 00:13:15,350 Here we have our R minus F that we want to 230 00:13:15,350 --> 00:13:17,700 be setting to 0. 231 00:13:17,700 --> 00:13:21,130 Or we want to rather, solve for the displacements t plus 232 00:13:21,130 --> 00:13:27,400 delta t U, which are the solution because this t plus 233 00:13:27,400 --> 00:13:32,060 delta t u displacements would give us this force vector. 234 00:13:32,060 --> 00:13:36,610 And this force vector equilibrates the R vector. 235 00:13:36,610 --> 00:13:41,580 Meaning f of U, this little f of U, is equal to 0. 236 00:13:41,580 --> 00:13:44,990 And that, of course, is what we want to achieve. 237 00:13:44,990 --> 00:13:47,130 Notice that we have here, of course, now 238 00:13:47,130 --> 00:13:48,720 n degrees of freedom. 239 00:13:48,720 --> 00:13:51,660 Therefore, n equations to solve. 240 00:13:51,660 --> 00:13:55,290 Previously, we only looked at one equation. 241 00:13:55,290 --> 00:14:01,070 If we apply the same principle as before to this little f, 242 00:14:01,070 --> 00:14:06,840 namely we expand f as a Taylor series about t plus delta t U 243 00:14:06,840 --> 00:14:15,000 i minus 1, then we get this equation on 244 00:14:15,000 --> 00:14:17,770 the right-hand side. 245 00:14:17,770 --> 00:14:20,300 And of course, on the left-hand side, we have f of t 246 00:14:20,300 --> 00:14:22,120 plus delta t Ui. 247 00:14:22,120 --> 00:14:26,660 Notice they are higher terms, which we will neglect just as 248 00:14:26,660 --> 00:14:28,630 we have done earlier for the single equation. 249 00:14:28,630 --> 00:14:31,650 Because we are looking for Taylor series expansion about 250 00:14:31,650 --> 00:14:34,790 t plus delta t U i minus 1. 251 00:14:34,790 --> 00:14:41,560 If we neglect this part here, we obtain directly this 252 00:14:41,560 --> 00:14:44,010 equation here. 253 00:14:44,010 --> 00:14:47,020 0 on the left-hand side, because once again, we are 254 00:14:47,020 --> 00:14:54,440 looking for the displacements values for which f is 0. 255 00:14:54,440 --> 00:14:56,770 So therefore, we set this deliberately equal to 0. 256 00:14:56,770 --> 00:15:00,750 And on the right-hand side, we have f t plus delta t U i 257 00:15:00,750 --> 00:15:04,540 minus 1 plus partial f with respect to U times the delta 258 00:15:04,540 --> 00:15:07,190 U. This delta U is unknown, of course, 259 00:15:07,190 --> 00:15:08,420 that we want to calculate. 260 00:15:08,420 --> 00:15:12,910 And this delta U is nothing else then the difference 261 00:15:12,910 --> 00:15:16,610 between the U value in iteration i and the U value in 262 00:15:16,610 --> 00:15:19,140 iteration i minus 1. 263 00:15:19,140 --> 00:15:22,670 Notice that this f is evaluated at t plus 264 00:15:22,670 --> 00:15:24,350 delta t U i minus 1. 265 00:15:27,510 --> 00:15:29,530 Let us look now at this equation. 266 00:15:29,530 --> 00:15:36,620 And if we write it a little bit more out, we find that on 267 00:15:36,620 --> 00:15:40,300 the left-hand side, we have, of course, a vector of 0's. 268 00:15:40,300 --> 00:15:43,990 On the right-hand side, we have a vector of f i 269 00:15:43,990 --> 00:15:47,070 components, f1 to fn. 270 00:15:47,070 --> 00:15:52,600 All of which are evaluated at t plus delta t U i minus 1. 271 00:15:52,600 --> 00:15:55,940 And then we have here a matrix, which is a square 272 00:15:55,940 --> 00:16:02,560 matrix in which the individual elements are partials of f i 273 00:16:02,560 --> 00:16:04,240 with respect to Uj. 274 00:16:04,240 --> 00:16:09,020 For example, f1, partial f1 with respect to U1 here. 275 00:16:09,020 --> 00:16:13,990 Partial fn with respect to U n here, et cetera. 276 00:16:13,990 --> 00:16:18,130 This matrix here will give us a tangent stiffness matrix. 277 00:16:18,130 --> 00:16:21,120 Notice this matrix is multiplied by the vector of 278 00:16:21,120 --> 00:16:24,090 incremental nodal point displacements corresponding to 279 00:16:24,090 --> 00:16:25,340 iteration i. 280 00:16:28,840 --> 00:16:37,040 If we remember that f, the little f, at t plus delta t U 281 00:16:37,040 --> 00:16:41,450 i minus 1 is equal to R minus 2 plus delta t F i minus 1. 282 00:16:41,450 --> 00:16:44,980 Capital F here, little f there, remember? 283 00:16:44,980 --> 00:16:50,230 Then, the partial of little f with respect to U at t plus 284 00:16:50,230 --> 00:16:53,733 delta t U i minus 1 is given via this equation. 285 00:16:53,733 --> 00:16:58,720 Now, please recognize that partial of R with respect to U 286 00:16:58,720 --> 00:17:03,870 is equal to 0 if the loads are deformation-independent. 287 00:17:03,870 --> 00:17:06,770 Of course, if the loads depend on the deformations, in other 288 00:17:06,770 --> 00:17:10,670 words, if the loads depend on the U displacements, then this 289 00:17:10,670 --> 00:17:15,010 would not be 0, and we would have to actually put a term in 290 00:17:15,010 --> 00:17:17,979 here, carry that term along in the solution of 291 00:17:17,979 --> 00:17:20,069 the nonlinear equations. 292 00:17:20,069 --> 00:17:24,050 But for the moment, we neglect this. 293 00:17:24,050 --> 00:17:27,210 We assume that the loads are deformation-independent. 294 00:17:27,210 --> 00:17:30,230 And then, this term is equal to 0. 295 00:17:30,230 --> 00:17:33,770 This here is now giving us the tangent stiffness matrix. 296 00:17:33,770 --> 00:17:37,430 And it is written down here. 297 00:17:37,430 --> 00:17:39,870 The tangent stiffness matrix is nothing else than the 298 00:17:39,870 --> 00:17:45,440 partials of capital F with respect to the U's, to the 299 00:17:45,440 --> 00:17:46,690 displacements. 300 00:17:48,520 --> 00:17:51,990 The final result then is given on this view graph. 301 00:17:51,990 --> 00:17:55,540 We substitute all the information that I shared with 302 00:17:55,540 --> 00:18:00,790 you into the Taylor series expansion around t plus delta 303 00:18:00,790 --> 00:18:02,470 t U i minus 1. 304 00:18:02,470 --> 00:18:05,110 We get directly this equation here. 305 00:18:05,110 --> 00:18:10,300 Notice tangent stiffness matrix corresponding to time t 306 00:18:10,300 --> 00:18:14,360 plus delta t and iteration i minus 1. 307 00:18:14,360 --> 00:18:20,580 Delta U i corresponding to iteration i we have to be very 308 00:18:20,580 --> 00:18:23,330 clear about this, that this is an increment in displacement 309 00:18:23,330 --> 00:18:28,460 from time t plus delta t iteration i minus 1 to 310 00:18:28,460 --> 00:18:30,220 iteration i. 311 00:18:30,220 --> 00:18:33,380 And of course, the nodal point force vector. 312 00:18:33,380 --> 00:18:36,020 Externally applied nodal point forces go in here. 313 00:18:36,020 --> 00:18:38,470 And this is the nodal point force vector corresponding to 314 00:18:38,470 --> 00:18:43,090 the internal element stresses at time t plus delta t and at 315 00:18:43,090 --> 00:18:45,090 the end of iteration i minus 1. 316 00:18:45,090 --> 00:18:49,100 Evaluated differently in the total Lagrangian formulation 317 00:18:49,100 --> 00:18:52,480 from the way we're evaluating it in the updated Lagrangian 318 00:18:52,480 --> 00:18:53,070 formulation. 319 00:18:53,070 --> 00:18:55,690 We talked about this vector abundantly in 320 00:18:55,690 --> 00:18:57,330 the previous lectures. 321 00:18:57,330 --> 00:19:00,510 We talked about this vector, of course, also in the 322 00:19:00,510 --> 00:19:01,630 previous lectures. 323 00:19:01,630 --> 00:19:05,970 We had in the previous lectures here at tK, a 324 00:19:05,970 --> 00:19:10,260 constant tangent stiffness matrix, which was set up at 325 00:19:10,260 --> 00:19:14,210 the beginning of this whole interactive process, and was 326 00:19:14,210 --> 00:19:16,030 never updated. 327 00:19:16,030 --> 00:19:19,850 Well, here we now update that stiffness matrix. 328 00:19:19,850 --> 00:19:23,650 And I think if you look through the information that I 329 00:19:23,650 --> 00:19:26,080 just discussed with you, you recognize why 330 00:19:26,080 --> 00:19:27,360 we're updating it. 331 00:19:27,360 --> 00:19:29,840 We are doing so because we are always starting with the new 332 00:19:29,840 --> 00:19:37,120 Taylor series expansion about the point i minus 1. 333 00:19:37,120 --> 00:19:41,990 This set of equations if of course solved for delta Ui. 334 00:19:41,990 --> 00:19:44,070 We add delta Ui to what we had already in terms of 335 00:19:44,070 --> 00:19:48,020 displacements to get our new estimate for the nodal point 336 00:19:48,020 --> 00:19:50,160 displacements. 337 00:19:50,160 --> 00:19:54,450 It is important to realize that the K matrix, which we 338 00:19:54,450 --> 00:19:57,440 are using in the solution process, is symmetric. 339 00:19:57,440 --> 00:20:01,190 Because first of all, we use symmetric stress and strain 340 00:20:01,190 --> 00:20:03,900 measures in our governing equation. 341 00:20:03,900 --> 00:20:06,990 Remember that when we apply the principle of virtual work, 342 00:20:06,990 --> 00:20:10,760 we use Cauchy stresses and infinitesimally 343 00:20:10,760 --> 00:20:12,250 small virtual strains. 344 00:20:12,250 --> 00:20:16,850 Now notice that both of these tenses are symmetric tenses. 345 00:20:16,850 --> 00:20:19,980 Of course, these are the tenses we're using in the UL, 346 00:20:19,980 --> 00:20:22,280 in the updated Lagrangian formulation. 347 00:20:22,280 --> 00:20:25,180 In the total Lagrangian formulation, we use the second 348 00:20:25,180 --> 00:20:28,340 Piola-Kirchhoff stress in Green-Lagrange strain. 349 00:20:28,340 --> 00:20:32,350 Both of these measures are again, symmetric measures. 350 00:20:32,350 --> 00:20:34,650 Both tenses are symmetric tenses. 351 00:20:34,650 --> 00:20:38,880 If we had introduced for a formulation non-symmetric 352 00:20:38,880 --> 00:20:42,210 tenses, for example a non-symmetric stress tenser, 353 00:20:42,210 --> 00:20:44,960 and of course, an energy conjugate non-symmetric strain 354 00:20:44,960 --> 00:20:48,780 tenser, we would have obtained a non-symmetric K matrix, 355 00:20:48,780 --> 00:20:53,500 which would be much more difficult to deal with. 356 00:20:53,500 --> 00:20:57,680 Much less cost effective in the solution process. 357 00:20:57,680 --> 00:21:01,470 Also, please recognize that we interpolated the real 358 00:21:01,470 --> 00:21:05,810 displacements and the virtual displacements with exactly the 359 00:21:05,810 --> 00:21:07,820 same functions. 360 00:21:07,820 --> 00:21:11,020 Whereas this point here is a continuum mechanics point, 361 00:21:11,020 --> 00:21:13,680 this is really a finite element point. 362 00:21:13,680 --> 00:21:17,590 If we had used different interpolation functions for 363 00:21:17,590 --> 00:21:21,930 the real displacements, then for the virtual displacements, 364 00:21:21,930 --> 00:21:26,300 or vice-versa, then we would have not necessarily obtained 365 00:21:26,300 --> 00:21:27,860 a symmetric K matrix. 366 00:21:27,860 --> 00:21:32,110 Of course, the most natural procedure is to use the same 367 00:21:32,110 --> 00:21:35,460 kind of interpolation functions for the virtual 368 00:21:35,460 --> 00:21:37,940 displacements as we used for the real displacements, and 369 00:21:37,940 --> 00:21:39,340 that's what we did. 370 00:21:39,340 --> 00:21:42,460 Finally, we also assumed that the loading was 371 00:21:42,460 --> 00:21:43,710 deformation-independent. 372 00:21:45,880 --> 00:21:49,000 If we have deformation-dependent loading, 373 00:21:49,000 --> 00:21:52,680 then if you go more back to the earlier view graph, then 374 00:21:52,680 --> 00:21:57,020 of course, this right-hand side vector R here would 375 00:21:57,020 --> 00:21:59,750 depend on the displacements. 376 00:21:59,750 --> 00:22:02,630 And there are basically two different approaches 377 00:22:02,630 --> 00:22:04,310 that one can take. 378 00:22:04,310 --> 00:22:08,480 In the first approach, one simply updates this vector 379 00:22:08,480 --> 00:22:11,770 with the iteration, or with taking into account the 380 00:22:11,770 --> 00:22:15,910 current or last calculated volume and surface areas for 381 00:22:15,910 --> 00:22:16,920 the elements. 382 00:22:16,920 --> 00:22:21,430 So we would simply put here an i minus 1 up there. 383 00:22:21,430 --> 00:22:24,660 The left-hand side matrix, we leave unchanged. 384 00:22:24,660 --> 00:22:27,630 And if that converges fast, certainly it's a very 385 00:22:27,630 --> 00:22:29,540 effective approach to use. 386 00:22:29,540 --> 00:22:32,950 We use that approach abundantly. 387 00:22:32,950 --> 00:22:36,490 However, another approach is to actually take the 388 00:22:36,490 --> 00:22:39,980 differentiation of R with respect to the U's, the way we 389 00:22:39,980 --> 00:22:41,890 have been looking at it earlier on an earlier view 390 00:22:41,890 --> 00:22:45,260 graph, and then you get some components that you 391 00:22:45,260 --> 00:22:47,110 add to the K matrix. 392 00:22:47,110 --> 00:22:51,200 And that K matrix may then out to be non-symmetric. 393 00:22:51,200 --> 00:22:54,400 However, if the loading is deformation-independent, then 394 00:22:54,400 --> 00:22:57,910 the differentiation of R with respect to the displacements U 395 00:22:57,910 --> 00:23:01,770 is equal to 0, and there's no component from that 396 00:23:01,770 --> 00:23:04,250 differentiation coming into the K matrix. 397 00:23:04,250 --> 00:23:07,550 And our K, of course, provided these are also 398 00:23:07,550 --> 00:23:10,850 satisfied, is symmetric. 399 00:23:10,850 --> 00:23:13,840 The iterative scheme that we just discussed is referred to 400 00:23:13,840 --> 00:23:16,780 as the full Newton-Raphson method. 401 00:23:16,780 --> 00:23:19,630 Full because we are setting up a new K matrix at the 402 00:23:19,630 --> 00:23:22,340 beginning of each iteration. 403 00:23:22,340 --> 00:23:26,120 The full Newton-Raphson method shows mathematically quadratic 404 00:23:26,120 --> 00:23:28,650 convergence the way we discussed it a bit earlier in 405 00:23:28,650 --> 00:23:31,640 the lecture. 406 00:23:31,640 --> 00:23:36,400 And that, of course, is always the case provided you are 407 00:23:36,400 --> 00:23:38,050 close enough to the root. 408 00:23:38,050 --> 00:23:41,550 This quadratic convergence only holds provided you are 409 00:23:41,550 --> 00:23:45,780 close enough to the root when you solve your equations. 410 00:23:45,780 --> 00:23:49,420 In finite element analysis, it is also important to recognize 411 00:23:49,420 --> 00:23:52,410 that a number of requirements must be fulfilled. 412 00:23:52,410 --> 00:23:57,410 For example, in elasto-plastic analysis, the stresses and 413 00:23:57,410 --> 00:23:59,020 strains must be properly-- 414 00:23:59,020 --> 00:24:03,520 plastic strains must be properly updated. 415 00:24:03,520 --> 00:24:07,050 And similarly, the rotations in a shell analysis must be 416 00:24:07,050 --> 00:24:07,870 properly updated. 417 00:24:07,870 --> 00:24:11,380 So it is not necessarily the case that you automatically 418 00:24:11,380 --> 00:24:14,780 get quadratic convergence when you do finite element 419 00:24:14,780 --> 00:24:18,250 announces with a full Newton-Raphson method. 420 00:24:18,250 --> 00:24:22,120 It is very important to also, on the level of updating the 421 00:24:22,120 --> 00:24:27,800 stresses and the rotations, really do things properly in 422 00:24:27,800 --> 00:24:33,040 quotes in order to obtain the full quadratic convergence of 423 00:24:33,040 --> 00:24:36,900 the Newton-Raphson method. 424 00:24:36,900 --> 00:24:39,580 We can depict the iteration process in 425 00:24:39,580 --> 00:24:41,350 two equivalent ways. 426 00:24:41,350 --> 00:24:45,300 The first way is shown here, the left. 427 00:24:45,300 --> 00:24:48,520 And it's really the way we've been discussing just now the 428 00:24:48,520 --> 00:24:54,280 solution of f equals R minus capital F. Where we want to 429 00:24:54,280 --> 00:24:58,150 solve for the root, the 0 of this equation. 430 00:24:58,150 --> 00:25:02,980 Notice here we have in red the little f depicted. 431 00:25:02,980 --> 00:25:07,280 Notice that at this point here, we have t plus delta t 432 00:25:07,280 --> 00:25:10,650 capital F i minus 1. 433 00:25:10,650 --> 00:25:14,180 And t plus delta t R is this value. 434 00:25:14,180 --> 00:25:18,130 Now, as we get closer to the root, which is of course, the 435 00:25:18,130 --> 00:25:23,800 point where the red curve crosses this U-axis, as we get 436 00:25:23,800 --> 00:25:28,620 closer to the root, this capital F will get closer and 437 00:25:28,620 --> 00:25:33,140 closer to the R. And that is, of course, when the little f 438 00:25:33,140 --> 00:25:34,310 is equal to 0. 439 00:25:34,310 --> 00:25:36,080 That's all we are showing here. 440 00:25:36,080 --> 00:25:41,830 Notice that we are setting up a slope f prime at the point 441 00:25:41,830 --> 00:25:44,330 corresponding to the i minus first iteration. 442 00:25:44,330 --> 00:25:49,180 And that slope brings us into this point, which will be the 443 00:25:49,180 --> 00:25:51,240 next point for-- 444 00:25:51,240 --> 00:25:55,330 or the next starting point for the next iteration. 445 00:25:55,330 --> 00:25:57,930 I should say, the point of starting 446 00:25:57,930 --> 00:25:59,600 with the next iteration. 447 00:25:59,600 --> 00:26:02,890 Now, this is one way of looking at 448 00:26:02,890 --> 00:26:04,420 the iterative process. 449 00:26:04,420 --> 00:26:06,960 We can also look at it as shown here. 450 00:26:06,960 --> 00:26:10,830 Notice we have here the R plotted now horizontally. 451 00:26:10,830 --> 00:26:13,240 The displacement vertically. 452 00:26:13,240 --> 00:26:19,680 Notice the R at a particular time t plus delta t is shown 453 00:26:19,680 --> 00:26:22,160 by this dashed line. 454 00:26:22,160 --> 00:26:27,110 And we really want to find this particular point here and 455 00:26:27,110 --> 00:26:30,970 the corresponding U displacement of course. 456 00:26:30,970 --> 00:26:34,380 At this point, the little f is 0 and the corresponding U 457 00:26:34,380 --> 00:26:36,250 displacement is down here. 458 00:26:36,250 --> 00:26:40,240 Notice that at the iteration i minus 1, it's the end of 459 00:26:40,240 --> 00:26:44,570 iteration i minus 1, we have obtained this point here. 460 00:26:44,570 --> 00:26:47,850 The U displacement corresponding to that point is 461 00:26:47,850 --> 00:26:52,180 obtained by projecting down on the displacement axes. 462 00:26:52,180 --> 00:26:57,800 And this slope here, the blue line, gives us a tangent 463 00:26:57,800 --> 00:27:00,350 stiffness matrix slope. 464 00:27:00,350 --> 00:27:02,700 These are two equivalent ways. 465 00:27:02,700 --> 00:27:06,070 This is here more like a force deflection curve, and we will 466 00:27:06,070 --> 00:27:10,270 use this representation now abundantly. 467 00:27:10,270 --> 00:27:13,500 But keep in mind, it's really the same thing if you look at 468 00:27:13,500 --> 00:27:15,300 it closely. 469 00:27:15,300 --> 00:27:17,890 Well, modifications to the iterative 470 00:27:17,890 --> 00:27:21,410 scheme are the following. 471 00:27:21,410 --> 00:27:25,100 If we leave the stiffness matrix constant throughout a 472 00:27:25,100 --> 00:27:26,670 complete solution process-- 473 00:27:26,670 --> 00:27:30,070 in other words, tau is equal to 0 here. 474 00:27:30,070 --> 00:27:31,390 It's never updated. 475 00:27:31,390 --> 00:27:34,570 We talk about the initial stress method. 476 00:27:34,570 --> 00:27:39,610 If tau is equal to t, where t, of course, corresponds to a 477 00:27:39,610 --> 00:27:44,670 particular solution step and meaning that the K matrix is 478 00:27:44,670 --> 00:27:48,620 constant, but it is updated as the beginning 479 00:27:48,620 --> 00:27:50,010 of a solution step. 480 00:27:50,010 --> 00:27:53,160 Then we talk about the modified Newton method. 481 00:27:53,160 --> 00:27:55,860 Modified Newton-Raphson method. 482 00:27:55,860 --> 00:27:59,670 Or, we may also find it effective to update the K 483 00:27:59,670 --> 00:28:04,040 matrix only at particular solution steps, at certain 484 00:28:04,040 --> 00:28:06,570 times only. 485 00:28:06,570 --> 00:28:09,590 We note that the initial stress method and the modified 486 00:28:09,590 --> 00:28:13,810 Newton methods are, of course, much less expensive than the 487 00:28:13,810 --> 00:28:16,646 full Newton method per iteration, however. 488 00:28:16,646 --> 00:28:21,070 We should add that many more iterations are necessary to 489 00:28:21,070 --> 00:28:27,150 achieve the same accuracy if we don't set up a new K matrix 490 00:28:27,150 --> 00:28:29,210 in every iteration. 491 00:28:29,210 --> 00:28:32,030 The initial stress method and modified Newton iteration 492 00:28:32,030 --> 00:28:36,240 technique do not exhibit quadratic convergence because 493 00:28:36,240 --> 00:28:38,770 to obtain quadratic convergence, you need to set 494 00:28:38,770 --> 00:28:41,750 up a new K matrix in each iteration. 495 00:28:41,750 --> 00:28:44,920 And of course, you have to be close enough to the root, to 496 00:28:44,920 --> 00:28:46,860 the solution that you're looking for. 497 00:28:46,860 --> 00:28:51,200 Let us now look at a simple example. 498 00:28:51,200 --> 00:28:55,250 An example where we have just one degree of freedom. 499 00:28:55,250 --> 00:28:57,790 And where we deal with two load steps. 500 00:28:57,790 --> 00:29:01,580 Here we show the force applied to the 501 00:29:01,580 --> 00:29:04,210 problem, or to the structure. 502 00:29:04,210 --> 00:29:07,460 In the first load step, we apply 1R. 503 00:29:07,460 --> 00:29:10,530 And then in the second load step, 2R. 504 00:29:10,530 --> 00:29:14,570 Notice horizontally here, I'm plotting displacements. 505 00:29:14,570 --> 00:29:17,460 And the force displacement curve is 506 00:29:17,460 --> 00:29:20,130 shown by the red line. 507 00:29:20,130 --> 00:29:26,020 Now, to obtain the solution for this displacement, U1, and 508 00:29:26,020 --> 00:29:31,340 that displacement U2, we, as we discussed, set up a K 509 00:29:31,340 --> 00:29:34,600 matrix, which corresponds to the slope, a 510 00:29:34,600 --> 00:29:36,940 slope on the red curve. 511 00:29:36,940 --> 00:29:40,490 And we iterate towards the sought solutions. 512 00:29:40,490 --> 00:29:44,680 Here, we have the iterative process using the initial 513 00:29:44,680 --> 00:29:45,950 stress method. 514 00:29:45,950 --> 00:29:49,540 In other words, where tau is equal to 0 in our governing 515 00:29:49,540 --> 00:29:52,680 finite element equations. 516 00:29:52,680 --> 00:29:58,380 This means that we are setting up a K matrix initially, and 517 00:29:58,380 --> 00:30:04,960 that K matrix is depicted here by the slope, the blue slope 518 00:30:04,960 --> 00:30:06,280 that you're seeing. 519 00:30:06,280 --> 00:30:12,380 Now with this slope and this load applied, we get the 520 00:30:12,380 --> 00:30:14,390 displacement shown down here. 521 00:30:14,390 --> 00:30:19,125 Notice the left superscript means load step 1. 522 00:30:19,125 --> 00:30:24,200 The right superscript means solution after iteration 1. 523 00:30:24,200 --> 00:30:28,980 And this is the solution obtained from this triangle. 524 00:30:28,980 --> 00:30:33,990 Notice there's a blue line up there, going up there, and 525 00:30:33,990 --> 00:30:39,650 where that blue line crosses the dashed line, we pick up a 526 00:30:39,650 --> 00:30:41,730 vertical fine, black line. 527 00:30:41,730 --> 00:30:45,890 And that vertical black line cuts through the displacement 528 00:30:45,890 --> 00:30:48,570 axes at this particular value. 529 00:30:48,570 --> 00:30:51,600 So this is the solution of our first iteration. 530 00:30:51,600 --> 00:30:55,740 Now, having calculated this displacement value, we go up 531 00:30:55,740 --> 00:31:00,145 vertically and get to the red line, the red curve. 532 00:31:00,145 --> 00:31:04,710 And that red curve corresponds, of course, to the 533 00:31:04,710 --> 00:31:06,560 internal forces. 534 00:31:06,560 --> 00:31:09,340 To the internal force I should say, corresponding to this 535 00:31:09,340 --> 00:31:11,320 displacement. 536 00:31:11,320 --> 00:31:16,200 There's a red dot right on that red curve. 537 00:31:16,200 --> 00:31:20,690 And at that red dot now, we've set up a K matrix again. 538 00:31:20,690 --> 00:31:24,440 Now notice, in the initial stress method, we keep the K 539 00:31:24,440 --> 00:31:26,000 matrix constant. 540 00:31:26,000 --> 00:31:32,410 This means that the slope here of this blue line that goes 541 00:31:32,410 --> 00:31:37,830 through this red point is the same as the slope that we had 542 00:31:37,830 --> 00:31:40,060 earlier right here. 543 00:31:40,060 --> 00:31:44,670 With that slope and the out of balance load-- and let's look 544 00:31:44,670 --> 00:31:46,330 now very closely here-- 545 00:31:46,330 --> 00:31:50,190 that corresponds to the distance between this red 546 00:31:50,190 --> 00:31:56,330 point and that dashed line with this out of balance load 547 00:31:56,330 --> 00:32:00,200 and this blue slope, we get an increment in displacement 548 00:32:00,200 --> 00:32:02,730 shown right here. 549 00:32:02,730 --> 00:32:05,260 And that increment in displacement added to the 550 00:32:05,260 --> 00:32:11,780 previous displacement gives us this value down here. 551 00:32:11,780 --> 00:32:16,240 And in fact, this value is very close to the correct 552 00:32:16,240 --> 00:32:20,380 solution, the exact solution, which is the solution 553 00:32:20,380 --> 00:32:24,750 corresponding to the dashed vertical line here. 554 00:32:24,750 --> 00:32:27,180 Notice that the dashed vertical line is the exact 555 00:32:27,180 --> 00:32:30,160 solution to that dashed horizontal line. 556 00:32:30,160 --> 00:32:34,920 And our vertical black line, which we are calculating, is 557 00:32:34,920 --> 00:32:37,970 virtually on that dashed black line. 558 00:32:37,970 --> 00:32:40,900 So we can accept now, convergence. 559 00:32:40,900 --> 00:32:42,960 Let's now go into the next load step. 560 00:32:42,960 --> 00:32:46,330 In the next load step, we want to, again, 561 00:32:46,330 --> 00:32:47,620 deal with a K matrix. 562 00:32:47,620 --> 00:32:50,640 And in the initial stress method, we use the same K 563 00:32:50,640 --> 00:32:52,970 matrix throughout the solution. 564 00:32:52,970 --> 00:32:56,980 Meaning that the slope now, which we are laying onto the 565 00:32:56,980 --> 00:33:02,860 red curve, that slope being this blue line, is the same 566 00:33:02,860 --> 00:33:08,470 slope that we used before in the first two iterations. 567 00:33:08,470 --> 00:33:13,750 With this slope and the out of balance load corresponding to 568 00:33:13,750 --> 00:33:17,720 the distance between this horizontal line, the dashed 569 00:33:17,720 --> 00:33:21,780 horizontal line, and that dashed horizontal line because 570 00:33:21,780 --> 00:33:25,040 we conversed in the first load step. 571 00:33:25,040 --> 00:33:27,450 With that out of balance load, we get an incremental 572 00:33:27,450 --> 00:33:34,610 displacement shown from here to there. 573 00:33:34,610 --> 00:33:37,220 And that incremental displacement then, added to 574 00:33:37,220 --> 00:33:39,990 the previous displacement, gives us this value in 575 00:33:39,990 --> 00:33:40,720 displacement. 576 00:33:40,720 --> 00:33:46,390 Notice time t plus delta t or load step t plus delta t, 577 00:33:46,390 --> 00:33:47,600 which is here. 578 00:33:47,600 --> 00:33:52,620 Load step 2 and iteration 1, end of iteration 1. 579 00:33:52,620 --> 00:33:54,960 This is now our current displacement. 580 00:33:54,960 --> 00:33:58,670 With this current displacement, once again, 581 00:33:58,670 --> 00:34:05,760 which comes from this triangle here, with this current 582 00:34:05,760 --> 00:34:10,610 displacement, we go vertically up and get to that red point, 583 00:34:10,610 --> 00:34:13,940 which lies on the force displacement curve-- 584 00:34:13,940 --> 00:34:17,500 internal force displacement curve more accurately. 585 00:34:17,500 --> 00:34:24,510 And now we have this red point, and we set up, or we 586 00:34:24,510 --> 00:34:27,170 use now, again, a K matrix. 587 00:34:27,170 --> 00:34:29,110 Of course, in the initial stress method we use the 588 00:34:29,110 --> 00:34:30,110 constant K matrix. 589 00:34:30,110 --> 00:34:32,949 So this slope is the same as that slope. 590 00:34:32,949 --> 00:34:35,320 And now watch closely. 591 00:34:35,320 --> 00:34:42,960 From this triangle here, which corresponds to this out of 592 00:34:42,960 --> 00:34:45,130 balance load, which is the same as that 593 00:34:45,130 --> 00:34:47,820 out of balance load. 594 00:34:47,820 --> 00:34:50,580 Out of balance load obtained by putting here 595 00:34:50,580 --> 00:34:53,469 a horizontal line. 596 00:34:53,469 --> 00:34:57,100 Notice and the difference between the dashed line up 597 00:34:57,100 --> 00:35:00,965 there and the horizontal line given by my pointer is the out 598 00:35:00,965 --> 00:35:02,270 of balance load. 599 00:35:02,270 --> 00:35:05,820 Which of course, is the same as the distance between this 600 00:35:05,820 --> 00:35:08,430 red point and that dashed line. 601 00:35:08,430 --> 00:35:12,580 This out of balance load with this slope gives us another 602 00:35:12,580 --> 00:35:15,190 increment in displacement, which brings us 603 00:35:15,190 --> 00:35:17,930 to this value here. 604 00:35:17,930 --> 00:35:22,660 We repeat that process and iteratively, we get closer and 605 00:35:22,660 --> 00:35:26,340 closer to the correct solution, the exact solution, 606 00:35:26,340 --> 00:35:30,050 which is the dashed vertical line here for the 607 00:35:30,050 --> 00:35:31,700 displacement. 608 00:35:31,700 --> 00:35:36,040 Because at that dashed vertical line, we are getting 609 00:35:36,040 --> 00:35:38,660 to that red curve. 610 00:35:38,660 --> 00:35:41,070 Now, this is the initial stress method. 611 00:35:41,070 --> 00:35:46,030 And once again, the point here is that we kept the same K 612 00:35:46,030 --> 00:35:48,900 matrix throughout the solution process. 613 00:35:48,900 --> 00:35:55,710 Surely, if you now were to look at this iterative scheme 614 00:35:55,710 --> 00:35:59,120 to obtain a faster solution, you would want to 615 00:35:59,120 --> 00:36:00,770 set up a new K matrix. 616 00:36:00,770 --> 00:36:04,490 In other words, once you have this displacement value 617 00:36:04,490 --> 00:36:12,670 calculated, you want to update the slope corresponding to the 618 00:36:12,670 --> 00:36:16,180 slope of the red curve. 619 00:36:16,180 --> 00:36:18,980 And if you do so-- 620 00:36:18,980 --> 00:36:22,470 in other words, you update this blue slope corresponding 621 00:36:22,470 --> 00:36:26,980 to the red dots that you have obtained in 622 00:36:26,980 --> 00:36:28,640 your iterative scheme. 623 00:36:28,640 --> 00:36:31,300 These red dots, of course, being different from what 624 00:36:31,300 --> 00:36:32,740 you're seeing here now. 625 00:36:32,740 --> 00:36:37,440 Then you would have the full Newton-Raphson method. 626 00:36:37,440 --> 00:36:42,750 If you only update the K matrix or the blue tangent 627 00:36:42,750 --> 00:36:47,230 that we are seeing here whenever you have converged to 628 00:36:47,230 --> 00:36:50,680 an equilibrium point. 629 00:36:50,680 --> 00:36:54,870 Meaning when you have converged corresponding for 630 00:36:54,870 --> 00:36:58,710 the first load step in this particular case. 631 00:36:58,710 --> 00:37:01,120 If you only update the stiffness matrix once you have 632 00:37:01,120 --> 00:37:03,560 converged to this load-- 633 00:37:03,560 --> 00:37:06,690 for this load step, then you would have the modified 634 00:37:06,690 --> 00:37:09,060 Newton-Raphson method. 635 00:37:09,060 --> 00:37:13,420 I think if you look at this picture with these 636 00:37:13,420 --> 00:37:16,700 possibilities, you will realize that surely, the full 637 00:37:16,700 --> 00:37:21,080 Newton-Raphson method with updating the slope after each 638 00:37:21,080 --> 00:37:23,970 iteration will converge fastest. 639 00:37:23,970 --> 00:37:26,690 And the modified Newton-Raphson method will 640 00:37:26,690 --> 00:37:29,150 converge a little slower than the full Newton-Raphson 641 00:37:29,150 --> 00:37:34,260 method, but still faster than the initial stress method. 642 00:37:34,260 --> 00:37:38,210 Well, this then brings us to the end of the discussion of 643 00:37:38,210 --> 00:37:39,080 this example. 644 00:37:39,080 --> 00:37:42,800 I like to now introduce you to another scheme that can be 645 00:37:42,800 --> 00:37:46,430 very important when we solve nonlinear equations, finite 646 00:37:46,430 --> 00:37:47,770 element equations. 647 00:37:47,770 --> 00:37:52,290 And that scheme is the scheme of line searches. 648 00:37:52,290 --> 00:37:55,480 The basic equations that we are looking at that we want to 649 00:37:55,480 --> 00:37:57,380 solve are once more depicted here. 650 00:37:57,380 --> 00:37:59,990 Tangent stiffness matrix, incremental displacement 651 00:37:59,990 --> 00:38:03,350 vector, externally applied loads, nodal point forces 652 00:38:03,350 --> 00:38:05,790 corresponding to the internal element stresses at time t 653 00:38:05,790 --> 00:38:10,140 plus delta t, and at the end of iteration i minus 1, 654 00:38:10,140 --> 00:38:12,640 beginning of iteration i, of course. 655 00:38:12,640 --> 00:38:15,070 Notice I don't have an i on here, I just want to denote 656 00:38:15,070 --> 00:38:18,220 this for the moment as an incremental nodal point 657 00:38:18,220 --> 00:38:22,200 displacement vector, which carries, however, a bar. 658 00:38:22,200 --> 00:38:29,560 Let us consider forming Fi using this equation where beta 659 00:38:29,560 --> 00:38:31,800 is an unknown. 660 00:38:31,800 --> 00:38:35,460 We want to choose beta now judiciously. 661 00:38:35,460 --> 00:38:40,380 We want to choose beta such that, in some sense, R minus 662 00:38:40,380 --> 00:38:44,580 Fi is small. 663 00:38:44,580 --> 00:38:48,515 And we have to, of course, look at the mechanism that we 664 00:38:48,515 --> 00:38:53,560 use to make R minus F small in some sense. 665 00:38:53,560 --> 00:38:59,040 Well, as a side note, we recognize that when this 666 00:38:59,040 --> 00:39:06,180 equation here is 0 for all possible U's, then clearly 667 00:39:06,180 --> 00:39:08,640 this must be 0. 668 00:39:08,640 --> 00:39:14,380 In other words, if we were to say to ourselves, let's make 669 00:39:14,380 --> 00:39:17,410 this equation 0 for all possible U's. 670 00:39:17,410 --> 00:39:22,770 Then really, we would have solved this equal to 0. 671 00:39:22,770 --> 00:39:25,640 The reason, of course, for it is that we could choose here 672 00:39:25,640 --> 00:39:31,310 for U a vector, which carries 0's everywhere except a 1 in 673 00:39:31,310 --> 00:39:32,610 one location. 674 00:39:32,610 --> 00:39:36,880 And if that is a particular row, then corresponding to 675 00:39:36,880 --> 00:39:44,170 that row, surely this vector here, this R minus F 676 00:39:44,170 --> 00:39:46,900 corresponding to that row must be equal to 0. 677 00:39:46,900 --> 00:39:50,390 So if we were simply to use here the identity matrix, or 678 00:39:50,390 --> 00:39:54,660 if we were to construct the identity matrix by allowing n 679 00:39:54,660 --> 00:39:57,730 U's that are linearly independent, of course, 680 00:39:57,730 --> 00:40:03,480 because they all make up such vector but with a 1 in a 681 00:40:03,480 --> 00:40:04,500 different location. 682 00:40:04,500 --> 00:40:08,950 As a matter of fact, we want to use U1 with a 1 here, 0's 683 00:40:08,950 --> 00:40:09,850 everywhere. 684 00:40:09,850 --> 00:40:14,650 U2, 1 here, 0 here, everywhere else 0, and so on. 685 00:40:14,650 --> 00:40:16,870 To construct an identity matrix here. 686 00:40:16,870 --> 00:40:20,500 Then with that identity matrix, clearly R minus F 687 00:40:20,500 --> 00:40:23,090 would be 0 everywhere. 688 00:40:23,090 --> 00:40:27,250 Of course, that is not a viable scheme really. 689 00:40:27,250 --> 00:40:34,170 But what we want to do is to choose a particular vector 690 00:40:34,170 --> 00:40:40,430 here so as to make R minus F still 0 in some sense. 691 00:40:40,430 --> 00:40:44,760 And the vector that is quite effectively chosen is the 692 00:40:44,760 --> 00:40:49,460 vector that we obtain from the solution of K delta U bar 693 00:40:49,460 --> 00:40:51,430 equals R minus F. 694 00:40:51,430 --> 00:40:58,070 If we use that vector here, basically projecting the R, 695 00:40:58,070 --> 00:41:01,650 capital R, minus F onto that vector. 696 00:41:01,650 --> 00:41:06,770 And we look for this to be 0. 697 00:41:06,770 --> 00:41:12,350 Then we have so to say, searched in the direction of 698 00:41:12,350 --> 00:41:20,210 delta U bar and made delta U bar transposed times this 699 00:41:20,210 --> 00:41:22,360 vector here equal to 0. 700 00:41:22,360 --> 00:41:23,970 How does this scheme work? 701 00:41:23,970 --> 00:41:29,240 Well, we add beta times delta U bar to the value that we had 702 00:41:29,240 --> 00:41:31,600 already from the previous iteration. 703 00:41:31,600 --> 00:41:39,340 And we search along beta such that this top here is much 704 00:41:39,340 --> 00:41:41,360 smaller than the bottom. 705 00:41:41,360 --> 00:41:45,650 In other words, STOL shall be a convergence tolerance. 706 00:41:45,650 --> 00:41:47,940 Notice the way it works. 707 00:41:47,940 --> 00:41:49,630 We choose a value of beta. 708 00:41:49,630 --> 00:41:53,550 With beta known, we add this quantity to this one. 709 00:41:53,550 --> 00:41:55,470 We get a value here. 710 00:41:55,470 --> 00:42:00,470 We take this value, calculate Fi corresponding to that 711 00:42:00,470 --> 00:42:02,510 value, see whether the convergence 712 00:42:02,510 --> 00:42:03,520 tolerance is satisfied. 713 00:42:03,520 --> 00:42:06,090 If not, we have to choose a new beta. 714 00:42:06,090 --> 00:42:13,590 And like that, we choose new betas until we have satisfied 715 00:42:13,590 --> 00:42:15,360 this convergence tolerance. 716 00:42:15,360 --> 00:42:17,320 That's what we call line search. 717 00:42:17,320 --> 00:42:21,700 We are searching along the line given by delta U bar 718 00:42:21,700 --> 00:42:25,110 until this is satisfied. 719 00:42:25,110 --> 00:42:29,210 It's a very important scheme, and this scheme can be 720 00:42:29,210 --> 00:42:32,590 combined with modified Newton iteration, with the initial 721 00:42:32,590 --> 00:42:35,590 stress method, and with full Newton iteration. 722 00:42:35,590 --> 00:42:40,680 And it adds a lot to the stability of the solution, of 723 00:42:40,680 --> 00:42:43,660 the overall solution or the nonlinear equations. 724 00:42:43,660 --> 00:42:47,470 We will find, or discuss applications just now. 725 00:42:47,470 --> 00:42:52,890 But at this point, we need to stop for a minute because we 726 00:42:52,890 --> 00:42:54,390 need to go onto a new tape. 727 00:42:58,160 --> 00:43:00,970 Finally, I would like to discuss with you a method that 728 00:43:00,970 --> 00:43:03,540 we also find to be very effective in nonlinear finite 729 00:43:03,540 --> 00:43:05,050 element computations. 730 00:43:05,050 --> 00:43:07,590 Namely, the BFGS method. 731 00:43:07,590 --> 00:43:08,835 BFGS stands for Broyden-Fletcher 732 00:43:08,835 --> 00:43:10,085 -Goldfarb-Shanno. 733 00:43:12,200 --> 00:43:15,930 And in this method, we proceed as follows. 734 00:43:15,930 --> 00:43:20,850 We calculate or define a delta i as shown 735 00:43:20,850 --> 00:43:23,070 here in this equation. 736 00:43:23,070 --> 00:43:28,560 And we define a gamma i as shown in this equation. 737 00:43:28,560 --> 00:43:32,200 We want to calculate the coefficient matrix such that 738 00:43:32,200 --> 00:43:34,110 this equation is here satisfied. 739 00:43:34,110 --> 00:43:37,720 In other words, given delta i and gamma i, we want to use 740 00:43:37,720 --> 00:43:41,040 the coefficient matrix that satisfies this equation. 741 00:43:41,040 --> 00:43:45,680 And that's basically what's being done in the BFGS scheme. 742 00:43:45,680 --> 00:43:50,610 Pictorially then, we can show for a single degree of freedom 743 00:43:50,610 --> 00:43:55,260 system, say the F plotted along this axis. 744 00:43:55,260 --> 00:43:58,730 The U displacements plotted along this axis. 745 00:43:58,730 --> 00:44:01,770 Notice here in green we show delta i. 746 00:44:01,770 --> 00:44:04,520 And we show here in green gamma i. 747 00:44:04,520 --> 00:44:08,890 And notice that the matrix that we want to use is 748 00:44:08,890 --> 00:44:11,760 indicated by this blue slope. 749 00:44:11,760 --> 00:44:16,560 Which is, in other words, given by delta i and gamma i. 750 00:44:16,560 --> 00:44:21,300 Well, the BFGS method is an iterative algorithm, which 751 00:44:21,300 --> 00:44:26,090 produces successive approximations to efficient 752 00:44:26,090 --> 00:44:27,980 stiffness matrix. 753 00:44:27,980 --> 00:44:30,410 Actually, we actually work with the inverse of that 754 00:44:30,410 --> 00:44:34,450 stiffness matrix as I will elaborate upon just now. 755 00:44:34,450 --> 00:44:36,770 It's really a compromise between the full Newton 756 00:44:36,770 --> 00:44:38,600 iteration method and the modified 757 00:44:38,600 --> 00:44:40,590 Newton iteration method. 758 00:44:40,590 --> 00:44:43,350 It shows very excellent convergence characteristics in 759 00:44:43,350 --> 00:44:45,950 the solution of many types of problems. 760 00:44:45,950 --> 00:44:48,680 Let's look at the individual steps that we use when we 761 00:44:48,680 --> 00:44:50,900 apply the BFGS method. 762 00:44:50,900 --> 00:44:54,530 We, first of all, calculate the direction of the 763 00:44:54,530 --> 00:44:56,160 displacement increment. 764 00:44:56,160 --> 00:45:03,090 Here we have delta U bar i is equal to K inverse, the 765 00:45:03,090 --> 00:45:06,120 inverse of the K matrix, corresponding to 766 00:45:06,120 --> 00:45:08,050 iteration i minus 1. 767 00:45:08,050 --> 00:45:12,150 And here we have the out of balance load vector. 768 00:45:12,150 --> 00:45:15,170 Notice once again we do not calculate actually an inverse 769 00:45:15,170 --> 00:45:16,080 of a matrix. 770 00:45:16,080 --> 00:45:20,730 We calculate the LDL transpose factorization of a matrix, and 771 00:45:20,730 --> 00:45:25,690 then we, as I will just now show, update that inverse in 772 00:45:25,690 --> 00:45:31,445 the BFGS iteration by specific matrices of rank 2. 773 00:45:31,445 --> 00:45:33,610 We get to that just now. 774 00:45:33,610 --> 00:45:39,360 Well, having calculated this delta U bar i, we use that 775 00:45:39,360 --> 00:45:44,850 delta U bar i in the line search to obtain a better 776 00:45:44,850 --> 00:45:49,650 displacement vector, a displacement vector 777 00:45:49,650 --> 00:45:51,550 corresponding to iteration i. 778 00:45:51,550 --> 00:45:56,450 A better value then just adding the delta U bar i to it 779 00:45:56,450 --> 00:45:57,910 by adjusting beta. 780 00:45:57,910 --> 00:46:02,960 In other words, we choose beta in such way as to satisfy this 781 00:46:02,960 --> 00:46:04,250 equation here. 782 00:46:04,250 --> 00:46:06,060 This is, of course, the equation that we are now 783 00:46:06,060 --> 00:46:11,680 familiar with as being equation that shows that 784 00:46:11,680 --> 00:46:16,300 convergence has been reached in a line search. 785 00:46:16,300 --> 00:46:22,650 Having now calculated the t plus delta t Ui and the 786 00:46:22,650 --> 00:46:28,390 corresponding t plus delta t Fi, we can calculate the delta 787 00:46:28,390 --> 00:46:30,660 i and gamma i. 788 00:46:30,660 --> 00:46:33,930 And therefore, apply the BFGS method. 789 00:46:33,930 --> 00:46:36,630 And that's now step three, calculation of 790 00:46:36,630 --> 00:46:38,760 the new secant matrix. 791 00:46:38,760 --> 00:46:46,420 K inverse of iteration i is equal to Ai transposed k 792 00:46:46,420 --> 00:46:50,100 inverse of iteration i minus 1 times Ai. 793 00:46:50,100 --> 00:46:55,393 This A matrix is obtained as shown here. 794 00:46:55,393 --> 00:46:57,910 v wi transposed. 795 00:46:57,910 --> 00:47:04,390 The vi and wi's are given as a function of the delta i, gamma 796 00:47:04,390 --> 00:47:10,450 i, and K i minus 1 value. 797 00:47:10,450 --> 00:47:15,120 You should look at the textbook to see the details of 798 00:47:15,120 --> 00:47:18,590 actually evaluating the vi and wi. 799 00:47:18,590 --> 00:47:23,430 But once you have vi, wi, you can calculate A. And clearly, 800 00:47:23,430 --> 00:47:29,100 by this equation then, you get a new stiffness matrix. 801 00:47:29,100 --> 00:47:33,670 It is important to note that in this iterative scheme, only 802 00:47:33,670 --> 00:47:38,550 vector products are needed to obtain vi and wi. 803 00:47:38,550 --> 00:47:41,680 And it's also then important to note that only vector 804 00:47:41,680 --> 00:47:45,520 products are used to calculate delta U bar i. 805 00:47:45,520 --> 00:47:50,040 If this would not be the case, particularly here. 806 00:47:50,040 --> 00:47:51,540 Then of course, the iterative scheme 807 00:47:51,540 --> 00:47:53,470 would be very expensive. 808 00:47:53,470 --> 00:47:59,960 Let me show you a bit of the details why we only need to 809 00:47:59,960 --> 00:48:03,230 use vector products to get the solution. 810 00:48:03,230 --> 00:48:09,230 Here, we have one typical step of the iteration. 811 00:48:09,230 --> 00:48:13,730 Notice on the left-hand side, we have delta U bar i. 812 00:48:13,730 --> 00:48:15,890 The value of displacement increment 813 00:48:15,890 --> 00:48:17,880 that we want to calculate. 814 00:48:17,880 --> 00:48:22,460 On the right-hand side, we have i plus wi minus 1 vi 815 00:48:22,460 --> 00:48:23,990 minus 1 transposed. 816 00:48:23,990 --> 00:48:28,890 In other words, the Ai minus 1 transposed. 817 00:48:28,890 --> 00:48:33,860 And we would have such matrices going on here until 818 00:48:33,860 --> 00:48:37,670 we are coming to the A1 transposed matrix. 819 00:48:37,670 --> 00:48:40,830 Then we have the inverse of a stiffness matrix. 820 00:48:40,830 --> 00:48:43,860 Here, corresponding to time tau. 821 00:48:43,860 --> 00:48:50,590 And then we have the Ai matrices as shown here. 822 00:48:50,590 --> 00:48:57,590 Now, notice that this total amount then is multiplied by 823 00:48:57,590 --> 00:49:02,690 the R minus F i minus 1, the out of balance load vector. 824 00:49:02,690 --> 00:49:06,370 Let's now go through the computations that are required 825 00:49:06,370 --> 00:49:10,500 to get delta U bar i. 826 00:49:10,500 --> 00:49:14,850 Here we have delta R. This delta R, the out of balance 827 00:49:14,850 --> 00:49:17,030 load vector, which I just call delta R, 828 00:49:17,030 --> 00:49:19,110 multiplies this matrix. 829 00:49:19,110 --> 00:49:22,010 But if you look at this multiplication, we find that 830 00:49:22,010 --> 00:49:29,430 this delta R times w i minus 1 transposed is simply a number. 831 00:49:29,430 --> 00:49:31,460 So it's one vector multiplication 832 00:49:31,460 --> 00:49:33,020 that gives us a number. 833 00:49:33,020 --> 00:49:39,110 This number multiplied by v i minus 1 is just the 834 00:49:39,110 --> 00:49:41,100 multiplication of a number times a vector. 835 00:49:41,100 --> 00:49:42,930 It's very simple. 836 00:49:42,930 --> 00:49:48,230 And notice, of course, this delta R i minus 1 is also 837 00:49:48,230 --> 00:49:49,850 multiplying the identity matrix. 838 00:49:49,850 --> 00:49:53,910 So in this step here, we have really nothing else then a 839 00:49:53,910 --> 00:49:55,080 vector multiplication. 840 00:49:55,080 --> 00:49:58,040 The only vector multiplication that, in fact, is there, is 841 00:49:58,040 --> 00:50:02,030 this value, this vector here, times that vector there, 842 00:50:02,030 --> 00:50:03,510 transposed. 843 00:50:03,510 --> 00:50:04,610 That gives us a number. 844 00:50:04,610 --> 00:50:07,120 This number times a vector [UNINTELLIGIBLE] 845 00:50:07,120 --> 00:50:08,040 operation. 846 00:50:08,040 --> 00:50:11,780 And of course, this one here times the i matrix is just 847 00:50:11,780 --> 00:50:13,500 this vector by itself. 848 00:50:13,500 --> 00:50:16,690 So in this step, what we obtain is 849 00:50:16,690 --> 00:50:19,950 simply another vector. 850 00:50:19,950 --> 00:50:24,240 We take this vector and repeat the process. 851 00:50:24,240 --> 00:50:28,500 And since the repetition is just what we discussed 852 00:50:28,500 --> 00:50:34,470 already, we find that we only need vector multiplications to 853 00:50:34,470 --> 00:50:38,790 roll up all of these multiplications up to here. 854 00:50:38,790 --> 00:50:43,045 Then we have a vector multiplied by K inverse. 855 00:50:43,045 --> 00:50:47,380 Well, that involves only vector multiplications again, 856 00:50:47,380 --> 00:50:49,810 because we have already the LDL transposed 857 00:50:49,810 --> 00:50:51,460 factorization of tau k. 858 00:50:54,140 --> 00:50:57,080 These vector multiplications result again, in a vector. 859 00:50:57,080 --> 00:51:00,880 And that vector multiplied by these matrices here in the 860 00:51:00,880 --> 00:51:04,330 same way as we have proceeded here, means again, only vector 861 00:51:04,330 --> 00:51:05,930 multiplications. 862 00:51:05,930 --> 00:51:10,440 So this very important step here only requires vector 863 00:51:10,440 --> 00:51:11,140 multiplications. 864 00:51:11,140 --> 00:51:16,520 And that makes this whole process really efficient. 865 00:51:16,520 --> 00:51:20,740 In summary then, we have the following solution procedures 866 00:51:20,740 --> 00:51:23,130 that we feel are very effective. 867 00:51:23,130 --> 00:51:27,020 The mortified Newton-Raphson iteration with line searches. 868 00:51:27,020 --> 00:51:29,540 Here we have the basic modified 869 00:51:29,540 --> 00:51:31,190 Newton iteration equation. 870 00:51:31,190 --> 00:51:34,110 And here we have the line search equation. 871 00:51:34,110 --> 00:51:38,930 We, of course, after each such solution, evaluate an 872 00:51:38,930 --> 00:51:39,810 efficient beta. 873 00:51:39,810 --> 00:51:43,600 A beta that makes U-- 874 00:51:43,600 --> 00:51:48,660 t plus delta U i a good vector when measured on the line 875 00:51:48,660 --> 00:51:52,010 search criteria the way I discussed it. 876 00:51:52,010 --> 00:51:56,210 And these two equations then, summarize really the modified 877 00:51:56,210 --> 00:51:58,475 Newton iteration with the line search. 878 00:51:58,475 --> 00:52:01,750 The BFGS method is another effective scheme. 879 00:52:01,750 --> 00:52:05,210 We use the BFGS method always with line searches. 880 00:52:05,210 --> 00:52:07,950 And then, as a third option, the full Newton-Raphson 881 00:52:07,950 --> 00:52:10,740 iteration with or without line searches. 882 00:52:10,740 --> 00:52:13,220 The full Newton-Raphson iteration with line searches 883 00:52:13,220 --> 00:52:14,890 is, of course, most powerful. 884 00:52:14,890 --> 00:52:18,210 But it is also most expensive per iteration. 885 00:52:18,210 --> 00:52:21,020 We should point out that these techniques that I discussed so 886 00:52:21,020 --> 00:52:24,760 far cannot be applied for post-buckling analysis. 887 00:52:24,760 --> 00:52:29,160 We will consider the solution of the response after buckling 888 00:52:29,160 --> 00:52:32,520 has been occurring, after the item load level has been 889 00:52:32,520 --> 00:52:36,020 reached in a forthcoming lecture. 890 00:52:36,020 --> 00:52:41,610 The next view graphs summarize these schemes: the modified 891 00:52:41,610 --> 00:52:44,450 Newton-Raphson, the BFGS method, and the full 892 00:52:44,450 --> 00:52:46,840 Newton-Raphson method. 893 00:52:46,840 --> 00:52:50,030 And I'd like to just very quickly go through the whole 894 00:52:50,030 --> 00:52:54,070 solution process for each of these methods. 895 00:52:54,070 --> 00:52:57,700 In the modified Newton iteration technique, we 896 00:52:57,700 --> 00:53:03,580 initialize our displacements corresponding to the beginning 897 00:53:03,580 --> 00:53:09,610 of the first iteration, the end of the 0 iteration as tU. 898 00:53:09,610 --> 00:53:12,920 And similarly, we initialize our force vector corresponding 899 00:53:12,920 --> 00:53:14,470 to the internal element stresses. 900 00:53:14,470 --> 00:53:19,010 We set i is equal to 1, the iteration counter equal to 1. 901 00:53:19,010 --> 00:53:22,610 We calculate a K matrix at the beginning of the load step, 902 00:53:22,610 --> 00:53:25,510 and we keep that K matrix constant. 903 00:53:25,510 --> 00:53:29,450 We then go into this step here where we calculate the out of 904 00:53:29,450 --> 00:53:30,700 balance loads. 905 00:53:33,210 --> 00:53:36,790 Calculate a new displacement vector. 906 00:53:36,790 --> 00:53:41,240 And in this box here, we measure whether we have 907 00:53:41,240 --> 00:53:42,160 converged already. 908 00:53:42,160 --> 00:53:45,390 Of course, this convergence would only be measured after 909 00:53:45,390 --> 00:53:49,560 we have gone say, at least twice through 910 00:53:49,560 --> 00:53:51,070 the iteration process. 911 00:53:51,070 --> 00:53:54,650 Or rather, we have calculated this increment in displacement 912 00:53:54,650 --> 00:53:56,020 vector twice. 913 00:53:56,020 --> 00:53:59,560 Because at the beginning corresponding to i equal to 1, 914 00:53:59,560 --> 00:54:03,400 we have now this fairly large out of balance load. 915 00:54:03,400 --> 00:54:05,870 We have just incremented the load vector. 916 00:54:05,870 --> 00:54:10,880 Therefore, this right inside should be non-zero. 917 00:54:10,880 --> 00:54:13,580 And we would get an increment in displacements. 918 00:54:13,580 --> 00:54:17,500 And to measure how large that increment in displacement is, 919 00:54:17,500 --> 00:54:20,350 we really want to go through this whole cycle at 920 00:54:20,350 --> 00:54:22,170 least one more time. 921 00:54:22,170 --> 00:54:25,610 We then, having calculated this incremental displacement 922 00:54:25,610 --> 00:54:30,620 vector with the bar on there, we perform a line search the 923 00:54:30,620 --> 00:54:34,230 way we discussed it to update our displacements 924 00:54:34,230 --> 00:54:36,720 corresponding to the iteration i. 925 00:54:39,960 --> 00:54:42,300 Of course, we only have come so far once through it, so i 926 00:54:42,300 --> 00:54:43,520 is equal to 1 still. 927 00:54:43,520 --> 00:54:48,360 And this gives us the accepted value of displacements 928 00:54:48,360 --> 00:54:50,936 corresponding to the first iteration. 929 00:54:50,936 --> 00:54:55,560 We now increment our iteration counter, and go in to the 930 00:54:55,560 --> 00:54:57,370 second iteration. 931 00:54:57,370 --> 00:55:01,060 Notice now in the second iteration, we have R minus t 932 00:55:01,060 --> 00:55:03,670 plus delta t F1 here. 933 00:55:03,670 --> 00:55:06,250 We calculate delta U bar 2. 934 00:55:06,250 --> 00:55:10,040 And now we would certainly measure convergence. 935 00:55:10,040 --> 00:55:13,360 And I will just now talk about how we measure convergence. 936 00:55:13,360 --> 00:55:17,240 If we have not converged yet, we keep on cycling through 937 00:55:17,240 --> 00:55:22,080 here until we do actually get convergence. 938 00:55:22,080 --> 00:55:25,840 The distinguishing feature of the modified Newton iteration 939 00:55:25,840 --> 00:55:29,642 is that we are having a constant K matrix. 940 00:55:29,642 --> 00:55:32,900 That we do not update the K matrix is this 941 00:55:32,900 --> 00:55:34,250 iterative loop here. 942 00:55:34,250 --> 00:55:36,280 It remains constant. 943 00:55:36,280 --> 00:55:40,570 Remember in the initial stress technique, initial stress 944 00:55:40,570 --> 00:55:46,020 method, we would not even update the K matrix ever after 945 00:55:46,020 --> 00:55:49,980 the first time it has been calculated. 946 00:55:49,980 --> 00:55:52,990 But in the modified Newton iteration, we update it at the 947 00:55:52,990 --> 00:55:55,800 beginning of each load step. 948 00:55:55,800 --> 00:56:00,140 In the BFGS method, we proceed as follows. 949 00:56:00,140 --> 00:56:03,420 We initial, again, our displacements and our nodal 950 00:56:03,420 --> 00:56:05,350 point force vector corresponding to the internal 951 00:56:05,350 --> 00:56:06,570 element stresses. 952 00:56:06,570 --> 00:56:08,470 We calculate tK. 953 00:56:08,470 --> 00:56:11,630 We calculate our incremental displacement vector with the 954 00:56:11,630 --> 00:56:13,120 bar on top. 955 00:56:13,120 --> 00:56:18,670 We measure convergence if i is greater or equal to. 956 00:56:18,670 --> 00:56:20,700 There's no point as I just mentioned, measuring 957 00:56:20,700 --> 00:56:24,070 convergence when you go first time through here. 958 00:56:24,070 --> 00:56:29,020 And we perform then the line search having calculated our 959 00:56:29,020 --> 00:56:32,440 actual incremental in nodal point displacement that we 960 00:56:32,440 --> 00:56:35,410 want to add to the previous nodal point displacements to 961 00:56:35,410 --> 00:56:39,150 get to this displacement vector. 962 00:56:39,150 --> 00:56:42,040 We, of course, also have this force vector. 963 00:56:42,040 --> 00:56:46,960 We now can update our inverse matrix, our inverse secret 964 00:56:46,960 --> 00:56:50,740 matrix the way we discussed it. 965 00:56:50,740 --> 00:56:53,120 And we increment our counter. 966 00:56:53,120 --> 00:56:58,110 And we keep on cycling through here until we find in this box 967 00:56:58,110 --> 00:57:01,300 that we have converged. 968 00:57:01,300 --> 00:57:06,870 In the full Newton iteration with line searches, we once 969 00:57:06,870 --> 00:57:10,460 again, initialize our displacement vector, our force 970 00:57:10,460 --> 00:57:13,490 vector, our iteration counter. 971 00:57:13,490 --> 00:57:17,710 We now calculate a new K matrix corresponding to the 972 00:57:17,710 --> 00:57:20,720 last conditions on displacements and internal 973 00:57:20,720 --> 00:57:22,320 forces, of course. 974 00:57:22,320 --> 00:57:25,140 We enter here to calculate our incremental displacement 975 00:57:25,140 --> 00:57:29,460 vector with the K matrix that was just evaluated. 976 00:57:32,110 --> 00:57:36,740 And we measure convergence once again only when i is 977 00:57:36,740 --> 00:57:38,460 greater equal to. 978 00:57:38,460 --> 00:57:45,450 We perform a line search to determine beta in here. 979 00:57:45,450 --> 00:57:49,010 We therefore, update our displacement vector to the 980 00:57:49,010 --> 00:57:51,570 displacements that now correspond to the end of 981 00:57:51,570 --> 00:57:53,230 iteration i. 982 00:57:53,230 --> 00:57:55,250 We update our iteration counter. 983 00:57:55,250 --> 00:57:57,720 Or rather, increase our iteration counter. 984 00:57:57,720 --> 00:58:02,580 And we are calculating a new stiffness matrix corresponding 985 00:58:02,580 --> 00:58:05,200 to the last calculated displacement. 986 00:58:05,200 --> 00:58:08,180 This is, of course, a distinguishing feature of the 987 00:58:08,180 --> 00:58:10,660 full Newton-Raphson method that we are calculating a new 988 00:58:10,660 --> 00:58:13,700 K matrix at the beginning of each iteration. 989 00:58:16,210 --> 00:58:21,360 Well, these view graphs then, summarize the full process of 990 00:58:21,360 --> 00:58:24,770 iteration for the modified Newton, the BFGS, and the full 991 00:58:24,770 --> 00:58:25,440 Newton method. 992 00:58:25,440 --> 00:58:28,660 Notice that I have included the line search as the option. 993 00:58:28,660 --> 00:58:31,480 If you don't have line searching, of course, it would 994 00:58:31,480 --> 00:58:37,290 simply mean that you set beta equal to 1, and skip that 995 00:58:37,290 --> 00:58:39,680 particular step of line searching. 996 00:58:39,680 --> 00:58:43,530 But in many cases, line searching is quite important. 997 00:58:43,530 --> 00:58:47,440 And certainly, can be necessary in 998 00:58:47,440 --> 00:58:49,790 some types of problems. 999 00:58:49,790 --> 00:58:52,690 Let us now look at what happens in this box here. 1000 00:58:52,690 --> 00:58:55,260 How do we measure convergence? 1001 00:58:55,260 --> 00:58:58,850 In each of these iterative schemes, we had such a box. 1002 00:58:58,850 --> 00:59:02,730 And I like to look now more closely at how we measure 1003 00:59:02,730 --> 00:59:05,040 convergence. 1004 00:59:05,040 --> 00:59:09,090 The measure of convergence should, of course, tell us, 1005 00:59:09,090 --> 00:59:11,740 how well do we satisfy equilibrium? 1006 00:59:11,740 --> 00:59:16,070 And there are basically three items that can be used. 1007 00:59:16,070 --> 00:59:20,870 Namely energy, force, or moment, displacement, and 1008 00:59:20,870 --> 00:59:23,110 rotation, of course, correspondingly as well. 1009 00:59:23,110 --> 00:59:25,940 But these are the items that one can work with to measure 1010 00:59:25,940 --> 00:59:27,280 convergence. 1011 00:59:27,280 --> 00:59:31,510 If you use an energy convergence criteria, and we 1012 00:59:31,510 --> 00:59:33,800 are very fond of this one. 1013 00:59:33,800 --> 00:59:40,970 You might want to use this equation here delta U bar i 1014 00:59:40,970 --> 00:59:46,800 transposed times the out of balance load divided by delta 1015 00:59:46,800 --> 00:59:50,870 U bar 1 transposed times t plus delta t R minus tF. 1016 00:59:50,870 --> 00:59:54,500 Now notice, this is the increment in nodal point 1017 00:59:54,500 --> 00:59:58,020 forces, or the out of balance load, at the 1018 00:59:58,020 --> 01:00:00,990 beginning of the load step. 1019 01:00:00,990 --> 01:00:05,850 So this here is an energy corresponding to the beginning 1020 01:00:05,850 --> 01:00:07,090 of the load step. 1021 01:00:07,090 --> 01:00:11,460 And we are comparing here a similar quantity, but 1022 01:00:11,460 --> 01:00:13,760 corresponding to the current load. 1023 01:00:13,760 --> 01:00:16,140 Corresponding to the current iteration. 1024 01:00:16,140 --> 01:00:21,190 And if this quotient here is very small measured on that 1025 01:00:21,190 --> 01:00:25,250 convergence tolerance, then we say we have converged. 1026 01:00:25,250 --> 01:00:28,910 An important point is that in this convergence measure, we 1027 01:00:28,910 --> 01:00:33,400 enter with displacements and with out of balance loads. 1028 01:00:33,400 --> 01:00:36,450 Both enter in this convergence measure. 1029 01:00:36,450 --> 01:00:41,270 And please realize also that if the loads don't increase 1030 01:00:41,270 --> 01:00:46,960 from time t to t plus delta t, meaning that this here is 1031 01:00:46,960 --> 01:00:49,330 equal to 0-- 1032 01:00:49,330 --> 01:00:50,410 0 vector. 1033 01:00:50,410 --> 01:00:52,860 Then, of course, we get 0 down here. 1034 01:00:52,860 --> 01:00:55,020 And we could not apply this convergence criteria 1035 01:00:55,020 --> 01:00:56,170 indirectly. 1036 01:00:56,170 --> 01:00:59,900 Therefore, we would do is instead of 0 here, a 1037 01:00:59,900 --> 01:01:03,890 reasonable small number in the computer program. 1038 01:01:03,890 --> 01:01:06,760 Also note that we apply this test here 1039 01:01:06,760 --> 01:01:09,300 prior to line searching. 1040 01:01:09,300 --> 01:01:13,590 In line searching, after line searching of course, this top 1041 01:01:13,590 --> 01:01:15,580 part here is small. 1042 01:01:15,580 --> 01:01:20,070 And we, therefore, do not want to apply this scheme here or 1043 01:01:20,070 --> 01:01:22,290 this convergence measure after line searching. 1044 01:01:22,290 --> 01:01:25,530 We should do it before line searching. 1045 01:01:25,530 --> 01:01:31,040 To measure convergence on forces, we can use this 1046 01:01:31,040 --> 01:01:32,420 equation here. 1047 01:01:32,420 --> 01:01:35,620 Notice here we have the out of balance loads. 1048 01:01:35,620 --> 01:01:39,640 We take the Euclidean norm on the out of balance loads, and 1049 01:01:39,640 --> 01:01:43,710 we compare that with RNORM, input by the 1050 01:01:43,710 --> 01:01:45,970 user, by the analyst. 1051 01:01:45,970 --> 01:01:50,560 It is interesting, important to note that we introduce here 1052 01:01:50,560 --> 01:01:53,860 only the components corresponding to the 1053 01:01:53,860 --> 01:01:57,220 translational degrees of freedom. 1054 01:01:57,220 --> 01:02:01,540 And here, of course, this RNORM is also a force 1055 01:02:01,540 --> 01:02:03,340 corresponding, so to say, to a 1056 01:02:03,340 --> 01:02:04,840 translational degree of freedom. 1057 01:02:04,840 --> 01:02:09,430 If we have also rotational degrees of freedom, then we 1058 01:02:09,430 --> 01:02:14,500 would take instead of RNORM, RMNORM. 1059 01:02:14,500 --> 01:02:18,760 And we would extract from these vectors, the components 1060 01:02:18,760 --> 01:02:20,130 corresponding to the rotational 1061 01:02:20,130 --> 01:02:21,700 degrees of freedom only. 1062 01:02:21,700 --> 01:02:24,530 The point is that we want to have the same units up 1063 01:02:24,530 --> 01:02:27,450 here as down here. 1064 01:02:27,450 --> 01:02:31,790 Forces, forces, moments, moments. 1065 01:02:31,790 --> 01:02:36,440 And we don't want to mix components from these 1066 01:02:36,440 --> 01:02:38,280 quantities. 1067 01:02:38,280 --> 01:02:43,160 Of course, RMNORM and RNORM must be chosen by the analyst. 1068 01:02:43,160 --> 01:02:48,190 The are input to the analysis as defined in the input to the 1069 01:02:48,190 --> 01:02:49,480 computer program. 1070 01:02:49,480 --> 01:02:53,030 Typically, RTOL is 0.01. 1071 01:02:53,030 --> 01:02:58,540 RNORM is simply the maximum of the NORM of the 1072 01:02:58,540 --> 01:03:00,920 loads that are applied. 1073 01:03:00,920 --> 01:03:06,780 Notice this Euclidean norm is evaluated as shown down here 1074 01:03:06,780 --> 01:03:09,250 for a vector a. 1075 01:03:09,250 --> 01:03:12,580 The symbolism, two bars each side with a 2 means nothing 1076 01:03:12,580 --> 01:03:16,440 else then taking the squares of all the components of the 1077 01:03:16,440 --> 01:03:19,630 vector, adding those squares up. 1078 01:03:19,630 --> 01:03:24,970 And then, taking the square root out of that addition. 1079 01:03:24,970 --> 01:03:29,200 On displacements, we can also measure convergence. 1080 01:03:29,200 --> 01:03:32,820 And here is the equation that can be used. 1081 01:03:32,820 --> 01:03:34,860 The Euclidean norm on the incremental 1082 01:03:34,860 --> 01:03:37,540 displacement vector. 1083 01:03:37,540 --> 01:03:40,970 And that is compared with DNORM. 1084 01:03:40,970 --> 01:03:43,410 Again, DNORM. 1085 01:03:43,410 --> 01:03:47,870 Notice that this here is again, input by the analyst to 1086 01:03:47,870 --> 01:03:49,900 the computer program. 1087 01:03:49,900 --> 01:03:53,150 It is part of the input that is provided for the analysis. 1088 01:03:53,150 --> 01:03:55,280 And it's a reference displacement. 1089 01:03:55,280 --> 01:03:58,300 Notice then, of course, we should only include here the 1090 01:03:58,300 --> 01:04:01,670 displacements, translational displacements in other words. 1091 01:04:01,670 --> 01:04:04,370 If we have also, rotational degrees of freedom. 1092 01:04:04,370 --> 01:04:08,370 Then we would here include only the rotational degrees of 1093 01:04:08,370 --> 01:04:13,780 freedom components and use here DMNORM. 1094 01:04:13,780 --> 01:04:18,700 This way again, we compare one quantity with certain units 1095 01:04:18,700 --> 01:04:21,580 with a quantity that has the same units. 1096 01:04:21,580 --> 01:04:25,220 Once for translations and once for rotations. 1097 01:04:25,220 --> 01:04:28,440 And so when we have rotational degrees of freedom in the 1098 01:04:28,440 --> 01:04:32,230 analysis, we this way make sure really, that the 1099 01:04:32,230 --> 01:04:35,710 incremental displacements and incremental rotations are 1100 01:04:35,710 --> 01:04:38,140 small at convergence. 1101 01:04:38,140 --> 01:04:41,270 This then completes what I wanted to share with you in 1102 01:04:41,270 --> 01:04:42,010 this lecture. 1103 01:04:42,010 --> 01:04:44,060 Thank you very much for your attention.