Heat Flow via Crank–Nicolson An Improved Leap Frog Rubin H Landau Sally Haerer, Producer-Director Based on A Survey of Computational Physics by Landau, Páez, & Bordeianu with Support from the National Science Foundation Course: Computational Physics II 1 / 1
Problem: How Does a Bar Cool? 100 K Insulated Metallic Bar Touching Ice Aluminum bar, L = 1 m, w along x Insulated along length, ends in ice ( T = 0 C) Initially T = 100 C How does temperature vary in space and time? 2 / 1
CN Improved Algorithm: Split t Step ⇒ ∂ t FD → CD Problem: Making the First Step Knowing Only T ( x , t = 0 ) Split time step for ∂ t at t + ∆ t / 2 ∂ T � x , t + ∆ t � ≃ T ( x , t + ∆ t ) − T ( x , t ) + O (∆ t 2 ) (1) ∂ t 2 ∆ t Bad as FD for t = t + ∆ t , Good at t + ∆ t / 2 CD ∂ 2 x at t = t + ∆ t / 2 ∂ 2 T � � x , t + ∆ t 1 (2) ≃ 2 (∆ x ) 2 × ∂ x 2 2 [ T ( x − ∆ x , t + ∆ t ) + T ( x + ∆ x , t + ∆ t ) − 2 T ( x , t + ∆ t ) � + T ( x − ∆ x , t ) + T ( x + ∆ x , t ) − 2 T ( x , t ) + O (∆ x 2 ) 3 / 1
Split-Time Discrete Form of Heat Equation Derivatives in Terms of Differences ∂ 2 T ( x , t ) ∂ T ( x , t ) = K (3) ∂ t C ρ ∂ x 2 T i , j + 1 − T i , j = η 2 [ T i − 1 , j + 1 − 2 T i , j + 1 + T i + 1 , j + 1 + T i − 1 , j − 2 T i , j + T i + 1 , j ] (4) K ∆ t x = i ∆ x , t = j ∆ t , η = (5) C ρ ∆ x 2 • Future in terms of present (grid values only, yet mid time): � 2 � 2 � � − T i − 1 , j + 1 + η + 2 T i , j + 1 − T i + 1 , j + 1 = T i − 1 , j + η − 2 T i , j + T i + 1 , j (6) • Implicit solution: solve all lattice sites simultaneously • Leap frog = explicit , permits single time step 4 / 1
Arrange Discrete Heat Equation as Matrix Equation � 2 � 2 T 0 , j + 1 + T 0 , j + η − 2 � T 1 , j + T 2 , j � η + 2 − 1 T 1 , j + 1 � 2 � 2 � T 2 , j + 1 T 1 , j + η − 2 � T 2 , j + T 3 , j − 1 η + 2 − 1 = � 2 � 2 T 3 , j + 1 � − 1 η + 2 − 1 T 2 , j + η − 2 � T 3 , j + T 4 , j . ... ... ... . . . . . Advance t , repeat step LHS future T ’s = j + 1 C-N = ↑ precise, stable, work RHS present T ’s = j vonN stability: all ∆ t , ∆ x OK: RHS ends @ future, OK as BC IC: j = 0, hot bar, cold ends ξ ( k ) = 1 − 2 η sin 2 ( k ∆ x / 2 ) 1 + 2 η sin 2 ( k ∆ x / 2 ) Solve matrix → j = 1 all x =i 5 / 1
Special Solution of Tridiagonal Matrix Equations ⊙ x 1 b 1 d 1 c 1 0 0 0 · · · · · · · · · x 2 b 2 a 2 d 2 c 2 0 0 · · · · · · · · · x 3 b 3 0 a 3 d 3 c 3 0 · · · · · · · · · = ... ... · · · · · · · · · · · · · · · · · · · · · · · · 0 0 0 0 a N − 1 d N − 1 c N − 1 · · · x N − 1 b N − 1 0 0 0 0 0 a N d N · · · x N b N A i , j ⇒ N 2 words, access Libes good for stnd linear eq [ A ] x = b Tridiagonal ⇒ only 3 vectors Yet [ A ] = tridiagonal { d i } i = 1 , N , { c i } i = 1 , N , { a i } i = 1 , N ⇒ ∃ more robust, faster solution Single subscript ⇒ 3N − 2 6 / 1
Soltn: Coef matrix → upper triangular , diagonals =1 Start: divide 1st eqtn by d 1 , subtract a 2 × 1st eqtn: c 1 b 1 1 0 0 0 d 1 · · · · · · · · · x 1 d 1 d 2 − a 2 c 1 b 2 − a 2 b 1 x 2 = 0 c 2 0 0 · · · · · · · · · d 1 d 1 x 3 b 3 0 a 3 d 3 c 3 0 · · · · · · · · · Next: divide 2nd eqtn by 2nd diagonal element c 1 1 0 0 0 · · · · · · · · · b 1 d 1 d 1 x 1 c 2 0 1 0 0 b 1 c 1 · · · · · · b 2 − a 2 d 2 − a 2 x 2 = a 1 d 1 c 1 0 a 3 d 3 c 3 0 x 3 d 2 − a 2 · · · · · · · · · d 1 b 3 · · · · · · · · · · · · · · · · · · · · · · · · Repeat steps → upper triangular form 1 h 1 0 0 0 x 1 p 1 · · · 0 1 h 2 0 0 x 2 p 2 = · · · 0 0 1 h 3 0 x 3 p 3 · · · Back substitution ⇒ explicit solution: x i = p i − h i x i − 1 ; i = n − 1 , n − 2 , . . . , 1 , x N = p N 7 / 1
Crank–Nicolson Implementation HeatCNTridiag.py Solve C-N linear equations using libe, esp for tridiagonal 1 Check stability for ∆ x and ∆ t 2 Contoured surface plot of T ( x , t ) 3 Compare precision and speed: leap-frog vs 4 Crank–Nicolson Assume stable, very small ∆ t = accurate 5 8 / 1
Recommend
More recommend
