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

152 trang ngocly 1830

Download

Bạn đang xem 20 trang mẫu của tài liệu "Giáo trình Các phương pháp tính truyền nhiệt", để tải tài liệu gốc về máy bạn click vào nút DOWNLOAD ở trên

Tài liệu đính kèm:

giao_trinh_cac_phuong_phap_tinh_truyen_nhiet.pdf

Nội dung text: Giáo trình Các phương pháp tính truyền nhiệt

Đại học Đà Nẵng Tr−ờng Đại học bách khoa Khoa công nghệ nhiệt điện lạnh PGS, TS. Nguyễn Bốn Các ph−ơng pháp tính truyền nhiệt - Đà Nẵng - 2001 -
Ch−ơng 1: Mô hình bài toán dẫn nhiệt 1.1. Định luật Fourier 1.1.1. Thiết lập Tính nhiệt l−ợng δQ dẫn qua mặt dS ở cách 2 lớp phân tử khí có nhiệt độ z T1 λ λ T2 T1 > T2 một đoạn bằng quãng đ−ờng tự do trung bình λ . y * Vì T1 và T2 sai khác bé, nên coi mật độ phân tử no và vận tốc trung bình x O ωr các phân tử trong hai lớp nh− nhau. H1. Để chứng minh Do đó, trong thời gian dτ, số phân tử ở định luật Fourier T1 và T2 qua dS là nh− nhau, bằng: 1 d2n = n ω dS dτ 6 o * L−ợng động năng qua dS từ T1 và T2 là: 1 i d2E = E d2n = n ω dS dτ kT 1 1 6 o 2 1 1 i d2E = E d2n = n ω dS dτ kT 2 2 6 o 2 2 Trừ hai ph−ơng trình cho nhau, ta đ−ợc: 1 ik δ2Q = ( E - E )d2n = n ω dSdτ (T - T ) 1 2 6 o 2 1 2 ⎛ dT ⎞ Vì T - T = - ⎜ ⎟ . 2 λ nên 1 2 ⎝ dx ⎠ i dT δ2Q = - n k ϖ λ dS dτ 6 o dx i i R 1 à iR 1 Do n k = n = (n ) ( ) = ρc nên 6 o 6 o N 3 o N 2à 3 o 3
1 dT dT δ2Q = - ( ρc ω λ ) dS dτ = - λ dS dτ 3 o dx dx δ2Q ⎛ ∂T ⎞ hay = q = - λ ⎜ ⎟ dSdτ ⎝ ∂x ⎠ * Khi dS có vị trí bất kỳ thì q = - λ gradT hay dạng vectơ dòng nhiệt là qr = - λ grardT 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ơ: qr = - λ grardT 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ó: 1 1 p ⎛ 8kT ⎞ ⎛ kT ⎞ ⎛ ⎞ ⎜ ⎟ ⎜ ⎟ λ = ρ ω λ cv = ⎜ ⎟ ⎜ ⎟ ⎜ 2 ⎟ Cv 3 3 ⎝ RT ⎠ ⎝ πm ⎠ ⎝ 2.πd p ⎠ 3 2 c v k T = cho thấy: λ không phụ thuộc p, và λ↑ khi T↑ 3 d 2 π3m 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 z 1.2.1. Định nghĩa: Ph−ơng trình vi phân dẫn nhiệt là ρ dV C ph−ơng trình cân bằng nhiệt cho một q λ λ phân tố dv bên trong vật. qω qω V 1.2.2. Thiết lập y x Luật cân bằng nhiệt cho dV ∈ V là: O H2. CBN cho dV 4
[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: ∂t q dVdτ - divqr dVdτ = ρdV.c dτ v p ∂τ ∂t q 1 hay = v - div qr , trong đó dòng nhiệt qua dV là: ∂τ ρcp ρcp r r r r r q = q λ + q ω = - λgradt + ρωcpt, r r r do đó: divq = div (ρcp ωt- λgradt ), coi (ρ, cp) = const ta có : r r r divq = ρcp div (tω) - div (λgradt ) r r r r r r = ρcp (tdiv ω + ω gradt ) - λdiv (gradt )- gradt . gradλ r r r 2 r r = ρcp (tdiv ω + ω gradt ) - λ∇ t - gradt . gradλ Vậy ph−ơng trình có dạng: ∂t qv λ = - tdiv ωr - ωr .grardt + ∇2t + (grardt grard λ)/ρc ∂τ ρc p p cpρ ∂t ∂t dt dx dt dy dt dz dt do + ωr . grardt = + . + . + . = ∂τ ∂τ dx dτ dy dτ dz dτ dτ nên ph−ơng trình vi phân dẫn nhiệt sau khi đặt a = λ , sẽ là: ρCp dt q v 1 = a∇2t + + grardt grardt λ) - tdivωr , với: dτ ρc p ρcp r r r grardt . gradλ là tích vô h−ớng của 2 vectơ gradt và grad λ, ∇2t = ∆t là toán tử Laplace của nhiệt độ, có dạng: ⎧ ∂ 2t ∂ 2t ∂ 2t ⎪ 2 + 2 + 2 (trong(trontọag tọax, y độ, z) vuông góc (xyz)) ⎪∂x ∂y ∂z 2 2 2 2 ⎪ ∂ t 1 ∂t 1 ∂ t ∂ t ∇ t = + ⋅ + ⋅ + (trong r,ϕ, z) ⎨ 2 2 2 2 (trong tọa độ trụ (rϕz)) ⎪∂r r ∂r r ∂ϕ ∂z ⎪ ∂ 2t 2 ∂t ∂ 2t cosθ ∂t ∂ 2t ⎪ + + + + (trong r,θ,ϕ) 2 2 2 2 2 2 2 θϕ ⎩⎪∂r r ∂r r ∂θ r sinθ ∂θ r sin θ∂ϕ (trong tọa độ cầu (r )) 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: 5
∂t 2 q 1 r = a∇ t + v + grardt . gradλ ∂τ ρcp ρcp ∂t q * Vật rắn có λ = const ∀xyz ph−ơng trình là: = a∇2t + v ∂τ ρc p ∂t * Vật rắn có λ = const , ổn định nhiệt = 0, ph−ơng trình là: ∂τ q a∇2t + v = 0. Nếu không có nguồn nhiệt, q = 0, thì ∇2t = 0. λ v 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 ∀M (x, y, z) ∈ V grardt = f(M, τ, t) ∀τ ∈ ∆τ 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: 6
t = t (M1, τ), ∀M1∈ W1, ∀τ ∂t 2. ĐKB loại 2: Cho biết dòng nhiệt dẫn qua biên: q (M ,τ) = -λ , 2 ∂n ∂t −1 tức cho biết = q (M , τ), ∀M ∈ W , ∀τ. ∂n λ 2 2 2 ∂t Khi = q = 0 tức biên W đ−ợc cách nhiệt tuyệt đối hoặc là ∂n 2 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 : ∂t(M ) ∂t (M ) -λ 4 = λ 4 4 và t(M ) = t (M ) ∂n 4 ∂n 4 4 4 5. ĐKB loại 5: Cho biết luật cân bằng nhiệt trên biên W5 di động, t do có sự chuyển pha, trao đổi chất (khối l−ợng thay đổi) hoặc đang biến dạng: dx5 dτ dx ∂t ∂t' ∂t(M5 ) 5 ∂t' -λ -λ -λ = rcρ - λ' (M5), ∂x ∂n ∂n dτ ∂n dx5 -rcρ dx dτ xx với r = nhiệt chuyển pha; 5 c dτ 0 x5 H3. CBN trên biên W 5 = 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, 7
tuỳ theo cách cho ĐKB, sẽ có các đặc điểm hình học sau đây: Đ−ờng cong W Cách cho ĐKB ý nghĩa hình học t(M,τ) t w t(M) đi qua một điểm cố 1 tw = const Mo định Mo ∈W V x t ∂t w q = 0 t(M) đạt cực trị trên W cách 2 = 0 V ∂n β = 0 nhiệt x W ∂t Các tiếp tuyến của t(M) tại w = const β = const ∂n V W song song, góc β = const x W ∂t Các tiếp tuyến của t(M) tại w t f − t w R 3 = tf λ ∂n λ / α λ W3 qua điểm R( , tf) V α x α W ∂t λ ∂t w = 4 ow γ Vo t(M) liên tục, không khả vi 4 ∂n λ ∂x tại W4 và γ = const tW = t4W V x ∂t dx -λ w = r ρ 5 W di chuyển với tốc độ ∂n e dτ 5 5 dx 5 δt'w dx ω = 5 dτ - λ V dτ x δn 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 8
trình vi phân (t), gồm ph−ơng trình ' t vi phân DN và các ph−ơng trình mô ∂ n ' ∂ λ t (M, ) − tả các ĐKĐT nh− sau: w1 τ W5 ∂t λ W1 − ∂n ρ, c, λ, q dx q v rcf ∂t v τ 2 ∂w2 d = a∇ t + và q(M,τ ) -1 ∂x ∂τ ρc (t) = M W2 ∂t 2 q −λ ∂t các ph−ơng trình mô tả các ĐKĐT. ∇ v ∂τ = a t + ∂n −λ ∂to ρc o t w ∂n ∂ x Mục đích chính của truyền −λ ∂ W4 W3 f] t nhiệt là tìm các ph−ơng pháp giải - w 2 [t hệ (t) để tìm hàm phân bố t(x,y,z,τ) α thoả mãn hệ (t). H5. Mô hình 1 bài toán DN 9
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 λ 2 q v −α tτ = a∇ t + , tx (L, τ) = [t(L, τ) - tf] là không TN ρc λ 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 n (tức ph−ơng trình vi phân và các ĐKB thuần nhất), thì t = ∑Ci t i cũng là i=1 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 10
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 ∂ 2 t ⎧ t − t f ⎪ = a 2 (TN) (FT) θ = ∂τ ∂x ⎪ t − t ⎪ ⎪ o f ⎪t(x,0) = t (DKD) ⎨ o ⎪ x (t) t (0,τ) = 0 (TN) (W ) Đổi biến ⎨X = ⎪ x o δ ⎪ ⎪ − α ⎪ aτ ⎪t x (δ,τ) = [t(δ,τ) − t f ] (0TN) (W3 ) F = ⎩ λ ⎪ 2 ⎩ δ αδ và đặt B = λ ∂t ∂t ∂θ ∂F a ∂θ thì do = . . = (t - t ) . ∂τ ∂θ ∂F ∂τ o f δ 2 ∂F ∂t ∂t ∂θ ∂X t − t ∂θ = . . = o f . ∂x ∂θ ∂X ∂x δ ∂X ∂ 2 t ∂ ∂t ∂ ∂θ ∂X ∂ 2θ t o − t f t o − t f 2 = ( ) = . ( . ) = ∂x ∂x ∂x ∂X δ ∂X ∂x δ2 ∂X 2 t − t − α t (δ, τ) = o f θ (1, F) = [t (δ, τ) - t ] có dạng TN là x δ x λ f −αδ θ (l, F) = θ [1, F] = Bθ(1,F) x λ 2 2 ∂t a ∂θ ∂ t t o − t f ∂ θ = (t - t ) = a = a . có dạng đơn giản ∂τ o f δ 2 ∂F ∂x 2 δ 2 ∂X 2 2 hơn là ∂θ = ∂ θ . Khi đó bài toán (t) đ−ợc chuyển đổi thành bài toán ∂F ∂X 2 không thứ nguyên (θ) t−ơng đ−ơng, có dạng chuẩn hoá là: 11
⎧∂θ ∂ 2θ ⎪ = 2 ⎪∂F ∂X ⎪θ (X ,0) = 1 (θ) ⎨ ⎪θ (0, F) = 0 (TN) ⎪ x ⎩⎪θ x (1, F) = −Bθ (1, F) (TN) 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 tại x = δ ra môi tr−ờng tf, α. Tìm tr−ờng t (x, τ) 2. Mô hình TH: t t = at o ⎧ τ xx a ⎪t(x,0) = t λ α ⎪ o q = 0 tf (t) t(x, τ ) ⎨t x (0,τ) = 0 ⎪ − α W2 W3 x ⎪t (δ,τ) = [t(δ,τ) − t ] ⎩⎪ x λ f O δ Bằng cách đổi biến: H6. Bài toán (2.2.2) 12
t − t f x aτ αδ θ = , X = , F = 2 , B = sẽ thu đ−ợc hệ ph−ơng trình t o − t f δ δ λ (θ) t−ơng đ−ơng, ở dạng chuẩn hoá: ⎧θ F = θ xx ⎪ ⎪θ (x,0) = 1 (θ) ⎨ ⎪θ x (0, F) = 0 ⎪ ⎩θ x (1, F) = −Bθ (1, F) 3. Tách biến bằng cách tìm nghiệm dạng θ(X,F) = X(x) F(F). X"(X) F"(F) Thay vào θ =θ có X(x) F'(F) = X"(X) F(F) hay = = -k2 F xx X(X) F(F) (do 2 hàm độc lập), chuyển thành 2 ph−ơng trình vi phân th−ờng: 2 ⎪⎧X"(X) +k X(X) = 0 → X(X) = c1 sin kX + c2 cos kX ⎨ 2 −k2F ⎩⎪F'(F) +k F(F) = 0 → F(F) = e −k2F Nghiệm tổng quát là θ(X,F) = (c1sin kX + c2coskX)e 4. Xác định các hằng số theo ĐKĐT −k2F θx(0,F) = 0 → (kc1cos0 + (-kc2sin0)e = 0 → −k2F c1 = 0 và θ (X,F) = c2 coskX e −k2F cotgk θx(1,F) = (-kc2sin0)e = 2 cosk k -Bθ (1,F)= -Bc coske−k F → B 2 sin k k π2π3π4π k = cotgk = , ph−ơng trình này có O k kk B 1 kk2 3 4 5 vô số nghiệm ki, i = 1 ữ n. Các k nghiệm riêng thoả mãn ĐKB có H7. Giải ph−ơng trình cotg k = B −k2F i dạng: θi(X,F) = c2coskiX.e , ∞ −k2F i nghiệm hợp là θ(X,F) = ∑ci cosk i Xe i=1 ∞ - Điều kiện đầu θ(X,0) = 1 → ∑cicoskiX=1 → coski X∑ci cosk i X = 1 1 sin k 1 cosk XdX i cos k X c cos k XdX coskiX → ∫ i = = ∫ i ∑ i i = 0 k i 0 13
1 4sin k 2 2k + sin 2k i cos k XdX i i ci ∫ i = ci → ci = 0 4ki 2k i + sin 2k i Vậy nghiệm bài toán là: ∞ sin k −k2F i i θ(X,F) = 4 ∑ cos(kiX)e i=1 2k i + sin 2k i * Đồ thị θ(X,F) và t(x, τ) có dạng: t θ F=0 1 F = 0 1 t τ = 0 2 o 3 2 4 3 4 5 6 5 τ = ∞ R tf F = ∞ x x O 1 O δ 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à 14
t(0, τ) = 2to. Tìm t(x, τ) * Mô hình TH: t ⎧t τ = at xx ⎪ W' ⎪t(x,0) = t o 1 (t)⎨ t(0,τ) = 2t ⎪ o 2t ⎪t(δ,τ) = t o ⎩ o a,λ ⎧ t − t o ⎪θ = to t W1 ⎪ o t(x, τ ) ⎪ x x Chuẩn hoá bằng cách đặt ⎨X = ⎪ δ δ ⎪ aτ O F = ⎪ 2 ⎩ δ H10. Bài toán (2.3.3) bài toán (t) trở thành dạng chuẩn hoá (θ) nh− sau: ⎧θ=θFxx ⎪ ⎪θ=(X,0) 0 (θ) ⎨ . Ta sẽ giải bài toán (θ) không thuần nhất ⎪θ=(0,F) 1 ⎩⎪θ=(1, F) 0 (0TN) này bằng ph−ơng pháp NROĐ 2. Tìm nghiệm riêng θ của bài toán ổn định: ⎧θ = 0 = θxx → θ = c1X + c 2 ⎪ ( θ ) ⎨θ(1) = 0 = c1 + c 2 → θ =1 − X ⎪ ⎩θ(0) =1= c2 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 F = θF = θxx = v xx − θxx = v xx ⎪ ⎪v(x,0) = θ(X,0) − θ(X) = X −1 (v) ⎨ ⎪v(1,F) = θ(1,F) − θ(1) = 0 − 0 = 0 ⎪ ⎩v(0,F) = θ(0,F) − θ(0) =1−1= 0 (TN) 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à: 15
−k2F v(X,F) = X(X)F(F) = (c1sinkx + c2coskx)e −k2F - Theo ĐKB: v(0,F) = 0 → c2 e = 0 → c2 = 0 −k2F → v(X,F) = c1sinkXe −k2F Theo v(1,F) = 0 ⇒ c1sinke = 0 → sin k = 0 → k = nπ ∞ (nπ)2 F → v(X,F) = ∑c n sin(nπX) e n=1 ∞ - Theo ĐKĐ: v(X,0) = X-1 → X - 1 = ∑cn sin(nπX)→ x=1 1 1 ∞ 1 cn − 2 ∫(X −1)sin(nπX)dX = ∫sin(nπX) ∑cn sin(nπX) dX → - = → cn = → 0 0 n=1 nπ 2 nπ 2 sin(nπX) 2 nghiệm ph−ơng trình (v) là: v(X,F) = - ∑ e (nπ) F . Do đó, π n nghiệm bài toán (θ) đã cho là: θ(X,F) = θ (X) + σ(X,F) 2 ∞ sin(nπX) θ(X,F) = (1-X) - ∑ exp (-n2π2F) π n=1 n * Phân bố nhiệt độ θ(X,F) và t(x,τ) có dạng: θ t 1 θ 2to = t = 1 2 - to - x t o x / F 2 δ = 2 ∞ 1 1 t τ= 0 x o x O F = 0 1 O δ 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. 16
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 khi m ≠ n ∫φn (X)φm (X)dX = ⎨ 0 ⎩c khi m = n ∞ 3. Biểu diễn nghiệm bài toán (θ) ở dạng θ(X,F) = ∑A n (F)φn (X) n=1 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 ∞ 1 θ(X,F)φ (X)dX φ (X)φ (X)dX ∫ m = ∑A n (F) ∫ n m = cAn(F) tức có quan 0 n=1 0 1 1 θ(x,F)φ (X)dX hệ An(F) = ∫ n c 0 d 4. Lập hệ ph−ơng trình th−ờng của A (F) bằng cách tính A (F), n dF n 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) = ∑φn (X)A n (F) n=1 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ó t δ, a, λ, t(x,0) = to, tx (δ, τ) = 0 và tx (0, t t τ) = - o . q = λ o δ δ t (t = - t o ) q = 0 Tìm t(x, τ) x δ tx = 0 t o = t( x ,0) a,λ x O δ 17
⎧t τ = at xx ⎪t(x,0) = t ⎪ o * Mô hình TH: (t) ⎨ t t (0,τ) = − o ⎪ x δ ⎪ ⎩⎪t x (δ,τ) = 0 t − t x aτ chuẩn hoá với θ = o , X = , F = , sẽ có: δ 2 t o δ ⎧θ F = θ xx ⎪ θ x (1, F) = 0 (TN) ⎪ ⎨ δ (θ) θ (0, F) = t (0,τ ) = −1 (0TN) ⎪ x t x ⎪ o ⎩⎪θ (X ,0) = 0 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 F = v xx ⎪ ⎪v x (1,F) = 0 (v) ⎨ (TN) v (0,F) = 0 ⎪ x ⎩⎪v(X,0) = 0 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 ⎧v x (0,F) = 0 → X x (0) = 0 = c1 → X(x) = c2 coskX ⎨ ⎩v x (1,F) = 0 → X x (1) = 0 = −kc2 sin k → k = nπ Do đó có X(X) = cncos (nπX) và hàm số riêng là φn(X) = cos(nπX). ∞ ∞ 3) Để θ(X,F) = ∑A n (F)φn (X) = ∑A n (F)cos(nπX) là nghiệm bài toán n=1 n=1 (θ) 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]: 18
⎧A (F), khi n = 0 1 1 o 2 ⎪ ∫θ(X,F)cos(nπX)dX = An(F) ∫cos (nπX)dX = ⎨1 0 0 A (F),∀n ≠ 0 ⎩⎪2 n Do đó, An(F) phải xác định theo θ(X,F) bởi quan hệ: ⎧ 1 A o (F) = ∫θ(X,F)dX ⎪ 0 (An) ⎨ 1 ⎪ A n (F) = 2∫θ(X,F)cos(nπX)dX,∀n ≠ 0 ⎩ 0 4. Lập ph−ơng trình vi phân th−ờng cho An(F) bằng cách tính d A (F) theo hệ (A ): dF n n dA (F) 1 1 - Khi n=0, o = θ dX = θ dX =θ 1 = θ (1,F) - θ (0,F) = 0 - ∫ F ∫ xx x 0 x x dF 0 0 (-1) = 1 → Ao(F) = F + c1 1 Điều kiện đầu cho Ao(0) = ∫θ(X,0)dX = 0 = c1 ⇒ Ao(F) = F 0 1 1 dA n (F) - khi ∀n ≠ 0, có: = 2 ∫θF cos(nπX)dX = 2 ∫θxx cos(nπX)dX , dF 0 0 1 (phân đoạn tích phân) = 2 {[θ cos(nπX)]|1 + nπ θ sin(nπX )dX }= x 0 ∫ x 0 1 2{1+2π [θsin(nπX) |1 - nπ θcos(nπX)dX]} 0 ∫ 0 A (F) = 2{1-n2π2 n } → ph−ơng trình vi phân cho A (F) là: 2 n 2 2 2 2 A'n = 2 - n π An →A'n +(n π )An = 2 có nghiệm tổng quát An(F) = 2 1 −(nπ)2 F 2 + c1 e . Điều kiện ban đầu cho An(0) = 2 ∫θ(X,0)cos(nπX)dX → (nπ) 0 2 2 2 2 −(nπ)2 F + c1 = 0 → c1 = - , do đó: An(F) = - e (nπ)2 n 2 π2 n 2π2 n 2 π2 5. Vậy nghiệm bài toán (θ) đã cho là: θ(X,F) = Ao(F) + ∞ 2 ∞ cos(nπX) 2 ∞ cos(nπX) ∑A n (F)cos(nπX) , tức: θ(X,F) = F + 2 ∑ 2 - 2 ∑ 2 . n=1 π n=1 n π n=1 n 19
2 2 ∞ 2cos(nπX) 1 2 1 exp(-n π F) hay, do tổng ∑ 2 2 = X - X + , có: n=1 π n 2 3 1 2 1 2 ∞ cos(nπX) 2 2 θ(X,F) = F + ( X - X + ) - 2 ∑ 2 exp(-n π F) 2 3 π n=1 n * Phân bố θ(X,F) và t(x,τ) có dạng: θ t q=0 1 3 q=0 3 2 2 1 1 to τ x F=0 x =0 O 1 O 1δ 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: ∞ cos(nπx / δ) 2 2 1 2 1 3 aτ 2 n π a t(x,τ) = to( 2 x - x+ ) + to[ 2 - 2 ∑ 2 exp (- 2 τ)] 2δ δ 4 δ π n=1 n δ 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. z 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 t(ρϕ , ,0) = g(ρϕ , ) * Phát biểu BT: Cho trụ l=∞ có a, t t1 t a (R, ϕ,τ) = t1 và ĐKĐ bất kỳ t(ρ,ϕ,0) = ϕ R ρ g(ρ,ϕ). Tìm tr−ờng nhiệt độ t(ρ,ϕ,τ). O H16. Bài toán trụ tổng quát 20
⎧ 1 1 t = a(t + t + t ) ⎪ τ ρρ ρ ρ ρ2 ϕϕ ⎪ t(R,ϕ,τ) = t * Mô hình TH: (t) ⎨ 1 ⎪t(ρ,ϕ,0) = g(ρ,ϕ) ⎪ ⎩⎪ t − t ρ a Chuẩn hoá bằng cách đặt: θ = 1 , r = , ϕ' = ϕ, F = τ , ta có: R 2 t1 R ⎧ 1 1 ⎪θF = θrr + θr + θϕϕ r r 2 ⎪ (θ) ⎨θ(1,ϕ,F) = 0 ⎪ g(ρ,ϕ) − t ⎪θ(r,ϕ,0) = 1 = f (r,ϕ) ⎩⎪ t1 * 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) 1 1 F' Wrr 1 Wr 1 Wϕϕ 2 WF' = W F + W F + WϕϕF → = + + = -k rr r r r 2 F W r W r 2 W 2 ⎧F ' + k 2 F = 0 → F(F) = e −k F → ⎪ ⎨ 1 1 2 ⎪Wrr + Wr + Wϕϕ + k W = 0 ⎩ r r 2 2. Tách biến r, ϕ bằng cách đặt 1 1 W(r, ϕ) = R(r)φ(ϕ) → R"φ + R'φ + Rφ" + k2 Rφ = 0 → r r 2 φ" - = r2 R" + r R' + k2r2 = n2 → φ R R ⎪⎧φ"+n 2φ = 0 → φ(ϕ) = Asin nϕ + B cos nϕ ⎨ 2 2 2 2 là ph−ơng trình Bessel cấp n, có ⎩⎪r R"+rR'+(k r − n )R = 0 , nghiệm là: R(r) = cJ (kr) + DY (kr) = cJ (kr) do lim Y (kr) = ∞ n n n r→0 n Trong đó Jn(kr) là hàm Bessel cấp n loại 1, có dạng: i kr n+2i ∞ (−1) ( ) 2 ∞ −x s−1 Jn(kr) = ,Với Γ(s) = e x dx là hàm giai thừa ∑ ∫0 i=1 i! Γ(n + i +1) 21
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 = ∑∑J n (kmnr) (Amnsin nϕ + Bnm cos m0=1n= 2 nϕ) e −kmn F - Anm và Bnm xác định theo điều kiện đầu: ∞ ∞ θ (r, ϕ, 0) = f(r, ϕ) = ∑∑J n (kmnr) (Amnsin nϕ + Bnm cos nϕ). m0=1n= 2π 2π f (r,ϕ)sin 2 nϕdϕ Nhân 2 vế với sinnϕdϕ và lấy ∫(~)dϕ có ∫ = 0 10 4244 4344 đặt = F(r,n) ∞ 2π 2 ∑J n (kmn'r) Amn' ∫sin n'ϕdϕ. Lại nhân hai vế với rJn(kmnr) dr và lấy m=1 10 4243 = π 1 1 1 2 rJ (k r)dr 1 2 ∫(~)dr có ∫ rJ n (k mn r)F(r,n)dr = πAmn' ∫ n mn = πAmn' [ J n+1 (k mn )] 0 0 0 2 2 1 2π Do đó có: Amn = 2 ∫∫rJ n (k mn r) f (r,ϕ) sin nϕ dϕ dr 0 0 πJ n+1 (k mn ) T−ơng tự xác định các hằng số Bmn: 1 1 2π rJ (k r) f (r,ϕ) - khi n = 0 → Bmo = 2 ∫ o mo ∫ dϕ dr nJ1 (kmo ) 0 0 1 2 2π rJ (k r) - ∀n ≠ 0 → Bmn 2 ∫ n mn ∫f (r,ϕ) cosϕ dϕ dx πJ n+1(kmn ) 0 0 * Ví dụ: Cho f(r,ϕ) = 1, tức t(ρ,ϕ,0) = 2t1, có: ⎧Amn = 0, ∀n ⎪ ⎪ ⎧0,∀n ≠ 0 ⎨ ⎪ 1 Bmn = ⎨ 1 ⎪ rJ (k r)dr, khi n = 0 có: ⎪ ⎪πJ 2 (k ) ∫ o mn ⎩ ⎩ 1 mn 0 22
2 2 1 → Bmo = 2 2 [kmorJ1(kmo)] 0 = kmo J1 (kmo ) k mo J1 (k mo ) Lúc này, nghiệm bài toán (θ) là: ∞ 2 2 θ(r,F) = ∑ Jo(kmr)exp (-k n F) m=1 km J1(km ) 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, z t(r,z,0) = ti, tw = to đồng nhất trên biên W ∈ V. Tìm t(r,z,τ) L t * Mô hình TH: o ⎧ 1 1 t t + t + t = t t i o ⎪ rr r r zz a τ t ⎪ ro r ⎨t(ro , z,τ ) = t(r,0,τ ) = t(r, L,τ ) = to , (0TN) (t) O ⎪t(r, z,0) = t ⎪ i H17. BT trụ hữu hạn ⎩ Đổi biến T = t - to, (t) trở thành: ⎧ 1 1 T + T + T = T ⎪ rr r r zz a τ ⎪ (T) ⎨T(ro , z, τ) = T(r,0, τ) = 0 (TN) ⎪T(r, z,0) = (t − t ) ⎪ i o ⎩ * 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 1 1 1 1 DN có dạng: R Z+ R Z + RZ = (R Z+RZ ) → (R + R - R ) Z rr r r zz a τ τ rr r r a τ 23
1 + (Z - Z )R = 0. Ph−ơng trình này đ−ợc thoả mãn khi đồng thời có: zz a τ ⎧ 1 1 R + R − R = 0 ⎪ rr r r a τ Đây là ph−ơng trình vi phân cho 2 bài ⎨ 1 toán 1 chiều ⎪Z − Z = 0 ⎩⎪ zz a τ - 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,τ) → Z(0,τ) = Z(L,τ) = 0 T(r,L,τ) = 0 = R(r,τ) Z(L,τ) } - Đ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à: ⎧ 1 1 ⎧ 1 R + R = R Z = Z ⎪ rr r n a τ ⎪ zz a τ ⎪ ⎪ ⎨R(ro ,τ) = 0 ⎨Z(0,τ ) = Z(L,τ ) = 0 (R) ⎪R(r,0) = (t − t ) và (Z) ⎪Z(z,0) = 1 ⎪ i o ⎪ ⎩ ⎩ Nghiệm của bài toán (R) và (Z) là: 2 aτ −(kmro ) r 2 ∞ o (với km là nghiệm ph−ơng J o (km r)e R(r,τ) = 2(ti - to)∑ trình Jo (kmro) = 0), m=1 (km ro )J o (km ro ) Z aτ sin(2n +1)π −(2n+1)2 π2 4 ∞ 2 Z(Z,τ) = ∑ L e L π n=0 (2n +1) Do đó, nghiệm của bài toán 2 chiều là T(r,Z,τ) = R,Z: Z π 2 ∞ ∞ J (k r)sin[(2n +1)π 2 2 o m −[km +(2n+1) ]aτ 8 L L2 T(r,z,τ) = (ti-to)∑∑ e π m0=1 n= k m ro J1 (k m ro )(2n +1) 2.5.3. Định lý giao nghiệm 2.5.3.1. Định lý 1: 24
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. 25
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 2π t(τ) = t1cosωτ, trong đó ω = , với τo là chu kỳ dao động, [s]. τo 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. ∂t - Khi tr−ờng t(x, y, z, τ) dao động, cả vectơ grardt = nr và dòng o ∂n ∂t nhiệt q = -λ cũng là những hàm tuần hoàn theo thời gian τ. ∂n 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 ∂t bài toán DĐN đ−ợc cho bởi một ph−ơng trình vi phân dạng = 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 26
⎧t f (τ) = t1 cos(ωτ) + t o ⎪ trong đó chứa to nh− nhiệt độ ban đầu. ⎨ α t (0,τ) = − [t (τ) − t(0,τ)] ⎩⎪ x λ f 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 2π có dao động t(0,τ) = to + t1cos τ = to + t1 cos (ωτ), sẽ t−ơng ứng mô t τo t + t cos ωτ hình sau: o1 to ⎧∂t ∂ 2 t = a t τ = atxx ⎪ 2 a, λ x ⎪∂τ ∂x ⎪ t(x, τ ) (t) ⎨t(0,τ) = t o + t1 cos(ωτ) ⎪ t(∞,τ) = lim t(x,τ) = t o ⎪ x→∞ ⎪ ⎩ H18. BT dao động nhiệ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 27
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 = aT ⎪ τ xx (T) ⎨T(0,τ) = t1 cos(ωτ) ⎪T(∞,τ) = lim T(x,τ) = 0 ⎩ x→∞ 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 = aS ⎪ τ xx (S) ⎨S(0,τ) = t1 sin(ωτ) ⎪S(∞,τ) = limS(x,τ) = 0 ⎩ x→∞ 2. Nhân (S) với i rồi cộng (T) để đ−ợc (W) = (T) + i(S) ⎧Wτ = aWxx ⎪ iωτ (W) ⎨W(0,τ) = t1[cos(ωτ) + isin(ωτ)] = t1e ⎪W(∞,τ) = lim W(x,τ) = 0 ⎩ x→∞ 3. Tìm nghiệm bài toán phức bằng cách tách biến: iωτ - Tìm nghiệm ph−ơng trình Wτ = aWxx ở dạng W(x,τ) = X(x)e , iω phải có: X(x)iωeiωτ = aX(x) eiωτ → X"(x) - X(x) = 0 → X(x) = a iω iω iω iω −x x −x x a a a a iωτ C1 e + C2 e do đó W(x,τ) = [C1 e + C2 e ]e iωτ ⎪⎧W(∞,τ) = lim W(x,τ) = 0 = c2∞e → c2 = 0 - Theo ĐKB ⎨ x→∞ iωτ iωτ ⎩⎪W(0,τ) = c1e = t1e → c1 = t1 28
iω Vậy nghiệm phức là W(x,τ) = t exp (-x ) eiωτ 1 a 4. Tìm dạng đại số của W(x,τ): 1 Do có i =i1/2 = (eiπ/2)1/2 = eiπ/4 = (1 + i) nên W(x,τ) có dạng: 2 ω ω 1 ω ω −x . i ωτ −x .(1+i) ωτ −x i(ωτ−x ) i a i W(x,τ) = t1 e a e = t1 e 2 e = t1 e 2a e 2a → ω −x ω ω W(x,τ) = t1 e 2a [cos(ωτ - x ) + isin(ωτ - x ) 2a 2a Vậy nghiệm bài toán (T) là: ω −x ω 2a T(x,τ) = ReW(x,τ) = t1 e cos(ωτ - x ) 2a 3.3.3. Khảo sát sóng nhiệt Đồ thị (t - x) của tr−ờng nhiệt độ t ω −x t(x, τ) = t + t exp (-x ω/ 2a ) cos (ωτ - 2a o 1 B(x)= t1 e to + t1 x ω/ 2a ) có dạng nh− H19 0 0 1 1 t 4 1. Dao động nhiệt có biên độ tắt τo to ω 2 −x τo τo 2a dần theo chiều sâu x, B(x) = t1 e . to - t1 B(x o ) t0-t1 Trị số xo , tại đó có = 1%, gọi là λx B(0) x o khoảng tác dụng của sóng ⇒ 1% = xo ω xo H19. Sóng nhiệt trong vật bán VH e 2a . ln100 aτo Vậy: - Khoảng cách tác dụng xo = = ln100 hay ω π 2a xo = 2,6 aτo → Chu kỳ τo↑ thì xo↑. 2 Ví dụ: Với đất có a = 0,00273m /h nên xo (τo = 24h) = 0,665m xo (τo = 365.24h) = 12,7m 2. B−ớc sóng và tốc độ truyền sóng nhiệt - B−ớc sóng λx là khoảng cách giữa 2 đỉnh sóng liên tiếp, xác định 29
theo ph−ơng trình ωτo cos (ωτo - λx ω/ 2a ) = 1 → λx = = 2 πaτo → Do đó λx↑ ω/ 2a khi τo↑. λ x πa - Tốc độ truyền sóng v = = 2 → Do đó v↓ khi τo↑. τo τo 3. Giá trị gradt: theo phân bố nhiệt độ: π − x aτo 2πτ π t(x,τ) = to + t1 e cos ( -x ), ta có : τo aτo π − x ∂t π aτo 2πτ π 2πτ π gradt = = t1 e [sin( - x ) - cos( -x )] ∂x aτo τo aτo τo aτo → gradt có dạng một dao động tắt dần theo độ sâu x. Ví dụ: π 2πτ 2πτ 2π π 2πτ * gradt|x=0 = -t1 (cos - sin ) = -t1 sin( - ) aτo τo τo aτo 4 τo 4. Dòng nhiệt tại x lúc τ là: π − x ∂t π aτo 2πτ π q(x,τ) = -λ → q(x,τ) = λt1 e [sin( -x ) - ∂x aτo τo aτo 2πτ π cos( - x )] cũng là một dao động tắt dần theo x, ví dụ: τo aτo 2π π 2πτ 2 * q(0,τ) = λt1 sin( - ), [W/m ], aτo 4 τo π 2πτ 1 3 7 q(0,τ) max khi sin ( - ) = ±1, ứng với τ = - τo, τo, τo, 4 τo 8 8 8 - L−ợng nhiệt tích trong vật sau 1/2 chu kỳ là: τo / 2 2π τo 2τ 2 Qτo/2 = ∫q(0,τ)dτ = λt1 . = λt1 o , [J/m ] 0 aτo π aπ 30
3.4. Dao động nhiệt không ổn định trong vách mỏng The problem of Transient Conduction Thought a Plane Wall When Oscillation of the Environment Temperature. by Nguyen Bon Da nang University, 1995. In this theme will be introduced the solution of the problem of nonstacionary conduction in a plane wall, with non-symetric 3rd class boundary conditions and general initial condition, in the environment temperature is oscillating. The introduced general solution may be applied to determine the solutions of special problems and to thermal calculation for a thin homegeneous wall or a thin layer of the finite anybody, when the environment temperature is constant or oscillating. Bài này sẽ đ−a ra lời giải của bài toán truyền nhiệt không ổn định tổng quát qua vách phẳng với điều kiện đầu tuỳ ý, trong môi tr−ờng dao động nhiệt. Nghiệm thu đ−ợc sẽ là hàm phân bố nhiệt độ trong vách theo tọa độ và thời gian, phụ thuộc vào 10 thông số cho tr−ớc tuỳ ý của các điều kiện đơn trị. Kết quả đ−a ra có thể ứng dụng để tính nhiệt khi nung nóng hay làm lạnh các vách mỏng có phân bố nhiệt độ ban đầu bất kỳ, khi mặt vách tiếp xúc môi tr−ờng có nhiệt độ không đổi hoặc dao động. bài toán truyền nhiệt dao động không ổn định trong vách mỏng 3.4.1. Đặt vấn đề Việc tính toán truyền nhiệt trong vách mỏng, là vách có chiều dày nhỏ hơn bán kính cong, th−ờng đ−ợc quy về vách phẳng. Bài toán dẫn nhiệt ổn định trong vách phẳng đã đ−ợc khảo sát đầy đủ với các điều kiện biên, kể cả khi nhiệt độ mặt vách dao động. Bài toán không ổn định trong vách đã đ−ợc khảo sát chủ yếu chỉ với biên 31
loại 1 hoặc 2 biên loại 3 đối xứng, trong đó nhiệt độ môi tr−ờng không đổi. Đề tài này có mục đích tìm lời giải của bài toán truyền nhiệt dao động không ổn định trong vách mỏng có 2 biên loại 3 không đối xứng, trong đó nhiệt độ môi tr−ờng có thể dao động theo thời gian, với điều kiện đầu cho tuỳ ý. Lời giải của bài toán tổng quát này sẽ có nhiều ý nghĩa trong truyền nhiệt, cả về lý thuyết lẫn ứng dụng. Nó có thể đ−ợc áp dụng để tính nhiệt cho các vỏ mỏng hoặc lớp mỏng gần biên các vật hữu hạn, khi nhiệt độ môi tr−ờng không đổi hoặc thay đổi tuần hoàn. 3.4.2. Phát biểu bài toán: 3.4.2.1. Phát biểu truyền nhiệt: Tìm phân bố nhiệt độ t (x, τ) trong vách phẳng rộng vô hạn dày δ có a, λ không đổi và nhiệt độ đầu t (x, 0) = to(x), khi cho mặt x = 0 tiếp xúc môi tr−ờng có nhiệt độ dao động chu kỳ τo theo luật: ⎛ τ ⎞ tf1(τ) = tf1 + t1 cos ⎜2π ⎟ , với hệ số toả nhiệt phức tạp α1, mặt x = ⎝ τo ⎠ δ tiếp xúc môi tr−ờng nhiệt độ tf2 với hệ số α2. 3.4.2.2. Mô hình toán học: Tìm trong miền x ∈ [0, δ] lúc τ ∈(0, ∞) hàm phân bố nhiệt độ t(x,τ) thoả mãn hệ ph−ơng trình vi phân (t) cho bởi: ⎧t τ = at xx ⎪ ⎪ − α1 ⎡ τ ⎤ t x (0, τ) = ⎢t f 1 + t1 cos 2π − t(0, τ)⎥ ⎪ λ τ (t)⎨ ⎣ o ⎦ ⎪ − α t (δ, τ) = 1 []t(δ, τ) − t ⎪ x λ f 2 ⎪ ⎩t(x,0) = t o (x) 32
t ao,λ tf1 a α1 1 t1 τo x) t o( tf2 a2 α2 o δ x H.20. Bài toán t(x,τ) tổng quát 3.4.2.3. Dạng chuẩn hoá của hệ (t): t − t Đ−a về dạng không thứ nguyên bằng phép đổi biến θ = f 2 , t f1 − t f 2 x aτ α δ α δ X = , F = với các thông số biên b = 1 , b = 2 , δ δ2 1 λ 2 λ t1 aτo t ô (x) − t f 2 B = , Fo = 2 và điều kiện đầu θo(X) = , bài toán t f1 − t f 2 δ t f1 − t f 2 (t) sẽ có dạng chuẩn hoá θ nh− sau: ⎧θF = θXX , ∀x ∈(0,1), ∀F ∈(0,∞) ⎪ ⎡ ⎛ F ⎞ ⎤ ⎪θ = − ⎜ π ⎟ + − θ ⎪ X (0,F) b1 ⎢Bcos⎜2 ⎟ 1 (0,F)⎥ (θ)⎨ ⎣ ⎝ Fo ⎠ ⎦ ⎪ ⎪θX (1,F) = −b2θ(1,F) ⎪ ⎩θ(X,0) = θo (X) Đây là hệ ph−ơng trình đạo hàm riêng cấp 2 của θ(X, F) có điều kiện biên phi tuyến không thuần nhất tại x = 0. 3.4.3. Phân tích bài toán (θ): Nghiệm của bài toán (θ) sẽ đ−ợc tìm ở dạng θ(X,F) = θ1(X,F) + θ2(X,F), trong đó θ1(X,F) là nghiệm riêng bài toán dẫn nhiệt không ổn định khi môi tr−ờng không dao động, θ2(X,F) là nghiệm bài toán dao động riêng ổn định, có tâm dao động cố định. Chọn tr−ớc bài toán không ổn định, thoả mãn hệ ph−ơng trình (θ1) nh− sau: 33
⎧θ1F = θ1XX ⎪ ⎪θ1X (0,F) = b1 []θ1 (0,F) −1 (θ1 )⎨ ⎪θ1X (1,F) = −b2θ1 (1,F) ⎪ ⎩θ1 (X,0) = θo (X) Khi đó θ2(X,F) thoả mãn hệ ph−ơng trình (θ2) = (θ) - (θ1), đ−ợc xác định nh− là hiệu các ph−ơng trình t−ơng ứng của hệ (θ) và (θ1), có dạng sau: ⎧θ2F = θ2XX ⎪ ⎡ ⎛ F ⎞⎤ ⎪θ = θ − ⎜ π ⎟ ⎪ 2X (0,F) b1 ⎢ 2 (0,F) Bcos⎜2 ⎟⎥ (θ2 )⎨ ⎣ ⎝ Fo ⎠⎦ ⎪ θ = − θ ⎪ 2X (1,F) b2 2 (1,F) ⎪ ⎩θ2 (X,0) = 0 (Không cần xét) Nghiệm bài toán không ổn định (θ1) sẽ đ−ợc tìm ở dạng θ1(X,F) = θ10(X) - θ11(X,F) trong đó θ10(X) là nghiệm riêng ổn định, thoả mãn điều kiện ổn định đ−ợc chọn là: ⎧θ10F = 0 = θ10XX ⎪ (θ10 )⎨θ10X (0) = b1 []θ10 (0) −1 ⎪ ⎩θ10X (1) = −b2θ10 (1) Còn θ11(X,F) là nghiệm không ổn định, phải thoả mãn hệ ph−ơng trình (θ11) = (θ10) - (θ1) có dạng sau: ⎧θ11F = θ11XX ⎪ ⎪θ11X (0,F) = b1 []θ11 (0,F) (θ11 )⎨ θ θ (1,F) = −b θ (1,F) R1 ⎪ 11X 2 11 1 ⎪ ⎩θ11 (X,0) = θ10 (X) − θ0 (X) θ1 0( x) a , λ 0 1 11 θ θ 2 θ 0 1 1 θ θ 1 2 θ θ θ x 0 x 1 R2 H.21 Phân tích bài toán θ = θ10 - θ11 + θ2 34
Tóm lại, nghiệm bài toán (θ) sẽ có dạng θ(X,F) = θ10(X) - θ11(X,F) + θ2(X,F), trong đó, θ10(X) là nghiệm riêng ổn định của (θ10), mô tả phân bố nhiệt độ tâm dao động khi ổn định, θ11(X,F) là nghiệm bài toán ch−a ổn định (θ11), mô tả độ lệch nhiệt độ tâm dao động so với ổn định, có thể đ−ợc gọi là độ không ổn định, còn θ2(X,F) là nghiệm của (θ2), mô tả dao động riêng của nhiệt độ quanh tâm dao động. 3.4.4. Nghiệm riêng ổn định: Phân bố nhiệt độ khi ổn định với θ10(x) đ−ợc xác định nh− nghiệm của hệ (θ10), sẽ có dạng: ⎧ ⎪ θ (X) = C X + C ⎪ 10 1 2 ⎪ − b1b2 b1 + b1b2 (1− X) ⎨C1 = hay θ10 (X) = ⎪ b1 + b2 + b1b2 b1 + b2 + b1b2 ⎪ b1 + b1b2 ⎪C 2 = ⎩ b1 + b2 + b1b2 3.4.5. Nghiệm riêng không ổn định: 3.4.5.1. Nghiệm tổng quát của hệ (θ11) đ−ợc tìm ở dạng θ11(X,F) = u(X).v(F). Để thoả mãn θ11F = θ11XX phải có u(X).v'(F) = u''(X). v(F) u'' (X) v ' (F) hay = Vì X và F là 2 biến độc lập, nên đẳng thức trên chỉ u(X) v(F) xảy ra khi cùng bằng 1 hằng số, đ−ợc chọn là 1 số thực âm (-k2). Khi u'' v ' ⎧u''+k 2 u = 0 = 2 đó = -k t−ơng đ−ơng 2 ph−ơng trình ⎨ 2 , có nghiệm là u v ⎩v'+k v = 0 2 ⎧v(F) = Ae−k F ⎨ và nghiệm tổng quát của θ11F = θ11XX sẽ có dạng ⎩u(X) = C sin kX + cos kX −k2F θ11(X,F) = A(CsinkX + coskX)e 35
3.4.5.2. Xác định các hằng số k và C: Theo 2 điều kiện biên loại 3 của hệ (θ11) sẽ tìm đ−ợc: ⎧ (b + b )k tgk = 1 2 ⎪ k 2 − b b ⎨ 1 2 b ⎪C = 1 ⎩⎪ k Có vô số nghiệm số ki, ∀i = 1 ữ ∞; và Ci = b1/ki, nên nghiệm của (θ11) sẽ là hợp tất cả các nghiệm riêng ứng với mỗi ki, ∀i = 1 ữ ∞; ∞ ⎛ b ⎞ 2 ⎜ 1 + ⎟ −ki F θ11(X,F) = ∑A i ⎜ sin k i X cosk i X⎟e i=1 ⎝ k i ⎠ y 2 b b 1 1 b b tgk (b1 + b2 )k 2 = = K -=b1 b2 K K k k k=- K 0 K1k1 k K2k2 k3kK3 /2 /2 /2 π π π K= k K=3 K=5 k k H.22. Các nghiệm số ki của ph−ơng trình đặc tr−ng 3.4.5.3. Xác định các hằng số Ai: Theo điều kiện đầu θ11(X,0) = θ10(X) - θo(X) = C1X + C2-θo(X) = ∞ ⎛ b ⎞ ⎜ 1 ⎟ = ∑A i ⎜ sin k i X + cos k i X⎟ i=1 ⎝ k i ⎠ ⎛ b ⎞ ⎜ 1 + ⎟ Sau khi nhân 2 vế với ⎜ sin k i X cos k i X⎟ rồi tích phân theo dX ⎝ k i ⎠ trong khoảng X ∈ [0, 1] nhờ tính trực giao của các hàm riêng, sẽ nhận đ−ợc: 36
1 ⎛ b ⎞ + − θ ⎜ 1 + ⎟ ∫[]C1X C 2 o (X) ⎜ sin k i X cos k i X⎟dX 0 ⎝ k i ⎠ Ai = 2 1 ⎛ b ⎞ ⎜ 1 + ⎟ ∫⎜ sin k i X cos k i X⎟ dX 0 ⎝ k i ⎠ t o − t f1 Nếu cho điều kiện đầu t (x, 0) = to = Const thì θo(X) = = t f1 − t f 2 θo, thay C1, C2 bởi các biểu thức ở mục 4, sẽ tìm đ−ợc: 2 2 b1 (k i − b1b2 )sin k i + k1b1 b2 cosk i + θo k i (b1 + b2 + b1b2 )(b1 cosk i − k i sin k i − b1 ) Ai= 4 2 2 2 2 2 (b1 + b2 + b1b2 )[(k i − b1 )sin2k i + 4b1k i sin k i + 2k i (k i + b1 ) Nghiệm bài toán (θ1) = (θ10) - (θ11) có dạng ∞ b + b b (1− X) ⎛ b ⎞ 2 1 1 2 ⎜ 1 ⎟ −ki F θ1(X,F) = - ∑ A i ⎜ sin k i X + cos k i X⎟e b1 + b2 + b1b2 i=1 ⎝ k i ⎠ Nghiệm này mô tả phân bố nhiệt độ không ổn định của vách trong môi tr−ờng tĩnh, có nhiệt độ bất biến. Đó chính là đ−ờng cong qua các tâm dao động nhiệt. 3.4.6. Nghiệm riêng dao động: 3.4.6.1. Biến đổi qua hệ toán tử phức: Bài toán dao động (θ2) sẽ đ−ợc giải bằng ph−ơng pháp toán tử phức. Nội dung ph−ơng pháp này, là chuyển bài toán thực (θ2) thành bài toán phức (W) = (θ2) + i(S), sau đó tìm nghiệm phức W và xác định nghiệm thực (θ2) nh− là phần thực của nghiệm phức, θ2 = Re(W). Bài toán ảo (S) đ−ợc lập bằng cách thay hàm θ2(X,F) và ⎛ F ⎞ ⎛ F ⎞ ⎜ ⎟ ⎜ ⎟ cos⎜2π ⎟ trong bài toán thực (θ2) bởi hàm ảo S(X,F) và sin⎜2π ⎟ , sẽ ⎝ Fo ⎠ ⎝ Fo ⎠ có dạng: ⎧S F = S XX ⎪ ⎪ ⎡ F ⎤ (S)⎨S X (0,F) = b1 ⎢S(0,F) − Bsin 2π ⎥ ⎪ ⎣ Fo ⎦ ⎪ ⎩S X (1,F) = −b 2S(1,F) 37
Sau khi nhân đơn vị ảo i với hệ (S) rồi cộng với hệ (θ2), theo phép biến đổi (W) = (θ2) + i(S), sẽ có bài toán phức dạng: ⎧WF = WXX ⎪ F ⎡ ⎤ ⎡ i2π ⎤ ⎪ F F Fo (W)⎨WX (0,F) = b1 ⎢W(0,F) − B(cos2π + isin2π )⎥ = b1 ⎢W(0,F) − Be ⎥ ⎪ ⎣ Fo F ⎦ ⎣⎢ ⎦⎥ ⎪ ⎩WX (1,F) = −b2 W(1,F) 3.4.6.2. Tìm nghiệm toán tử phức: Nghiệm tổng quát của ph−ơng trình WF = WXX đ−ợc tìm ở dạng F i2π F iωF 2π tách biến W(X, F) = V(X).e o = V(X).e , với ω = đ−ợc gọi là Fo tần số không thứ nguyên, sẽ có dạng: −X iω X iω iωF W(X,F) = (C1 e + C2 e )e Chọn các hằng số phức C1, C2 thoả mãn 2 điều kiện biên loại 3 trong hệ (W), sau khi đặt u = iω , sẽ thu đ−ợc. ⎧ 2u Bb1.e (u + b2 ) iϕ1 ⎪C1 = 2u = Bb1.r1e ⎪ e (u + b )(u + b ) − (u − b )(u − b ) ⎨ 1 2 1 2 Bb .(u − b ) ⎪C = 1 2 = Bb .r eiϕ2 ⎪ 2 2u 1 2 ⎩ e (u + b1 )(u + b 2 ) − (u − b1 )(u − b2 ) C1 C 2 Trong đó r1, r2 và ϕ1, ϕ2 là modul và argument của số phức , Bb1 Bb1 ⎧ ⎛ C ⎞ ⎛ C ⎞ r mod⎜ 1 ⎟, arg⎜ 1 ⎟ ⎪ 1 = ⎜ ⎟ ϕ1 = ⎜ ⎟ ⎪ ⎝ Bb1 ⎠ ⎝ Bb1 ⎠ ⎨ ⎛ C ⎞ ⎛ C ⎞ ⎪r = mod⎜ 2 ⎟,ϕ = arg⎜ 2 ⎟ ⎪ 2 ⎜ ⎟ 2 ⎜ ⎟ ⎩ ⎝ Bb1 ⎠ ⎝ Bb1 ⎠ Do đó, nghiệm phức của hệ (W) sẽ là hàm phức có dạng: i(ϕl + ωF) −X iω i(ϕl + ωF) X iω W(X.F) = Bb1[r1e e + r2 e e ] 3.4.6.3. Tìm nghiệm thực θ2: 38
Nghiệm của hệ (θ2) sẽ là phần thực của nghiệm phức W(X,F) vừa ω tìm đ−ợc. Để ý rằng u = iω = (1 + i) , nghiệm W(X,F) sẽ có dạng: 2 ω ω ω ω ⎡ −X i(ϕ +ωF−X ) X i(ϕ +ωF+X ) ⎤ 2 1 2 2 2 2 W(X,F) = Bb1 ⎢r1e e + r2e e ⎥ ⎣⎢ ⎦⎥ Do đó, nghiệm riêng dao động của hệ (θ2), đ−ợc tìm theo θ2 = ReW, sẽ là: ω ω ⎡ −X X ⎤ 2 ω 2 ω θ2(X,F) = Bb1 ⎢r1e cos(ϕ1 + ωF − X ) + r2e cos(ϕ2 + ωF + X )⎥ ⎣⎢ 2 2 ⎦⎥ 2πδ2 2π 2u 2ω Khi chu kỳ τo nhỏ, ω = = sẽ khá lớn và mod(e ) = e rất Fo a.τo lớn, nên có thể coi r2 = 0. Lúc đó biên độ dao động A(X) = ω −X 2 Bb1r1 e giảm khi X tăng, nên dao động tắt dần trong vật. Nếu gọi độ sâu thấm nhiệt là Xt, ở đó có A(Xt) = 1%, thì X = 2 ln100. A(0) t ω 3.4.7. Kết luận: 3.4.7.1. Nghiệm của bài toán tổng quát về truyền nhiệt dao động không ổn định trong vách mỏng sẽ có dạng: ∞ b + b b (1− X) ⎛ b ⎞ 2 1 1 2 ⎜ 1 + ⎟ −ki F θ1(X,F) = - ∑A i ⎜ sin k i X cos k i X⎟e b1 + b2 + b1b2 i=1 ⎝ k i ⎠ ω ω ⎡ −X X ⎤ 2 ω 2 ω + Bb1 ⎢r1e cos(ϕ1 + ωF − X ) + r2e cos(ϕ2 + ωF + X )⎥ hay ⎣⎢ 2 2 ⎦⎥ λ (t − t )(δ + − x) f1 f2 α 2 t(x,τ) = tf2 + − (t f − t f ) . ⎛ λ λ ⎞ 1 2 ⎜ + δ + ⎟ ⎝ α1 α 2 ⎠ ∞ ⎛ α δ x x ⎞ ⎛ a.τ ⎞ A ⎜ 1 sin k + cos k ⎟exp − k 2 ∑ i ⎜ i i ⎟ ⎜ i 2 ⎟ + i=1 ⎝ λk i δ δ ⎠ ⎝ δ ⎠ π π α δ ⎡ −X τ π X τ π ⎤ + t 1 ⎢r e aτo cos(ϕ + 2π − x ) + r e aτo cos(ϕ + 2π + x )⎥ 1 λ 1 1 τ a.τ 2 2 τ a.τ ⎣⎢ o o o o ⎦⎥ 39
Đây là hàm tổng quát mô tả phân bố nhiệt độ trong vách, phụ thuộc vào x, τ và 10 thông số cho tuỳ ý của các điều kiện đơn trị: t = τ (x, τ; δ, a, λ, to(x), tf1,ti, τo, α1, tf2, α2). 3.4.7.2. Nếu cho các thông số lấy những giá trị đặc biệt, chẳng hạn cho δ → ∞, (α1 hoặc α1) → (0 hay → ∞), to(x) = (tf1, tf2 hay hàm bất kỳ), tf1 = tf2 > to(x), thì sẽ tìm đ−ợc nghiệm các bài toán có các điều kiện đơn trị khác nhau, nh− là các tr−ờng hợp đặc biệt của bài toán tổng quát. Một vài phân bố t (x, τ) đặc biệt đ−ợc minh hoạ ở hình (H.23) t t t t t01 t = x) τ 8 t 0( = ) τ (x τ = t 0 8 τ τ τ 8 τ τ tf2 τ = t02 8 τ = to x t x o o 8 x o o o xo δ δ δ δ x o δ x 1) τ0 = 0 2) α2 = 0 3) α2 = ∞ 4) t0(x) cắt t(x,∞) 5) α2 = ∞, tf1 > t0 > tf2 H.23. Vài phân bố t(x, τ) đặc biệt tk (τ ) t (x,τ ) t + t đ tmax o 1 to 24h to tmin to - t1 o o o x x 0,25 0,5 0,75 1 τ n N x H.24. Dao động nhiệt độ môi tr−ờng τk(τ) và trong vật τđ(x, τ) theo chu kỳ 24h và 1 năm. 3.4.7.3. Các kết quả nêu trên có thể đ−ợc ứng dụng để khảo sát tính toán quá trình truyền nhiệt không ổn định trong vách hoặc trong lớp biên mỏng của các vật hữu hạn, trong môi tr−ờng có nhiệt độ bất biến hoặc dao động. Ví dụ có thể áp dụng để khảo sát dao động phổ biến của nhiệt độ môi tr−ờng không khí và trong lớp mỏng gần bề mặt của mọi vật trên địa cầu, d−ới tác dụng tuần hoàn của bức xạ Mặt trời theo chu kỳ 24 giờ và chu kỳ năm (Hình 24) 40
Ch−ơng 4: ph−ơng pháp toán tử laplace 4.1. Nội dung ph−ơng pháp toán tử Laplace 4.1.1. Ph−ơng pháp toán tử. Ph−ơng pháp toán tử là ph−ơng pháp giải một hệ ph−ơng trình vi phân bằng cách: - Chuyển hệ ph−ơng trình vi phân thành hệ ph−ơng trình toán tử nhờ một phép biến đổi nào đó, ví dụ phép biến đổi Laplace. Hệ ph−ơng trình vi phân của toán tử này có bậc thấp hơn ph−ơng trình vi phân ban đầu (th−ờng là hệ ph−ơng trình vi phân th−ờng hoặc ph−ơng trình đại số). - Tìm nghiệm toán tử và dùng phép biến đổi ng−ợc xác định nghiệm gốc của bài toán đã cho. 4.1.2. Phép biến đổi Laplace thuận 1. Định nghĩa: Phép biến đổi Laplace thuận là phép biến đổi hàm số f(t) thành hàm fˆ (s), theo quan hệ: ∞ fˆ (s) = ∫e −st f (t)dt t=0 Hàm fˆ (s) gọi là ảnh (hoặc toán tử) Laplace của hàm gốc f(t), ký hiệu là L [f(t)] = fˆ (s) ∞ (a−s)t at −st at e ∞ 1 Ví dụ: L [e ] = ∫e e dt = | t=0 = t=0 a − s s − a 2. Tính chất của phép biến đổi Laplace: Từ định nghĩa trên suy ra các tính chất: ˆ L [c1f(t) + c2g(t)] = c1 f (s) + c2 gˆ (s) ∞ ∞ ' −st ' -st ∞ −st L [f(t)] = ∫e f(t)dt = [f(t)e ] t=0 - s ∫f (t)e dt 0 0 = -f(0) + sfˆ (s) L [f"(t)] = -f'(0) + s[-f(0) + sfˆ (s)] = -f'(0) - sf(0) + s2 fˆ (s) 41
3. Biến đổi Laplace của hàm nhiều biến: - ảnh Laplace của hàm f(x,y) có thể nhận đ−ợc bằng cách tích phân theo x hoặc theo y: ∞ −sx ˆ Lx[f(x,y)] = ∫e f(x,y)dx = f (s;y) x=0 ∞ −sy ˆ Ly[f(x,y)] = ∫e f(x,y)dy = f (x;s) y=0 - ảnh của các đạo hàm của f(x,y) cũng t−ơng tự ∂ ∞ ∂ ∞ e −sx -sx ∞ f (x,y)e−sx dx Lx[ f(x,y)] = ∫ f(x,y)dx] = f(x,y)e | x=0 + s ∫ ∂x x=0 ∂x x=0 = -f(0,y) + sfˆ (s;y) ∞ ∂ −sx ∂ d ˆ Lx[ f(x,y)] = ∫e f(x,y)dx = f (s;y) ∂y x=0 ∂y dy ∂ 2f (x, y) L [ ] = -f (0,y) - sf(0,y) + s2 fˆ (s;y) x ∂x 2 x 2 d 2 L [ ∂ f (x, y) ] = fˆ (x,s) v.v. y ∂x 2 dx 2 4.1.3. Phép biến đổi Laplace ng−ợc: Là phép biến đổi ảnh fˆ (s) cho tr−ớc ra hàm gốc f(t) của nó. - Phép tính này có thể tìm f(t) theo công thức f(t) = 1 x+iω lim est fˆ(s)ds trong đó s là biến phức có cận giữa 2 số phức liên ω→∞ ∫ 2πi s=x−iω hợp s = x - iω đến s = x + iω với Res = x > 0 - Hoặc có thể phân tích fˆ (s) ra các phân thức hữu tỷ dạng fˆ (s) = A As + B ∑c s k + + và dùng bảng tìm các gốc t−ơng ứng k ∑ (s − a) k ∑ (s 2 + a 2 ) k rồi cộng lại. - Tính chất phép biến đổi Laplace ng−ợc (định lý tính chập): Nếu biết f(t), g(t) là gốc của ảnh fˆ(s), gˆ(s) thì gốc của tích các ảnh ∞ hˆ(s)= fˆ(s) gˆ(s) là h(t) = ∫f (σ)g(t − σ)dσ σ=0 4.1.4. Các b−ớc của ph−ơng pháp Laplace giải một hệ ph−ơng 42
trình vi phân 1. Chuyển hệ ph−ơng trình (T) đã cho thành hệ ph−ơng trình toán tử (Tˆ ) bằng cách nhân cả hai vế của các ph−ơng trình (T) với e-sF rồi tích phân theo F từ 0 đến ∞. Hệ ph−ơng trình (Tˆ ) là hệ ph−ơng trình vi phân th−ờng của toạ độ, khi T = T(x,τ). 2. Tìm nghiệm toán tử Tˆ (x;s) của hệ (Tˆ ). 3. Biến đổi Laplace ng−ợc để tìm hàm gốc T(x,τ) nh− là nghiệm của ph−ơng trình đã cho. Phạm vi sử dụng ph−ơng pháp này chủ yếu phụ thuộc vào khả năng thực hiện phép biến đổi ng−ợc Laplace. 4.2. Ph−ơng pháp toán tử cho bài toán vách phẳng biên W1: Để làm ví dụ minh hoạ ph−ơng pháp toán tử, ta giải bài toán vách phẳng biên loại 1 cho ở mục (2.3.3): ⎧θF = θxx ⎪ F ⎪θ(X,0) = 0 (θ)⎨ ⎪θ(0,F) =1 ⎩⎪θ(1,F) = 0 θ(0,F)=1 θ (1,F)=0 θ = θ 1. Chuyển hệ (θ) thành hệ Fxx θ (x,0)=0 toán tử (θˆ ) bằng biến đổi Laplace x theo thời gian F. O 1 H21. Bài toán (2.3.3) ∞ ∞ θ e −sFdF θ e −sF dF - Ph−ơng trình θF = θxx có dạng: ∫ F = ∫ xx → F=0 F =0 ∞ 2 −sF ∞ −sF d ˆ ˆ [θ (X , F)e ]F =0 + s ∫e θ(X,F)dF = 2 θ (X,s) hay [0-θ(X,0)] + sθ (X,s) = F=0 dx d 2 θˆ (X,s) dx 2 d 2 → ph−ơng trình vi phân toán tử là θˆ (X,s) - sθˆ (X,s) = θ (X,0)= 0 dx 2 (do điều kiện đầu θ(X,0) = 0) - Các ĐKB có dạng: 43
∞ ∞ 1 θ(0,F) = 1 → ∫θ(0,F)e−sFdF = ∫ e −sF dF → θˆ (0,s) = . F=0 0 s ∞ θ(1,F) = 0 → ∫θ(1,F)e−sFdF = 0 → θˆ (1;s) = 0 F=0 ⎧θˆ"−sθˆ = 0 ⎪ ⎪ 1 Do đó, hệ ph−ơng trình toán tử có dạng (θˆ )⎨θˆ(0,s) = ⎪ s ˆ ⎩⎪θ(1,s) = 0 2. Tìm nghiệm θˆ (X,s): Ph−ơng trình θˆ "(X,s) - sθˆ (X,s) = 0 có nghiệm tổng quát là θˆ (X,s) = Ash(x s ) + Bch(x s ) 1 θˆ (0,s) = = B s − ch s θˆ (1,s) = 0 = Ash s + Bch s → A = . ssh s Nghiệm toán tử là: sh s ch(x s) − ch s sh(x s) θˆ (X,s) = s.sh s 1 sh[(1 − X ) s] θˆ (X;s) = . s sh s 3. Tìm nghiệm gốc bằng phép biến đổi ng−ợc: 1 Theo bảng fˆ(s) f(X) có 1 g(X,F) và s ∞ sh[(1− X) s] 2 2 -2π ∑(−1)n nsin[nπ(1− X)] . e −n π F f(X,F) sh s n=1 Theo định lý tính chập (4.1.3) ảnh: θˆ (X,s) = fˆ(X;s)gˆ(X;s) có gốc là: F θ(X,F) = ∫f (X,σ)g(X,F − σ)dσ σ=0 F ∞ 2 2 = ∫ {−2π∑(−1)2 nsin[nπ(1− X)]e−n π σ}.dσ σ=0 1 ⎛ −n2π 2 F ⎞ = -2π∑(-1)nnsin[nπ(1-X)] . ⎜ e −1⎟ . Vì sin [nπ(1-X)] ⎜ 2 2 ⎟ ⎝ − n π ⎠ = sin nπ cos nπX - cosnπsinnπX = 0 - (-1)n sin nπX, và vì 2 ∞ 1 2 ∞ sin(nπX) ∑ sin(nπX) = 1-X nên θ(X,F)= (1-X) - ∑ exp (-n2π2F). π n=1 n π n=1 n Khi ổn định, F → ∞ có θ(X) = 1-X 44
4.3. Ph−ơng pháp toán tử tìm (x,f) trong vật bán vô hạn 4.3.1. Phát biểu bài toán: τ Cho vật bán vô hạn x ≥ 0 có nhiệt độ đầu phân bố đều t(x,0) = t và biên loại 1 tại x = 0 là t(0, τ) o t1 t(x→∞ ) tτ = atxx = to = t1, ∀τ > 0. Tìm t(x, τ) trong vật x O t = to H22. BT vật bán vô hạn 4.3.2. Mô hình BT: ⎧t = at ⎧T = aT ⎪ τ xx ⎪ τ xx (t) ⎨t(x,0) = t o = limt(x,τ) Đổi biến T=t-to có (T)⎨T(x,0) = 0 = limT(x,τ) ⎪ x→∞ ⎪ x→∞ ⎩t(0,τ) = t1 ⎩T(0,τ) = t1 − t o = To 4.3.3. Giải bằng ph−ơng pháp toán tử: 1. Chuyển hệ (T) thành hệ (Tˆ ): 2 d ˆ s ˆ 1 Tτ = aT → X(x,s) - T(x,s) = T(x,0) xx dx 2 a a limT(x,τ) = 0 →lim Tˆ (x,s) = 0 x→∞ x→∞ ∞ ∞ T -sτ -sτ ˆ o T(0,τ) = To → ∫T(0,τ)e dτ = ∫To e dτ → T (0,s) = . Vậy hệ τ=0 0 s d 2 s (Tˆ ) có dạng: Tˆ (x,s) - Tˆ (x,s) = 0 dx 2 a (Tˆ ) lim Tˆ (x,s) = 0 x→∞ T Tˆ (0,s) = 0 s 2. Tìm nghiệm Tˆ (x,s): s s s Tˆ " - Tˆ = 0 → Tˆ (x,s) = A exp (-x ) + B exp (x ) a a a lim Tˆ (x,s) = 0 = B x→∞ 45
T Tˆ (0,s) = o = A + B = A s T s Vậy nghiệm toán tử là Tˆ (x,s) = o exp(-x ) s a T s 3. Xác định nghiệm gốc T(x,τ): Tra bảng có o exp (-x ) s a x To.erfc ( ). Vậy nghiệm bài toán đã cho là T0(x,τ) = To erfc 2 aτ x T(x,τ) x ( ) hay = θ(x,τ) = erfc ( ), trong đó, erfc(x) là hàm sai 2 aτ To 2 aτ số bù Gauss, Gauss' complementary errow function, xác định bởi: ∞ 2 2 erfc(x) = ∫ e−δ dσ π σ =x Hàm erfc(x) liên hệ với hàm sai số Gauss, Gauss’ errow function, n 2n+1 x ∞ 2 2 2 (−1) x là erf(x) = ∫ e−δ δσ = ∑ theo biểu thức erf(x) + erfc(x) π σ=0 π n=0 n!(2n +1) = 1 n 2n+1 1 ∞ (−1) x Do đó có nghiệm θ(x,τ) = 1 - ∑ n πaτ n=0 n!(2n +1)(4aτ) *Có thể chứng minh rằng θ ∞ 2 −σ2 erf(x) + erfc(x) = ∫e dσ = 1 1 x π σ=0 θ= erfc ( ) ∞ 2 2 aτ nh− sau: Gọi I = ∫e−x dx ta có: 0 ∞ ∞ 2 2 2 2 2 −x −y −(x +y ) I = ∫e dx ∫e dy = ∫∫e dxdy O 12x 0 0 (y,x)∈[0,∞) 2 aτ 2 2 = lime−(x +y ) dxdy . r ∫∫ H23. Phân bố θ(x,τ) →∞ D ⎧x = rcosϕ Đổi biến ⎨ → y ⎩y = rsin ϕ ϕ / 2 r 2 −r2 π −1 −r2 r I = limdϕ e rdr = lim [ ( e )| 0 ] r→∞ ∫ ∫ r→∞ 0 0 2 2 D π π = . Do đó I = ϕ 4 2 O r→∞ x → erf(x) + erfc(x) = 1 H24. Để tính I 46
Ch−ơng 5: ph−ơng pháp SAI PHÂN HữU HạN Ph−ơng pháp sai phân hữu hạn, finite deference method (FDM) và ph−ơng pháp phần tử hữu hạn finite element method (FEM) là những ph−ơng pháp xấp xỉ. FDM là phép xấp xỉ vi phân và FEM là phép xấp xỉ tích phân, với độ chính xác tuỳ ý. Cả hai ph−ơng pháp đều cho phép nhận đ−ợc nghiệm của bài toán ở dạng số. 5.1. Nội dung và các b−ớc áp dụng ph−ơng pháp sai phân hữu hạn 5.1.1. Nội dung FDM T− t−ởng của ph−ơng pháp sai phân hữu hạn là việc thay một cách gần đúng tất cả các vi phân, nguyên là những số gia "vô cùng bé", bằng các số gia "bé hữu hạn", gọi là các sai phân. Nhờ đó có thể chuyển gần đúng mọi ph−ơng trình vi phân thành hệ ph−ơng trình đại số. Hệ ph−ơng trình đại số này, mặc dù số ẩn số có thể rất lớn, luôn giải đ−ợc bằng máy tính. 5.1.2. Các b−ớc áp dụng FDM Khi áp dụng FDM để giải một hệ ph−ơng trình vi phân bất kỳ ta th−ờng tiến hành theo các b−ớc sau: 1. Chia vật khảo sát ra n phần tử bằng cách chia các toạ độ t−ơng ứng ra các phần ∆x, ∆y, ∆z. Nhiệt độ của mỗi phần tử đ−ợc coi là nhiệt độ tại điểm nút của nó. Điểm nút của phần tử bên trong vật là tâm của khối chữ nhật ∆x∆y∆z. Nút của phần tử trên biên nằm giữa mặt biên. Những phần tử trên cạnh hoặc góc của vật có thể có 2 hoặc 3 mặt biên khác loại, th−ờng chọn nút của phần tử loại này tại giữa cạnh hoặc góc các mặt biên. 2. Chuyển hệ ph−ơng trình đạo hàm riêng thành hệ n ph−ơng trình vi phân th−ờng cấp 1. Mỗi ph−ơng trình cấp 1 của hệ này là biểu thức mô tả quan hệ giữa đạo hàm theo thời gian của nhiệt độ nút (ijk), là dtijk với nhiệt độ của 6 nút lân cận nó (nút trong); hoặc với các thông dτ số cho tr−ớc của biên (nút biên). 47
Mỗi ph−ơng trình cấp 1 này có thể nhận đ−ợc bằng cách thay tất cả các đạo hàm theo toạ độ bằng biểu thức sai phân t−ơng ứng hoặc bằng cách viết ph−ơng trình cân bằng nhiệt độ mỗi nút. 3. Chuyển hệ n ph−ơng trình vi phân th−ờng cấp 1 của các đạo dtijk hàm thành hệ n ph−ơng trình đại số bằng cách chia thời gian ra dτ các khoảng ∆τ bằng nhau và dùng các phép xấp xỉ đạo hàm theo thời gian . 4. Lập dạng ma trận cho hệ ph−ơng trình đại số và giải nó nhờ máy tính. Tất cả các b−ớc trên có thể lập trình cho máy tính thực hiện. 5.1.3. Phạm vi sử dụng FDM: - FDM có thể giải mọi bài toán về truyền nhiệt (hoặc khoa học kỹ thuật nói chung) ổn định và không ổn định 1,2 hoặc 3 chiều, với các điều kiện vật lý bất kỳ, điều kiện biên bất kỳ, có nguồn hoặc không có nguồn nhiệt. Nói chung, cho mọi tr−ờng hợp, với các thông số của điều kiện đơn trị (λ,ρ,α,tw,tf, q, qv) = f(x,y,z,τ) đều có thể dùng FDM. - Với vật có biên dạng không quy tắc phải tốn công lớn khi lập cân bằng nhiệt cho các nút biên, khi đó FDM tỏ ra kém hiệu quả. - Khi hệ số dẫn nhiệt λ = λ(t) hoặc khi điều kiện biên phi tuyến tính (ví dụ khi cho biết nhiệt độ môi tr−ờng Tf và dòng nhiệt bức xạ đi 4 4 qua biên theo luật q = εσ0 (T − Tf ),thì việc áp dụng FDM không có gì trở ngại. Những tr−ờng hợp trên ta th−ờng dùng FDM 5.2. Dạng sai phân của các đạo hàm theo toạ độ 5.2.1. Phép sai phân toán học dt 5.2.1.1. Sai phân của dx 48
t Theo định nghĩa đạo hàm của ti+1 hàm t(x) tại xi, ta có: dt t(x + ∆x) − t(xi ) i ti xi = lim , dx ∆x→0 ∆x t Khi chọn ∆x đủ nhỏ, có thể thay i-1 đạo hàm bởi các dạng sai phân sau: ∆x ∆x x x x dt t i+1 − t i O i-1 i i+1 x xi+0 = dx ∆x dt dt t − t H25. Để dẫn ra dạng sai phân của = i i −1 dx dx xi−0 ∆x 5.2.1.2. Sai phân của đạo hàm cấp 2 txx Đạo hàm cấp 2 có thể nhận đ−ợc từ khai triển Taylor: dt 1 d2t 1 d3t t = t + ∆x + ∆2 + ∆3 + 0(∆4 x) i+1 i dx 2! dx 2 x 3! dx3 x + 2 3 dt 1 d t 2 1 d t 3 4 ti−1 = ti − ∆x + ∆ x − ∆ x + 0(∆ x) dx 2! dx2 3! dx3 2 d t 2 4 ti +1 + ti−1 = 2ti + ∆ x + 0(∆ x), dx2 4 Bỏ qua vô cùng bé bậc 4 là 0 (∆x ) , ta có: d 2t 1 = (t − 2t + t ) dx2 i ∆2 x i −1 i i+1 ⎛ dt dt ⎞ d 2 t ⎜ − ⎟ / ∆x = Biểu thức này cũng bằng ⎜ xi+0 xi−0 ⎟ 2 i ⎝ dx dx ⎠ dx T−ơng tự nh− trên ta có các biểu thức sai phân của đạo hàm t theo y và z. Ví dụ: Với mọi nút (ije) ∈ V của phần tử trong V, ph−ơng trình vi phân DN (hay ph−ơng trình cân bằng nhiệt) có dạng: ∂t ∂ 2t ∂ 2t ∂2t = a( 2 + 2 + 2 ) đ−ợc thay bằng: ∂τ ∂x ∂y ∂z 49
dt a a ije = a 2 (ti-1 - 2ti + ti+1)jι 2 (tj-1 - 2tj + tj+1) iι 2 (te-1 - 2te + te+1) ij dτ ∆ x ∆ y ∆ z 5.2.2. Phép sai phân vật lý: Phép sai phân vật lý là cách thay z ph−ơng trình vi phân DN trong mỗi ij-1 l ∆ ijl+1 )jl y (i + 1 phần tử nút bằng ph−ơng trình cân a λ ijl ∆z )jl i - 1 ( ij+1 l bằng nhiệt: [độ tăng nội năng của ∆x ∆ y x z ∆ nút]=[tổng các dòng nhiệt vào nút]- x ∆ ijl [tổng các dòng nhiệt ra khỏi nút] hay -1 V O dạng cụ thể là: dtije λ ρc∆x∆y∆z = [ (t - t ) ι∆x∆z dτ ∆x i-1 i j y λ ∂t + (tj-1 - tj)ie∆x∆z + H26. Để sai phân ft = a∇2t ∆y ∂τ λ λ λ λ (t - t ) ij∆x∆y] - [ (t - t ) ∆y∆z + (t - t ) ie∆x∆z + (t - ∆z e-1 e ∆x i i+1 je ∆y j j+1 ∆z e te+1) ij∆x∆y] Phép sai phân này th−ờng dùng cho các nút biên - Nếu nút (ijl) ∈W1 thì trị số tijl đã cho tr−ớc x - Nếu (ijl) ∈ W20 (cách nhiệt) thì cho qi,i+1 (hoặc qi-1,i) bằng 0, hay coi ti-1= ti+1 tức là biên đối xứng (ĐX). y - Nếu (ijl) ∈ W2 có q ≠ 0 thì thay qj,,j+1 (hoặc qj-1,,j) bằng q hoặc q(x,y,z,τ) z - Nếu (ijl) ∈ W3 với α, tf đã cho thì thay qι,ι+1 (hay qι-1, ι) bằng α (tijι-tf) - Nếu (ijl) nằm trên cạnh hay góc, nơi có 2 hay 3 ĐKB phân biệt, thì thay đồng thời các dòng nhiệt t−ơng ứng. Chú ý: Các phần tử trên biên th−ờng chọn kích th−ớc theo ph−ơng vuông góc biên bằng ∆/2, do đó, thể tích phần tử trong ph−ơng trình 50
1 1 1 3 5 7 cân bằng nhiệt có thể bằng , , hoặc , , thể tích phân tử trong 2 4 8 4 8 8 (∆x∆y∆z) 5.3. Các ph−ơng pháp xấp xỉ đạo hàm theo thời gian dt i Để xấp xỉ ta chia thời gian t i dτ t i(τ ) ra các khoảng ∆τ bằng nhau và ký tik tik+1 I I hiệu t , là nhiệt độ tại nút i lúc τ = tik+1 C i k k C k∆τ. Phép xấp xỉ này cho phép tik+1 tính nhiệt độ tik+1 của nút i lúc tik+1 E E dt i τk+1= (k+1)∆τ theo tik, ∆τ và dτ τ bằng 1 trong 4 cách sau đây: ∆τ τk=k τk+1 5.3.1. Ph−ơng pháp Euler: O ∆τ dt i dt i Tính tik+1 theo giá trị đạo hàm dτ H27.Để xấp xỉ dτ tại thời điểm τk: dt t ,k+1= t + i ∆τ (E) i ik dτ k Ph−ơng pháp này sẽ cho một hệ ph−ơng trình đại số dạng t−ờng minh. Ví dụ: Ph−ơng trình cân bằng nhiệt cho nút trong 1 chiều, sau khi a∆τ đặt p = 2 , sẽ có dạng: ∆ x a tik+1 = tik + 2 (ti-1 - 2ti + ti+1)k∆τ = [pti-1 + (1-2p)ti + pti+1]k ∆ x dt 53.2. Ph−ơng pháp ẩn (Implicit): Tính t theo i và ∆τ bằng ik+1 dτ k+1 công thức: dt t = t , k + i ∆τ (I) ik+1 i dτ k+1 Ph−ơng pháp này sẽ cho một hệ ph−ơng trình đại số thuần ẩn - Ví dụ: Ph−ơng trình cân bằng nhiệt của nút trong 1 chiều sẽ có 51
dạng: [-pti-1+(1-2p)ti-pti+1]k+1 = tik 5.3.3. Ph−ơng pháp Crank-Nicolson: dt Tính t theo trị trung bình của i tại lúc τ và τ bởi công thức ik+1 dτ k k+1 1 ⎡dt dt ⎤ i + i ∆τ tik+1 = tik + ⎢ k k+1 ⎥ 2 ⎣ dτ dτ ⎦ Ph−ơng trình này sẽ cho một hệ ph−ơng trình đại số dạng thuần ẩn. Ví dụ ph−ơng trình cho các nút trong có dạng: −1 1 1 ⎡dt i dt i ⎤ [ pti−1 + (1+ p)ti − pti+1]k+1 = ⎢ k + k+1 ⎥∆τ (C) 2 2 2 ⎣ dτ dτ ⎦ dt i 5.3.4. Ph−ơng pháp tổng quát: Tính t ,k+1 theo k và i dτ dt i với thông số tự chọn ϕ, 0 ≤ ϕ ≤ 1, theo công thức dτ k+1 dt dt t ,k+1 = t ,k + [(1- ϕ) i + ϕ i ]∆τ (T) i i dτ k dτ k+1 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: ∆2 ff Euler: ∆τ ≤ 2a(1+ B) α∆ ∆2 Với B = ff Cr-Nic:∆τ ≤ λ a(1+ B) Ph−ơng pháp ẩn nghiệm số luôn ổn định, không hạn chế ∆τ 52
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, n a x ∀i =1ữn là ∑ ij j =bi, ∀i =1ữn hay : j=1 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 ⎡a11 a12 a1n ⎤ ⎡b1 ⎤ ⎡x1 ⎤ ⎢a a a ⎥ ⎢ ⎥ ⎢ ⎥ 21 22 2n x b 2 ⎢ ⎥ ⎢ 2 ⎥ ⎢ ⎥ ⎢ ⎥ = . . . ⎢ . ⎥ ⎢ . ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ x b ⎣a n1 a 22 a nn ⎦ ⎣ n ⎦ ⎣ n ⎦ Nếu định thức det (A) ≠ 0 hệ Ax = b có nghiệm: det(D ) x = i ∀i =1ữn, trong đó D là ma trận vuông nhận đ−ợc i det(A) i 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 53
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 t − q t(x,0) = f(x) t (0,τ) = W2 W3 x λ (τ) q α −α 1 2 3 4 tf t (L, τ) = [t(L,τ)-ι ] ∆ x λ f ∆ 2 ∆ ∆ 2 x t(x,0) = f(x) ii+1 O L ta sẽ giải hệ (t) bằng FDM theo H28. FDM cho bài toán một các b−ớc đã nêu. chiều tổng quát L 1) Chia vách ra các phần ∆ = (bằng nhau), đặc tr−ng cách bởi 4 3 phần tử nh− trên. dt 2) Chuyển hệ (t) về hệ ph−ơng trình vi phân của i bằng cách dτ viết ph−ơng trình cân bằng nhiệt cho các nút (SP vật lý): ∆dt1 λ - Với nút 1 ∈ W : ρC = q − (t1 − t2), → 2 v 2dτ ∆ dt 2a 2a gọi a = λ/ρC → 1 = q − (t + t ) v dτ λ∆ ∆2 1 2 λ λ dt i - Với các nút i ∈ V: (t i−1 − t i ) − (t i − t i+1 ) = ρC v ∆ → ∆ ∆ dτ dt a i = ( t - 2t + t ), ∀i= 2,3 dτ ∆2 i-1 i i+1 ∆ dt ∆ -Với nút 4 ∈ W : ρC 4 = (t − t ) −α (t - t ) → 3 v 2 dτ 2 3 4 4 f dt a α∆ 2a α∆ 4 = [2t -2(1+ )t ]+ t dτ ∆2 3 λ 4 ∆2 λ f 54
α∆ Đặt = B ta có hệ ph−ơng trình vi phân cấp 1 là: λ dt a a q∆ i = (-2t + 2t ) + . i = 1 dτ ∆2 i i+1 ∆2 λ dt a (t ) i = (t - 2t + t ) i = 2,3 i dτ ∆2 i-1 i i+1 dt a 2a i = [t - 2(1 + B)t ] + Bt i = 4 dτ ∆2 i-1 i ∆2 f 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: dt t , = t ,k + i ]∆τ i k+1 i dτ k a a∆τ - ∀i ∈ V có: t , = t ,k + (t - 2t + t ) ∆τ, đặt = F → i k+1 i ∆2 i-1 i i+1 k ∆2 ti,k+1= [Fti-1+(1-2F)ti+ Fti+1]k ,(i = 2,3) a a q∆ - ∀i ∈ W : t , = t ,k + [ (-2t + 2t ) + . ]∆τ → 2 i k+1 i ∆2 i i+1 ∆2 λ 1 ∆ t , = F[( − 2)t + 2t ] + 2Fq ,(i = 1) i k+1 F i i+1 k λ ⎧ a 2a ⎫ - ∀i ∈ W3: ti, k+1 = tik + ⎨ []2t i−1 − 2(1+ B)t i + Bt f ⎬ ∆τ → ⎩∆2 ∆2 ⎭ 1 t,k+1 = F[2t +( -2-2B)t ] + 2FBt , (i=4) i i-1 F i k f Vậy hệ (ti) chuyển thành hệ ph−ơng trình đại số sau: 1 ∆ t ,k+1 = F[( -2)t +2t ] +2Fq 1 F 1 2 k λ 1 t ,k+1 = F[(t + ( -2)t +t ] 2 1 F 2 3 k 1 t ,k+1 = F[(t + ( -2)t +t ] 3 2 F 3 4 k 1 t ,k+1 = F[(2t + ( -2-2B)t ] +2FBt 4 3 F 4 k f 55
4) Lập dạng ma trận cho hệ ph−ơng trình đại số trên: 1 ∆ t1 ( − 2) 2 t1 q F λ 1 t 1 ( − 2) 1 t +2F 0 2 = F 2 1 t 1 ( − 2) 1 t 0 3 F 3 1 t 2 ( − 2 − 2B) t k Bt 4 F 4 f k+1 Các chú ý và nhận xét: ∆2 L2 Với ph−ơng pháp Euler, phải chọn ∆τ ≤ = , 2a(1+ B) 18a(1+ B) do ∆ = L 3 - Nếu xấp xỉ theo ph−ơng pháp Crank-Nicolson, hệ ph−ơng trình là: 1 ( + 1) -1 t F 1 1 1 1 - (+ 1) - t 2 F 2 2 1 1 1 = - (+ 1) - t 2 F 2 3 1 -1 ( +1+ B) t F 4 k+1 1 ∆ ( - 1) 1 t q F 1 λ 1 1 1 (- 1) t2 0 2 F 2 1 1 1 (- 1) t +2. 0 2 F 2 3 1 1 ( -1- B) t Bt F 4 f k 56
- Nếu xấp xỉ theo ph−ơng pháp thuần ẩn, ph−ơng trình là: 1 ∆ ( + 2) -2 t t q F 1 1 λ 1 -1 ( + 2) -1 F t2 = t2 +2F 0 F 1 -1 ( + 2) -1 t t 0 F 3 3 1 1 ( +2+2 B) t t Bt F 4 k+1 4 k f 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 [t]0 = [f(0), f(∆), f(2∆), f(3∆)] 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 = ,t q đều nằm trong ma trận hệ số tự do,nên các đại λ f 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à: ∆2 L2 ∆τ ≤ = tức có thể chọn gấp 2 b−ớc ∆τ của CN a(1+ B) 9a(1+ B) E 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: 57
t q=0 J W20 369 12 λ = λ (x, y, τ ) α ( y ,τ) hoặc λ = λ (t) q(y,τ ) 2 5118 q = q (x, y, τ ) fj v v t f (τ ) 1 4 7 10 W2 W3 x Oo W1 i I t(x,0, τ ) = t (x, τ ) o t(x,y,0) = f(x,y) 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ệ: ∂t ∂ 2t ∂ 2t q = a( + ) + v ∂τ ∂x2 ∂y 2 ρc t(x,0,τ) = t1(x,τ) (W1) ty(x,j,τ) = 0 (W20) (t q(y,τ) ) tx(0,y,τ) = (W2) − λ α(y,τ ) t (I,y,τ) = [ t(I, y,τ ) - t (y,τ)] (W ) x − λ f 3 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ạ độ: dt a a - ∀(ij) ∈ V, q = 0 có : ij = (t -2t +t ) + (t -2t +t ) v dτ 2 i-1 i i+1 j 2 j-1 j j+1 i ∆ x ∆ y dt (t − 2t + t )j (t j−1 − 2t j + t j+1 )i q Nếu q ≠ 0: ij = a[ i−1 i i+1 + ]+ v v dτ 2 2 ρc ∆ x ∆ y dt a a - ∀(ij) ∈W : ij = (t -2t +t ) + (2t -2t ) 20 dτ 2 i-1 i i+1 j 2 j-1 j i ∆ x ∆ y 58
- ∀(ij) ∈W2: ∆x dtij tij − ti t − t t − t ∆x ρc ∆y = (q-λ +1 j )∆y + (λ ij−1 ij − λ ij ij+1 ) 2 dτ ∆x ∆y ∆y 2 - ∀(ij) ∈W3: ∆x dtij t i − t ij t − t t ij − t ij ∆x ρc ∆y =[λ −ij ) − α(t − t )]∆y +[λ ij−1 ij − λ +1 )] 2 dτ ∆x ij f ∆y ∆y 2 ∀(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ó : ∆x∆y dtij t i−ij − t ij ∆y (t − t ) ∆x ρc =[λ -α(t -t )] + λ ij−1 ij 4 dτ ∆x ij f 2 ∆y 2 ∆x∆y dtij t ij − t i+ij ∆y (t − t ) ∆x Với (0,J) tại nút 3 có: ρc =(q − λ ) + λ ij−1 ij 4 dτ ∆x 2 ∆y 2 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 = a∆τ , B = α(y,τ )∆ = B(y,τ), ∆2 λ(y,τ ) nếu dùng xấp xỉ Euler, hệ ph−ơng trình đại số có dạng: (V): t = F[t +t +( 1 -4)t +t +t ] ijk+1 i-ij ij-1 F ij ij+1 i+ij k (W ): t = F[t +t +( 1 -4)t +t ] 20 ijk+1 i-ij ij-1 F ij i+ij k (t 1 2∆F (W ): t = F[t +( -4)t +2t +t ] + q ) 2 ijk+1 ij-1 F ij i+ij ij+1 k λ k (W ): t = F[2t +t +( 1 -4-2B)t +t ] +2FBt 3 ijk+1 i-ij ij-1 F ij ij+1 k f Nếu vật có nguồn nhiệt qv = qv x,y,τ) thì vế phải ph−ơng trình của q (i∆, j∆,k∆τ) mọi nút cần cộng thêm với ∆2.F v λ(i∆, j∆,k∆τ) 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∆τ : ⎡2F∆ ⎤ ⎡∆2 F ⎤ q q [t]k+1 = = Ak[Ft]k + [Ft1]k + ⎢ ⎥ + [2FtfB]k+ ⎢ v ⎥ ⎣ λ ⎦ k ⎣ λ ⎦ k Nếu ký hiệu p = 1 , F = a∆τ , B = α∆ ,C = p - 4 - 2B(∆,τ), F ∆2 λ 1 59
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 F t7 t8 1 1 p 1 1 t8 t9 1 2 p 1 t9 t10 2 C2 1 t10 t11 2 1 C2 1 t11 t12 t12 k+1 2 2 C3 k 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 2F∆ 0 0 2 F qV(∆,3∆,τ) +F + +2Ft + ∆ λ f λ t1(2∆,τ) 0 0 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 60
- 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 F∆ F∆2 t = A Ft + [Ft + 2 q + 2Ft B + q ] k+1 k k 1 λ f λ v 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 1 a∆τ 1 - 4 - 2B > 0 tức F = < F ∆2 2(2 + B) ∆2 Vì B = B (x,y,τ) nên cần lấy B = maxB và có ∆τ < ijk 2a(2 + max B) 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 α∆ ∆2 ∆x = ∆y = 10-2m thì B = = 0,125, ∆τ < = 23,5 s λ 2a(2 + B) a∆τ 1 → chọn ∆τ = 20s thì F = = 0,2 và p = = 5 ∆2 F 61
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 dt cho mọi nút) để nhận đ−ợc hệ n ph−ơng trình vi phân của ije : dτ α(x,y, τ ) ,τ ) q= , z α(x,y, τ ) 0 α( x Z α(x,y, τ ) τ ) ,y, f(x t α(x,y, τ ) z ,τ ) ( x , L α α y ( y , l+1 z ,τ ) qv (x,y,z, τ ) l +1 to (x,y,z) j λ j ( x , y , -1 z , τ l j-1 ) o q(x,y, τ ) i-1 q t i w(x , z ,τ ) i+1 I x H.30 Sai phân bài toán t(x,y,z,τ) tổng quát dtije a a a ∀ij ∈V: = 2 (ti-1-2ti+ti+1)je+ 2 (tj-1-2tj+tj+1)ie+ 2 (te-1-2te+te+1)ij dτ ∆ x ∆ y ∆ z ∀ij ∈Wn: ví dụ biên là mặt có W3 thì: λ ∆z dtije λ ∆z 2 ∆z ρc∆x∆y = ∆ (ti-1-2ti+ti+1)je.∆y + ∆2 (tj-1-2tj+tj+1)ie∆x 2 + 2 dτ x 2 y λ [ (te-1-te)ij-α(te-tf)ij]∆x∆y ,.v.v ∆ z ∀ij ∈ góc giao 3 mặt, ví dụ W20(x) x W2(y) x W3(z): 62
∆x∆y∆z dtije λ ∆y∆z ⎡ λ ⎤ ∆x∆z ρc = (2ti-1-2ti)je + ⎢ (t j−1 − t j )ie − q⎥ 8 dτ 2 4 ∆y 4 ∆ x ⎣ ⎦ je ⎡ λ + (t -t )-α(t -t ) ⎤ ∆x∆y ⎢∆z e-1 e e f ⎥ 4 ⎣ ⎦ij - áp dụng phép xấp xỉ thời gian, ví dụ ph−ơng pháp Euler, có hệ 1 t = F[t + t + t + ( − 6)t ijek+1 i-ije ij-ie ije-1 F ije + tije+1+ tij+ie+ti+ijl]k (∀ijl∈V) 1 t = F[t + t +2t + ( − 6-2B)t ijek+1 i-ije ij-ie ije-1 F ije + 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. ∆2 α∆ Để nghiệm ổn định, cần chọn ∆τ < 2a(3 + ) λ - Biểu diễn ma trận, hệ ph−ơng trình đại số có dạng: 2 [t] = {A[t]F+F[t ]+[2F ∆ q]+[2Ft B]+[ ∆ F q ]} k+1 w1 λ f λ v 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: dtijl λ ∆r λ ∆r ρc(r ∆ϕ.∆r∆z) =[ (t -t )je(r - )∆ϕ∆z- (t -t )je(r + )∆ϕ∆z i dτ ∆r i-1 i i 2 ∆r i i+1 i 2 63
λ λ λ λ +[ (tj-1-tj)ie- (tj-tj+1)ie∆r∆z]+[ (te-1-te)if - ( te- te-1)]ri∆ϕ∆z ri ∆ϕ ri ∆ϕ ∆z ∆r Khi chọn ∆r = ∆z = ∆, ph−ơng trình trên có dạng: a ∆ ∆ 2 t’ije= 2 [(1- )ti-ije+( ) tij-1e+tije-1- ∆ 2ri ri ∆ϕ 2 ∆ ∆ 2 ∆ 2(2+ 2 2 )tije+ tije+1+ ( ) tij+1e+ (1+ )ti+ije ri ∆ ϕ ri ∆ϕ 2ri * 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à : 2 a ⎡ ∆ ∆ 2 ∆ t’iJe = 2 ⎢ (1- ti-iJe+2( ) tiJ-1e- 2(2+ 2 2 )tiJe ∆ ⎣ 2ri ri ∆ϕ ri ∆ ϕ ∆ ⎤ tije+1+(1+ )t1+1Je ⎥ 2ri ⎦ * ∀(Rje) trên biên r = R loại 3, ph−ơng trình cân bằng nhiệt là: ρc (R- ∆r )∆ϕ ∆r ∆z dtRje = 4 2 dτ [ λ (t -T .(R ∆r )∆ϕ∆z-α(t -t )R∆ϕ∆z] + [ λ (t -t ) ∆r R-1 R)je 2 Rje f R∆ϕ j-1 j Re λ ∆r λ λ ∆r ∆r - R∆ϕ (tj-tj+1)Re] )∆z + [ (te-1-te)Rj- (te-te+1)Rj](R- )∆ϕ . 2 ∆z ∆z 4 2 α∆ Khi chọn ∆r = ∆z = ∆ và đặt B = λ ph−ơng trình có dạng: ⎧ ∆ a R − 2 2 t' = ⎪ 2 ∆ ⎡ ∆ Rje 2 ⎨2 t R -ije + t Rj-1e + t Rje-1 − 2 ⎢ 2 ∆ ∆ ∆ 2 ∆ 2 ⎪ R − R(R − )∆ ϕ ⎣ R(R − )∆ ϕ ⎩ 4 4 4 ∆ RB ∆2 2aRB - ⎤ + + t R je+1 + t Rj+1e ⎫ + t j . ∆ + ∆ ⎥ tRje ∆ 2 ⎬ 2 ∆ R(R − ) R(R − ) ⎦ R(R − )∆ ϕ ∆ (R − ) 4 4 4 ⎭ 4 * 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 64
nhiệt sẽ kết hợp cả 3 điều kiện biên dtije 3. áp dụng các ph−ơng pháp xấp xỉ , sẽ đ−ợc hệ ph−ơng trình dτ đại số của tije, k+1, với n ẩn số. Z α α α Zl+1 ∆z Zl ∆ z ϕ j+1 ∆ϕ Zl-1 ϕ j+1 ϕj q ∆r r o ri ri ri+1 R H31. Sai phân bài toán trụ t(π,ϕ,z,τ) Dùng xấp xỉ Euler: t = t + dtije ⎢ ∆τ, ijek+1 ijek dτ k đặt ∆ = δi, ∆ = δ, ∆ϕ = γ và F = a∆τ , ta có hệ ph−ơng trình đại số: ri R ∆2 2 ⎧ δi δi 2 1 δ tijek+1 = F ⎨(1− ) ti-1je+ ( ) tij-1e+tije-1+( -4-2 )tije ⎩ 2 γ F γ 2 + t +( δi ) t + (1+ δi )t ⎫ ∀(ije)∈ V ije+1 γ ij+1e 2 i+1je ⎬ ⎭k ⎧ δ2 1 2 ⎪ 2 − δ − 4 − (t) tijek+1 = F ⎨ tR-ije+ 2 δ tRj-1e+tRje-1+[ F δ . δi γ (1− ) (1− ) ⎪1− 4 ⎩⎪ 4 4 δ 2 δ δ2 ⎫ 2BF ( 2 - -B)]tRje+tRje+1+ tRj+1e ⎬ + tf, ∀(ije)∈ W3 γ 4 2 δ δ γ (1− ) ⎭k (1− ) 4 4 và các ph−ơng trình có nút trên biên, cạnh, góc, khác. 2 Dạng ma trận:[t] = F{A[t]+[t ]+[ 2∆ q]+[2Bt ]+[ ∆ q ]} nh− trên k+1 w1 a f λ v k với A là ma trận hệ (nxn), các ma trận khác là ma trận cột (nx1). 65
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: n m 1. Ph−ơng trình F(T , Tx ) = 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: 4 4 -λTx (W,τ) = εwσ0[T (W,τ)- Tf ] Đó 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à: 66
Tτ = aTxx t T(x,0) = Ti δT = 0 (t) Tf δx − εσ s a 0 4 4 λ Tx(0,τ)= [ Tf -T (0,τ)] q = ε∂ (T 4 - T 4 ) 0 1 2 3 4 λ f o S T (0,τ) = 0 x x Nhận xét: bài toán này có cùng O L một mô hình TH với bài toán biên H32. Bài toán biên phi tuyến thuần nhất, khi vách dày 2δ và TĐ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 ∆ = L , ví dụ ∆ = L x n x 4 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: ∆x dT λ ρcS o = εσ ( T 4 - T 4 )S - (T -T )S → 2 2dτ 0 f 0 ∆x 0 1 dTo 2εσ0 4 4 2λ = ( Tf - T0 ) - 2 (T0-T1) 2dτ ρc∆x ρc∆ x T Nếu chuẩn hoá bằng cách đặt θ = Ti x ∆x X = → ∆X = L L λτ aτ F = = cρL2 L2 ∂T ∂T ∂θ ∂F T a ∂θ Thì do = . . = i . nên: ∂τ ∂θ ∂F ∂τ L2 ∂F dθ L2 dT 2εσ L2 2λL2 (T − T ) 0 0 0 4 4 0 4 hay = . = (Tf − T0 ) − dF aT dτ λ λ 2 i Ta(L∆x) (L∆x) aT a i a i 67
dθ 3 3 0 1 ⎡ εσ0 .LTi 3 ⎤ εσ0 .LTi 4 = 2 ⎢ -2(1+ ∆Xθ0 )θ0+2θ1 ⎥ +2 .θf dF (∆X) ⎣ λ ⎦ λ∆X εσ .LT3 Nếu đặt R= 0 i là đại l−ợng không thứ nguyên, ta có: λ 1 2R 4 θ = [-2(1+R∆X.θ3 )θ +2θ ]+ θ , ∈ W 0F (∆X) 2 0 0 1 ∆X f 6 1 θ = [θ - 2θ +θ ] iF (∆X) 2 0 1 2 1 (θF) θ = [θ - 2θ +θ ] ∈V 2F (∆X) 2 1 2 3 1 θ = [θ - 2θ +θ ] 3F (∆X) 2 2 3 4 1 θ = [2θ - 2θ ] ∈W 4F (∆X) 2 3 4 20 ∆F 3. Xấp xỉ Euler θ = θ +θ ∆F, đặt = p ta có hệ ph−ơng k+1 k Fk (∆X) 2 trình đại số, viết ở dạng ma trận nh− sau: 3 4 θ0,k+1= [1-2p(1+R∆X θ0 ]θ0,k+2pθ1,k+2pR∆X θf (θ ) i θ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ó: 3 2p 4 θ0 [1-2p(1+R∆X θ0 ] θ0 2pR∆X θf p (1-2p) p 0 θ1 θ1 p (1-2p) p 0 θ2 = θ2 + p (1-2p) p 0 θ3 θ3 2p (1-2p) 0 θ4 k+1 k θ4 k Thay điều kiện đầu lúc τ = 0 (tức F = 0) theo θ(X,0) = 1 vào 68
3 [θ]k=0 tính phần tử bức xạ của ma trận [1-2p(1+R∆Xθ0k ] 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 T(x,0) = Ti (t) Tf1 − εσ s a 0 4 4 λ ε∂ 4 4 Tx(0,τ) = [Tf1 -T (0,τ)] q = (T f1 - T o ) 0 1 2 3 4 α λ Tf2 α T (L,τ) = [T(L,τ)-T ] x − λ f2 x O δ Bài toán này đ−ợc giải t−ơng tự nh− trên, chỉ khác ph−ơng H33. Bài toán truyền nhiệt KOD trình cân bằng nhiệt cho nút phi tuyến tính biên W2 là: ∆x λ α∆x αL pc T = (T -T )-α[T -T ] hay với B = = ∆X → 2 4τ ∆x 3 4 4 f2 λ λ 1 2B θ = [2T -2(1-B)T ]+ T F (∆X) 2 3 4 (∆X) 2 f2 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à: 3 4 θ0 [1-2p(1+R∆X θ0 ] 2p θ0 2pR∆X θf θ1 p (1-2p) p θ1 0 θ2 = p (1-2p) p θ2 + 0 θ3 p (1-2p) p θ3 0 2p [1-2p(1-B) θ4 k+1 k θ4 k 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 69
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à : 3 min{[1-2p(1-B)],[1-2p(1+ R∆ θ0 ]} > 0 3 tức là: 1-2p(1+ R∆Xθ0 > 0 → ∆F 1-2 (1+R∆X θ3 ) > 0 tức phải chọn; (∆X) 2 0 2 (∆X) Tf1 ∆F < 3 , với max θ0 = 2[1+ R∆X(max θ0 )] Ti ∆ 2 x ,[s] Nghĩa là cần chọn ∆τ < L2 ∆xεσ 2a(1+ 0 T 3 ) λ f1 Khi giải hệ ph−ơng trình trên, phải th−ờng xuyên tính lại số hạng 3 đầu [1-2p(1+R∆ θ 0 )]k theo θ0, k sau mỗi b−ớc. 70
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: Cd [t] = -(K+H)[t] + (h+q) dτ 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 1 ⎤ ⎢ ⎥ [t] = ⎢ ⎥ là ma trận cột (Mx1) của các nhiệt độ nút. ⎣⎢t M ⎦⎥ 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: 71
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: t 2 I[t(x,y,z)] =∫ f(x,y,z,t,tx,ty,tz) dV + ∫ (qt + α )dw v w 2 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. t t(x,y,τ =const) tζ 1 2 k fi f(e) y yk yj S yi 3 i e k O x i x j W 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ó (e) dạng mặt phẳng qua 3 điểm ti, tj, tk), và xác định phân bố t nh− hàm (e) tuyến tính của các biến, có dạng f(x,y,ti,tj, tk) = t (xem H.32) 72
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 E tử (e): I = ∑ I (e) e=1 Cho thấy I = I(t1, t2, , tM) là phiếm hàm của M hàm nhiệt độ nút t (τ) và điều kiện cực tiểu δI = 0 trở thành dI = 0, i d[t] t E e ⎡ 1 ⎤ dI dI với [t] = ⎢ ⎥ tức = ∑ = 0 ⎢ ⎥ d[t] e=1 d[t] ⎢ ⎥ ⎣t M ⎦ 2.4. Dùng phép biến đổi ma trận để đ−a điều kiện trên về dạng: dI = K[t] + C[&t ] + H[t] - h - q = 0 d[t] dt Đó là hệ ph−ơng trình đại số của t và &t = dτ 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 73
du d 2 u cực tiểu tại x là: = 0 và > 0 dx dx 2 - Điều kiện cần và đủ để u = u(x,y) đạt cực tiểu tại (x,y) là: ⎪⎧u x = 0, u y = 0 ⎨ 2 ⎩⎪u xx > 0, (u xx u yy − u xy ) > 0 - 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: ∂u = 0, ∀i = 1 ữ n. Sau đây sẽ dùng ký hiệu "A ∆ B", có nghĩa là " B ∂x i đ−ợc định nghĩa là A". ⎡ ⎤ ⎢x1 ⎥ ⎢ ⎥ Định nghĩa một ma trận cột [x] ∆ x 2 , điều kiện trên có thể viết ở ⎢ ⎥ ⎢ ⎥ ⎣x n ⎦ ⎡ ∂u ⎤ ⎡0 ⎤ ⎢∂x ⎥ ⎢ ⎥ ⎢ 1 ⎥ ⎢ ⎥ ⎢ ∂u ⎥ ⎢0 ⎥ ∂u ∂u ⎢ ⎥ dạng: = 0hay cụ thể , ∆ ⎢∂x ⎥ = ∂[x] ∂[x] ⎢ 2 ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ∂u ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ ⎣⎢∂xn ⎦⎥ ⎣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 74
t ∆x = xj - xi , hàm F(x,y) là mặt phẳng e F(x) trong phần tử tam giác 2 chiều ∆ , hàm 1 (Fi + fj) 2 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− (1) (2) (e) (E) x tổng các giá trị trung bình của tích Oo x1 xi xj xM E ∆x phân trên mỗi phần tử I = ∑ I (e) e=1 H.33. Để xấp xỉ tích phân 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. x E E M (e) Fi + Fj I = F(x)dx là I = ∑ I = ∑ (xj-xi) ∫x1 e=1 e=1 2 với Fi = F(xi), Fj = F(xj), (Fi+Fj)/2 là trị trung bình của F(x) trên e phần tử (e) có kích th−ớc ∆ x = xj-xi ∆ Nếu chọn xij = xj - xi = ∆, M = 4, E = 5 thì ta có: x M ∆ F(x)dx = (F + F ) e ∆ ∫x ∑ i j = (F1+2F2 + 2F3 + 2F4 + F5) 1 2 2 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. 75
Ví dụ: + Công mà hệ thực hiện trong một quá trình p = p(v) bất kỳ giữa v 2 p(v) hai trạng thái (p1, v2), (p2,v2) là phiếm hàm l = l[p(v)] = ∫v dv 1 + Nhiệt l−ợng do hệ trao đổi trong quá trình bất kỳ T = T(s) giữa s 2 T(s) hai trạng thái (T1,s1), (T2, s2) là phiếm hàm q = q[T(s)] = ∫s ds. 1 + Độ dài đ−ờng cong y(x) bất kỳ qua 2 điểm (x0,y0) (x1,y1) là x 2 1+ y'2 dx phiếm hàm l = l[y(x)] = ∫x 1 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. yt - 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: y(x) Tìm đ−ờng cong y(x) nối hai điểm O x xo S x1 (x0,y0) (x1,y1) cho tr−ớc, để khi quay z 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 H.34. Bài toán tìm mặt cực tiểu y (1) là tìm cực tiểu của phiếm x hàm:S[y(x)] = 1 1+ y'2 dx ∫x 0 v y(x) ds 2. Tìm đ−ờng đoản thời: g dx Tìm đ−ờng cong y(x) nối hai điểm x A(x0,y0), B (x1,y1) cho tr−ớc để chất O điểm lăn không ma sát trên nó từ A đến xo x1 B tốn ít thời gian nhất. H35. BT đ−ờng đoản thời Theo định luật bảo toàn năng l−ợng, tìm thấy quan hệ τ = τ[y(x)] 76
và bài toán biến phân là tìm cực tiểu phiếm hàm nh− sau: 2 2 2 m dx + dy 2 1 ⎛ dy ⎞ mgy = ( ) → dτ 1+ ⎜ ⎟ dx → 2 dτ 2gy ⎝ dx ⎠ 2 1 x1 1+ y' τ = τ[y(x)] = dx 2g ∫x0 y 3. Tìm đ−ờng trắc địa: Tìm đ−ờng cong f(x,y,z) nối hai z điểm A,B trên mặt cong ϕ(x,y,z) = 0 cho tr−ớc, có độ dài ngắn nhất. Theo ϕ(x,y,z)=0 ý 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 B x l= l[y(x),z(x)] = 2 1+ y'2 +Z'2 dx ∫x1 A y ⎧y = y(x) trong đó ⎨ là các hàm xác định O ⎩z = z(x) x dạng tham số của đ−ờng cong f(x,y,z)∈ϕ(x,y,z) = 0. H.36. BT tìm đ−ờng ⎧y = y(x) 6.3.3. Biến phân của phiếm trắc địa ⎨ z = z(x) 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 = u (x) - u(x) trong u u(x) đó u (x) thay đổi tuỳ ý trong lớp hàm δu u(x,ε) nào đó, gần đ−ờng cong u(x). u(x) T−ơng tự định nghĩa: δu(x,y) = u (x, y) - u(x,y) x - ý nghĩa hình học: Ox0 x1 u(x) thay đổi trong lớp đ−ờng cong H.37. Biến phân δu(x) qua hai điểm biên (x0,y0) (x1,y1) phụ thuộc thông số ε nhỏ nào đó: 77
u u(x) ∈ u(x;ε) = u(x) + ε[u (x) -u(x) u(x,y) δu u(x,y) = u(x) + εδu. u (x, y) thay đổi trong lớp mặt cong qua s y chu tuyến S, phụ thông số ε nhỏ nào đó: D u (x, y) ∈ u(x,y;ε) và u(x,y;ε) = u(x,y) + x ε[ u(x, y) -u(x,y)] = u(x,y) + εδu(x,y). H.38. Biến phân δu(x,y) Biến phân δu là sự sai khác nhỏ giữa hai hàm u cùng chung biên nào đó. 6.3.3.2. Biến phân của phiếm hàm: - Định nghĩa: Biến phân δI của phiếm hàm I[u(x)] là đạo hàm của phiếm hàm I[u(x) + εδu(x)] theo tham số ε, tính tại ε = 0. d δI[u(x)] = I[u(x)] + εδu(x)]⏐ε dε =0 T−ơng tự có các định nghĩa d δI[u(x,y)] = I[u(x,y)] + εδu(x,y)]⏐ε ,và dε =0 d δI[u (x), ,u (x)] = I[(u (x) + εδu ), ,(u (x) + εδu )]⏐ε 1 n dε 1 1 n n =0 - ý nghĩa: Biến phân δI là phần chính bậc 1 (tuyến tính) của số gia phiếm hàm ∆I = I[u + δu] - I[u] khi max [δu] → 0 δI mô tả sai số các tích phân (phiếm hàm) suy từ một định luật bảo toàn nào đó, so với định luật bảo toàn năng l−ợng. 6.3.4. Định lý Euler -Lagrange x 6.3.4.1. Điều kiện cực tiểu phiếm hàm I[u(x)] = 1 F(x,u,u')dx ∫x0 - Cần tìm hàm u = u(x) qua 2 điểm biên cho tr−ớc [x0, u0 = u(x0)] và [x1, u1 = u(x1)] làm cực tiểu phiếm hàm có dạng: x 1 F I[u(x)] = ∫x (x,u(x), u'(x))dx. 0 - Tập hợp các đ−ờng cong qua 2 điểm biên có thể mô tả bằng hàm 78
u u(x) u(x,ε) của x và thông số ε ∈ [0,1] bởi: δu u(x,ε) = u(x) + ε[ u (x) - u(x)] u(x) = u(x) + εδu(x) với u (x) là đ−ờng cong gần x O x0 x1 u(x) và δu(x) là biến phân của hàm H.39. Để cực tiểu u(x). I[u(x)] - Với đ−ờng cong u(x;ε) thì phiếm hàm I[u(x,ε)] là một hàm chỉ của ε, do đó điều kiện cần để cực tiểu phiếm hàm của hàm u(x), ứng với ε = 0, là: d I[u(x,ε)] = 0 tức là δI = 0 dε ε=0 d x1 → δI = 0 tức là δI = F(x,u(x;ε), u'(x;ε))dx|ε=0 dε ∫x0 x1 ⎡∂F ∂u(x,ε) ∂F ∂u'(x,ε)⎤ . + . dx = ∫x ⎢ ⎥ 0 ⎣ ∂u ∂ε ∂u' ∂ε ⎦ε=0 x 1 (Fu.δu + Fu'δu')dx = 0 = ∫x 0 Phân đoạn tích phân sau, ta có: x1 x1 d x1 ( Fu ' )δudx δI = ∫ Fu δu dx + Fu' δu | x - 0 ∫x 0 dx x0 x1 d x1 = (Fu − Fu ')δudx = 0 (do δu | x = 0 - 0 = 0) ∫x0 dx 0 d Vì δu (x) ≠ 0 ∀x ∈ (x , x ) nên Fu - Fu' = 0 0 1 dx → Suy ra định lý Euler - Lagrange: - Định lý E - L1: x1 Để cực tiểu phiếm hàm I[à(x)] = ∫x F(x, u, u')dx, hàm u(x) cần 0 79
∂F d ∂F thoả mãn điều kiện sau: − ( ) = 0 ∂u dx ∂u ' 6.3.4.2. Ph−ơng pháp biên phân: - Ph−ơng pháp biến phân là ph−ơng pháp sử dụng điều kiện cực tiểu để tìm đ−ờng cong cực trị. d - Điều kiện cực tiểu I[y(x)] theo EL là: F - F = 0 → suy ra: y dx y' Fy- Fy'x- Fy'y.y'-Fy'y'.y'' = 0 ,(E) Ph−ơng trình (E) gọi là ph−ơng trình Euler. - Ph−ơng pháp tìm đ−ờng cong cực trị là: 1) Tích phân ph−ơng trình Euler thu đ−ợc hàm y(x,c1, c2) 2) Xác định c1, c2 theo hai điều kiện biên y0 = y(x0), y1 = y(x1) Ví dụ: Khi hàm F chỉ phụ thuộc y và y': F=F(y,y') → ph−ơng trình (E) có dạng Fy- Fy'y.y'-Fy'y'.y'' = 0 → Nhân hai vế với y' có: 2 y'(Fy- Fy'y.y'-Fy'y'.y'') = Fy.y'-Fy'yy' -Fy'y'y'y'' = = d (F-F .y')= 0, (dạng đạo hàm toàn phần theo x) dx y' Tích phân đầu theo x có F- Fy'y' = C1 là một ph−ơng trình vi phân không chứa x, giải đ−ợc bằng cách tách biến hoặc đổi biến, ví dụ * Tìm mặt cực tiểu (bài toán 1) : x - c2 y = c1 ch c1 Mặt tròn xoay qua (o,y0), (x1,y1) y cực tiểu khi phiếm hàm S[y(x)] cực y tiểu. Do dS = 2πyds = 1 yo ds 2 2 dx x 2πy dx + dy nên suy ra O x1 x 1 y 1 + y'2 dx z S[y(x)] = 2π ∫x = 0 . H.40. Mặt Catenoid Với F = y 1+ y'2 = F (y,y'), ta có tích phân đầu: 80
2 2 yy' y F- y'Fy' = y 1+ y' - = = C1 1+ y'2 1+ y'2 dy C1 shtdt Đổi biến y' = sht thì y = C cht và dx = = = C dt 1 y' sht 1 → x = C1t + C2 x = C t + C ⎧ 1 2 x − C2 Khử t từ ⎨ có y = C1ch là đ−ờng dây xích, quay ⎩y = C1cht C1 quanh ox tạo ra mặt Catenoid, với C1, C2 tìm theo điều kiện biên yo = y(0) , y1 = y(x1). * Tìm đ−ờng đoản thời: (bài toán 2) Thời gian chất điểm lăn trên đ−ờng cong y(x) từ A→B là: 0 ≡A 2 x1 1 x1 1+ y' + dx τ = τ[y(x)] = ∫x=0 ds 2g y v = dy dτ g ds 1+ y'2 d Với F = F(y,y') = x y y 1 B ta có tích phân đầu: y 2 2 1+ y' y' H.41. Đ−ờng đoản thời F- y'Fy' = - y y(1+ y'2 ) ⎧x = C(t − sin t) Cycloid ⎨ 1 ⎩y = C(1− cos t) = = C y(1+ y'2 ) 2 → y(1 + y' ) = C1 là ph−ơng trình vi phân cấp 1 phi tuyến. C1 C1 Đổi biến y' = cotg t có y = 2 = 2 1+ y' 1+ cot g t C1 = C sin2t = (1-cos2t) và : 1 2 dy 2C1 sint cost.dt 2 dx = y' = cot gt = 2C1sin t dt 81
= C (1- cos2t)dt → x = C (t- 1 sin2t) + C 1 1 2 2 C1 hay x = (2t - sin2t) + C 2 2 Đ−ờng đoản thời có dạng họ đ−ờng cycloid ⎧ C1 ⎫ ⎪x = (2t − sin 2t) + C 2 ⎪ ⎪ 2 ⎪ C1 ⎨ ⎬ C qua A(0,0) → C2 = 0 đặt C = , ⎪y = 1 (1 − cos 2t) ⎪ 2 ⎩⎪ 2 ⎭⎪ 2t = ϕ ta có đ−ờng đoản thời là đ−ờng cycloid, mô tả bởi ph−ơng trình dạng tham số sau: ⎧x = C(ϕ − sin ϕ) y y ⎨ hay x = c[arcos(1- ) - sin(arcos(1- ))] ⎩y = C(1 − cos ϕ) c c với C = bán kính vòng tròn lăn để vẽ ra Cycloid, xác định theo y1 = y(x1). 6.3.4.3. Điều kiện cực tiểu phiếm hàm dạng I[u(x,y)] F I[u(x,y)] = ∫∫ (x,y,u,ux,uy)dxdy D - Cần tìm hàm u(x,y) trong miền D có trị số cho tr−ớc trên biên W ∈ D, làm cực tiểu phiếm hàm ∂u ∂u I[u(x,y)] = F(x,y,u , )dxdy ∫∫ ∂x ∂y D - ứng với các mặt thay đổi trong họ mặt cong u(x,y,ε) = u(x,y) + ε[ u (x,y) - u(x,y) ] = u(x,y) + εδu(x,y) thì phiếm hàm I[u(x,y,ε)] chỉ phụ thuộc ε, đ−ợc cực tiểu khi: d I[u(x,y,ε)] = δI = 0, tức là khi: dε ε=0 82
⎡ ∂u ∂u ∂u ⎤ Fu + Fu x + Fuy y δI = ∫∫ ⎢ x ⎥ dxdy D ⎣ ∂ε ∂ε ∂ε ⎦ε=0 F δà + F δu + F δu )dxdy = 0 U = ∫∫ u ux x uy y u D u(x,yU(x,y) δUu u(x,yU(x,y) ⎧ ∂ ∂Fu x ⎪ (Fu x δu) = δu + Fu x δu x ⎪∂x ∂x - Vì có: ⎨ ∂Fu s ⎪ ∂ y (Fu y δu) = δu + Fu y δu y y ⎩⎪∂y ∂y 0 D (Fu δu + F δu ) nên ∫∫ x x uy y dxdy = W D x ⎡ ∂ ∂ ⎤ (Fu δu) + (Fu δu) H.42. Để cực tiểu = ∫∫ ⎢ x y ⎥ dxdy- D ⎣∂x ∂y ⎦ phiếm hàm I[à(x,y)] ∂Fu x ∂Fu y ∫∫( + )δudxdy D ∂x ∂y ∂Fu ∂Fu y = - ∫∫( x + )δudxdy D ∂x ∂y ∂ ∂ [ (Fu δu) + (Fu δu)]dxdy (do theo công thức Green có: ∫∫ x y D ∂x ∂y (Fu δu − Fu δu)dx = 0 = ∫ x y vì δu|W = 0) w - Vì vậy, điều kiện cực tiểu phiếm hàm là: ∂Fu ∂Fu δI = (Fu − x − y )δudxdy ∫∫ = 0, do δu|D ≠ 0 suy ra: D ∂x ∂y ∂F ∂ ∂F ∂ ∂F → − ( ) − ( ) = 0 ∂u ∂x ∂u x ∂y ∂u y Định lý Euler-Lagrange N02: F Để cực tiểu phiếm hàm: I[u(x,y)] = ∫∫ (x,y,u,ux,uy,)dxdy, hàm D 83
∂F ∂ ∂F ∂ ∂F u(x,y) cần thoả mãn điều kiện: − ( ) − ( ) = 0 ∂u ∂x ∂u x ∂y ∂u y 6.3.4.4. Điều kiện cực tiểu phiếm hàm dạng I[u(x,y,z)]. Các bài toán truyền nhiệt 3 chiều không ổn định với biên loại 2 và 3 th−ờng t−ơng ứng với phiếm hàm dạng: t 2 F (qt + α )ds I[u(x,y,z)] = ∫∫∫ (x,y,z,t,tx,ty,tz)dV+ ∫ D S 2 - Khi cực tiểu phiếm hàm I, cần có δI = 0, tức là: ⎡∂F ∂F ∂F ∂F ⎤ δI = ⎢ δt + δtx + δt y + δtz ⎥ dV+ (qδt +αtδt)ds = 0 ∫ ∂t ∂t ∂t ∂t ∫ V ⎣⎢ x y z ⎦⎥ S ∂t ∂ - Theo δt = δ ( ) = (δt) , v.v., ta có: x ∂x ∂x δI= z ⎡∂F ∂F ∂ ∂F ∂ ∂F ∂ ⎤ δt + . (δt) + . (δt) + . (δt) d → ∫ ⎢ ⎥ ny0 V ⎣⎢ ∂t ∂t x ∂x ∂t y ∂y ∂t z ∂z ⎦⎥ S V+ ∫δt(q + αt)ds = 0 S V - Sau khi tích phân từng phần (per partes) ta y 0 có: x ∂F ∂ ∂F ∂ ∂F . (δt)dV = δt.l dS - δt .( )dV H.43. Miền tích phân V,S ∫ ∫ x ∫ V ∂t x ∂x S ∂tx V ∂x ∂t x của bài toán 3 chiều ∂F ∂ ∂F ∂F ∂ δt.l δt .( )dV ∫ . (δt)dV = ∫ y dS - ∫ V ∂t y ∂y S ∂t y V ∂y ∂t y ∂F ∂ ∂F ∂ ∂F . (δt)dV = δt.l dS - δt .( )dV ∫ ∫ z ∫ V ∂t z ∂z S ∂tz V ∂z ∂t z với lx, ly, lz là cosin chỉ ph−ơng của pháp tuyến n phía ngoài V tại M ∈ S. Do đó, điều kiện δI = 0 sẽ có dạng: 84
⎡∂F ∂ ∂F ∂ ∂F ∂ ∂F ⎤ δI = δt⎢ − ( ) − ( ) − ( )⎥dV ∫ ∂t ∂x ∂t ∂y ∂t ∂z ∂t V ⎣⎢ x y z ⎦⎥ ∂F ∂F ∂F δt l +l + l + q + αt + ∫ [ x y z ]dS = 0 D ∂t x ∂t y ∂t z Do δt|v ≠ 0 → biểu thức trong dấu móc [∼]v = 0 → suy ra định lý E-L3: Định lý Euler -Lagrange N03: Để cực tiểu phiếm hàm I[t(x,y,z)] cho theo công thức (3), hàm t(x,y,z) cần thoả mãn điều kiện: ⎧∂F ∂ ∂F ∂ ∂F ∂ ∂F ⎪ − ( ) − ( ) − ( ) = 0,∀(x, y, z) ∈V ⎪ ∂t ∂x ∂tx ∂y ∂t y ∂z ∂tz ⎨ ∂F ∂F ∂F ⎪l + l + l + q + αt = 0,∀(x, y, z) ∈ S ⎪ x y z ⎩ ∂tx ∂t y ∂tz τ 6.4. Ví dụ minh hoạ các b−ớc áp dụng FEM 6.4.1. Bài toán biên cô lập Ta sẽ dùng FEM giải bài toán DN t = to tx = 0 ∂t ∂ 2t x không ổn định một chiều, biên cô lập ρ = λ c ∂τ ∂x 2 W1+ W20 cho bởi mô hình. ⎧ ∂t ∂ 2 t t(x,0) = ti ⎪ρc = λ 2 ∂τ ∂x ⎪ O x x (t) t(0,τ) = t , t (L,τ) = 0 ⎨ 0 x H.44.Bài toán (t) ⎪t(x,0) = t ⎪ i ⎩⎪ Ta sẽ thực hiện các b−ớc đã nêu ở (6.1.2.) 6.4.2. Phát biểu biến phân: ( Variational Statement) - Xác định phiếm hàm I[t(x,τ)] t−ơng ứng bài toán (t) bằng cách so sánh ph−ơng trình vi phân dẫn nhiệt và ph−ơng trình Euler- ∂F ∂ ∂F Lagrange để tìm hàm F: − ( ) = 0 (ph−ơng trình E-L) ∂t ∂x ∂t x 85
∂t ∂ ∂t ρc − (λ ) = 0 (ph−ơng trình vi phân DN) ∂τ ∂x ∂x ∂F ∂t ∂F Phép so sánh cho thấy = ρc và = λt x ∂t ∂τ ∂t x Tích phân 2 ph−ơng trình trên theo t và tx, coi t và tx là biến độc lập, ta có: 2 ∂t ∂ 1 ∂t F = ρc dt = ρc tdt = ρc + f (t x ) ∫ ∂τ ∂τ∫ 2 ∂τ ∂t 1 2 F = λ dt x = λ t x dt x = g(t) + λt x ∫ ∂x ∫ 2 So sánh hai tích phân cho thấy hàm tích phân là: 2 1 ⎡ ∂t ∂t 2 ⎤ F= ⎢ρc + λ( ) ⎥ 2 ⎣ ∂τ ∂x ⎦ Bài toán biến phân t−ơng ứng với bài toán (t) là tìm hàm t = t(x,τ) sao cho tại thời điểm bất kỳ làm làm cực tiểu phiếm hàm L 2 1 ⎡ ∂t 2 ∂t ⎤ I = ∫ ⎢λ( ) + ρc ⎥ dx 2 x=0 ⎣ ∂x ∂τ ⎦ 6.4.3.Phát biểu phần tử hữu hạn (Finite Element Formulation) Đầu tiên chia phiếm hàm I thành tổng hai tích phân: 2 1 L ∂t 1 L ∂t λ( ) 2 dx ρc dx I = Iλ+Ic với Iλ = ∫ , Ic = ∫ 2 0 ∂x 2 0 ∂τ 6.4.3.1. Phân hoạch các phần tử hữu t i hạn t1 Tiếp theo chia miền tích phân [0,L] ra E t2 te phần tử nhờ M = E + 1 điểm nút ti tj tM Tại thời điểm τ bất kỳ, trong mỗi phần tử e, O (1) (2) (e) (E) x e x x ta coi nhiệt độ t thay đổi tuyến tính theo x từ ti xx11 x22 xii xjj xMM = t(xi) đến tj = t(xj), tuy nhiên ti và tj ch−a biết. H45. FE formulation Toàn bộ t− t−ởng của FEM là tìm các nhiệt độ (t1, t2, ,ti,tj, ,tM) 86
sao cho phiếm hàm I nhỏ nhất. Phiếm hàm I thực chất là sai số của phép xấp xỉ hàm t(x) bởi tập giá trị (t1, t2, ,ti,tj, ,tM) tại các phần tử hữu hạn, so với hàm t(x) chính xác ch−a biết thỏa mãn ∂t ∂ 2 t ph−ơng trình cân bằng nhiệt ρc − λ = 0 . ∂τ ∂x 2 Tại thời điểm bất kỳ, phiếm hàm I là hàm của M nhiệt độ nút I = I(t1, t2, , tM) nên điều kiện cần để cực tiểu I là sự triệt tiêu của các đạo hàm I lần l−ợt theo mỗi nhiệt độ nút ti. ∂I ∂I ∂I c = λ + = 0 , ∀i = 1ữ M ∂t i ∂ti ∂t i ⎡t1 ⎤ ⎢t ⎥ ⎢ 2 ⎥ Định nghĩa ma trận [t] ⎢ ⎥ (1), thì điều kiện cực tiểu I là: ⎢ ⎥ ⎣t M ⎦ ⎡ ∂I ⎤ ⎢∂t ⎥ ⎢ 1 ⎥ ⎢ ∂I ⎥ ⎢ ⎥ dI dI λ dIc dI ∂t = + = 0, trong đó định nghĩa ⎢ 2 ⎥ d[t] d[t] d[t] d[t] ⎢. . . ⎥ ⎢ ⎥ ⎢ ∂I ⎥ ⎢ ⎥ ⎣∂t M ⎦ Biểu diễn I nh− tổng các Ie tại mỗi phần tử e; E E E e e e I = ∑ I = ∑ I λ + ∑ I c với e=1 e=1 e=1 2 x x e 1 j ∂t 2 e 1 j ∂t Iλ = λ( ) dx và Ic = ρc dx 2 ∫xi ∂x 2 ∫xi ∂τ Chú ý: Các chỉ số e hiểu là "của phần tử e", không phải số mũ. e ⎡t i ⎤ Định nghĩa ma trận (2x1) nhiệt độ nút của phân tử e: [t] ⎢ ⎥ (2) ⎣t j ⎦ 87
⎡0 0⎤1 ⎢ ⎥ ⎢. . ⎥. ⎢1 0 ⎥i e ⎢ ⎥ và ma trận (Mx2) định vị phần tử e là: D ⎢0 1⎥j (3) ta có: ⎢. . ⎥ ⎢ ⎥. ⎢0 0⎥ ⎣ ⎦M ⎡0 0⎤ ⎡ 0 ⎤ ⎢ ⎥ e ⎢ e ⎥ . . ⎡∂Iλ ⎤ E e ⎢∂Iλ ⎥ ⎢ ⎥ E e E ⎢ ⎥ dI λ dIλ ⎢ ∂t ⎥ E ⎢1 0 ⎥ ∂t e dIλ ⎢ i ⎥ . ⎢ i ⎥ D = ∑ = ∑⎢ e ⎥ = ∑ ⎢ ⎥ e = ∑ e d[t] d[t] ∂Iλ 0 1 ⎢ ⎥ d[t] e=1 e=1 ⎢ ⎥ e=1 ⎢ ⎥ ∂Iλ e=1 ⎢ ∂t j ⎥ ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ . . ⎢ ∂t j ⎥ ⎣⎢ 0 ⎦⎥ ⎢ ⎥ ⎣ ⎦ ⎣⎢0 0⎦⎥ E dI e T−ơng tự có: dIc = De c d[t] ∑ d[t]e e=1 Nh− vậy, để cực tiểu I, phải xác định [t] sao cho: E e E e dI dI e dI c λ D e λ D = ∑ e + ∑ e = 0 d[t] e=1 d[t] e=1 d[t] 6.4.3.2. Mô tẩ hàm te: Giả thiết tại thời điểm bất kỳ, nhiệt độ trong phần tử e thay đổi nh− một hàm tuyến tính của ti, tj và x (tức là đoạn thẳng qua xiti, xjtj), e e e e ⎡ϕ1 ⎤ T e T cho bởi: t = ϕ1 + ϕ2 x= [1x] ⎢ ⎥ p [ϕ] , ở đây gọi ρ [1 x] (4), ⎣ϕ2 ⎦ e e các hệ số ϕ1 ,ϕ2 xác định theo: e e e t i = ϕ1 + ϕ2 x i ⎫ ⎪ e ⎡1 x i ⎤⎡ϕ1 ⎤ e e e e ⎬ hay dạng ma trận là: [t] = ⎢ ⎥ P [ϕ] t = ϕ + ϕ x 1 x ⎢ϕ ⎥ j 1 2 j ⎭⎪ ⎣ j ⎦⎣ 2 ⎦ e ⎡1 x i ⎤ Với P ⎢ ⎥ (5). Từ đó suy ra: ⎣1 x j ⎦ x x 1 ⎡ j i ⎤ [ϕ]e = Pe-1[t]e = Re[t]e với Re Pe-1 = ⎢ ⎥ (6) là ma trận x − x j i ⎣⎢−1 1⎦⎥ nghịch đảo của Pe E e dI λ e dIλ 6.4.3.3. Tính = ∑ D e : d[t] e1= d[t] 88
e e 2 x j x e 1 ∂t 2 λ j ⎡ ∂ T e e ⎤ Ta có: Iλ = λ( ) dx = ( p R [t] ) dx ∫x ∫x ⎢ ⎥ 2 i ∂x 2 i ⎣∂x ⎦ e x 2 λ j T e e T ∂ = []p x R [t] ) dx với px [][]1 x = 0 1 (7) 2 ∫xi ∂x - Ta sẽ áp dụng quy tắc đạo hàm của hàm vô h−ớng (tức ma trận (1x1) nh− là tích hai vectơ A [t]e ) theo ma trận [t]e : (1x2) (2x1) ⎡ ∂ ⎤ (a t + a t ) e e ⎢ 1 i 2 j ⎥ d(A[t] ) d ⎡t i ⎤ ∂t i a = ) ⎢ ⎥ ⎡ 1 ⎤ T e e ([a1a2] ⎢t ⎥ = ⎢ ∂ ⎥ = ⎢ ⎥ = A d[t] d[t] ⎣ j ⎦ ⎣a2 ⎦ ⎢ (a1t i + a 2 t j ) ⎥ ⎣∂tj ⎦ Chú ý: Ký hiệu AT hiểu là chuyển vị của ma trận A, và quy tắc chuyển vị tích (AB)T = BTAT (không giao hoán): e e 2 e x x dIλ λ j d T e e λ j T e e T e T (px R [t] ) 2( px R [t] )( px R ) - Tính e = ∫x e dx= ∫x dx d[t] 2 i d[t] 2 i x T x T e j T e e e e j e T e e e e λ (p R [t] )(R p ) λ (R px )( px R [t] ) = ∫x x x dx= ∫x dx, vì R , [t] , i i px không phụ thuộc vào x nên: ký hiệu xij = (xj - xi) ta có: e ⎡x −1⎤ dI x j T e 1 j ⎡0⎤ λ e e T e e λ x . e = λ dx.R px px R [t] = ij ⎢ ⎥⎢ ⎥ d[t] ∫xi x − x 1 1 ij ⎣⎢ i1 ⎦⎥⎣ ⎦ e 1 −1 1 ⎡x j − x i ⎤ λ ⎡ ⎤ e e e .[01] [t]e = [t] = K [t] ⎢ ⎥ x ⎢−1 1⎥ x ij ⎣−1 1 ⎦ ij ⎣ ⎦ λe ⎡1 −1⎤ Với Ke ⎢ ⎥ (8) (đối xứng) là ma trận dẫn nhiệt của phần x ij ⎣−1 1⎦ ⎡t1 ⎤ ⎢ ⎥ ⎡t ⎤ 0 .10 0 t T tử e. - Chú ý rằng [t]e = i = ⎡ ⎤ ⎢ i ⎥ = De [t] ta có: ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ t j ⎣0 . 01 0⎦ t j ⎣ ⎦ ⎢ ⎥ ⎣⎢tM ⎦⎥ E e E dI dI T λ = De λ e e e ∑ e = ∑ D K D [t] = K[t] d[t] e=1 d[t] e=1 89
Với K là ma trận vuông (MxM) đối xứng, gọi là ma trận dẫn nhiệt: ⎡00⎤ ⎢ ⎥ ⎢ ⎥ E ⎢10⎥ ⎡i j ⎤ T K K ⎡0 .10 0⎤ E e e e E ⎡ 11 12 ⎤ ⎢ ⎥ K D K D ⎢ ⎥ ⎢ ⎥ ⎢ ⎥ = K K i ∑ ∑⎢01⎥ K K 0 . 01 0 ∑⎢ 11 12 ⎥ e=1 e1= ⎣ 21 22 ⎦ ⎣ ⎦ e=1 ⎢ ⎥ ⎣⎢K21 K 22 ⎦⎥ j ⎢ ⎥ ⎢ ⎥ ⎣⎢00⎦⎥ e Nếu λ = const, ∀e E ⎡i j ⎤ λ ⎢ ⎥ → K = ∑ 1 −1 i ⎢ ⎥ xij = xj - xi = ∆x = const ∆x e1= ⎣⎢−1 1⎦⎥j E e dIc e dIc 6.4.3.4. Tính = ∑ D e ta có: d[t] e=1 d[t] e 2 e x e 1 j ∂(t ) (pc) d x j T e e 2 Ic = dx {p R [t] } dx ∫x = ∫x → 2 i ∂τ 2 dτ i e e xj xj dI (C)dρ e T c = 2(pTe R [t] e )(p TeT R ) dx (ρC) RppR[t]dxeTee e ∫ = ∫ = d[t] 2 dτ xi dτ xi xj d ⎪⎧ eTeeT =(ρC)e ⎨R∫ pp dxR [t] , ở đây: dτ ⎩⎪ xi ⎡ 1 ⎤ (x−− x ) (x22 x ) xj xj xj ⎢ ji j i⎥ ⎡⎤1 ⎡1x⎤ 2 T dx ⎢ ⎥ ppdx = ⎢⎥[]1.x dx = ⎢ 2 ⎥ = = ∫ ∫ x ∫ 1122 33 xi xi ⎣⎦ xi ⎣xx⎦ ⎢ (x−− x ) (x x ) ⎥ ⎣⎢23ji ji⎦⎥ ⎡⎤1 1(xx)+ ⎢⎥2 ji ⎢⎥e = xij 11 . Vì chỉ [t] phụ thuộc thời gian τ ⎢⎥(x+++ x ) (x22 x x x ) ⎣⎦⎢⎥23ji j jii 90