Pseudo-Derived the Divergence Theorem

Lecture 5 5.1 Administration • Hand out SV proof. (Jim – Do you mean “SVD proof ”). • Do problem 3 by hand; can verify with MATLAB. Be sure to normalize your eigenvectors. 5.2 5.2.1 Shear strain versus rigid rotation (continuation) Removing mathematical impact of shear strain The last class was ended by showing how a proper choice of axes could remove all apparent shear strain (see figure ??). Here, ∂γ ∂u
, = ∂t ∂y
∂α ∂v
, = ∂t ∂x
where ∂v
∂u
=−

∂x ∂y (5.1) The presence of the negative sign is due to the fact that ∂α ∂t is negative. Thus, the shear strain rate is

 1 ∂u ∂v 1 ∂u ∂u = =0 (5.2) + − 2 ∂y ∂x 2 ∂y ∂y And the rigid rotation rate is

 1 ∂u ∂u ∂u 1 ∂v ∂u = − = − − 2 ∂x ∂y 2 ∂y ∂y ∂y Recall that the 1 2 (5.3) keeps one from double counting. Figure 5.1: (fig:Lec5AxisEvolveOverall) A proper choice of axes can remove all apparent shear strain. 1 Figure 5.2: (fig:Lec5RemoveRotation1) Initial axes such that ∆γ = ∆α. Figure 5.3: (fig:Lec5RemoveRotation2) Initial axes such that ∆γ = ∆α, allowing for the presence of vorticity. 5.2.2 Removing mathematical impact of rigid rotation Can you also remove the mathematical impact of rigid rotation? One would have to pick initial axes so that ∆γ = ∆α (see figure ??). If this can be done, then γ̇ = α̇ ⇒ ∂v ∂u = ∂x ∂y (5.4) Thus, the shear rate is 1 2
∂u ∂v + ∂y ∂x  = 1 2
∂u ∂u + ∂y ∂y  = ∂u ∂y (5.5) And the rigid rotation rate is 1 2
∂v ∂u − ∂x ∂y  1 = 2
∂v ∂v − ∂x ∂x  =0 (5.6) But a perpendicular axis only evolves such that γ̇ = α̇ when there is no vorticity. When there is vorticity, to get γ̇ = α̇ with respect to x − y axis, you could not have a perpendicular set of initial axis as seen in figure ??. The angles would then have to be measured with respect to funky initial axes instead of a more conventional x − y axis. 5.3 Recap • Linear strain rate ≡ • Shear strain rate ≡ ∂u ∂v ∂w ∂x , ∂y , ∂z . 1 2 • Rigid rotation rate ≡  1 2 ∂u ∂y  + ∂v ∂x ∂v ∂x −  . ∂u ∂y  . Note that all these quantities are with respect too perturbations. They say how some small perturbation is stretched and rotated. They say nothing about the solid body translation of the blob moving through space. 2 Figure 5.4: (fig:Lec5Lagranigantranslation) Solid body translation. Figure 5.5: (fig:Lec5Stretching) Example ??. Parcel stretched by 2 in x-direction and contracts by 2 in y-direction. We are in the Lagrangian world, so the solid body translation can be obtained through the ∼ blob’s velocity. To first order: v ≈ ∆x ∆t ⇒ ∆x = v∆t. Think of this as the translation of the center of the blob. The linear and shear strain rates, and the rigid rotation rate tell us about how the shape of the blob changes. That is, the deformation of the blob. 5.3.1 Example A fluid blob has velocity u = 2xi − 2yj. Describe the deformation. Linear strain rate (dilation): ∂u ∂v = 2, = −2 ∂x ∂y (5.7) Thus, there is a stretching in the x-direction of 2 units per unit length per second. There is a contraction in the y-direction of 2 units per unit length per second (see figure ??). Shear strain rate: 1 2
∂u ∂v + ∂y ∂x  = 1 (0 + 0) = 0 2 (5.8) = 1 (0 − 0) = 0 2 (5.9) Rigid rotation rate: 1 2 5.4
∂v ∂u − ∂x ∂y  Velocity Gradient Tensor We can state all this stuff about dilation, deformation and rotation a bit more clearly. 3 Figure 5.6: (fig:Lec5VelGradTen) The velocity gradient tensor, G, acting on a small perturbation, δx, yields the velocity experienced by the perturbation. 5.4.1 Definition  G ≡ velocity gradient tensor = ∂u ∂x ∂v ∂x ∂u ∂y ∂v ∂y  (5.10) If I have some perturbation vector  δx = δx δy  (5.11) Then δu = Gδx  ∂u = ∂x δx ∂v ∂x δx ∂u ∂y δy ∂v ∂y δy  (5.12) (5.13) The velocity gradient tensor, G, provides the velocity experienced by a small perturbation, δx. This is shown schematically in figure ??. 5.4.2 Shear strain and rigid rotation rates in terms of G Define GT ≡ transpose of G (5.14) 1 1 1 1 G + GT + G − GT 2 2 2 2 (5.15) Can express G as G= Now 1 1 G + GT 2 2 1 1 1 r ≡ rotation rate tensor = G − GT 2 2 2 e ≡ strain rate tensor = (5.16) (5.17) (answer to Q4). With the definition of the velocity gradient tensor (equation ??), e (through equation ??) becomes   ⎤ ⎡ ∂u ∂v 1 ∂u + ∂x 2 ∂y ∂x ⎦  e=⎣  (5.18) 1 ∂u ∂v ∂v + 2 ∂y ∂x ∂y 4 Note that this is a symmetric matrix, the diagonal terms are the linear strain rate (dilation), and that the off-diagonal terms are the shear strain rate. Similarly, 12 r (through equations ?? and ??) becomes   ⎤  ⎡  1 ∂v ∂u − 0 − 0 −ω 2 ∂x ∂y 1 ⎦=  r=⎣  (5.19) 1 ∂v 2 ω 0 0 − ∂u 2 ∂x ∂y Here we have an antisymmetric matrix, with the off-diagonal terms the rigid rotation rate. Thus, the impact of linear and shear strain rates come from e δx (5.20) The impact of rigid rotation rate comes from This can also be written as 5.4.3 1 2 1 r δx 2 (5.21) u = x2 yi + 2xy 2 j (5.22) (ω × δx). Example A fluid blob has the velocity Calculate normal and shearing rate of strain at x = −2, y = −1. Determine if the fluid is rotating at that location. Linear and shearing rates of strain:  G=  1 1 e = G + GT = 2 2 2xy 2y 2 x2 4xy 1 2 (4) 1 2 (2) + 12 (4) + 12 (4)   4 8  = 4 2 1 2 (4) 1 2 (8) + 12 (2) + 12 (8) (5.23)   = 4 3 3 8  (5.24) e tells us we have a linear strain rate, or dilation, of 4 units per unit length per second in the x-direction and 8 units per unit length per second in the y-direction. We also have a shear strain rate of 3 rad/sec. Rigid rotation rate:  1 1 1 r = G − GT = 2 2 2 1 2r 5.4.4 1 2 (4) 1 2 (2) − 21 (4) − 21 (4) 1 2 (4) 1 2 (8) − 12 (2) − 12 (8)   = 0 −1 1 0   1 = 2 0 ω −ω 0  (5.25) tells us we have a rotation rate of -1 rad/sec. Examples with pictures See figures ?? through ?? for examples. Jim - Here, Lec5PicExpl1a corresponds to expl 1.1 from your notes, Lec5PicExpl1b corresponds to expl 1.2, and so on and so forth. Let me know if there is any confusion. 5 Figure 5.7: (fig:Lec5PicExpl1a through fig:Lec5PicExpl1f) First picture example of G. Here G is the same as in example ??. The upper left panel shows δx on the unit circle (black circle). The upper right panel shows G δx (the green ellipse) along with δx (now blue). The middle left panel adds e δx (the red ellipse). The middle right panel adds 12 r δx = 12 (ω × δx) (black). The lower left panel show just 12 r δx = 12 (ω × δx). Finally, the lower right panel shows e δx + 12 r δx (cyan). Here the addition of the two components is clearly seen. 6 Figure 5.8: (fig:Lec5PicExpl2a through fig:Lec5PicExpl2f) Second picture example of G. Here there is quite a bit of rotation. G = [−1, 6; 1, −1], e = [−1, 3.5; 3.5, −1], and 12 r = [0, 2.5; −2.5, 0] (Jim – I could not wrestle LATEX into putting arrays in the caption of a figure so I just adopted MATLAB syntax). The upper left panel shows δx (blue). The upper right panel adds G δx (green ellipse). The middle left panel adds e δx (the red ellipse). The middle right panel adds 12 r δx = 12 (ω × δx) (black). The lower left panel show just 12 r δx = 12 (ω × δx) (black) along with δx (blue). Finally, the lower right panel shows e δx + 12 r δx (cyan). Again, the addition of the two components is clearly seen. 7 Figure 5.9: (fig:Lec5PicExpl3a through fig:Lec5PicExpl3c) Third picture example of G. Here there is no rotation. G = [−1, 6; 6, −1], e = [−1, 6; 6, −1], and 12 r = [0, 0; 0, 0]. The upper left panel shows δx (blue). The upper right panel adds G δx (green ellipse). The lower panel shows e δx (the red ellipse) which is the same as G δx. 8 Figure 5.10: (fig:Lec5PicExpl4a through fig:Lec5PicExpl4e) Second picture example of G. Here there is no shear strain. G = [−3, 6; −6, −1], e = [−3, 0; 0, −1], and 12 r = [0, 6; −6, 0]. The upper left panel shows δx (blue). The upper right panel adds G δx (green ellipse). The middle left panel adds e δx (the red ellipse). Note that there is only stretching and not shearing here. The middle right panel adds 21 r δx = 12 (ω × δx) (black). Finally, the lower panel shows e δx + 12 r δx (cyan). Again, the addition of the two components is clearly seen. 9 5.4.5 Rock example (Show pictures of rock) Now a real world example. This rock was taken from a fault zone somewhere in the world. My friends on the single-digit floors of the Green building tell me that as you walk around near the fault you see all this amalgomated rock that has spherical stones embedded in it. I suppose, millions of years ago there was a stream that eroded the stones into spheres, and then, over time, the combination of stones and mud turned to rock. A fault occurred in the area, and applied a huge shear stress to the rock. The result was to turn spherical stones into ellipses. You can see the dominant direction of shear. 5.4.6 Linear uncertainty propagator I also want to give a real world example for the folks on the double digit floors of the Green building – The atmosphere and ocean folks. First, we have to look at the velocity gradient tensor in a different way. G is defined in such a way that δu = G δx (5.26) For a given small displacement vector, you obtain the associated velocity change. Now consider G δx(0) ∆t = δx(∆t) (5.27) Instead of the rate of change, I find the change itself after some ∆t. Define M = G∆t ≡ linear uncertainty propagator (5.28) All the decomposition I did with G, I can do with M 1 1 1 1 M + MT + M − MT 2 2 2 2 1 = ẽ + r̃ 2 M = (5.29) (5.30) Where ẽ ≡ strain tensor, 12 r̃ ≡ rotation tensor. If a circle gets converted to an ellipse (see figure ??), what we really want to know are the principal axes of the ellipse, the initial axes that evolve into the principal axes, and the (Jim This got cut off in the copying, bottom of page 13). How do we determine these principal axes? The first thing we want is the major axis of the ellipse. A growth factor can be obtained for an arbitrary initial vector through σ2 = δx(∆t)T δx(∆t) δx(0)T δx(0) (5.31) The numerator is the squared length of the final-time vector. The denominator is the squared length of the initial-time vector. To get the major axis of the final time ellipse we need to maximize σ 2 : the ratio of the squared length of the major axis of the final-time ellipse to the squared length of the initial-time vector is the largest value σ 2 can take. We know that δx(∆t) = M δx(0), so we can substitute this into equation ?? yielding σ2 = δx(0)T MT M δx(0) δx(0)T δx(0) 10 (5.32) Figure 5.11: (fig:Lec5CircleConversion, fig:Lec5StrictCircleConversion) The left panel is a circle changing to an ellipse after being acted on by M. Note that there is no translation so strictly speaking, the right panel show the initial circle and the final ellipse in their actual final positions. Figure 5.12: (fig:Lec5vEvolution) eig(MT M) gives V information. (Jim – You might want some more comments here). The short answer is that σ 2 is maximized for δx(0) when δx(0)σ 2 = MT M δx(0) (5.33) I’ve handed around a simple proof that was done by two students in 12.800 last year – your TA was one of them. Note that MT M is symmetric positive definate – This means that the eigevectors are orthogonal. What about the directions at final-time (in a manner more interesting than δx(∆t) = M δx(0)!). Take δx(0)σ 2 = MT M δx(0) and pre-multiply both sides by M ⇒ M δx(0)σ 2 δx(∆t)σ 2 = MMT M δx(0) = MMT δx(∆t) (5.34) The directions of the ellipse’s principal axes are given by the eigenvectors of MMT . Note that the eigenvalues of MT M are the same as those for MMT . This is all seen schematically in figures ?? and ?? and again in another manner in figures ?? through ??. There is an alternative approach that can provide the same information: singular value decomT T position (SVD). I’m not going to prove √ this, but SVD(A) = UΣV where U = eig(AA ), V = eig(AT A), and Σ us diagonal and = D where D ≡ eigval(AAT ). If I relate A’s to M’s: • SVD(M) = UΣVT • U = eig(MMT ) ⇒ directions of principle axes of final time ellipse. • V = eig(MT M) ⇒ directions of initial-time circle that evolves into the principle axes of the ellipse. 11 Figure 5.13: (fig:Lec5uEvolution) eig(MMT ) gives V information. (Jim – You might want some more comments here). Figure 5.14: (fig:Lec5MInitial) M = [−0.4945, −.5752; −2.2910, 1.2716]. Blue is δx(0) and red is δx(∆t) = M δx(0). Figure 5.15: (fig:Lec5MFirstStep) M is the same from figure ??. Blue is again δx(0), green is M δx(0), and red is MT M δx(0). The eigenvectors of red give V, and the eigenvalues give Σ2 . Figure 5.16: (fig:Lec5MFinalStep) M is the same from figure ??. Blue is again δx(0), red is MT δx(∆t), and green is MMT δx(0). The eigenvectors of green give U, and the eigenvalues give Σ2 . 12 Figure 5.17: (fig:Lec5SVD) Initial axes (blue) and principle axes (red) along with their growth factors. • Σ are the factors by which the δx(0) grow. This is all seen in figure ??. 5.5 Reading for class 6 KC01: 4.1 - 4.4 KFF: V 13