Giáo trình Các phương pháp tính truyền nhiệt

nhau. Do đó, trong thời gian dτ, số phân tử ở T1 và T2 qua dS là nh− nhau, bằng: z x T2T1 λλ y O H1. Để chứng minh định luật Fourier d2n = 6 1 no ω dS dτ * L−ợng động năng qua dS từ T1 và T2 là: d2E1 = E 1 d 2n = 6 1 no ω dS dτ 2 i kT1 d2E2 = E 2 d 2n = 6 1 no ω dS dτ 2 i kT2 Trừ hai ph−ơng trình cho nhau, ta đ−ợc: δ2Q = ( E 1 - E 2)d2n = 6 1 no ω dSdτ 2 ik (T1 - T2) Vì T1 - T2 = - ⎟⎠ ⎞⎜⎝ ⎛ dx dT . 2 λ nên δ2Q = - 6 i no k ϖ λ dx dT dS dτ Do 6 i no k = 6 i no N R = 3 1 (no N à ) ( à2 iR ) = 3 1 ρco nên 4 δ2Q = - ( 3 1 ρco ω λ ) dx dT dS dτ = - λ dx dT dS dτ hay τ δ dSd Q2 = q = - λ ⎟⎠ ⎞⎜⎝ ⎛ ∂ ∂ x T * Khi dS có vị trí bất kỳ thì q = - λ gradT hay dạng vectơ dòng nhiệt là qr = - λ dTagrr 1.1.2. Phát biểu: Vectơ dòng nhiệt tỷ lệ thuận với gradient nhiệt độ: Biểu thức vectơ: q r = - λ dTagrr Dạng vô h−ớng: q = - λgradT, [W/m2]; δQ = - λgradT.dS, [W] 1.1.3. Hệ số dẫn nhiệt Hệ số dẫn nhiệt là hệ số của định luật Fourier: λ = |q/gradT| [W/mK] Theo chứng minh trên ta có: λ = 3 1 ρ ω λ cv = 3 1 ⎟⎠ ⎞⎜⎝ ⎛ RT p ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ πm kT8 ⎟⎟⎠ ⎞ ⎜⎜⎝ ⎛ π pd.2 kT 2 Cv = 3 2 m Tk d c 3 3 2 v π cho thấy: λ không phụ thuộc p, và λ↑ khi T↑ hoặc cv↑ hoặc đ−ờng kính d cùng khối l−ợng phân tử m giảm. Định luật Fourier đúng cho mọi chất rắn, lỏng, khí. 1.2. Ph−ơng trình vi phân dẫn nhiệt 1.2.1. Định nghĩa: Ph−ơng trình vi phân dẫn nhiệt là ph−ơng trình cân bằng nhiệt cho một phân tố dv bên trong vật. 1.2.2. Thiết lập Luật cân bằng nhiệt cho dV ∈ V là: H2. CBN cho dV z x y qλ qω qω λ ρ C dV V O 5 [L−ợng nhiệt phát sinh trong dV] - [Thông l−ợng nhiệt qua dV]= [Biến thiên entanpy của dV] Cho tr−ớc (qv, ρ, cp, λ) ∈ dV, có thể viết ph−ơng trình trên ở dạng: qvdVdτ - divqr dVdτ = ρdV.cp τ∂ ∂t dτ hay τ∂ ∂t = p v c q ρ - pc 1 ρ div q r , trong đó dòng nhiệt qua dV là: qr = qr λ + q r ω = - λ dtagrr + ρωr cpt, do đó: divqr = div (ρcp ωr t- λ dtagrr ), coi (ρ, cp) = const ta có : divqr = ρcp div (tωr ) - div (λ dtagrr ) = ρcp (tdiv ωr + ωr dtagrr ) - λdiv ( dtagrr )- dtagrr . λdagrr = ρcp (tdiv ωr + ωr dtagrr ) - λ∇2t - dtagrr . λdagrr Vậy ph−ơng trình có dạng: τ∂ ∂t = pc qv ρ - tdiv ω r - ωr . dtagrr + ρ λ pc ∇2t + ( dtagrr dagrr λ)/ρcp do τ∂ ∂t + ωr . dtagrr = τ∂ ∂t + dx dt . τd dx + dy dt . τd dy + dz dt . τd dz = τd dt nên ph−ơng trình vi phân dẫn nhiệt sau khi đặt a = Cpρ λ , sẽ là: τd dt = a∇2t + p v c q ρ + pc 1 ρ dtagr r dtagr r λ) - tdivωr , với: là tích vô h−ớng của 2 vectơ và ∇2t = ∆t là toán tử Laplace của nhiệt độ, có dạng: ∇2t = ⎪⎪ ⎪⎪ ⎩ ⎪⎪ ⎪⎪ ⎨ ⎧ ∂ ∂+∂ ∂+∂ ∂+∂ ∂+∂ ∂ ∂ ∂+∂ ∂⋅+∂ ∂⋅+∂ ∂ ∂ ∂+∂ ∂+∂ ∂ ),,( sinsin cos2 ),,(11 ),,( 222 2 222 2 2 2 2 2 2 2 22 2 2 2 2 2 2 2 ϕθϕθθθ θ θ ϕϕ rtrong r tt rr t r t rr t zrtrong z tt rr t rr t zyxtọatrong z t y t x t 1.2.3. Các dạng đặc biệt của ph−ơng trình vi phân dẫn nhiệt * Với vật rắn, ωr = 0, ph−ơng trình có dạng: dtagrr . λdagrr dtagrr dagrr λ, (trong tọa độ vuông góc (xyz)) (trong tọa độ trụ (rϕz)) (trong tọa độ cầu (rθϕ)) 6 τ∂ ∂t = a∇2t + p v c q ρ + pc 1 ρ dtagr r . λdagrr * Vật rắn có λ = const ∀xyz ph−ơng trình là: τ∂ ∂t = a∇2t + p v c q ρ * Vật rắn có λ = const , ổn định nhiệt τ∂ ∂t = 0, ph−ơng trình là: a∇2t + λ vq = 0. Nếu không có nguồn nhiệt, qv = 0, thì ∇2t = 0. 1.3. Các điều kiện đơn trị (ĐKĐT) 1.3.1. Định nghĩa: ĐKĐT là những điều kiện cho tr−ớc nhằm xác định duy nhất nghiệm của một hệ ph−ơng trình. 1.3.2. Phân loại các ĐTĐT: Theo nội dung, các ĐKĐT đ−ợc phân ra 4 loại sau: 1. Điều kiện hình học: Cho biết mọi thông số hình học đủ để xác định hình dạng, kích th−ớc, vị trí của hệ. 2. Điều kiện vật lý: Cho biết luật phân bố các thông số vật lý theo nhiệt độ t tại ∀M∈ hệ; tức cho luật xác định (ρ, cp, λ, a...) = f(t, M∈V). 3. Điều kiện ban đầu: Cho biết luật phân bố nhiệt độ lúc τ = 0 tại mọi điểm M ∈ hệ, tức cho biết t = t(x, y, z, τ = 0), ∀(x, y, z) ∈ V. 4. Điều kiện biên: Cho biết luật phân bố nhiệt độ hoặc luật cân bằng nhiệt tại mọi điểm trên biên W, ở mọi thời điểm τ, tức cho biết: t = t(M, τ) hoặc dtagrr = f(M, τ, t) ∀M (x, y, z) ∈ V ∀τ ∈ ∆τ xét 1.3.3. Các loại điều kiện biên (ĐKB) Tại mỗi miền Wi của mặt biên kín W = ∑Wi, tuỳ theo cách phân bố t hoặc cách trao đổi nhiệt, ta có thể cho biết các loại ĐKB sau đây: 1. ĐKB loại 1: Cho biết luật phân bố nhiệt độ t tại mọi điểm M1 ∈ W1 ở mọi thời điểm: 7 t = t (M1, τ), ∀M1∈ W1, ∀τ 2. ĐKB loại 2: Cho biết dòng nhiệt dẫn qua biên: q (M2,τ) = -λ n t ∂ ∂ , tức cho biết n t ∂ ∂ = λ −1 q (M2, τ), ∀M2 ∈ W2, ∀τ. Khi n t ∂ ∂ = q = 0 tức biên W2 đ−ợc cách nhiệt tuyệt đối hoặc là biên đối xứng, lúc này t đạt cực trị tại W2, và đ−ờng cong t(M) có tiếp tuyến nằm ngang. 3. ĐKB loại 3: Cho biết biên W3 tiếp xúc chất lỏng có tf, α và toả nhiệt ra chất lỏng theo luật: -λtn (M3, τ) = α[t(M3, τ) - tf], tức cho biết gradt (M3) = [tf - t(M3)]/(λ/α), ∀M3 ∈ W3, ∀τ. 4. ĐKB loại 4: Cho biết luật CBN khi biên W4 tiếp xúc vật rắn khác, có nhiệt độ t4 và λ4, tại M4 ∈ W4, ph−ơng trình cân bằng nhiệt có dạng : -λ n )M(t 4 ∂ ∂ = λ4 n )M(t 44 ∂ ∂ và t(M4) = t4 (M4) 5. ĐKB loại 5: Cho biết luật cân bằng nhiệt trên biên W5 di động, t x H3. CBN trên biên W5 do có sự chuyển pha, trao đổi chất (khối l−ợng thay đổi) hoặc đang biến dạng: -λ n )M(t 5 ∂ ∂ = rcρ τd dx 5 - λ' n 't ∂ ∂ (M5), với rc = nhiệt chuyển pha; τd dx5 = vận tốc biên W5; ρ: khối l−ợng riêng pha mới. 1.3.4. ý nghĩa hình học của các loại ĐKB Dạng đ−ờng cong phân bố nhiệt độ t(x, y, z, τ) tại lân cận biên W, -λ t x ∂ ∂ -λ t ' n ∂ ∂ -rcρ 5dxdτ 5dx dτ x 0 x5 8 tuỳ theo cách cho ĐKB, sẽ có các đặc điểm hình học sau đây: W Cách cho ĐKB Đ−ờng cong t(M,τ) ý nghĩa hình học 1 tw = const x w M V t o t(M) đi qua một điểm cố định Mo ∈W 2 n t w ∂ ∂ = 0 x q = 0 V t = 0β t(M) đạt cực trị trên W cách nhiệt n t w ∂ ∂ = const x β V Các tiếp tuyến của t(M) tại W song song, góc β = const 3 n t w ∂ ∂ = αλ − / tt wf xV tf R λα Các tiếp tuyến của t(M) tại W3 qua điểm R( α λ , tf) 4 n t w ∂ ∂ = λ λ 4 x t ow ∂ ∂ tW = t4W xV Voγ t(M) liên tục, không khả vi tại W4 và γ = const 5 -λ n t w ∂ ∂ = re ρ τd dx5 - λ n 't w δ δ xV W5 di chuyển với tốc độ ω = τd dx5 H4. Minh hoạ ý nghĩa hình học các ĐKB 1.4. Mô hình một bài toán dẫn nhiệt Mô hình toán học của một bài toán dẫn nhiệt là một hệ ph−ơng W = const W W 5dx dτ 9 trình vi phân (t), gồm ph−ơng trình vi phân DN và các ph−ơng trình mô tả các ĐKĐT nh− sau: (t) = τ∂ ∂t = a∇2t + c vq ρ và các ph−ơng trình mô tả các ĐKĐT. Mục đích chính của truyền nhiệt là tìm các ph−ơng pháp giải hệ (t) để tìm hàm phân bố t(x,y,z,τ) thoả mãn hệ (t). ρ, λ,c, qv −λo to n −λ t n W4 W3 W2 W1 W5t (M, )w1 τ −λ t' n' cf dx τdr −λ t n w x-1q(M, )τ 2 t = a t +∇2 qvρc α [tw - tf] 2 −λ t x w M ∂ ∂ ∂∂τ ∂∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ H5. Mô hình 1 bài toán DN 10 Ch−ơng 2: các Ph−ơng pháp giải tích 2.1. phép chuẩn hoá và định lý hợp nghiệm: 2.1.1. Nội dung cơ bản của các ph−ơng pháp giải tích ý t−ởng của Fourier là chuyển ph−ơng trình đạo hàm riêng tuyến tính thành một số ph−ơng trình vi phân th−ờng t−ơng đ−ơng, bằng cách tách biến, tìm nghiệm riêng ổn định và biến thiên hằng số. Các cách trên đ−ợc sử dụng tuỳ thuộc tính thuần nhất hay không thuần nhất của ph−ơng trình dẫn nhiệt và ph−ơng trình vi phân mô tả các điều kiện biên. 2.1.2. Ph−ơng trình vi phân thuần nhất và không TN - Định nghĩa: Ph−ơng trình vi phân F(t, tx, txx) = 0 đ−ợc gọi là thuần nhất khi: nếu t là nghiệm của ph−ơng trình thì ct, ∀c =const, cũng là nghiệm của F(t, tx, txx) = 0. - Ví dụ: tτ = atxx, tx(0,τ) = - λ α t(0,τ) là TN tτ = a∇2t + c vq ρ , tx (L, τ) = −α λ [t(L, τ) - tf] là không TN Nhận xét: Ph−ơng trình truyền nhiệt không chứa số hạng tự do, nh− qv và tf, là ph−ơng trình thuần nhất. 2.1.3. Nguyên lý hợp nghiệm Nếu các ti,∀i = 1ữn, là nghiệm riêng của bài toán biên thuần nhất (tức ph−ơng trình vi phân và các ĐKB thuần nhất), thì t = ∑ = n 1i ii tC cũng là nghiệm của bài toán TN đó, ∀Ci = const 2.1.4. Phép chuẩn hoá - Định nghĩa: Phép chuẩn hoá một hệ ph−ơng trình là cách đổi các biến và thông số có thứ nguyên thành các biến và thông số không thứ nguyên. - Lợi ích của phép chuẩn hoá là đơn giản hệ ph−ơng trình và cách 11 giải, khiến cho nghiệm có tính tổng quát, không phụ thuộc các đại l−ợng có thứ nguyên, và trong vài tr−ờng hợp, có thể thuần nhát hoá các điều kiện biên không thuần nhất. - Ví dụ: Bài toán làm nguội tấm phẳng với 2 biên Wo/W3 có mô hình: (t) ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ −τδλ α−=τδ =τ = ∂ ∂=τ∂ ∂ )W()TN0(]t),(t[),(t )W()TN(0),0(t )DKD(t)0,x(t )FT()TN( x tat 3fx ox o 2 2 Đổi biến ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ δ τ= δ= − −=θ 2 fo f aF xX tt tt và đặt B = λ αδ thì do τ∂ ∂t = θ∂ ∂t . F∂ θ∂ . τ∂ ∂F = (to - tf) 2 a δ . F∂ θ∂ x t ∂ ∂ = θ∂ ∂t . X∂ ∂θ . x X ∂ ∂ = δ − fo tt . X∂ θ∂ 2 2 x t ∂ ∂ = x∂ ∂ ( x t ∂ ∂ ) = X∂ ∂ . ( δ − fo tt . X∂ θ∂ ) x X ∂ ∂ = 2 fo tt δ − 2 2 X∂ ∂ θ tx (δ, τ) = δ − fo tt θx (1, F) = λ α− [t (δ, τ) - tf] có dạng TN là θx (l, F) = λ αδ− θ [1, F] = Bθ(1,F) τ∂ ∂t = (to- tf) 2 a δ F∂ θ∂ = a 2 2 x t ∂ ∂ = a 2 fo tt δ − . 2 2 X∂ ∂ θ có dạng đơn giản hơn là F∂ ∂θ = 2 2 X∂ ∂ θ . Khi đó bài toán (t) đ−ợc chuyển đổi thành bài toán không thứ nguyên (θ) t−ơng đ−ơng, có dạng chuẩn hoá là: 12 (θ) ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ −= = = ∂ ∂=∂ ∂ )(),1(),1( )(0),0( 1)0,( 2 2 TNFBF TNF X XF x x θθ θ θ θθ Bài toán (θ) có hai điều kiện biên ở dạng thuần nhất. 2.2. Ph−ơng pháp tách biến Fourier 2.2.1. Nội dung ph−ơng pháp Fourier Là tìm nghiệm ở dạng tách biến, nh− là tích của một hàm của tọa độ với một hàm của thời gian. Nhờ đó có thể chuyển một ph−ơng trình đạo hàm riêng thành hệ hai ph−ơng trình vi phân th−ờng t−ơng đ−ơng. Ph−ơng pháp này th−ờng dùng để giải các hệ ph−ơng trình thuần nhất. 2.2.2. Cách giải các bài toán thuần nhất Các bài toán thuần nhất có thể giải bằng ph−ơng pháp tách biến Fourier theo các b−ớc: tách biến ph−ơng trình vi phân DN tìm nghiệm tổng quát, xác định các nghiệm riêng theo các ĐKĐT, hợp nghiệm. Đó là các b−ớc của ph−ơng pháp tách biến. 2.2.3. Ví dụ: Bài toán làm nguội tấm phẳng biên (W2+W3) 1. Phát biểu bài toán: Cho vách phẳng có δ, a, λ, to = t(x,0) cách nhiệt tại x = 0, toả nhiệt tại x = δ ra môi tr−ờng tf, α. Tìm tr−ờng t (x, τ) 2. Mô hình TH: (t) ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ −τδλ α−=τδ =τ = =τ ]t),(t[),(t 0),0(t t)0,x(t att fx x o xx Bằng cách đổi biến: O t x a t(x, ) λ τq = 0 W2 δ W3 α tf to H6. Bài toán (2.2.2) 13 θ = fo f tt tt − − , X = δ x , F = 2 a δ τ , B = λ αδ sẽ thu đ−ợc hệ ph−ơng trình (θ) t−ơng đ−ơng, ở dạng chuẩn hoá: (θ) ⎪⎪⎩ ⎪⎪⎨ ⎧ −= = = = ),1(),1( 0),0( 1)0,( FBF F x x x xxF θθ θ θ θθ 3. Tách biến bằng cách tìm nghiệm dạng θ(X,F) = X(x) F(F). Thay vào θF=θxx có X(x) F'(F) = X"(X) F(F) hay )X(X )X("X = )F(F )F("F = -k2 (do 2 hàm độc lập), chuyển thành 2 ph−ơng trình vi phân th−ờng: ⎪⎩ ⎪⎨⎧ =→=+ +=→=+ − Fk2 21 2 2 e)F(F0)F(Fk)F('F kXcosckXsinc)X(X0)X(Xk)X("X Nghiệm tổng quát là θ(X,F) = (c1sin kX + c2coskX) F2ke− 4. Xác định các hằng số theo ĐKĐT θx(0,F) = 0 → (kc1cos0 + (-kc2sin0) F2ke− = 0 → c1 = 0 và θ (X,F) = c2 coskX F2ke− θx(1,F) = (-kc2sin0) F2ke− = -Bθ (1,F)= -Bc2 cosk F2ke− → ksin kcos = cotgk = B k , ph−ơng trình này có vô số nghiệm ki, i = 1 ữ n. Các nghiệm riêng thoả mãn ĐKB có dạng: θi(X,F) = c2coskiX. F2ike− , O cotgk k k k k k k1 2 3 4 5 π 2π 3π 4π k B H7. Giải ph−ơng trình cotg k = B k nghiệm hợp là θ(X,F) = ∑∞ = − 1i F2ik ii Xekcosc - Điều kiện đầu θ(X,0) = 1 → ∑cicoskiX=1 → coski X∑∞ 1 ii Xkcosc = coskiX → ∫1 0 iXdXkcos = i i k ksin = ∫ ∑1 0 iii dXXkcoscXkcos = 14 ci ∫ 1 0 i 2 XdXkcos = ci i ii k kk 4 2sin2 + → ci = ii i k2sink2 ksin4 + Vậy nghiệm bài toán là: θ(X,F) = 4 ∑∞ = +1i ii i k2sink2 ksin cos(kiX) F2ike− * Đồ thị θ(X,F) và t(x, τ) có dạng: O xF =∞ 6 5 4 3 2 1 1 F=0 1 O t x tf δ to 5 4 3 2 = 0τ =τ ∞ R H8. Phân bố θ(X,F) H9. Phân bố t(x, τ) 2.3. Ph−ơng pháp nghiệm riêng ổn định 2.3.1. Phạm vi sử dụng ph−ơng pháp NROĐ Để giải các bài toán không thuần nhất có nghiệm riêng khi ổn định, tức là khi tτ = θF = 0 2.3.2. Nội dung ph−ơng pháp NROĐ Gồm các b−ớc sau: 1. Tìm nghiệm riêng ổn định θ (x) của bài toán (θ), ứng với lúc ổn định, theo ph−ơng trình θ F = 0 = θxx 2. Thay (v = θ - θ ) vào bài toán (θ) để lập bài toán (v), sẽ đ−ợc bài toán (v) thuần nhất. 3. Tìm nghiệm v của bài toán (v) bằng ph−ơng pháp tách biến, sau đó lập nghiệm của bài toán (θ) đã cho là θ = θ + v 2.3.3. Ví dụ: Bài toán gia nhiệt vách phẳng biên (W1) 1. Phát biểu BT: Cho vách phẳng có δ, a, λ, t(x,0) = to = t(δ, τ) và F = 0 θ 15 t(0, τ) = 2to. Tìm t(x, τ) * Mô hình TH: (t) ⎪⎪⎩ ⎪⎪⎨ ⎧ =τδ =τ = =τ o o o xx t),(t t2),0(t t)0,x(t att Chuẩn hoá bằng cách đặt ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ δ τ= δ= −=θ 2 o o aF xX t tt O t x a, t(x, ) λ τ to δ W1 2to W'1 H10. Bài toán (2.3.3) bài toán (t) trở thành dạng chuẩn hoá (θ) nh− sau: (θ) F xx (X,0) 0 (0,F) 1 (1,F) 0 (0TN) θ = θ⎧⎪θ =⎪⎨θ =⎪⎪θ =⎩ . Ta sẽ giải bài toán (θ) không thuần nhất này bằng ph−ơng pháp NROĐ 2. Tìm nghiệm riêng θ của bài toán ổn định: ( θ ) ⎪⎩ ⎪⎨ ⎧ ==θ +==θ +=θ→θ==θ 2 21 21xx c1)0( cc0)1( cXc0 → θ X1 −= 3. Thay v(X,F) = θ(X,F) - θ (X) = θ(X,F) + X - 1 vào (θ): bài toán (θ) trở thành bài toán (v) thuần nhất nh− sau: (v) ⎪⎪⎩ ⎪⎪⎨ ⎧ =−=θ−θ= =−=θ−θ= −=θ−θ= =θ−=θ=θ= )TN(011)0()F,0()F,0(v 000)1()F,1()F,1(v 1X)X()0,X()0,x(v vvv xxxxxxxxFF 4. Tìm nghiệm bài toán (v) bằng ph−ơng pháp tách biến, t−ơng tự nh− bài toán 2.2.2: - Tách biến ph−ơng trình vF = vxx có nghiệm tổng quát là: 16 v(X,F) = X(X)F(F) = (c1sinkx + c2coskx) F2k e − - Theo ĐKB: v(0,F) = 0 → c2 F 2k e − = 0 → c2 = 0 → v(X,F) = c1sinkX F 2k e − Theo v(1,F) = 0 ⇒ c1sink F 2k e − = 0 → sin k = 0 → k = nπ → v(X,F) = ∑ π∞ =1n n )Xnsin(c F 2)n(e π - Theo ĐKĐ: v(X,0) = X-1 → X - 1 = ∑ π∞ =1x n )Xnsin(c → ∫ π− 1 0 dX)Xnsin()1X( = ∫ π 1 0 )Xnsin( dX)Xnsin(c 1n n π∑ ∞ = → - πn 1 = 2 cn → cn = π − n 2 → nghiệm ph−ơng trình (v) là: v(X,F) = - π 2 ∑ π n )Xnsin( F2)n(e π . Do đó, nghiệm bài toán (θ) đã cho là: θ(X,F) = θ (X) + σ(X,F) θ(X,F) = (1-X) - π 2 ∑ π∞ =1n n )Xnsin( exp (-n2π2F) * Phân bố nhiệt độ θ(X,F) và t(x,τ) có dạng: O 1 x θ θ = 1- x F = ∞ F = 0 1 1 2 O t x δ t t = 2t - t = 0 o τ 2to o o 1 2 δ x/ H11. Phân bố θ(X,F) H12. Phân bố t(x, τ) 2.4. Ph−ơng pháp biến thiên hằng số 2.4.1. Phạm vi sử dụng: Ph−ơng pháp biến thiên hằng số thời gian An(F) đ−ợc sử dụng khi: - Bài toán (θ) không tồn tại nghiệm riêng ổn định - hoặc có nghiệm riêng ổn định θ nh−ng không tìm đ−ợc - Bài toán với vật có nguồn nhiệt trong, hoặc đ−ợc gia nhiệt bằng điện. 17 2.4.2. Nội dung ph−ơng pháp BTHS Gồm các b−ớc sau: 1. Lập bài toán (v) thuần nhất, bằng cách cho bằng 0 tất cả các ĐKB không thuần nhất trong bài toán (θ). 2. Tách biến v(X,F) = X(X).F(F) và tìm X(X) thoả mãn các ĐK biên thuần nhất, sẽ đ−ợc các nghiệm riêng dạng Xn(X) = cφn(X), trong đó φn(X) = f(n,X) là hàm số riêng, thoả mãn điều kiện trực giao: ∫ φφ 1 0 mn dX)X()X( = ⎩⎨ ⎧ = ≠ nmkhic nmkhi0 3. Biểu diễn nghiệm bài toán (θ) ở dạng θ(X,F) = ∑ φ∞ =1n nn )X()F(A và biến thiên hằng số thời gian An(F), tức tìm biểu thức xác định An(F) nhờ điều kiện trực giao của φn(X): ∫ φθ 1 0 m dX)X()F,X( = ∑∞=1n n )F(A ∫ φφ 1 0 mn dX)X()X( = cAn(F) tức có quan hệ An(F) = c 1 ∫ φθ 1 0 n dX)X()F,x( 4. Lập hệ ph−ơng trình th−ờng của An(F) bằng cách tính dF d An(F), tìm nghiệm An(F) thoả mãn điều kiện ban đầu. 5. Viết nghiệm bài toán (θ) ở dạng θ(X,F) = )F(A)X( n 1n n∑φ ∞ = 2.4.3. Bài toán tấm phẳng biên (W2 + W20) 1. Phát biểu: Cho vách phẳng có δ, a, λ, t(x,0) = to, tx (δ, τ) = 0 và tx (0, τ) = - δ ot . Tìm t(x, τ) O t xa, t( ,0) λ q = t =o δ q = 0t toλ δ (t = - )toδx x = 0 x t 18 * Mô hình TH: (t) ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ =τδ δ−=τ = =τ 0),(t t),0(t t)0,x(t att x o x o xx chuẩn hoá với θ = o o t tt − , X = δ x , F = 2 a δ τ , sẽ có: (θ) ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ = −== = = 0)0,( )0(1),0(),0( )(0),1( X TNt t F TNF x o x x xxF θ τδθ θ θθ 2. Giải ph−ơng trình 0TN (θ) bằng ph−ơng pháp BTHS: 1) Lập bài toán (v) thuần nhất từ (θ): (v) ⎪⎪⎩ ⎪⎪⎨ ⎧ = = = = 0)0,X(v 0)F,0(v 0)F,1(v vv x x xxF (TN) 2) Tìm nghiệm riêng bài toán biên, vx (1,F) = vx(0,F) = 0, bằng cách tách biến v(X,F) = X(x)F(F) có X(x) = c1 sin kX + c2 cos kX ⎩⎨ ⎧ π=→−==→= =→==→= nkksinkc0)1(X0)F,1(v kXcosc)x(Xc0)0(X0)F,0(v 2xx 21xx Do đó có X(X) = cncos (nπX) và hàm số riêng là φn(X) = cos(nπX). 3) Để θ(X,F) = )X()F(A n 1n n φ∑ ∞ = = )Xncos()F(A 1n n π∑ ∞ = là nghiệm bài toán (θ) thì hằng thời gian An(F) phải xác định theo điều kiện trực giao của hàm riêng φn(X)= cos (nπX), bằng cách nhân ph−ơng trình với cos(nπX)dX rồi tích phân trong khoảng X ∈ [0,1]: 19 dX)Xncos()F,X( 1 0 ∫ πθ = An(F) dX)Xn(cos 1 0 2∫ π = ⎪⎩ ⎪⎨ ⎧ ≠∀ = 0n),F(A 2 1 0nkhi),F(A n o Do đó, An(F) phải xác định theo θ(X,F) bởi quan hệ: (An) ⎪⎩ ⎪⎨ ⎧ ≠∀∫ πθ= ∫θ= 0n,dX)Xncos()F,X(2)F(A dX)F,X()F(A 1 0 n 1 0 o 4. Lập ph−ơng trình vi phân th−ờng cho An(F) bằng cách tính dF d An(F) theo hệ (An): - Khi n=0, dF )F(dAo = dX 1 0 F∫ θ = dX 1 0 xx∫θ =θx 10 = θx(1,F) - θx(0,F) = 0 - (-1) = 1 → Ao(F) = F + c1 Điều kiện đầu cho Ao(0) = dX)0,X( 1 0 ∫ θ = 0 = c1 ⇒ Ao(F) = F - khi ∀n ≠ 0, có: dF )F(dAn = 2 dX)Xncos( 1 0 F π∫θ = 2 dX)Xncos( 1 0 xx π∫θ , (phân đoạn tích phân) = 2 { 10x |)]Xncos([ πθ + nπ ∫1 0 )sin( dXXnx πθ }= 2{1+2π 10|)Xnsin([ πθ - nπ ∫ πθ1 0 ]}dX)Xncos( = 2{1-n2π2 2 )F(An } → ph−ơng trình vi phân cho An(F) là: A'n = 2 - n 2π2An →A'n +(n2π2)An = 2 có nghiệm tổng quát An(F) = 2)n( 2 π + c1 F2)n(e π− . Điều kiện ban đầu cho An(0) = 2 dX)Xncos()0,X( 1 0 ∫ πθ → 2)n( 2 π + c1 = 0 → c1 = - 22n 2 π , do đó: An(F) = 22n 2 π - 22n 2 π F2)n(e π− 5. Vậy nghiệm bài toán (θ) đã cho là: θ(X,F) = Ao(F) + ∑ π∞ =1n n )Xncos()F(A , tức: θ(X,F) = F + 22π ∑ π∞ =1n 2n )Xncos( - 2 2 π ∑ π∞ =1n 2n )Xncos( . 20 exp(-n2π2F) hay, do tổng ∑ π π∞ =1n 22n )Xncos(2 = 2 1 X2 - X + 3 1 , có: θ(X,F) = F + ( 2 1 X2 - X + 3 1 ) - 2 2 π ∑ π∞ =1n 2n )Xncos( exp(-n2π2F) * Phân bố θ(X,F) và t(x,τ) có dạng: O 1 x θ q=0 F=0 1 2 31 O 1 x q=0 =0 1 2 3 to τ t H14. Phân bố θ(X,F) H15. Phân bố t(x,τ) Tr−ờng nhiệt độ trong vách tăng vô hạn, có dạng: t(x,τ) = to( 22 1 δ x 2- δ 1 x+ 4 3 ) + to[ 2 a δ τ - 22π ∑ ∞ = δπ 1n 2n )/xncos( exp (- 2 22 an δ π τ)] 2.5. Ph−ơng pháp Fourier cho bài toán không ổn định nhiều chiều Các bài toán nhiều chiều không ổn định có thể giải bằng ph−ơng pháp tách biến lặp, hoặc ph−ơng pháp quy về nhiều bài toán không ổn định một chiều. 2.5.1. Ph−ơng pháp tách biến lặp 2.5.1.1. Nội dung ph−ơng pháp tách biến lặp gồm các b−ớc: 1. Tách riêng biến thời gian tìm hàm thời gian F(F) 2. Lần l−ợt tách các biến toạ độ và tìm các nghiệm riêng theo từng toạ độ. 3. Xác định các hằng số theo ĐKĐT và biểu diễn nghiệm bài toán ở dạng tích các nghiệm thu đ−ợc. 2.5.1.2. Ví dụ: Bài toán trụ vô hạn biên W1 với điều kiện đầu tổng quát * Phát biểu BT: Cho trụ l=∞ có a, t (R, ϕ,τ) = t1 và ĐKĐ bất kỳ t(ρ,ϕ,0) = g(ρ,ϕ). Tìm tr−ờng nhiệt độ t(ρ,ϕ,τ). z = g( , ) O ρϕ ρ ϕt( , ,0) ρϕ t a t1 R H16. Bài toán trụ tổng quát δ 21 * Mô hình TH: (t) ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ϕρ=ϕρ =τϕ ρ+ρ+= ϕϕρρρτ ),(g)0,,(t t),,R(t )t1t1t(at 1 2 Chuẩn hoá bằng cách đặt: θ = 1 1 t tt − , r = R ρ , ϕ' = ϕ, F = 2R aτ , ta có: (θ) ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ ϕ=−ϕρ=ϕθ =ϕθ θ+θ+θ=θ ϕϕ ),r(f t t),(g)0,,r( 0)F,,1( r 1 r 1 1 1 2rrrF * Giải bằng ph−ơng pháp tách biến lặp: 1. Tách biến thời gian bằng cách đặt θ(r,ϕ,F) = W(r,ϕ)F(F) WF' = WrrF + r 1 WrF + 2r 1 WϕϕF → F 'F = W Wrr + r 1 W Wr + 2r 1 W Wϕϕ = -k2 → ⎪⎩ ⎪⎨ ⎧ =+++ =→=+ − 011 )(0 2 2 2' 2 WkW r W r W eFFFkF rrr Fk ϕϕ 2. Tách biến r, ϕ bằng cách đặt W(r, ϕ) = R(r)φ(ϕ) → R"φ + r 1 R'φ + 2r 1 Rφ" + k2 Rφ = 0 → - φ φ" = r2 R "R + r R 'R + k2r2 = n2 → ⎪⎩ ⎪⎨⎧ =−++ +=→=+ ,0)('" cossin)(0" 2222 2 RnrkrRRr nBnAn ϕϕϕφφφ là ph−ơng trình Bessel cấp n, có nghiệm là: R(r) = cJn(kr) + DYn(kr) = cJn(kr) do 0rlim→ Yn(kr) = ∞ Trong đó Jn(kr) là hàm Bessel cấp n loại 1, có dạng: Jn(kr) = ∑∞ = + ++Γ − 1i i2ni )1in(!i ) 2 kr()1( ,Với dxxe)s( 0 1sx∫∞ −−=Γ là hàm giai thừa 22 đặt = F(r,n) Gama 3. Xác định các hằng số: -Xác định k theo ĐK biên: θ (1, ϕ, F) = 0 = R(1) φ (ϕ) F(F) → R(1) = 0 = cJn(k) → Jn(k) = 0 → có vô số nghiệm kmn ứng với giá trị cấp n của ph−ơng trình Jn(kmn) = 0. Do đó, θ(r, ϕ,τ) = RφF = ∑ ∑∞ = ∞ =1m 0n n J (kmnr) (Amnsin nϕ + Bnm cos nϕ) Fe mnk 2− - Anm và Bnm xác định theo điều kiện đầu: θ (r, ϕ, 0) = f(r, ϕ) = ∑ ∑∞ = ∞ =1m 0n n J (kmnr) (Amnsin nϕ + Bnm cos nϕ). Nhân 2 vế với sinnϕdϕ và lấy ∫π2 0 (~)dϕ có 444 344 21 ϕϕϕ∫π dnsin),r(f 22 0 = ∑∞ =1m n J (kmn'r) Amn' 43421 π= ∫ ϕϕ π2 0 2 d'nsin . Lại nhân hai vế với rJn(kmnr) dr và lấy ∫ 1 0 (~)dr có dr)n,r(F)rk(rJ 1 0 mnn∫ = πAmn' dr)rk(rJ 1 0 mn 2 n∫ = πAmn' [ )k(J21 mn2 1n+ ] Do đó có: Amn = )k(J 2 mn 2 1n+π )rk(rJ 1 0 mnn∫ ∫ ϕ π2 0 ),r(f sin nϕ dϕ dr T−ơng tự xác định các hằng số Bmn: - khi n = 0 → Bmo = )( 1 2 1 moknJ )( 1 0 rkrJ moo∫ ∫ ϕπ2 0 ),r(f dϕ dr - ∀n ≠ 0 → Bmn )( 2 2 1 mnn kJ +π )( 1 0 rkrJ mnn∫ ∫ ϕπ2 0 ),r(f cosϕ dϕ dx * Ví dụ: Cho f(r,ϕ) = 1, tức t(ρ,ϕ,0) = 2t1, có: ⎪⎪⎩ ⎪⎪⎨ ⎧ ⎪⎩ ⎪⎨ ⎧ = ≠∀ = ∀= ∫1 0 2 1 0,)( )( 1 0,0 ,0 nkhidrrkrJ kJ n B nA mno mn mn mn π có: 23 → Bmo = )( 2 2 1 2 momo kJk [kmorJ1(kmo)] 10 = )k(Jk 2 mo1mo Lúc này, nghiệm bài toán (θ) là: θ(r,F) = ∑∞ =1 1 )( 2 m mm kJk Jo(kmr)exp (-k 2n F) 2.5.2. Ph−ơng pháp quy về các bài toán 1 chiều 2.5.2.1. Nội dung: ý t−ởng của ph−ơng pháp này là quy đổi bài toán không ổn định nhiều chiều thành ra tích của các bài toán không ổn định một chiều. Phạm vi sử dụng: Để giải các bài toán nhiều chiều biên đồng nhất. 2.5.2.2. Ví dụ: Bài toán làm nguội thanh trụ W1. * Phát biểu BT: Cho trụ hữu hạn roxL có a, t(r,z,0) = ti, tw = to đồng nhất trên biên W ∈ V. Tìm t(r,z,τ) * Mô hình TH: (t) ⎪⎪⎩ ⎪⎪⎨ ⎧ = === =++ i oo zzrrr tzrt TNtLrtrtzrt t a tt r t )0,,( )0(,),,(),0,(),,( 11 τττ τ z O t to r ro t i toL H17. BT trụ hữu hạn Đổi biến T = t - to, (t) trở thành: (T) ⎪⎪⎩ ⎪⎪⎨ ⎧ −= =τ=τ =++ τ )tt()0,z,r(T )TN(0),0,r(T),z,r(T T a 1TT r 1T oi o zzrrr * Giải bằng ph−ơng pháp quy về 2 bài toán 1 chiều: - Tìm nghiệm dạng T(r,z,τ) = R(r,τ)Z(z,τ) ⇒ Ph−ơng trình vi phân DN có dạng: RrrZ+ r 1 RrZ + RZzz = a 1 (RτZ+RZτ) → (Rrr+ r 1 Rr - a 1 Rτ) Z 24 + (Zzz - a 1 Zτ)R = 0. Ph−ơng trình này đ−ợc thoả mãn khi đồng thời có: ⎪⎪⎩ ⎪⎪⎨ ⎧ =− =−+ τ τ 0Z a 1Z 0R a 1R r 1R zz rrr Đây là ph−ơng trình vi phân cho 2 bài toán 1 chiều - Các ĐK biên t−ơng ứng là: T(ro,z,τ) = 0 = R(ro,τ) Z(z,τ) → R(ro, τ) = 0 T(r,0,τ) = 0 = R(r,τ) Z(0,τ) T(r,L,τ) = 0 = R(r,τ) Z(L,τ) } → Z(0,τ) = Z(L,τ) = 0 - ĐK đầu là T(r,z,0) = (ti - to) = R(r,0) Z(z,0) Chọn R(r,0) = (ti - t1) và Z(Z,0) = 1, ta có 2 bài toán không ổn định 1 chiều là: (R) ⎪⎪⎩ ⎪⎪⎨ ⎧ −= =τ =+ τ )tt()0,r(R 0),r(R R a 1R r 1R oi o nrr và (Z) ⎪⎪⎩ ⎪⎪⎨ ⎧ = == = 1)0,( 0),(),0( 1 zZ LZZ Z a Zzz ττ τ Nghiệm của bài toán (R) và (Z) là: R(r,τ) = 2(ti - to)∑∞= − 1 )( )()( )( 2 2 m omoom r ark mo rkJrk erkJ o om τ (với km là nghiệm ph−ơng trình Jo (kmro) = 0), Z(Z,τ) = π 4 ∑ + π+∞ =0n )1n2( L Z)1n2sin( 2L a22)1n2( e τπ+− Do đó, nghiệm của bài toán 2 chiều là T(r,Z,τ) = R,Z: T(r,z,τ) = π 8 (ti-to)∑∑∞ = ∞ = + π+ 1m 0n om1om mo )1n2)(rk(Jrk L Z)1n2sin[()rk(J τπ a L nkm e ])12([ 2 2 22 ++− 2.5.3. Định lý giao nghiệm 2.5.3.1. Định lý 1: 25 Nếu X(x, τ) là nghiệm của ph−ơng trình Xτ = aXxx Y(y, τ) là nghiệm của ph−ơng trình Yτ = aYyy Z(z, τ) là nghiệm của ph−ơng trình Zτ = aZzz thì θ(x,y,z,τ) = X(x,τ)Y(y,τ) Z(z,τ) là nghiệm của ph−ơng trình θτ = a(θxx + θyy + θzz) 2.5.3.2. Chứng minh: θτ = XτY Z+ XYτZ + XYZτ = aXxxYZ + XaYyyZ + XYaZzz = a(θxx + θyy + θzz) 2.5.3.3. Định lý 2: Nếu X, Y, Z nói trên đồng thời thoả mãn cùng một điều kiện biên tuyến tính thuần nhất fw(θ, θx,θy,θz) = 0, thì θ = XYZ cũng thoả mãn điều kiện biên đó. Chứng minh theo định nghĩa ĐKB thuần nhất. 2.5.3.4. ứng dụng các định lý giao nghiệm: Có thể tìm nghiệm bài toán không ổn định ở các vật thể hữu hạn, thoả mãn ph−ơng trình θτ = a∇2θ và cùng một ĐKB thuần nhất, nh− là tích các nghiệm của bài toán một chiều t−ơng ứng, ví dụ với vật V có dạng hộp, trụ hữu hạn, đới cầu v.v. 26 Ch−ơng 3: ph−ơng pháp toán tử phức và các bài toán dao động nhiệt 3.1. Bài toán dao động nhiệt 3.1.1. Khái niệm dao động nhiệt - Dao động nhiệt là hiện t−ợng nhiệt độ của vật thay đổi tuần hoàn theo thời gian. - Nếu một vật đ−ợc gia nhiệt và làm lạnh tuần hoàn thì trong vật xuất hiện giao động nhiệt. - Khi đó tr−ờng nhiệt độ t (trong vật và của môi tr−ờng) có dạng một hàm tuần hoàn theo thời gian τ, tức là t = t(τ). Hàm tuần hoàn này luôn có thể coi nh− tổng các hàm điều hoà dạng t(τ) = t1sin ωτ hoặc t(τ) = t1cosωτ, trong đó ω = o 2 τ π , với τo là chu kỳ dao động, [s]. Ví dụ: Dao động nhiệt do sự quay của địa cầu theo ngày đêm có τo = 24h, theo năm τo = 365,24922 ngày, hoặc do chu kỳ cấp nhiệt của động cơ đốt trong có τo = 0,05s. - Khi tr−ờng t(x, y, z, τ) dao động, cả vectơ dtagrr = onr n t ∂ ∂ và dòng nhiệt q = -λ n t ∂ ∂ cũng là những hàm tuần hoàn theo thời gian τ. 3.1.2. Mô hình một bài toán dao động nhiệt - Cũng nh− một bài toán không ổn định tổng quát, mô hình một bài toán DĐN đ−ợc cho bởi một ph−ơng trình vi phân dạng τ∂ ∂t = a∇2t và các điều kiện đơn trị, trong đó nhiệt độ hoặc dòng nhiệt q tại các biên đ−ợc cho nh− một hàm tuần hoàn theo τ. - Điều kiện ban đầu th−ờng coi là phân bố đều, dạng t(x, y, z, τo = 0) = to = const. Khi đó nhiệt độ trên biên W1 hoặc nhiệt độ môi tr−ờng gần biên W3 đ−ợc cho ở dạng: t(0,τ) = t1cos (ωτ) + to hoặc 27 ⎪⎩ ⎪⎨ ⎧ τ−τλ α−=τ +ωτ=τ )],0(t)(t[),0(t t)cos(t)(t fx o1f trong đó chứa to nh− nhiệt độ ban đầu. Do đó, điều kiện ban đầu không cần xét riêng trong mô hình bài toán DĐN. - Ví dụ: Bài toán dao động nhiệt tìm t (x, τ) trong vật nửa vô hạn 0 < x < ∞, có nhiệt độ đầu phân bố đều t(x,0) = to, khi nhiệt độ biên W1 có dao động t(0,τ) = to + t1cos ττ π o 2 = to + t1 cos (ωτ), sẽ t−ơng ứng mô hình sau: (t) ⎪⎪ ⎪ ⎩ ⎪⎪ ⎪ ⎨ ⎧ =τ=τ∞ ωτ+=τ ∂ ∂=τ∂ ∂ ∞→ ox 1o 2 2 t),x(tlim),(t )cos(tt),0(t x t a t Trong đó to là nhiệt độ ban đầu phân bố đều, không đổi, và t1 là biên độ của dao động nhiệt. 3.2. Ph−ơng pháp toán tử phức hay tổ hợp phức (Complex Combination): 3.2.1. Nội dung ph−ơng pháp toán tử phức (TTP): Dùng phép biến đổi toán tử phức (sẽ giải thích sau đây) để chuyển bài toán thực (T) thành bài toán phức dạng (W) = (T) + i(S) với (S) là bài toán ảo, tìm nghiệm bài toán phức bằng cách tách biến ở dạng W(x,τ) = X(x)eiωτ, sau đó xác định nghiệm bài toán thực (T) nh− phần thực của nghiệm phức T(x,τ) = ReW(x,τ). Ta sẽ giải thích nội dung này bằng các b−ớc sau: 3.2.2. Các b−ớc của ph−ơng pháp toán tử phức: 1. Lập bài toán ảo (S) theo mô hình bài toán thực (T), trong đó thay cos(ωτ) bởi sin(ωτ) hoặc ng−ợc lại, và thay T bởi hàm số ảo S. 2. Tổ hợp bài toán thực (T) và bài toán ảo (S) thành bài toán phức (W) theo công thức của phép biến đổi toán tử phức: (W) = (T) + i(S), với i là đơn vị ảo. Trong bài toán phức, các ĐKB phức sẽ có dạng H18. BT dao động nhiệt x t to t + t coso 1 ωτ t = atxxτ a, λ t(x, )τ 28 cos(ωτ) + i sin(ωτ) = eiωτ theo công thức Euler. 3. Giải bài toán phức bằng cách tách biến, tức tìm nghiệm phức dạng W(x,τ) = X(x)eiωτ. Tr−ờng hợp này sẽ tìm đ−ợc X(x) nh− nghiệm một ph−ơng trình vi phân th−ờng. 4. Biểu diễn nghiệm phức ra dạng đại số W(x,τ)= T(x,τ) + iS(x,τ), và thu đ−ợc nghiệm bài toán thực T(x,τ) = ReW(x,τ). Ta sẽ minh hoạ các b−ớc bằng ví dụ sau. 3.3. Bài toán dao động nhiệt trong vật bán vô hạn 3.3.1. Phát biểu và mô hình (Nh− mục 3.1.2) Đổi biến số T = t(x,τ) - to, bài toán (t) trở thành (T) ⎪⎩ ⎪⎨ ⎧ =τ=τ∞ ωτ=τ = ∞→ τ 0),x(Tlim),(T )cos(t),0(T aTT x 1 xx 3.3.2. Giải bằng ph−ơng pháp THP: 1. Lập bài toán ảo (S) t−ơng ứng với bài toán (T) (S) ⎪⎩ ⎪⎨ ⎧ =τ=τ∞ ωτ=τ = ∞→ τ 0),x(Slim),(S )sin(t),0(S aSS x 1 xx 2. Nhân (S) với i rồi cộng (T) để đ−ợc (W) = (T) + i(S) (W) ⎪⎩ ⎪⎨ ⎧ =τ=τ∞ =ωτ+ωτ=τ = ∞→ ωτ τ 0),x(Wlim),(W et)]sin(i)[cos(t),0(W aWW x i 11 xx 3. Tìm nghiệm bài toán phức bằng cách tách biến: - Tìm nghiệm ph−ơng trình Wτ = aWxx ở dạng W(x,τ) = X(x)eiωτ, phải có: X(x)iωeiωτ = aX(x) eiωτ → X"(x) - a iω X(x) = 0 → X(x) = C1 a ix e ω− + C2 a ix e ω d...à τd dt i 1k+ với thông số tự chọn ϕ, 0 ≤ ϕ ≤ 1, theo công thức ti,k+1 = ti,k + [(1- ϕ) τd dt i k + ϕ τd dt i 1k+ ]∆τ (T) Ph−ơng pháp xấp xỉ Euler, Implicit, Crank-Nicolson là các dạng đặc biệt của ph−ơng pháp tổng quát, t−ơng ứng ϕ = 0, ϕ = 1, ϕ =1/2. Ví dụ cho nút trong :[ϕpti-1+ (1-2ϕp)ti+ϕpti]k+1 = {(1-ϕ)pti-1+ [1- 2(1-ϕ)pti +(1-ϕ)pti-1}]k Các chú ý: - Khi chọn b−ớc thời gian ∆τ càng nhỏ thì nghiệm số càng chính xác, nh−ng khối l−ợng tính toán càng lớn. - Khi chọn ∆τ phải đảm bảo điều kiện ổn định nghiệm, điều kiện này tuỳ thuộc ph−ơng pháp xấp xỉ thời gian đ−ợc sử dụng: ff Euler: ∆τ ≤ )B1(a2 2 + ∆ ff Cr-Nic:∆τ ≤ )B1(a 2 + ∆ Với B = λ ∆α Ph−ơng pháp ẩn nghiệm số luôn ổn định, không hạn chế ∆τ 53 Ph−ơng pháp Crank - Nicolson cho nghiệm chính xác nhất, b−ớc ∆τ lớn hơn ph−ơng pháp Euler, nh−ng hệ ph−ơng trình đại số dạng ẩn và phức tạp. 5.4. Các ph−ơng pháp giải hệ ph−ơng trình đại số tuyến tính: - Dạng tổng quát của hệ ph−ơng trình đại số tuyến tính n ẩn xi, ∀i =1ữn là ∑ = n 1j jijxa =bi, ∀i =1ữn hay : a11x1 + a12 x2 +...+ a1nxn = b1 a21x1 + a22 x2 +...+ a2nxn = b2 an1x1+an2x2+...+annxn =bn Hệ này có thể viết ở dạng ma trận Ax = b ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ a ... a a ... a ... a a a ... a a nn22n1 2n2221 1n1211 ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ x . x x n 2 1 ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ = n 2 1 b . b b Nếu định thức det (A) ≠ 0 hệ Ax = b có nghiệm: xi = )Adet( )Ddet( i ∀i =1ữn, trong đó Di là ma trận vuông nhận đ−ợc bằng cách thay cột ai ∈ A bởi cột bi. Để tính các định thức có thể dùng quy tắc Cramer, tuy nhiên không tiện lợi khi n > 4. - Khi n khá lớn có thể giải hệ trên bằng ph−ơng pháp khử của Gauss hay Jordan hoặc các ph−ơng pháp lặp (Jacob, Gauss- Seidel.v.v.) Nội dung các ph−ơng pháp này có thể xem trong các tài liệu về ph−ơng pháp tính. Các b−ớc giải theo các ph−ơng pháp trên thông th−ờng đã có các ch−ơng trình mẫu 54 5.5. FDM cho bài toán KOĐ một chiều tổng quát: Để minh hoạ các b−ớc áp dụng FDM, sẽ giải bài toán không ổn định một chiều qua vách tổng quát với biên W2 + W3 và điều kiện đầu tổng quát t(x,τ = 0) = f(x) tτ = atxx tx(0,τ) = λ q− tx(L, τ) = λ α− [t(L,τ)-ιf] t(x,0) = f(x) ta sẽ giải hệ (t) bằng FDM theo các b−ớc đã nêu. O t ∆ 2 ∆ ∆2 1 2 43 α ft ∆ t(x,0) = f(x) W2 W3 L x i i+1 q H28. FDM cho bài toán một chiều tổng quát 1) Chia vách ra các phần ∆ = 3 L (bằng nhau), đặc tr−ng cách bởi 4 phần tử nh− trên. 2) Chuyển hệ (t) về hệ ph−ơng trình vi phân của τd dt i bằng cách viết ph−ơng trình cân bằng nhiệt cho các nút (SP vật lý): - Với nút 1 ∈ W2: ρCv ),2(2 1 1 ttq d dt −∆−= ∆ λ τ → gọi a = λ/ρCv→ τd dt1 = )tt( a2qa2 212 +∆−∆λ - Với các nút i ∈ V: ∆ρ=−∆ λ−−∆ λ +− v1iii1i C)tt()tt( τd dt i → τd dt i = ( a 2∆ ti-1 - 2ti + ti+1), ∀i= 2,3 -Với nút 4 ∈ W3: ρCv 2 ∆ ατ −− ∆= )( 2 43 4 tt d dt (t4 - tf) → τd dt4 = 2 a ∆ [2t3-2(1+ λ ∆α )t4]+ λ ∆α ∆2 a2 tf (τ) 55 Đặt λ ∆α = B ta có hệ ph−ơng trình vi phân cấp 1 là: τd dt i = 2 a ∆ (-2ti+ 2ti+1) + 2 a ∆ . λ ∆q i = 1 (ti) τd dt i = 2 a ∆ (ti-1- 2ti + ti+1) i = 2,3 τd dt i = 2 a ∆ [ti-1- 2(1 + B)ti] + 2 a2 ∆ Btf i = 4 3) Chia thời gian ra các ∆τ bằng nhau và áp dụng ph−ơng trình xấp xỉ, ví dụ theo Euler, để chuyển hệ (t) về hệ ph−ơng trình đại số, theo ph−ơng trình: ti,k+1 = ti,k + τd dt i k ]∆τ - ∀i ∈ V có: ti,k+1 = ti,k + 2a∆ (ti-1- 2ti + ti+1)k∆τ, đặt F a 2 =∆ τ∆ → ti,k+1= [Fti-1+(1-2F)ti+ Fti+1]k ,(i = 2,3) - ∀i ∈ W2: ti,k+1 = ti,k + [ 2a∆ (-2ti+ 2ti+1) + 2 a ∆ . λ ∆q ]∆τ → ti,k+1 = F[( +− it)2F 1 2ti+1]k + 2Fq λ ∆ ,(i = 1) - ∀i ∈ W3: ti, k+1 = tik + [ ] ⎭⎬ ⎫ ⎩⎨ ⎧ ∆++−∆ − f2i1i2 Bt a2t)B1(2t2a ∆τ → ti,k+1 = F[2ti-1+( F 1 -2-2B)ti]k + 2FBtf, (i=4) Vậy hệ (ti) chuyển thành hệ ph−ơng trình đại số sau: t1,k+1 = F[( F 1 -2)t1+2t2]k+2Fq λ ∆ t2,k+1 = F[(t1+ ( F 1 -2)t2+t3]k t3,k+1 = F[(t2+ ( F 1 -2)t3+t4]k t4,k+1 = F[(2t3+ ( F 1 -2-2B)t4]k+2FBtf 56 4) Lập dạng ma trận cho hệ ph−ơng trình đại số trên: Các chú ý và nhận xét: Với ph−ơng pháp Euler, phải chọn ∆τ ≤ )B1(a18 L )B1(a2 22 +=+ ∆ , do ∆ = 3 L - Nếu xấp xỉ theo ph−ơng pháp Crank-Nicolson, hệ ph−ơng trình là: k+1 ( F 1 + 1) -1 t1 - 2 1 ( F 1 + 1) - 2 1 t2 - 2 1 ( F 1 + 1) - 2 1 t3 -1 ( F 1 +1+ B) t4 t1 )2F 1( − 2 t1 q λ ∆ t2 1 )2F 1( − 1 t2 +2F 0 t3 1 )2F 1( − 1 t3 0 t4 2 )B22F 1( −− t4 k Btf = k+1 ( F 1 - 1) 1 t1 q λ ∆ 2 1 ( F 1 - 1) 2 1 t2 0 2 1 ( F 1 - 1) 2 1 t3 +2. 0 1 ( F 1 -1- B) t4 Btf k = 57 ( F 1 + 2) -2 t1 t1 q λ ∆ -1 ( F 1 + 2) -1 F t2 = t2 +2F 0 -1 ( F 1 + 2) -1 t3 t3 0 1 ( F 1 +2+2 B) t4 k+1 t4 k Btf - Nếu xấp xỉ theo ph−ơng pháp thuần ẩn, ph−ơng trình là: Các nhận xét khác: 1. Các hệ ph−ơng trình đại số trên có thể giải bằng cách thay điều ban đầu t(x,0) = f(x) vào t(i,0) = f(i∆) của ma trận cột [ti]k=0: [t]0 = [f(0), f(∆), f(2∆), f(3∆)]T sau đó tìm ti tại k=1, tiếp theo dùng ma trận[t]1 tính [t]2 v.v. 2. Vì B = λ ∆α ,tf q đều nằm trong ma trận hệ số tự do,nên các đại l−ợng α, λ, tf q trong điều kiện biên có thể cho nh− những hàm của (x,y,z, τ), chẳng hạn: (α, λ, tf q) = f(x,τ) Khi giải ph−ơng trình loại này chỉ cần thay f(x,τ) bởi f(i∆, k∆τ) vào các vị trí t−ơng ứng trong các ma trận. 3. B−ớc thời gian của phép xấp xỉ Crank-Nicolson là: ∆τCN ≤ )B1(a9 L )B1(a 22 +=+ ∆ tức có thể chọn gấp 2 b−ớc ∆τE của ph−ơng pháp Euler. 5.6. FDM giải bài toán không ổn định 2 chiều tổng quát: 5.6.1. Phát biểu bài toán : Tìm tr−ờng nhiệt độ t(x,y,τ) trong vật hai chiều với điều kiện biên loại 1,2,3 và điều kiện đầu đ−ợc cho ở dạng tổng quát: 58 1 4 107 τft ( ) W2 W3 I x iW1 2 5 118 3 6 129 t(x,y,0) = f(x,y)t(x,0, ) = t (x, )ττ o q=0 W20J q(y, )τ f yα ( ,τ) (x, y, )λ = λ τ hoặc (t)λ = λ q = q (x, y, )τ v v O t H.29-MH bài toán KOĐ hai chiều tổng quát 5.6.2. Mô hình TH: Tìm t(x,y,τ) thoả mãn hệ: c q y t x tat vρτ +∂ ∂+∂ ∂=∂ ∂ )( 2 2 2 2 t(x,0,τ) = t1(x,τ) (W1) ty(x,j,τ) = 0 (W20) tx(0,y,τ) = λ− τ),y(q (W2) tx(I,y,τ) = [ ),,(),( τλ τα yIty− - tf(y,τ)] (W3) t(x,y,0) = f(x,y) (τo) 0 ≤ x ≤ I, 0 ≤ y ≤ J, (τ > 0) 5.6.3. Giải bằng ph−ơng pháp sai phân hữu hạn: 1) Chia I, J ra các khoản ∆x, ∆y, đánh số thứ tự các nút (theo ph−ơng có số nút ít nhất), không kể các nút ∈ W1 ( có t đã biết). 2) Sai phân theo toạ độ: - ∀(ij) ∈ V, qv = 0 có : τd dtij = 2 x a ∆ (ti-1-2ti+ti+1)j+ 2y a ∆ (tj-1-2tj+tj+1)i Nếu qv ≠ 0: τd dtij = a[ 2 y 1jj1j 2 x 1ii1i i)tt2t(j)tt2t( ∆ +−+∆ +− +−+− ]+ c qv ρ - ∀(ij) ∈W20: τd dtij = 2 x a ∆ (ti-1-2ti+ti+1)j+ 2y a ∆ (2tj-1-2tj)i (t ) j o 59 - ∀(ij) ∈W2: ρc 2 x∆ ∆y τd dtij = (q-λ 2 )() 111 x y tt y tt y x tt ijijijijjiij ∆ ∆ −−∆ −+∆∆ − +−+ λλ - ∀(ij) ∈W3: ρc 2 x∆ ∆y τd dtij=[λ 2 x)] y tt y tt [y)]tt() x tt 1ijijij1ij fij ijiji ∆ ∆ −λ−∆ −λ+∆−α−∆ − +−− ∀(ij) trên hai biên (góc) thì ph−ơng trình cân bằng nhiệt kết hợp điều kiện của cả hai biên. Ví dụ biên (IJ) ∈(W3+ W20), tại nút 12 có : ρc 4 yx∆∆ τd dtij=[λ x tt ijiji ∆ −− -α(tij-tf)] 2 )( 2 1 x y tty ijij ∆ ∆ −+∆ −λ Với (0,J) tại nút 3 có: ρc 4 yx∆∆ τd dtij= ) x tt q( ijiij ∆ −λ− + 2 )( 2 1 x y tty ijij ∆ ∆ −+∆ −λ 3) Xấp xỉ tijτ để lập hệ ph−ơng trình đại số. Chẳng hạn chọn: ∆x = ∆y = ∆ và đặt F = 2∆ ∆τa , B = ),( ),( τλ τα y y ∆ = B(y,τ), nếu dùng xấp xỉ Euler, hệ ph−ơng trình đại số có dạng: (V): tijk+1 = F[ti-ij+tij-1 +( F 1 -4)tij+tij+1+ti+ij]k (W20): tijk+1 = F[ti-ij+tij-1 +( F 1 -4)tij+ti+ij]k (W2): tijk+1 = F[tij-1+( F 1 -4)tij+2ti+ij+tij+1]k+ kq F λ ∆2 (W3): tijk+1 = F[2ti-ij+tij-1+( F 1 -4-2B)tij+tij+1]k+2FBtf Nếu vật có nguồn nhiệt qv = qv x,y,τ) thì vế phải ph−ơng trình của mọi nút cần cộng thêm với ∆2.F )k,j,i( )k,j,i(q v τ∆∆∆λ τ∆∆∆ 4) Lập ma trận đại số Hệ ph−ơng trình đại số của bài toán t(x,y,τ) tổng quát có thể viết d−ới dạng ma trận là [t]k+1 = Ak[t]k+[b]k, với [b]k là tổng của 4 ma trận cột, n hàng 1 cột có trị khác 0 tại các nút có quan hệ với ti, q, αtf, qv, tính tại thời điểm τk = k∆τ : [t]k+1 = = Ak[Ft]k + [Ft1]k + k qF2 ⎥⎦ ⎤⎢⎣ ⎡ λ ∆ + [2FtfB]k+ k v 2 qF ⎥⎦ ⎤⎢⎣ ⎡ λ ∆ Nếu ký hiệu p = F 1 , F = 2∆ ∆τa , B = λ α∆ ,C1 = p - 4 - 2B(∆,τ), (t ) 60 C2 = p - 4 - 2B(2∆τ), C3 = p - 4 - 2B(3∆τ) thì hệ ma trân cụ thể là: t1 p 2 1 t1 t2 1 p 1 2 t2 t3 2 p 2 t3 t4 1 p 1 1 t4 t5 1 1 p 1 1 t5 t6 1 2 p 1 t6 t7 1 p 1 1 t7 t8 = 1 1 p 1 1 F t8 + t9 1 2 p 1 t9 t10 2 C2 1 t10 t11 2 1 C2 1 t11 t12 k+1 2 2 C3 k t12 k t1(0,τ) q(∆,τ) 0 qV(0,∆,τ) 0 q(2∆,τ) 0 qV(0,2∆,τ) 0 q(3∆,τ) 0 qV(0,3∆,τ) t1(∆,τ) 0 0 qV(∆,∆,τ) 0 0 0 qV(∆,2∆,τ) 0 0 0 qV(∆,3∆,τ) +F t1(2∆,τ) + λ ∆F2 0 +2Ftf 0 + λ F2∆ qV(2∆,∆,τ) 0 0 0 qV(2∆,2∆,τ) 0 0 0 qV(2∆,3∆,τ) t1(3∆,τ) 0 B(∆,τ) qV(3∆,∆,τ) 0 0 B(2∆,τ) qV(3∆,2∆,τ) 0 K 0 K B(3∆,τ) K qV(3∆,3∆,τ) K 61 - Ma trận của hệ ph−ơng trình đại số có dạng viết gọn là: tk+1 = AkFtk + bk hay tk+1 = AkFtk + [Ft1 + 2 q F λ ∆ + 2FtfB + λ ∆2F qv]k Trong đó Ak là ma trận vuông (n x n) các ma trận khác đều là ma trận cột (n x 1) Ma trận hệ Ak có 5 đ−ờng chéo liên hệ 5 giá trị nhiệt độ các nút. Ak là ma trận hệ, t1 là ma trận biên W1, q là ma trận biên W2, B là ma trận biên W3, qv là ma trận nguồn nhiệt. Tất cả các ma trận này có thể tính theo toạ độ x, y và thời gian theo giá trị g(i∆, j∆, k∆τ) tại mỗi nút nếu cho (λ, α, tI, tf, q, qv hay ρ) = g(x, y, τ) - Để đảm bảo tính hội tụ của các nghiệm, phải chọn b−ớc thời gian ∆τ sao cho F 1 - 4 - 2B > 0 tức F = )B2(2 1a 2 +<∆ τ∆ Vì B = B (x,y,τ) nên cần lấy B = maxBijk và có ∆τ < )Bmax2(a2 2 + ∆ Ta có thể dùng các ph−ơng pháp xấp xỉ khác để chuyển sang hệ ph−ơng trình đại số dạng ẩn, để thu đ−ợc nghiệm chính xác hơn. Hệ ph−ơng trình đ−ợc giải bằng ma trận bằng cách thay điều kiện đầu t(x,y,0) = f(x,y) vào ma trận tk lúc k = 0, là tk=0 = [f(i∆x, j∆y)] = t0→ tính tk=1. Nếu các ma trận hệ số đều phụ thuộc thời gian thì sau mỗi b−ớc phải tính lại các phần tử của ma trận tại thời điểm τk = k∆τ Ví dụ: nếu cho a = 10-6m2/s, λ = 2 W/m2K, α = 25W/m2K, và nếu chọn ∆x = ∆y = 10-2m thì B = λ α∆ = 0,125, ∆τ < )2(2 2 Ba + ∆ = 23,5 s → chọn ∆τ = 20s thì 2,02 =∆ ∆= τaF và p = F 1 = 5 62 5.7. FDM giải bài toán không ổn định 3 chiều t(x,y,z,τ) 5.7.1. Trong tọa độ vuông góc xyz. - Chia các x, y, z trong vật ra ∆x, ∆y, ∆z, tạo ra n phần tử, đặc tr−ng bởi n nút (ije). Gọi tijek là nhiệt độ tại nút có toạ độ (xi = i∆x, yj = j∆y, ze= e∆z) lúc τk = k∆τ. Sai phân theo toạ độ (bằng cách viết ph−ơng trình cân bằng nhiệt cho mọi nút) để nhận đ−ợc hệ n ph−ơng trình vi phân của τd dtije : (x,y, )α τ (x,z , )α τ (x,y, )α τ (x,y, )α τ tf(x ,y, )τ (x,y, )α τ (x,y,z, ) λ τ α τ (x,z , ) α τ(y,z, ) q (x,y,z, ) τ q(x,y, ) τ t (x,y,z) t (x,z, ) τ i i-1 i+1 I o q y x l l-1 l+1 L Z q=0 j j+1 j-1 v o w H.30 Sai phân bài toán t(x,y,z,τ) tổng quát ∀ij ∈V: τd dtije = 2 x a ∆ (ti-1-2ti+ti+1)je+ 2y a ∆ (tj-1-2tj+tj+1)ie+ 2z a ∆ (te-1-2te+te+1)ij ∀ij ∈Wn: ví dụ biên là mặt có W3 thì: ρc∆x∆y 2 z∆ τd dtije = 2 x∆ λ (ti-1-2ti+ti+1)je.∆y 2 z∆ + 2 y∆ λ (tj-1-2tj+tj+1)ie∆x 2 z∆ + [ z∆ λ (te-1-te)ij-α(te-tf)ij]∆x∆y ,.v.v... ∀ij ∈ góc giao 3 mặt, ví dụ W20(x) x W2(y) x W3(z): 63 ρc 8 zyx ∆∆∆ τd dtije = 2 x∆ λ (2ti-1-2ti)je 4 zy∆∆ + ⎥⎦ ⎤⎢⎣ ⎡ −−∆ λ − qie)tt(y j1j je 4 zx∆∆ + ⎢⎣ ⎡ ∆z λ (te-1-te)-α(te-tf) ij ⎥⎦ ⎤ 4 yx∆∆ - áp dụng phép xấp xỉ thời gian, ví dụ ph−ơng pháp Euler, có hệ tijek+1 = F[ti-ije + tij-ie + tije-1 + ( 6F 1 − )tije + tije+1+ tij+ie+ti+ijl]k (∀ijl∈V) tijek+1 = F[ti-ije + tij-ie +2tije-1 + ( 6F 1 − -2B)tije + tij+ie+ ti +ije]k +(2FtfB)k (∀ijl ∈ mặt W3) và các ph−ơng trình cho các nút ∈ mặt, cạnh, góc khác. Để nghiệm ổn định, cần chọn ∆τ < )3(2 2 λ α∆+ ∆ a - Biểu diễn ma trận, hệ ph−ơng trình đại số có dạng: [t]k+1= {A[t]F+F[tw1]+[2F λ ∆ q]+[2FtfB]+[ λ F2∆ qv]}k Trong đó, ma trận hệ A là ma trận vuông n x n có 7 đ−ờng chéo, tất cả các hệ số λ, α, tw1, q, tf, qv tính theo toạ độ (i∆x, j∆y, e∆z) tại lúc τk = k∆τ tại mỗi b−ớc tính, bắt đầu từ điều kiện đầu, lúc τ = 0. 5.7.2. Ph−ơng pháp sai phân hữu hạn trong toạ độ trụ (r,ϕ,z) Bài toán 3 chiều tổng quát t(r, ϕ, z,τ ) với điều kiện đầu và các điều kiện biên tuỳ ý sẽ đ−ợc giải một cách t−ơng tự: 1. Chia vật ra n phần tử bởi các mặt cắt ri = i∆r, ϕj = j∆ϕ, ze = e∆z đánh số thứ tự các nút từ (1→ n), (nh− hình H31) 2. Viết ph−ơng trình cân bằng nhiệt cho ∀ phần tử lúc τk bất kỳ: * ∀ (ije)∈ V, ph−ơng trình cân bằng nhiệt có dạng: ρc(ri∆ϕ.∆r∆z) τd dtijl =[ r∆ λ (ti-1-ti)je(ri- 2 r∆ )∆ϕ∆z- r∆ λ (ti-ti+1)je(ri+ 2 r∆ )∆ϕ∆z 64 +[ ϕ λ ∆ir (tj-1-tj)ie- ϕ λ ∆ir (tj-tj+1)ie∆r∆z]+[ z∆ λ (te-1-te)if - r∆ λ ( te- te-1)]ri∆ϕ∆z Khi chọn ∆r = ∆z = ∆, ph−ơng trình trên có dạng: t’ije= 2 a ∆ [(1- ir2 ∆ )ti-ije+( ϕ∆ ∆ ir )2tij-1e+tije-1- 2(2+ ϕ22 2 ∆ ∆ ir )tije+ tije+1+ ( ϕ∆ ∆ ir )2tij+1e+ (1+ ir2 ∆ )ti+ije * Với các nút trên mặt biên W, ví dụ: ∀(ije) trên biên cách nhiệt W20 tại mặt ϕ = ϕJ có ph−ơng trình cân bằng nhiệt lúc chọn ∆r = ∆z = ∆ là : t’iJe = 2 a ∆ ⎢⎣ ⎡ (1- ir2 ∆ ti-iJe+2( ϕ∆ ∆ ir )2tiJ-1e- 2(2+ ϕ22 2 ∆ ∆ ir )tiJe tije+1+(1+ ir2 ∆ )t1+1Je ⎥⎦ ⎤ * ∀(Rje) trên biên r = R loại 3, ph−ơng trình cân bằng nhiệt là: ρc (R- 4 r∆ )∆ϕ 2 r∆ ∆z τd dtRje = [ r∆ λ (tR-1-TR)je.(R 2 r∆ )∆ϕ∆z-α(tRje-tf)R∆ϕ∆z] + [ ϕ∆ λ R (tj-1-tj)Re - ϕ∆ λ R (tj-tj+1)Re] 2 r∆ )∆z + [ z∆ λ (te-1-te)Rj- z∆ λ (te-te+1)Rj](R- 4 r∆ )∆ϕ 2 r∆ . Khi chọn ∆r = ∆z = ∆ và đặt B = λ ∆α ph−ơng trình có dạng: 2Rje a t' ∆= ije-Rt 4 R 2 R 2 ⎪⎩ ⎪⎨ ⎧ ∆− ∆− 2tt ) 4 R(R 1-Rje1e-Rj 2 2 −+ ϕ∆∆− ∆+ 2⎢⎣ ⎡ ϕ2 2 ) 4 ( ∆∆− ∆ RR - ) 4 R(R ∆− ∆ + ) 4 R(R RB ∆− + ⎥⎦ ⎤ tRje 1eRj2 2 1jeR t ) 4 R(R t ++ ϕ∆∆− ∆++ ⎭⎬ ⎫ j 2 t ) 4 ( 2 ∆−∆ + R aRB . * Với các nút trên cạnh, nh− giao tuyến 2 mặt biên, thì ph−ơng trình cân bằng nhiệt kết hợp cả hai điều kiện biên. * Với các nút tại góc, nơi giao điểm 3 mặt biên ph−ơng trình cân bằng 65 nhiệt sẽ kết hợp cả 3 điều kiện biên 3. áp dụng các ph−ơng pháp xấp xỉ τd dtije , sẽ đ−ợc hệ ph−ơng trình đại số của tije, k+1, với n ẩn số. o ri Rri ri+1 rq r∆ Zl-1 Zl Zl+1∆ ∆ z z α α α ∆ϕ ϕ j+1 ϕ j+1ϕj Z H31. Sai phân bài toán trụ t(π,ϕ,z,τ) Dùng xấp xỉ Euler: tijek+1 = tijek + τd dtije ⎢k∆τ, đặt ri ∆ = δi, R ∆ = δ, ∆ϕ = γ và F = 2a∆ τ∆ , ta có hệ ph−ơng trình đại số: tijek+1 = F ⎩⎨ ⎧ δ− ) 2 i1( ti-1je+ 2) i( γ δ tij-1e+tije-1+( F 1 -4-2 2 2 γ δ )tije + tije+1+( γ δi ) tij+1e+ (1+ 2 iδ )ti+1je k⎭⎬ ⎫ ∀(ije)∈ V tijek+1 = F ⎪⎪⎩ ⎪⎪⎨ ⎧ δ− δ− 4 i1 2 tR-ije+ ) 4 1(2 2 δ−γ δ tRj-1e+tRje-1+[ ) 4 1( 24 F 1 δ− −− . ( 2 2 γ δ - 4 δ -B)]tRje+tRje+1+ ) 4 1(2 2 δ−γ δ tRj+1e k⎭⎬ ⎫ + ) 4 1( BF2 δ− tf, ∀(ije)∈ W3 và các ph−ơng trình có nút trên biên, cạnh, góc, khác. Dạng ma trận:[t]k+1 = F{A[t]+[tw1]+[ a 2∆ q]+[2Btf]+[ λ ∆2 qv]}k nh− trên với A là ma trận hệ (nxn), các ma trận khác là ma trận cột (nx1). (t) 66 5.8. FDM cho bài toán biên phi tuyến: 5.8.1. Điều kiện biên phi tuyến tính: - Định nghĩa: 1. Ph−ơng trình F(Tn, mxT ) = 0 có n ≠ 1 hoặc m ≠ 1 gọi là ph−ơng trình phi tuyến tính 2. Điều kiện biên đ−ợc mô tả bởi một ph−ơng trình vi phân phi tuyến gọi là điều kiện biên phi tuyến - Ví dụ: điều kiện biên loại 3, khi mặt vách tiếp xúc chất khí hoặc chân không, trao đổi nhiệt với môi tr−ờng chủ yếu bằng bức xạ, xác định nhờ định luật Stefan-Boltzmann, thì ph−ơng trình cân bằng nhiệt trên biên có dạng: -λTx (W,τ) = εwσ0[T4(W,τ)- 4fT ] Đó là một ph−ơng trình vi phân (hay điều kiện biên) phi tuyến tính. Nếu W tiếp xúc môi tr−ờng chân không vô hạn ngoài vũ trụ thì coi Tf = 3.K hoặc gần đúng, coi Tf = O.K Chẳng hạn, tr−ờng hợp mặt W có nhiệt độ TW lớn, trao đổi nhiệt phức tạp với không khí hoặc chân không, hay vách TĐN phức tạp với sản phẩm cháy có Tf lớn, là các bài toán biên phi tuyến 5.8.2. Bài toán biên phi tuyến 1. Phát biểu bài toán: Cho vách phẳng rộng ∞, dày δ = L có (a, λ) không đổi (hoặc phụ thuộc (M,τ)), nhiệt độ ban đầu T(x,0) = Ti[K], mặt x = δ đ−ợc cách nhiệt , mặt x = 0 tiếp xúc chân không có nhiệt độ Tf và TĐN bức xạ ra môi tr−ờng này. Tìm phân bố t(x,τ) khi τ > 0. 2. Mô hình toán học của bài toán là: 67 Tτ = aTxx T(x,0) = Ti Tx(0,τ)= λ εσ− 0 [ 4fT -T4(0,τ)] Tx(0,τ) = 0 Nhận xét: bài toán này có cùng một mô hình TH với bài toán biên thuần nhất, khi vách dày 2δ và TĐN O t 1 2 43 L x q = (T - T ) 0ε∂ f4 o4 s aλ δTδx = 0Tf H32. Bài toán biên phi tuyến bức xạ với hai môi tr−ờng khí đồng chất, có cùng nhiệt độ Tf. 3. Giải bằng FDM. 1. Chia chia chiều dày δ = L ra các khoảng ∆x = n L , ví dụ ∆x = 4 L tạo ra 5 phần tử có 5 nút là i, (i = 0ữ4) 2. Sai phân theo x: ph−ơng trình cân bằng nhiệt là: - Trên biên bức xạ xi = 0: ρcS 2 x∆ τd2 dTo = εσ0( 4fT - 40T )S - x∆ λ (T0-T1)S → τd2 dTo = xc 2 0 ∆ρ εσ ( 4fT - 40T ) - 2 xc 2 ∆ρ λ (T0-T1) Nếu chuẩn hoá bằng cách đặt θ = Ti T X = L x → ∆X = L x∆ F = 2Lcρ λτ = 2L aτ Thì do τ∂ ∂T = θ∂ ∂T . F∂ θ∂ . τ∂ ∂F = 2 i L aT . F∂ θ∂ nên: −− ∆λ εσ=τ= θ )TT( )xL(aT a L2 d dT. aT L dF d 4 0 4 f i 2 00 i 2 0 i 2 40 2 aT)xL( a )TT(L2 ∆λ −λ (t) hay S 68 (θF) dF d 0θ = ⎢⎣ ⎡ ∆ 2)X( 1 -2(1+ 30 3 i0 XLT. θ∆λ εσ )θ0+2θ1 ⎥⎦ ⎤ +2 4f 3 i0 . X LT. θ∆λ εσ Nếu đặt R= λ εσ 3i0 LT. là đại l−ợng không thứ nguyên, ta có: θ0F= 2)X( 1 ∆ [-2(1+R∆X. 30θ )θ0+2θ1]+ X R2 ∆ 4 fθ , ∈ W6 θiF= 2)X( 1 ∆ [θ0- 2θ1+θ2] θ2F= 2)X( 1 ∆ [θ1- 2θ2+θ3] ∈V θ3F= 2)X( 1 ∆ [θ2- 2θ3+θ4] θ4F= 2)X( 1 ∆ [2θ3- 2θ4] ∈W20 3. Xấp xỉ Euler θk+1= θk+θFk∆F, đặt 2)X( F ∆ ∆ = p ta có hệ ph−ơng trình đại số, viết ở dạng ma trận nh− sau: θ0,k+1= [1-2p(1+R∆X 30θ ]θ0,k+2pθ1,k+2pR∆X 4fθ θi,k+1= pθi-1k+(1-2p)θik+pθ1+1,k với i = 1,2,3 θ4,k+1= 2pθ3,k+(1-2p)θ4,k ( Với biên cách nhiệt θX = 0, do đối xứng, coi θi-1 = θi+1) 4. Chuyển hệ (θi) sang dạng ma trận, có: θ0 [1-2p(1+R∆X 30θ ] 2p θ0 2pR∆X 4fθ θ1 p (1-2p) p θ1 0 θ2 = p (1-2p) p θ2 + 0 θ3 p (1-2p) p θ3 0 θ4 k+1 2p (1-2p) θ4 k 0 Thay điều kiện đầu lúc τ = 0 (tức F = 0) theo θ(X,0) = 1 vào (θi) k 69 [θ]k=0 tính phần tử bức xạ của ma trận [1-2p(1+R∆X 3k0θ ] theo θ0,0 = 1, ta tính đ−ợc [θ]k=1. Lặp lại chu trình này, tính đ−ợc [θ]k=2, ...v. v 5.8.3. Bài toán truyền nhiệt qua vách phi tuyến: Nếu thay ĐKB tại x = L của bài toán tại H.32 bởi điều kiện biên loại 3, ta có bài toán sau: Tτ = aTxx T(x,0) = Ti Tx(0,τ) = λ εσ− 0 [ 41fT -T4(0,τ)] Tx(L,τ) = λ− α [T(L,τ)-Tf2] Bài toán này đ−ợc giải t−ơng tự nh− trên, chỉ khác ph−ơng trình cân bằng nhiệt cho nút biên W2 là: O t 1 2 43 x q = (T - T ) 0ε∂ f14 o4 s aλ Tf1 δ α Tf2 H33. Bài toán truyền nhiệt KOD phi tuyến tính pc 2 x∆ T4τ = x∆ λ (T3-T4)-α[T4-Tf2] hay với B = λ ∆α x = λ αL ∆X → θF= 2)X( 1 ∆ [2T3-2(1-B)T4]+ 2)X( B2 ∆ Tf2 Xấp xỉ Euler dẫn tới ph−ơng trình: θ4k+1= 2pθ3+[1-2(1-B)]θ4+2pBθf2 Do đó, dạng ma trận của hệ ph−ơng trình đại số là: θ0 [1-2p(1+R∆X 30θ ] 2p θ0 2pR∆X 4fθ θ1 p (1-2p) p θ1 0 θ2 = p (1-2p) p θ2 + 0 θ3 p (1-2p) p θ3 0 θ4 k+1 2p θ4 k 2B2pθf2 Ph−ơng trình này cũng đ−ợc giải với điền kiện đầu θ (X,0) = 1 Bài toán này có thể áp dụng để tính nhiệt cho vách ống sinh hơi (t) kk [1-2p(1-B) 70 của lò hơi, khi nó nhận nhiệt BX từ buồng lửa và trao cho n−ớc trong lò. Việc chọn ∆τ thoả mãn điều kiện ổn định là : min{[1-2p(1-B)],[1-2p(1+ R∆ 30θ ]} > 0 tức là: 1-2p(1+ R∆X 30θ > 0 → 1-2 2)X( F ∆ ∆ (1+R∆X 30θ ) > 0 tức phải chọn; ∆F < )](maxXR1[2 )X( 3 0 2 θ∆+ ∆ , với max θ0 = i 1f T T Nghĩa là cần chọn ∆τ < ]s[, )TxL1(a2 31f0 2 2 x λ εσ∆+ ∆ Khi giải hệ ph−ơng trình trên, phải th−ờng xuyên tính lại số hạng đầu [1-2p(1+R∆ 30θ )]k theo θ0, k sau mỗi b−ớc. 71 Ch−ơng 6 ph−ơng pháp phần tử hữu hạn finite element method (FEM) 6.1. Nội dung và các b−ớc của ph−ơng pháp phân tử hữu hạn 6.1.1. Nội dung FEM. T− t−ởng của FEM là thay bài toán giải hệ ph−ơng trình vi phân (t) bằng một bài toán biến phân t−ơng ứng, tức là tìm hàm số t làm cực tiểu một phiếm hàm I t−ơng ứng bài toán (t). Bài toán biến phân đ−ợc giải gần đúng nhờ phép xấp xỉ tích phân bằng cách thay hàm t(x,y,z,τ) bởi một hệ M hàm thời gian tn(τ) tại các nút (đỉnh) của một số hữu hạn E phần tử tạo ra vật cần xét. Kết quả cho biết, để cực tiểu biến hàm I, hàm t phải tìm cần thoả mãn hệ M ph−ơng trình vi phân th−ờng cấp 1 nh− sau: C τd d [t] = -(K+H)[t] + (h+q) Trong đó C, K, H là các ma trận vuông (MxM) của các hệ số nhiệt dung C, dẫn nhiệt λ, toả nhiệt α, còn h, q là các ma trận cột (Mx1) của các giá trị nhiệt độ môi tr−ờng tf và dòng nhịêt q qua biên loại 2, [t] = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ M 1 t .. t là ma trận cột (Mx1) của các nhiệt độ nút. Nếu bài toán (t) là ổn định , [ ]t& = 0, hàm t đ−ợc xác định theo hệ ph−ơng trình đại số: (K+ H) [t] = h + q. Nếu bài toán (t) không ổn định, thì sau khi sử dụng phép sai phân thời gian, ta thu đ−ợc một hệ M ph−ơng trình đại số, cho phép xác định t theo điều kiện ban đầu. 6.1.2. Các b−ớc áp dụng FEM Để giải hệ ph−ơng trình vi phân (t) theo FEM, có thể tiến hành các b−ớc nh− sau: 72 1. Xác định phiếm hàm I t−ơng ứng bài toán (t). Tr−ờng hợp bài toán t với biên loại 1, 2, 3 tổng quát, phiếm hàm I sẽ có dạng tích phân sau: I[t(x,y,z)] = ∫ v f(x,y,z,t,tx,ty,tz) dV + ∫ w (qt + α 2 t 2 )dw Với W là biên của vật V, còn dạng hàm f xác định theo bài toán (t) và định lý Euler-Lagrange về điều kiện cực tiểu phiếm hàm. Điều kiện cực tiểu phiếm hàm I là biến phân δI = 0. O t e 3 ki S W t(x,y, =const)τ tζ k12 fi f(e) y yk yj yi xi xj xk x H. 32 Phân bố nhiệt độ t(e) trong phần tử e hai chiều 2. Mô tả điều kiện cực tiểu δI = 0 ở dạng hệ ph−ơng trình vi phân th−ờng cấp 1 của các nhiệt độ nút (hoặc hệ ph−ơng trình đại số khi (t) là bài toán ổn định). B−ớc này có thể chia ra các b−ớc nhỏ nh− sau: 2.1. Chia V (hoặc S, hoặc δ) ra một số hữu hạn phần tử có dạng khối tứ diện (hoặc tam giác, hoặc đoạn ∆x) bởi hệ thống M điểm nút (coi là đỉnh phần tử). Đánh số thứ tự và ghi địa chỉ nút theo toạ độ (i, j, k) của mỗi nút. 2.2. Giả thiết rằng: Tại một thời điểm τ bất kỳ, phân bố nhiệt độ t(e) trong phần tử e là một hàm tuyến tính của các nhiệt độ nút (tức có dạng mặt phẳng qua 3 điểm ti, tj, tk), và xác định phân bố t (e) nh− hàm tuyến tính của các biến, có dạng f(x,y,ti,tj, tk) = t (e) (xem H.32) 73 2.3. Xấp xỉ tích phân I nh− là tổng các tích phân I(e) trên mỗi phân tử (e): I = ∑ = E 1e )e(I Cho thấy I = I(t1, t2,..., tM) là phiếm hàm của M hàm nhiệt độ nút ti(τ) và điều kiện cực tiểu δI = 0 trở thành ]t[d dI = 0, với [t] = ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢ ⎣ ⎡ M 1 t .. t tức ]t[d dI = ∑ = E 1e ]t[d dIe = 0 2.4. Dùng phép biến đổi ma trận để đ−a điều kiện trên về dạng: ]t[d dI = K[t] + C[ t& ] + H[t] - h - q = 0 Đó là hệ ph−ơng trình đại số của t và t& = τd dt 3. Nếu bài toán không ổn định, [ t& ] ≠ 0, tiếp tục dùng phép sai phân thời gian, chuyển hệ ph−ơng trình vi phân theo dτ thành hệ ph−ơng trình đại số và giải theo điều kiện đầu. B−ớc 1, 2 và 3 có thể ch−ơng trình hoá cho máy tính thực hiện: 6.1.3. Phạm vi ứng dụng FEM: - Cũng nh− FDM,FEM có khả năng giải mọi bài toán biên bất kỳ, có điều kiện vật lý và điều kiện đầu cho tùy ý, với độ chính xác cao tuỳ ý. - FEM rất tiện lợi khi hệ vật có biên dạng không quy tắc, vì khi đó chỉ cần ghi điạ chỉ các nút biên, không cần tính thể tích và diện tích mặt các phần tử nh− trong FDM. - Khi biên di động (bài toán biên loại 5), cả FDM, FEM sẽ trở nên phức tạp. 6.2. Cực tiểu của hàm nhiều biến và phép xấp xỉ tích phân Cơ sở của FEM là điều kiện cực tiểu của hàm số và phiếm hàm cùng phép xấp xỉ tích phân. 6.2.1. Điều kiện cực tiểu hàm số u = u(x1, x2,...,xn) - Theo phép tính vi phân, điều kiện cần và đủ để hàm u = u(x) đạt 74 cực tiểu tại x là: dx du = 0 và 2 2 dx ud > 0 - Điều kiện cần và đủ để u = u(x,y) đạt cực tiểu tại (x,y) là: ⎪⎩ ⎪⎨ ⎧ >−> == 0)uuu(,0u 0u,0u xy 2 yyxxxx yx - Tổng quát, điều kiện cần để hàm n biến u = u(x1, x2,...,xn) đạt cực tiểu là tất cả các đạo hàm riêng u theo lần l−ợt các biến đều triệt tiêu, tức: ix u ∂ ∂ = 0, ∀i = 1 ữ n. Sau đây sẽ dùng ký hiệu "A ∆ B", có nghĩa là " B đ−ợc định nghĩa là A". Định nghĩa một ma trận cột [x] ∆ ⎥⎥ ⎥⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎣ ⎡ n .. 2 1 x x x , điều kiện trên có thể viết ở dạng: 0 ]x[ u =∂ ∂ hay cụ thể , ]x[ u ∂ ∂ ∆ ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ ∂ ∂ ∂ ∂ ∂ ∂ xn u ... x u x u 2 1 = ⎥⎥ ⎥⎥ ⎥⎥ ⎥⎥ ⎥ ⎦ ⎤ ⎢⎢ ⎢⎢ ⎢⎢ ⎢⎢ ⎢ ⎣ ⎡ 0 ... 0 0 Sau đây ta chỉ quan tâm tới điều kiện cần nói trên, còn điều kiện đủ để đạt cực tiểu đ−ợc xác định bởi nội dung của bài toán biến phân cụ thể. Trong kỹ thuật, việc cực tiểu phiếm hàm I[t], t−ơng ứng việc xác định hàm t sao cho sai số cực tiểu, so với hàm t xác định chính xác theo định luật bảo toàn năng l−ợng. 6.2.2. Phép xấp xỉ tích phân: - Để xấp xỉ một tích phân, ta chia miền tích phân ra các phần tử hữu hạn và giả thiết rằng, hàm tích phân thay đổi tuyến tính trong mỗi phần tử. Cụ thể, coi hàm F(x) là đ−ờng thẳng trong phần tử một chiều 75 ∆x = xj - xi , hàm F(x,y) là mặt phẳng trong phần tử tam giác 2 chiều ∆e, hàm F(x,y,z) là tuyến tính với các biến tọa độ (xi, yj, zk) tại 4 đỉnh của tứ diện. - Nhờ cách chia và giả thiết nói trên, ta có thể xác định tích phân I nh− tổng các giá trị trung bình của tích phân trên mỗi phần tử I = ∑ = E 1e )e(I - Ví dụ: Biểu thức xấp xỉ tích phân trên E = 5 phần tử ở hình H.33. I = ∫ M 1 x x dx)x(F là I = ∑ = E 1e )e(I = ∑ = +E 1e ji 2 FF (xj-xi) với Fi = F(xi), Fj = F(xj), (Fi+Fj)/2 là trị trung bình của F(x) trên phần tử (e) có kích th−ớc ex∆ = xj-xi Nếu chọn xij xj - xi = ∆, M = 4, E = 5 thì ta có: ∫ ∑ +∆=M 1 x x e ji )FF(2 dx)x(F = 2 ∆ (F1+2F2 + 2F3 + 2F4 + F5) 6.3. Lý thuyết biến phân (variation Theory) (Có thể tham khảo tài liệu: Variationoe istrislenie của L.E.Elgols). 6.3.1. Phiếm hàm - Khái niệm: Phiếm hàm là một đại l−ợng biến thiên mà trị số của nó phụ thuộc vào một hay một vài hàm số nào đó. Ký hiệu: I = I[t(x)], I = I[t(x,y)] hay I = I[u1(x), [u2(x)] v.v - Đinh nghĩa: đại l−ợng biến thiên v đ−ợc gọi là phiếm hàm phụ thuộc hàm số y(x), (y(x) gọi là đối thức), ký hiệu là ν = ν[y(x)], nếu ứng với mỗi hàm y(x)thuộc một lớp hàm nào đó, có một giá trị ν xác định. T−ơng tự định nghĩa phiếm hàm ν = ν [t(x,y)] hoặc ν = ν [u1(x),u2(x)], ...v.v - Phiếm hàm th−ờng là một tích phân của một hàm F nào đó trên một miền xác định cho tr−ớc. ∆ = O t (1) x 1 2 (Fi + fj) (2) (e) (E) x1 xi xj xM x∆ F(x) H.33. Để xấp xỉ tích phân I o 76 Ví dụ: + Công mà hệ thực hiện trong một quá trình p = p(v) bất kỳ giữa hai trạng thái (p1, v2), (p2,v2) là phiếm hàm l = l[p(v)] = ∫ 2 1 v v )v(p dv + Nhiệt l−ợng do hệ trao đổi trong quá trình bất kỳ T = T(s) giữa hai trạng thái (T1,s1), (T2, s2) là phiếm hàm q = q[T(s)] = ∫ 2 1 s s )s(T ds. + Độ dài đ−ờng cong y(x) bất kỳ qua 2 điểm (x0,y0) (x1,y1) là phiếm hàm l = l[y(x)] = dx'y1 2 1 x x 2∫ + 6.3.2. Nội dung của lý thuyết biến phân: - Lý thuyết biến phân là một ngành của toán học chuyên nghiên cứu các ph−ơng pháp tìm cực trị của các phiếm hàm. Nó cho phép tìm đ−ợc một hoặc một số hàm số làm cực tiểu một phiếm hàm đã cho. - Bài toán biến phân là bài toán tìm cực tiểu của một phiếm hàm cho tr−ớc, phụ thuộc một vài hàm số ch−a biết. Phép tính biến phân cho phép tìm đ−ợc các hàm số này. - Ví dụ: các bài toán biến phân tiêu biểu: 1. Tìm mặt cực tiểu: Tìm đ−ờng cong y(x) nối hai điểm (x0,y0) (x1,y1) cho tr−ớc, để khi quay quanh Ox tạo ra mặt có diện tích cực tiểu. Phát biểu biến phân của bài toán (1) là tìm cực tiểu của phiếm hàm:S[y(x)] = dx'y11 0 x x 2∫ + 2. Tìm đ−ờng đoản thời: Tìm đ−ờng cong y(x) nối hai điểm A(x0,y0), B (x1,y1) cho tr−ớc để chất điểm lăn không ma sát trên nó từ A đến B tốn ít thời gian nhất. Theo định luật bảo toàn năng l−ợng, tìm thấy quan hệ τ = τ[y(x)] H35. BT đ−ờng đoản thời O y x xo x1 y(x) ds dx g v H.34. Bài toán tìm mặt cực tiểu O t x y(x) S x1xo z y 77 và bài toán biến phân là tìm cực tiểu phiếm hàm nh− sau: mgy = →τ + 222 ) d dydx ( 2 m dx dx dy 1 gy2 1d 2 ⎟⎠ ⎞⎜⎝ ⎛+τ → τ = τ[y(x)] = dxy 'y1 g2 1 1x 0x 2∫ + 3. Tìm đ−ờng trắc địa: Tìm đ−ờng cong f(x,y,z) nối hai điểm A,B trên mặt cong ϕ(x,y,z) = 0 cho tr−ớc, có độ dài ngắn nhất. Theo ý nghĩa hình học, bài toán (3) theo biến phân là tìm cực tiểu phiếm hàm l= l[y(x),z(x)] = dx'Z'y12 1 x x 22∫ ++ trong đó ⎩⎨ ⎧ = = )x(zz )x(yy là các hàm xác định dạng tham số của đ−ờng cong f(x,y,z)∈ϕ(x,y,z) = 0. 6.3.3. Biến phân của phiếm hàm: 6.3.3.1. Biến phân của đối thức: - Định nghĩa: Biến phân δu của hàm u(x) (hay đối thức u(x)) là hiệu hai hàm số δu = )(xu - u(x) trong đó )(xu thay đổi tuỳ ý trong lớp hàm u(x,ε) nào đó, gần đ−ờng cong u(x). T−ơng tự định nghĩa: δu(x,y) = ),( yxu - u(x,y) - ý nghĩa hình học: )x(u thay đổi trong lớp đ−ờng cong qua hai điểm biên (x0,y0) (x1,y1) phụ thuộc thông số ε nhỏ nào đó: H.36. BT tìm đ−ờng trắc địa ⎩⎨ ⎧ = = )x(zz )x(y..., x 0, 0 (4) T , 0, x , 0 (5) x T , T , T const, x , 0 (6) T , x ∂ τ ∂ τ= ∀ ∂τ ∂ ∂ τ ∂ τ= ∀ ξ ∂τ ∂ = = ≥ ∀ ξ < < ∞ τ = τ = = ∂ ∞ τ = →∞ τ >∂ ξ τ = ξ τ = = ∀ = ξ τ > ∂ ξ τλ − λ∂ ( ) ( )22 2T , dWl , x , 0x d ⎧⎪⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪ ∂ ξ τ ξ= ρ ∀ = ξ τ >⎪ ∂ τ⎪⎩ 7.2.3. Giải bằng ph−ơng pháp Stefan * Theo kết quả của bài toán (4.3) về vật bán vô hạn, ta sẽ tìm nghiệm của ph−ơng trình (1) và (2) ở dạng sau: T1 (x, τ) = 1 1 1 xA B erf 2 a ⎛ ⎞+ ⎜ ⎟⎜ ⎟τ⎝ ⎠ T2(x,τ)= 2 2 2 xA B erf 2 a ⎛ ⎞+ ⎜ ⎟⎜ ⎟τ⎝ ⎠ ,ở đây erf(x) = ( )( ) 2 n 2n 1x n 00 1 x2 2e d n! 2n 1 +∞−δ =δ= −δ= +π π∑∫ (7) 131 là hàm sai số Gauss. Các hằng số A1B1A2B2 đ−ợc xác định theo các ĐK đơn trị nh− sau: * A1 xác định theo ĐKB (4): T1 (0, τ) = Tw = A1 A2 tìm theo giả thiết cho rằng T2 (∞, τ) = To T2 (∞, τ) = To = A2 + B2 → A2 = To - B2 Vậy nghiệm riêng của (1) + (4) và (2) + (5) là: T1 (x, τ) = w 1 1 xT B erf 2 a ⎛ ⎞+ ⎜ ⎟⎜ ⎟τ⎝ ⎠ và T2 (x, τ) = o 2 o 2 2 2 x xT B 1 erf T B erfc 2 a 2 a ⎡ ⎤⎛ ⎞ ⎛ ⎞− − = −⎢ ⎥⎜ ⎟ ⎜ ⎟⎜ ⎟ ⎜ ⎟τ τ⎢ ⎥⎝ ⎠ ⎝ ⎠⎣ ⎦ * B1 và B2 sẽ đ−ợc xác định theo ĐKB (6) nh− sau: T1 (ξ, τ) = T2 (ξ, τ) = Ts có dạng: w 1 1 T B erf 2 a ⎛ ⎞ξ+ ⎜ ⎟⎜ ⎟τ⎝ ⎠ = o 2 s 2 T B erfc T 2 a ⎛ ⎞ξ− =⎜ ⎟⎜ ⎟τ⎝ ⎠ Vì (B1, B2) = const ∀τ nên các đẳng thức trên chỉ thực hiện đ−ợc khi ξ = C τ , với C là 1 hằng số nào đó sẽ đ−ợc xác định. Do đó, ĐKB (6) sẽ là: w 1 1 CT B erf 2 a ⎛ ⎞+ ⎜ ⎟⎜ ⎟⎝ ⎠ = o 2 s 2 CT B erfc T 2 a ⎛ ⎞− =⎜ ⎟⎜ ⎟⎝ ⎠ Suy ra s w1 1 T TB Cerf 2 a −= ⎛ ⎞⎜ ⎟⎜ ⎟⎝ ⎠ và o s2 2 T TB Cerfc 2 a −= ⎛ ⎞⎜ ⎟⎜ ⎟⎝ ⎠ Vậy nghiệm riêng của [(1) + (4), (2) + (5)] x (6) là: 132 T1 (x, τ) = ( )s wW 1 1 T T xT erf 2 aCerf 2 a ⎛ ⎞−+ ⎜ ⎟⎜ ⎟⎛ ⎞ τ⎝ ⎠⎜ ⎟⎜ ⎟⎝ ⎠ T2 (x, τ) = ( )o so 2 2 T T xT erfc 2 aCerfc 2 a ⎛ ⎞−− ⎜ ⎟⎜ ⎟⎛ ⎞ τ⎝ ⎠⎜ ⎟⎜ ⎟⎝ ⎠ * C đ−ợc xác định theo ĐKB loại 5 (7) nh− sau: ( )1 1 T C , x ∂ τ τλ ∂ - ( )2 2 T C , x ∂ τ τλ ∂ = 2 CWl 2 ρ τ ở đây d d ξ τ = ( )d CCd 2τ =τ τ là vận tốc di động của biên, tức là vận tốc đóng băng. Các hàm sai số Gauss có dạng: erf(x) = ( ) ( ) 2 n 2n 1 x 0 n 0 1 x2 2e d n! 2n 1 +∞δ= −δ δ= = −δ = +π π ∑∫ , erfc(x) = ( ) ( ) 2 n 2n 1 x n 1 1 x2 2e d 1 erf (x) 1 n! 2n 1 +∞∞ −δ δ= = −δ = − = − +π π ∑∫ Đạo hàm của chúng là: ( )n 2n n 0 1 xd 2erf (x) dx n! ∞ = −= π ∑ = ( ) 2n2 x n 0 x2 2 e n! ∞ − = − =π π∑ 2xd d 2erfC(x) erf (x) e dx dx −= − = − π Do đó, ĐKB (7) là λ1T1x( )C ,τ τ - λ2T2x( )C ,τ τ = 2 CWl 2ρ τ sẽ ứng với ph−ơng trình sau: 0 133 ( ) 2 1s w 1 1 1 Cexp 4aT T . aCerf 2 a ⎛ ⎞−⎜ ⎟⎜ ⎟− ⎝ ⎠λ ⎛ ⎞ π τ⎜ ⎟⎜ ⎟⎝ ⎠ + ( ) 2 2o s 2 2 2 Cexp 4aT T . aCerfc 2 a ⎛ ⎞−⎜ ⎟⎜ ⎟− ⎝ ⎠λ ⎛ ⎞ π τ⎜ ⎟⎜ ⎟⎝ ⎠ = 2 CWl 2 ρ τ . Nếu đặt C = 1K2 a , tức K = 1 C 2 a ta có ph−ơng trình để xác định C nh− sau: ( ) ( ) 222 1o s2 1 1 s w 2 2 1 aexp Kexp K aT T a erf K T T a aerfc K a ⎛ ⎞−⎜ ⎟− ⎛ ⎞⎛ ⎞ −λ ⎝ ⎠+ ⎜ ⎟⎜ ⎟λ − ⎛ ⎞⎝ ⎠⎝ ⎠ ⎜ ⎟⎝ ⎠ = ( )2 1s w 1 lW a K T T π ρ − λ . Đặt ( )2 1 os w 1 lW a K T T ρ =− λ , ph−ơng trình có dạng: f(K) = ( )oK Kπ → Giải bằng đồ thị ta có K và tìm đ−ợc C = 1K2 a K. Hằng số Ko là 1 đại l−ợng không thứ nguyên, đ−ợc gọi là tiêu chuẩn (hoặc số) Koccivich * Chuyển về dạng không thứ nguyên bằng cách đặt Fox = 12 a x τ , Fox gọi là biến Fourier của toạ độ và thời gian, Ka = 2 1 a a , ta có nghiệm của bài toán đã nêu ở dạng không thứ nguyên nh− sau: ( )1 w 1 s w T x, T T T τ −θ = − = ( ) ( ) ox 1 ox 1erf 2 F F erf K ⎛ ⎞⎜ ⎟⎜ ⎟⎝ ⎠ = θ y o K y=f(K) y= π KoK K=c/2 1a H58. Để xác định K và C. 134 ( )o 2 2 o s T T x, T T − τθ = − = ( ) ( )a ox 2 oxa 1erfc 2 K F F erfc K K ⎛ ⎞⎜ ⎟⎜ ⎟⎝ ⎠ = θ 7.2.4. Tính gần đúng trong kỹ thuật: * Do các chuỗi của erf(x) và exp(x2) hội tụ rất nhanh khi n tăng, nên với độ chính xác cho phép của kỹ thuật, có thể chỉ cần lấy số hạng đầu của các chuỗi này (ứng với n = 0) khi tính toán, tức là coi: ( )erf x = ( )( ) n 2n 1 n 0 1 x2 n! 2n 1 +∞ = − +π ∑ =& 2 xπ ( ) ( )erfc x 1 erf x= − =& 21 x− π 2Cexp 4a ⎛ ⎞−⎜ ⎟⎜ ⎟⎝ ⎠ = n2 n 0 1 C n! 4a ∞ = ⎛ ⎞−⎜ ⎟⎜ ⎟⎝ ⎠ ∑ =& 1. Khi đó có: 1 Cerf 2 a ⎛ ⎞⎜ ⎟⎜ ⎟⎝ ⎠ =& 1 2 C. 2 aπ = 1 C aπ 2 Cerfc 2 a ⎛ ⎞⎜ ⎟⎜ ⎟⎝ ⎠ =& 1 - 2 C aπ Khi đó ph−ơng trình ĐKB loại 5 để xác định C sẽ có dạng: ( )s w 1 1 1 T T a C a − πλ π τ + ( )o s 2 2 2 T T C1 a a −λ ⎛ ⎞− π τ⎜ ⎟⎜ ⎟π⎝ ⎠ = 2 CWl 2 ρ τ hay C2 = ( ) ( )1 s w 1 o s 2 2 2 2 T T 2 T T C lw lw a C ⎛ ⎞λ − λ −+ ⎜ ⎟⎜ ⎟ρ ρ π −⎝ ⎠ * Xét tr−ờng hợp To = Ts, tức là khi nhiệt độ ban đầu của pha ẩm bằng nhiệt độ đóng băng. 135 Khi To = Ts ta có: C = ( )1 s w 2 2 T T l W λ − ρ Nếu pha ẩm (2) là n−ớc, có độ ẩm w = 1, thì C = ( ) 1 21 s w 2 2 T T l ⎡ ⎤λ −⎢ ⎥ρ⎣ ⎦ - Lúc này, tr−ờng nhiệt độ trong 2 pha có dạng: T1 (x, τ) =& Tw + (Ts - Tw) 1aC π 1 xerf 2 a ⎛ ⎞⎜ ⎟⎜ ⎟τ⎝ ⎠ hay T1 (x, τ) =& Tw + (Ts - Tw) ( )21 s w l W x. 2 T T ρ λ − τ → ( ) ( ) ( ) 2 1 w s w 1 2 o s lf W xT x, T T T . 2 T x, T T const ⎧ τ = + −⎪ λ τ⎨⎪ τ = = =⎩ - Vận tốc dịch chuyển biên, tức vận tốc đóng băng, là: d d ξ τ = C 2 τ = ( ) ( )1 s w 2 T T f 2l W. λ − = τρ τ , tổng quát d d ξ τ = K 1a τ , với K = 1 C 2 a = 1 s w 2 1 (T T ) 2l Wa λ − ρ . Vậy vận tốc đóng băng chỉ phụ thuộc τ, đồng biến theo λ1, Ts nghịch biến theo Tw, l, ρ2, W và τ. Vận tốc đóng băng tỷ lệ nghịch với τ , tức là khi τ tăng 4 lần thì vận tốc giảm 2 lần. Biên chuyển động chậm dần với gia tốc ξ'' = 2 2 d d ξ τ = ( )1 s w 2 3 T T1 2 2l W λ −− ρ τ , [m/s 2] Nhận xét: Gia tốc có trị âm, làm biên di chuyển chậm dần. Khi τ lớn, có thể coi gia tốc ξ'' = 0. ρ2 =& & 3 136 Lúc này biên di chuyển gần nh− đều, nh−ng rất chậm. 7.2.5. Tính độ dày lớp băng tại thời điểm τ * Tr−ờng hợp tổng quát, độ dày lớp băng tại thời điểm τ là x = ξ = C τ , với C = 12 a K , tức x = ξ = 12K a τ , [m] * Tr−ờng hợp To = Ts và tính gần đúng bậc 1 theo x, có C = ( )1 s w 2 2 T T l W λ −ρ nên độ dày lớp băng là x = ξ = ( )1 s w 2 2 T T l W λ − τρ , [m] * Nếu pha (2) là n−ớc, có W = 1, ở điều kiện To = Ts thì x = ξ = ( )1 s w 2 2 T T l λ − τρ , m 7.2.6. Tính thời gian đóng băng đến độ dày đã cho ξ = L. * Tr−ờng hợp tổng quát với lớp băng phẳng, rộng ∞, thời gian đạt tới độ dày ξ = L = C τ là τ = 22 2 2 1 1 L L L C 2 a K 4a K ⎛ ⎞⎛ ⎞ = =⎜ ⎟⎜ ⎟ ⎜ ⎟⎝ ⎠ ⎝ ⎠ , [s] * Tr−ờng hợp To = Ts và tính gần đúng bậc 1, có τ = ( ) 2 2 1 s w l WL 2 T T ρ λ − , [s] * Với n−ớc ở To = Ts thì thời gian để tạo lớp băng phẳng, dày L là (cho W = 1): τ = ( ) 2 2 o 1 s w 1 l L L. K T T 2 2a ρ =λ − τo ξ ξ' " c' 2 ξ = τ H59. Vận tốc và gia tốc của mặt băng x = ξ 3 " 4 c τ −=ξ 137 7.3. Bài toán đông lạnh các vật ẩm hữu hạn 7.3.1. Mục đích chủ yếu khi tính đông lạnh các vật ẩm hữu hạn là tính thời gian để nhiệt độ cực đại trong vật bằng 1 trị số cho tr−ớc. Thời gian đông lạnh τ gồm 2 giai đoạn: τ = τo + τ1, trong đó τo là thời gian để hoá rắn toàn bộ vật ẩm, có nhiệt độ tâm vật bằng Ts, còn τ1 là thời gian để nhiệt độ tâm vật giảm trừ Ts đến nhiệt độ Tk cho tr- −ớc, theo yêu cầu của công nghệ cấp đông Việc tính τ1 có thể dựa vào kết quả của bài toán dẫn nhiệt không ổn định trong vật rắn 1 pha. Sau đây ta sẽ tính τo theo ph−ơng pháp gần đúng. Phép tính gần đúng sẽ dựa trên các giả thiết sau: 7.3.2. Các giả thiết 1. Các vật ẩm hữu hạn có dạng đối xứng 2. Điều kiện biên ngoài vật có tính đối xứng, loại 1 3. Nhiệt độ ban đầu trong vật ẩm là đồng nhất, và bằng nhiệt độ hoá rắn: T2 (M, τ) = Ts 4. Trong lớp vật rắn tạo thành sau chuyển pha, phân bố nhiệt độ là tuyến tính đối với biên di động x = ξ 7.3.3. Tính thời gian làm đông τo 1. Đông đặc vật ẩm phẳng, rộng 2L, có To = Ts, có λ1, l, ρ2 hai biên ngoài có Tw = const < To đối xứng. Bài toán này có mô hình giống mô hình bài toán ở trên. Điều kiện biên loại 5 trên biên di động x = ξ là: λ1T1x(ξ,τ)-λ2T2x(ξ,τ)=lρ2 dd ξ τW o x T ξ L-L T T0 s TW TW H60. Làm đông vật phẳng do T2(x, τ) = Ts = const nên T2x(ξ, τ) = 0 138 H61. Làm đông vật trụ H62. Làm đông vật cầu ξ TW Ts t ro R ξ R r TWTso Do T1(x, τ) tuyến tính với x = ξ tại ∀τ nên T1x(ξ, τ) = s wT T−ξ Vậy ph−ơng trình cân bằng nhiệt trên biên W5 có dạng: s w 1 T T−λ ξ = 2 dl W d ξρ τ hay ( )1 s w 2 T T d d l W λ −ξ ξ = τρ Thời gian làm đông τo ứng với khi ξ = L nên có: L o dξ ξ∫ = ( )o 1 s wo 2 T T d l W τ λ − τρ∫ → ( )2 1 s w o 2 T TL 2 l W λ −= τρ Vậy τo = ( ) 2 2 2 o 1 s w 1 l W L L. K T T 2 2a ρ =λ − , với Ko = ( )2 11 s w l Wa T T ρ λ − 2. Đông đặc vật ẩm hình trụ giải t−ơng tự nh− trên, ta đ−ợc τo = ( ) 2 2 2 o 1 s w 1 l W R R. K T T 4 4a ρ =λ − 3. Bài toán làm đông vật ẩm hình cầu cho kết quả τo = ( ) 2 2 2 o 1 s w 1 l W R R. K T T 6 6a ϕ =λ − Các công thức trên khi tính cho khối chất lỏng hoàn toàn thì lấy W = 1 7.3.4. So sánh thời gian τo: - Nếu các vật phẳng, trụ, cầu có cùng độ dầy tức R = L thì ta có: τof = 2τot = 3τoc Với vật ẩm hình dạng bất kỳ, thời gian đóng băng τo tỷ lệ thuận với bình ph−ơng độ dầy của vật. Độ dầy của vật đ−ợc hiểu là khoảng cách trung bình giữa hai mặt đ−ợc làm lạnh của vật. Do đó, để giảm 139 thời gian thời gian đông kết, nên giảm độ dầy cuả vật ẩm. 7.4. Bài toán đông kết vật đúc 7.4.1. Phát biểu bài toán Khi tính đông kết vật đúc, th−ờng coi vùng kim loại lỏng có nhiệt độ phân bố đều, bằng nhiệt độ nóng chảy ts. Khi đó chỉ cần tìm độ dày lớp kim loại đông kết ξ = ξ(τ) và tốc độ biên ξ, tức tốc độ ngng kết d d ξ τ = f(τ) trên cơ sở giả thiết nh− ở mục (7.2.3), tức là coi tr−ờng nhiệt độ trong lớp đã hoá rắn là tuyến tính với x = ξ. Khi đó bài toán là: ( ) ( ) x x s s w x t t ' d' ' l x x d 0 0 t ' x , t const t tt x =ξ =ξ =ξ ∂ ∂ ξ⎧λ −λ = ϕ⎪ ∂ ∂ τ⎪⎪ξ τ = =⎪⎨ > ξ τ = =⎪⎪ −∂⎪ =∂ ξ⎪⎩ d dτξ ρλ t t t o s w x qs(τ) ' c λ, , ,ρc x ξ= H63. BT đông kết vật đúc 7.4.2. Tính ξ(τ) và tốc độ đông kết Do t' = const nên t ' x ∂ ∂ = 0. Ta có ph−ơng trình: ξ −λ ws tt = d'l d ξρ τ hay ( )s wd t t d'l λξ ξ = − τρ Tích phân ph−ơng trình có 21 2 ξ = ( )s wt t C'l λ − τ +ρ Theo ξ (τ = 0) = 0 = C. Vậy: ξ = ( ) [ ]s w2 t t , m'l λ − τρ Tốc độ đông kết là ξ' = d d ξ τ = ( )s wt t 2 ' λ − ρ τ , [m/s] ρ' 140 Nếu vật đúc dày 2L, 2 biên loại 1 đối xứng thì thời gian đông kết là: τ = ( ) 2 2 o s w 'l L L. K t t 2 2a ρ =λ − , [s] 7.5. Tính truyền nhiệt khi nóng chảy lớp bảo vệ vỏ phi thuyền có vận tốc lớn 7.5.1. Vấn đề bảo vệ nhiệt cho vỏ phi thuyền Khi bay vào khí quyển với vận tốc lớn, do ma sát với không khí, vỏ phi thuyền sẽ nhận 1 l−ợng nhiệt rất lớn. H64. Lớp bảo vệ vỏ tàu bằng vật liệu nóng chảy L−ợng nhiệt này tỷ lệ với lực cản của không khí F = 2k 1 K v S 2 ρ và vận tốc v của tàu, và bằng: Qo = 3 k 1 K v S 2 ρ , [W] hay qo = 3o k Q 1 K v S 2 = ρ , [W/m2] L−ợng nhiệt nhận vào có thể làm nhiệt độ vỏ tàu tăng rất cao, gây nguy hại cho cả con tàu. Do đó, ng−ơì ta phải tìm cách giải thoát l−ợng nhiệt này, bảo đảm cho nhiệt độ thành tàu không v−ợt quá 1 giá trị an toàn Tk nào đó. Giải pháp hiện nay là bọc vỏ tàu bằng 1 lớp vật liệu có nhiệt độ nóng chảy Ts không lớn hơn Tk nói trên, Ts < Tk . Nhiệt ma sát làm nóng chảy lớp vỏ này rồi thoát ra khí quyển. Việc thiết kế lớp bảo vệ nhiệt bao gồm việc chọn vật liệu thích hợp, xác định tr−ờng nhiệt độ trong lớp nhận nhiệt nóng chảy, tính vận Ts TK TO c lTS δ ρ λ v qo 141 tốc nóng chảy và xác định độ dày đủ an toàn cho chuyến bay. Sau mỗi chuyến bay, lớp bảo vệ sẽ bị nóng chảy rồi thoát cả nhiệt lẫn chất vào khí quyển, và ng−ời ta sẽ bọc lại cho lần bay tiếp theo. 7.5.2. Phát biểu bài toán lớp nóng chảy Tìm tr−ờng nhiệt độ T(y, τ) trong lớp vật liệu có các thông số vật lý (ρ, C, λ, l, Ts) cho tr−ớc, có biên nóng chảy cho bởi hệ ph−ơng trình vi phân sau: ( ) ( ) ( ) ( ) 2 2 o s o 0 0 T Ta y T y, 0 T , T TT 0 T 0, T Tq l ξ→∞ ξ= ξ= ⎧∂ ∂⎪ =∂τ⎪ ∂⎪ τ = = ξ→∞ τ =⎪⎪∂⎪ =⎨∂ξ⎪⎪ ξ = τ =⎪⎪ ∂ξ ∂⎪ = ρ − λ∂τ ∂ξ⎪⎩ 7.5.3. Xác định tr−ờng nhiệt độ trong lớp nóng chảy Gọi vận tốc di động biên nóng chảy ξ là W = d d ξ τ . Chuyển bài toán (T) sang hệ toạ độ động (ξ, τ) bằng cách đổi biến ξ = y - Wτ. Khi đó T T T. W∂ ∂ ∂ξ ∂= = −∂τ ∂ξ ∂τ ∂ξ và 2 2 2 2 T T y ∂ ∂=∂ ∂ξ nên ph−ơng trình vi phân Tτ = aTyy có dạng: 2 2 T TW a∂ ∂− =∂ξ ∂ξ hay Tξξ + W T a ξ = 0 → Nghiệm tổng quát là T(ξ) = A exp W a ⎛ ⎞− ξ⎜ ⎟⎝ ⎠+ B, với các hằng số A, B tìm theo ĐKB: T(ξ = 0) = Ts = A + B B = To T(ξ → ∞) = To = B A = Ts - To (W5) → H65. Bài toán biên nóng chảy Ts TK TO qO T c ,lρ λW W = τ d d ξτ = y - W δO ξ y τ 142 Vậy tr−ờng T có dạng: T(ξ) = (Ts - To) exp Wa ⎛ ⎞− ξ⎜ ⎟⎝ ⎠+ To Hay ở dạng không thứ nguyên θ(ξ) = o s o T T T T − − = exp W a ⎛ ⎞− ξ⎜ ⎟⎝ ⎠ 7.5.4. Xác định vận tốc nóng chảy Vận tốc nóng chảy W = d d ξ τ đ−ợc xác định theo ĐKB (W5) qo = 0 TlW ξ= ∂ρ − λ ∂ξ = ρlw + λ(Ts - To) W a hay, do c a ρ λ= nên: qo = ρlw + Cρw(Ts - To) = wρ [l + C(Ts - To)] Vậy vận tốc nóng chảy bằng W = d d ξ τ = ( )os o q l C T Tρ + −⎡ ⎤⎣ ⎦ , [m/s] Tr−ờng nhiệt độ trong lớp vỏ bảo vệ cho bởi: ( )( ) ( ) ( ) ( ) ( ) o s o o s o o s o q CT T T exp T l C T T T y, qy l C T T ⎧ ⎧ ⎫− ξ⎪ ⎪ξ = − +⎪ ⎨ ⎬λ + −⎡ ⎤⎪ ⎪ ⎪⎣ ⎦⎩ ⎭τ ⎨ τ⎪ ξ = −⎪ ρ + −⎡ ⎤⎣ ⎦⎩ ( ) ( ) ( ) ( )o os o os o s o q C qT y, T T exp y T l C T T l C T T ⎧ ⎫⎡ ⎤− τ⎪ ⎪τ = − − +⎢ ⎥⎨ ⎬λ + − ρ + −⎡ ⎤ ⎡ ⎤⎢ ⎥⎪ ⎪⎣ ⎦ ⎣ ⎦⎣ ⎦⎩ ⎭ 7.5.5. Tính l−ợng nhiệt dẫn vào vỏ tàu Mục đích của lớp bảo vệ là khử bỏ phần lớn nhiệt l−ợng sinh ra do ma sát. Phần nhiệt còn lại sẽ dẫn vào trong, làm tăng nội năng của lớp bảo vệ còn lại và dẫn tiếp vào thành tàu, phần nhiệt này bằng: với , hoặc cụ thể hơn, là: 143 qv = 0 T ξ= ∂−λ ∂ξ = ρCW (To - Ts) hay qv = ( ) ( )o s os o q C T T l C T T − + − , → ( ) ( )s ovo s o C T Tq q l C T T −= + − Công thức trên cho thấy nếu chọn vật liệu có nhiệt nóng chảy lớn, l ↑, thì dòng nhiệt thừa qv sẽ nhỏ. 7.5.6. Xác định chiều dày an toàn của lớp cách nhiệt nóng chảy Gọi thời gian con tàu cần bay trong khí quyển là τ. Để chuyến bay an toàn, chiều dày δ lớp cách nhiệt nóng chảy phải đ−ợc chọn sao cho δ > Wτ, hay δ = kWτ với k > 1 là hệ số dự phòng chọn tr−ớc. δ = k ( )os o q l C T T τ ρ + −⎡ ⎤⎣ ⎦ Nếu liên hệ với biểu thức của qo, ta có: δ = k ( ) [ ] 3 k s o K v , m 2 l C T T ρ τ ρ + −⎡ ⎤⎣ ⎦ Tóm lại, khi thiết kế lớp an toàn nhiệt cho vỏ tàu, phải chọn vật liệu có nhiệt độ nóng chảy Ts ≤ Tk, có nhiệt nóng chảy l lớn, và độ dày δ thoả mãn công thức nêu trên. 144 Mục lục Trang Ch−ơng 1: Mô hình bài toán dẫn nhiệt ........................3 1.1. Định luật Fourier ..............................................................................3 1.1.1. Thiết lập ....................................................................................3 1.1.2. Phát biểu ...................................................................................4 1.1.3. Hệ số dẫn nhiệt .........................................................................4 1.2. Ph−ơng trình vi phân dẫn nhiệt...........................................................4 1.2.1. Định nghĩa .................................................................................4 1.2.2. Thiết lập .....................................................................................4 1.2.3. Các dạng đặc biệt của ph−ơng trình vi phân dẫn nhiệt .........5 1.3. Các điều kiện đơn trị (ĐKĐT) ................................................................6 1.3.1. Định nghĩa .................................................................................6 1.3.2. Phân loại các ĐTĐT..................................................................6 1.3.3. Các loại điều kiện biên (ĐKB) .................................................6 1.3.4. ý nghĩa hình học của các loại ĐKB .........................................7 1.4. Mô hình một bài toán dẫn nhiệt...........................................................8 Ch−ơng 2: các Ph−ơng pháp giải tích .........................10 2.1. phép chuẩn hoá và định lý hợp nghiệm...............................................10 2.1.1. Nội dung cơ bản của các ph−ơng pháp giải tích ..................10 2.1.2. Ph−ơng trình vi phân thuần nhất và không TN ...................10 2.1.3. Nguyên lý hợp nghiệm ............................................................10 2.1.4. Phép chuẩn hoá .......................................................................11 2.2. ph−ơng pháp tách biến fourier .........................................................12 2.2.1. Nội dung ph−ơng pháp tách biến Fourier ............................12 2.2.2. Cách giải các bài toán thuần nhất .........................................12 145 2.2.3. Ví dụ: Bài toán làm nguội tấm phẳng biên (W2+W3)...........12 2.3. Ph−ơng pháp nghiệm riêng ổn định ...................................................14 2.3.1. Phạm vi sử dụng ph−ơng pháp NROĐ .................................14 2.3.2. Nội dung ph−ơng pháp NROĐ ..............................................14 2.3.3. Ví dụ: Bài toán gia nhiệt vách phẳng biên (W1) ..................14 2.4. Ph−ơng pháp biến thiên hằng số........................................................16 2.4.1. Phạm vi sử dụng ......................................................................16 2.4.2. Nội dung ph−ơng pháp BTHS................................................17 2.4.3. Bài toán tấm phẳng biên (W2 + W20) .....................................17 2.5. Ph−ơng pháp Fourier cho bài toán không ổn định nhiều chiều .............20 2.5.1. Ph−ơng pháp tách biến lặp .....................................................20 2.5.2. Ph−ơng pháp quy về các bài toán 1 chiều .............................23 2.5.3. Định lý giao nghiệm ................................................................25 Ch−ơng 3: ph−ơng pháp toán tử phức và các bài toán dao động nhiệt .............26 3.1. Bài toán dao động nhiệt ...................................................................26 3.1.1. Khái niệm dao động nhiệt ......................................................26 3.1.2. Mô hình một bài toán dao động nhiệt ...................................26 3.2. Ph−ơng pháp toán tử phức hay tổ hợp phức (Complex Combination) ....27 3.2.1. Nội dung ph−ơng pháp toán tử phức (TTP)..........................27 3.2.2. Các b−ớc của ph−ơng pháp toán tử phức..............................27 3.3. Bài toán dao động nhiệt trong vật bán vô hạn ...................................28 3.3.1. Phát biểu và mô hình (Nh− mục 3.1.2) ...................................28 3.3.2. Giải bằng ph−ơng pháp THP .................................................28 3.3.3. Khảo sát sóng nhiệt .................................................................29 3.4. Dao động nhiệt không ổn định trong vách mỏng ................................31 3.4.1. Đặt vấn đề ................................................................................31 146 3.4.2. Phát biểu bài toán ...................................................................32 3.4.3. Phân tích bài toán (θ)..............................................................33 3.4.4. Nghiệm riêng ổn định .............................................................35 3.4.5. Nghiệm riêng không ổn định..................................................35 3.4.6. Nghiệm riêng dao động...........................................................37 3.4.7. Kết luận....................................................................................39 Ch−ơng 4: ph−ơng pháp toán tử laplace.................41 4.1. Nội dung ph−ơng pháp toán tử Laplace..............................................41 4.1.1. Ph−ơng pháp toán tử ..............................................................41 4.1.2. Phép biến đổi Laplace thuận..................................................41 4.1.3. Phép biến đổi Laplace ng−ợc..................................................42 4.1.4. Các b−ớc của ph−ơng pháp Laplace giải một hệ ph−ơng trình vi phân .........................................43 4.2. Ph−ơng pháp toán tử cho bài toán vách phẳng biên W1.......................43 4.3. Ph−ơng pháp toán tử tìm (x,f) trong vật bán vô hạn ...........................45 4.3.1. Phát biểu bài toán ...................................................................45 4.3.2. Mô hình BT..............................................................................45 4.3.3. Giải bằng ph−ơng pháp toán tử .............................................45 Ch−ơng 5: ph−ơng pháp SAI PHÂN HữU HạN ................47 5.1. Nội dung và các b−ớc áp dụng ph−ơng pháp sai phân hữu hạn .............47 5.1.1. Nội dung FDM.........................................................................47 5.1.2. Các b−ớc áp dụng FDM..........................................................47 5.1.3. Phạm vi sử dụng FDM ............................................................48 5.2. Dạng sai phân của các đạo hàm theo toạ độ.......................................48 5.2.1. Phép sai phân toán học ...........................................................48 5.2.2. Phép sai phân vật lý ................................................................50 5.3. Các ph−ơng pháp xấp xỉ đạo hàm theo thời gian .................................51 147 5.3.1. Ph−ơng pháp Euler ................................................................51 5.3.2. Ph−ơng pháp ẩn (Implicit) .....................................................51 5.3.3. Ph−ơng pháp Crank-Nicolson................................................52 5.3.4. Ph−ơng pháp tổng quát ..........................................................52 5.4. Các ph−ơng pháp giải hệ ph−ơng trình đại số tuyến tính.....................53 5.5. FDM cho bài toán KOĐ một chiều tổng quát.........................................54 5.6. FDM giải bài toán không ổn định 2 chiều tổng quát ............................57 5.6.1. Phát biểu bài toán ..................................................................57 5.6.2. Mô hình TH ............................................................................58 5.6.3. Giải bằng ph−ơng pháp sai phân hữu hạn ...........................58 5.7. FDM giải bài toán không ổn định 3 chiều t(x,y,z,τ) ...............................62 5.7.1. Trong tọa độ vuông góc xyz ...................................................62 5.7.2. Ph−ơng pháp sai phân hữu hạn trong toạ độ trụ (r,ϕ,z)......63 5.8. FDM cho bài toán biên phi tuyến ........................................................66 5.8.1. Điều kiện biên phi tuyến tính .................................................66 5.8.2. Bài toán biên phi tuyến ..........................................................66 5.8.3. Bài toán truyền nhiệt qua vách phi tuyến .............................69 Ch−ơng 6: ph−ơng pháp phần tử hữu hạn, finite element method (FEM) ....................71 6.1. Nội dung và các b−ớc của ph−ơng pháp phân tử hữu hạn.....................71 6.1.1. Nội dung FEM .........................................................................71 6.1.2. Các b−ớc áp dụng FEM ..........................................................71 6.1.3. Phạm vi ứng dụng FEM .........................................................73 6.2. Cực tiểu của hàm nhiều biến và phép xấp xỉ tích phân .........................73 6.2.1. Điều kiện cực tiểu hàm số u = u(x1, x2,...,xn) ..........................73 6.2.2. Phép xấp xỉ tích phân .............................................................74 6.3. Lý thuyết biến phân (variation Theory) ..............................................75 148 6.3.1. Phiếm hàm ...............................................................................75 6.3.2. Nội dung của lý thuyết biến phân ..........................................76 6.3.3. Biến phân của phiếm hàm ......................................................77 6.3.4. Định lý Euler -Lagrange ........................................................78 6.4. Ví dụ minh hoạ các b−ớc áp dụng FEM ................................................85 6.4.1. Bài toán biên cô lập................................................................85 6.4.2. Phát biểu biến phân: ( Variational Statement).....................85 6.4.3.Phát biểu phần tử hữu hạn (Finite Element Formulation)...86 6.5. Bài toán biên TĐN W2 + W3..................................................................95 6.5.1. Phát biểu và mô hình ..............................................................95 6.5.2. Phát biểu biến phân ................................................................95 6.5.3. Phát biểu FEM ........................................................................96 6.5.4. Phát biểu sai phân (theo Euler) .............................................97 6.6. Ph−ơng pháp phần tử hữu hạn giải bài toán không ổn định 2 chiều t(x, y,τ ) với biên cô lập...................................99 6.6.1. Phát biểu mô hình ...................................................................99 6.6.2. Phát biểu biến phân ................................................................99 6.6.3. Phát biểu theo phần tử hữu hạn...........................................100 6.6.4. Phát biểu sai phân .................................................................106 6.6.5. Ví dụ áp dụng cụ thể .............................................................106 6.7. Ph−ơng pháp phần tử hữu hạn giải bài toán không ổn định t (x,y,τ) tổng quát .......................................................107 6.7.1. Phát biểu mô hình .................................................................107 6.7.2. Phát biểu biến phân ..............................................................108 6.7.3. Phát biểu phần tử hữu hạn ...................................................109 6.7.4. Tính đạo hàm theo [t] của Iλ và IC .....................................109 6.7.5. Tính dIg/d[t] ...........................................................................109 149 6.7.6. Tính dIα/d[t] ..........................................................................110 6.7.7. Tính dIq/d[t] ..........................................................................112 6.7.8. Phát biểu phần tử hữu hạn ...................................................113 6.7.9. Phát biểu sai phân để đ−a về hệ ph−ơng trình đại số ........113 6.8. Ph−ơng pháp phần tử hữu hạn giải bài toán có λ = λ(t) .....................114 6.9. Ph−ơng pháp phần tử hữu hạn giải bài toán biên phi tuyến ................117 6.10. Ph−ơng pháp phần tử hữu hạn giải bài toán 3 chiều t (x,y,z,τ) không ổn định tổng quát ....................................................118 6.10.1. Phát biểu bài toán ...............................................................119 6.10.2. Phát biểu biến phân ............................................................119 6.10.3. Phát biểu FEM ....................................................................121 6.10.4. Phát biểu sai phân ..............................................................125 Ch−ơng 7: Các bài toán biên di động .........................127 7.1. Mô tả bài toán biên di động ............................................................127 7.1.1. Khái niệm biên di động.........................................................127 7.1.2. Cân bằng nhiệt trên biên chuyển pha .................................127 7.1.3. Mô hình TH bài toán biên di động ......................................129 7.2. Bài toán biên hoá rắn.....................................................................129 7.2.1. Phát biểu bài toán đóng băng vùng đất −ớt ........................129 7.2.2. Phát biểu mô hình .................................................................130 7.2.3. Giải bằng ph−ơng pháp Stefan ............................................130 7.2.4. Tính gần đúng trong kỹ thuật ..............................................134 7.2.5. Tính độ dày lớp băng tại thời điểm τ ...................................136 7.2.6. Tính thời gian đóng băng đến độ dày đã cho ξ = L ...........136 7.3. Bài toán đông lạnh các vật ẩm hữu hạn ...........................................137 7.3.1. Mục đích chủ yếu khi tính đông lạnh ..................................137 150 7.3.2. Các giả thiết ...........................................................................137 7.3.3. Tính thời gian làm đông τo....................................................137 7.3.4. So sánh thời gian τo ...............................................................138 7.4. Bài toán đông kết vật đúc ...........................................................139 7.4.1. Phát biểu bài toán .................................................................139 7.4.2. Tính ξ(τ) và tốc độ đông kết .................................................139 7.5. Tính truyền nhiệt khi nóng chảy lớp bảo vệ vỏ phi thuyền có vận tốc lớn .........................................................140 7.5.1. Vấn đề bảo vệ nhiệt cho vỏ phi thuyền................................140 7.5.2. Phát biểu bài toán lớp nóng chảy........................................141 7.5.3. Xác định tr−ờng nhiệt độ trong lớp nóng chảy...................141 7.5.4. Xác định vận tốc nóng chảy .................................................142 7.5.5. Tính l−ợng nhiệt dẫn vào vỏ tàu ..........................................142 7.5.6. Xác định chiều dày an toàn của lớp cách nhiệt nóng chảy 143 151 152