ویکی‌جزوه/دانشکده:فنی و مهندسی/انتقال حرارت/حل عددی مسائل هدایت، پایا و گذرا

روش تفاضل محدود برای هدایت دایمی

گسسته سازی معادله گرما

$\begin{matrix} \frac{\partial T}{\partial x} = \frac{Δ T}{Δ x} \\ T_{i + 1} = T_{i} + {\frac{\partial T}{\partial X} |}_{i} (x_{i + 1} - x_{i}) + \dots {(Δ x)}^{2} \\ Δ T = T_{i + 1} - T_{i} = {\frac{\partial T}{\partial x} |}_{i} Δ x \end{matrix}$

.تقریب پیشرو:

$\begin{matrix} {\frac{\partial T}{\partial x} |}_{i} ≃ \frac{T_{i + 1} - T_{i}}{Δ x} \end{matrix}$

خطا از مرتبهٔ $Δ x$ است.

تقریب پسرو:

$\begin{matrix} {\frac{\partial T}{\partial x} |}_{i} ≃ \frac{T_{i} - T_{i + 1}}{Δ x} \end{matrix}$

خطااز مرتبهٔ $Δ x$ است.

$\begin{matrix} \int_{C . V .} ρ c \frac{\partial T}{\partial t} d V = \int_{C . V .} \nabla . (K \nabla T) d V \\ ρ c \frac{\partial T}{\partial t} V_{C . V .} = \int_{C . S .} K \nabla T . d \vec{s} = - \sum_{f} \vec{q^{″}} . \vec{s} \\ Δ x \to 0 \Rightarrow T = c t e . \\ ρ c \frac{\partial T}{\partial t} V_{C . V .} = - {({q^{″}}_{x} s)}_{e} - {(- {q^{″}}_{x} s)}_{w} = {(K \frac{d T}{d x} s)}_{e} - {(K \frac{d T}{d x} s)}_{w} \\ ρ c \frac{T_{i}^{n + 1} - T_{i}^{n}}{Δ t} A Δ x = K A {(\frac{d T}{d x})}_{e} - K A {(\frac{d T}{d x})}_{w} \\ {(\frac{d T}{d x})}_{e} = \frac{T_{i + 1} - T_{i}}{Δ x}, {(\frac{d T}{d x})}_{w} = \frac{T_{i} - T_{i - 1}}{Δ x} \\ \Rightarrow ρ c \frac{T_{i}^{n + 1} - T_{i}^{n}}{Δ t} = K \frac{T_{i + 1} - T_{i}}{Δ x} + K \frac{T_{i - 1} - T_{i}}{Δ x} \end{matrix}$

در اینجا برای حل مسایل از روش حجم محدود استفاده می‌کنیم که معادلهٔ بقای حجم را روی سلول می‌نویسیم.

درحالت دایمی که تولید نیز داشته باشیم:

$T_{i + 1} - 2 T_{i} + T_{i - 1} + \frac{{q^{″}}_{g e n}}{K} Δ x^{2} = 0$

حل عددی برای هدایت دایمی یک بعدی با روش حجم محدود

$\begin{matrix} K \frac{T_{i + 1} - T_{i}}{Δ x} + K \frac{T_{i - 1} - T_{i}}{Δ x} = 0 \\ i = 0 \Rightarrow K \frac{T_{1} - T_{0}}{Δ x} + h (T_{\infty} - T_{0}) = 0 \\ i = N \Rightarrow T = 100^{\circ} C \end{matrix}$

دستگاه معادلات

$\begin{matrix} A T = B \end{matrix}$

برای هر نقطه یا معادله یا مقدار آن را داریم.

n معادله و n مجهول داریم.

$\begin{matrix} A = (\begin{matrix} - (1 + h \frac{Δ x}{K}) & 1 & 0 & 0 & 0 & \dots \\ 1 & - 2 & 1 & 0 & 0 & \dots \\ 0 & 1 & - 2 & 1 & 0 & \dots \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ \dots & 0 & 0 & 1 & - 2 & 1 \\ \dots & 0 & 0 & 0 & 0 & 1 \end{matrix}) \\ T = (\begin{matrix} T_{0} \\ T_{1} \\ T_{2} \\ T_{3} \\ ⋮ \\ T_{N} \end{matrix}), B = (\begin{matrix} - h \frac{Δ x}{K} T_{\infty} \\ 0 \\ 0 \\ 0 \\ ⋮ \\ 0 \end{matrix}) \end{matrix}$

با حل دستگاه بالا دمای کل گره‌ها رو بدست می‌آوریم.

مثال‌ها

مثال1

با استفاده از روش حجم محدود (شکل ۱) معادله گرما بصورت زیر تبدیل خواهد شد:

$\begin{matrix} \frac{d (k \frac{d T}{d x})}{d x} + q_{g e n}^{'^{″}} = 0 \\ Δ x = \frac{L}{N} \\ q_{l} = - k \frac{\partial T}{\partial x})_{l} A = - k A \frac{T - T_{0}}{Δ x} \\ \frac{\partial T}{\partial x} ≃ \frac{T - T_{0}}{Δ x} \\ q_{r} = - k \frac{\partial T}{\partial x})_{r} A = - k A \frac{T_{1} - T_{2}}{Δ x} \end{matrix}$

از بقای انرژی برای هر سلول داریم:

قرارداد: از آنجا که جهت شار گرمایی معلوم نیست فرض را بر این می‌گذاریم که تمام شارهای گرمایی بصورت ورودی به سیستم باشند.

$\begin{matrix} \sum q = - k A \frac{T - T_{0}}{Δ x} - k A \frac{T_{1} - T_{2}}{Δ x} = k A \frac{T_{0} - 2 T_{1} + T_{2}}{Δ x} \\ q_{g e n} = q_{g e n}^{'^{″}} Δ x A \\ k \frac{T_{0} - 2 T_{1} + T_{2}}{Δ x^{2}} + q_{g e n_{1}}^{'^{″}} = 0 \\ k \frac{T_{1} - 2 T_{2} + T_{3}}{Δ x^{2}} + q_{g e n_{2}}^{'^{″}} = 0 \\ k \frac{T_{2} - 2 T_{3} + T_{4}}{Δ x^{2}} + q_{g e n_{3}}^{'^{″}} = 0 \\ T_{0} = T_{0} \\ T_{4} = T_{l} \\ - 2 T_{1} + T_{2} = - q_{g e n}^{'^{″}} \frac{Δ x^{2}}{k} - T_{0} \\ T_{1} - 2 T_{2} + T_{3} = - q_{g e n}^{'^{″}} \frac{Δ x^{2}}{k} \\ T_{2} - 2 T_{3} = - q_{g e n}^{'^{″}} \frac{Δ x^{2}}{k} - T_{l} \\ A T = B \\ A = [\begin{matrix} - 2 & 1 & 0 \\ 1 & - 2 & 1 \\ 0 & 1 & - 2 \end{matrix}] \\ T = [\begin{matrix} T_{1} \\ T_{2} \\ T_{3} \end{matrix}] \\ B = - \frac{Δ x^{2}}{k} [\begin{matrix} q_{g e n_{1}}^{'^{″}} \\ q_{g e n_{2}}^{'^{″}} \\ q_{g e n_{3}}^{'^{″}} \end{matrix}] + [\begin{matrix} T_{0} \\ 0 \\ T_{l} \end{matrix}] \end{matrix}$

برای بدست آوردن مقدار گرمای وارد یا خارج شده از دیوار، با نوشتن بقای انرژی برای قسمت هاشور خورده خواهیم داشت:

$\begin{matrix} q + k A \frac{T_{1} - T_{0}}{Δ x} + q_{g e n_{0}}^{'^{″}} \frac{Δ x}{2} A = 0 \\ q = - A [k \frac{T_{1} - T_{0}}{Δ x} + q_{g e n_{0}}^{'^{″}} \frac{Δ x}{2}] \end{matrix}$

مثال2

دمای نقاط نمایش داده شده در شکل ۲ را به روش عددی در دیوار مسطح زیر گسسته کنید.

نقطهٔ صفر:

$T_{0} = 520 R$

نقطهٔ یک:

$\frac{T_{0} - 2 T_{1} + T_{2}}{Δ x^{2}} = 0$

نقطهٔ دو:

$\frac{T_{1} - 2 T_{2} + T_{3}}{Δ x^{2}} = 0$

نقطهٔ سه:

$\begin{matrix} 0 = k A \frac{T_{2} - T_{3}}{Δ x} + α q_{s} - [ε σ ({T_{s u r r}}^{4} - {T_{3}}^{4}) = h_{r} A (T_{s u r r} - T_{3})] \\ h_{r} = ε σ ({T_{s u r r}}^{2} + {T_{3}}^{2}) (T_{s u r r} + T_{3}) \end{matrix}$

$\begin{matrix} T_{2} - (1 - \frac{h_{r} Δ x}{k}) T_{3} = - \frac{h_{r} Δ x}{k} - α {q_{s}}^{'^{″}} \\ [\begin{matrix} - 2 & 1 & 0 \\ 1 & - 2 & 1 \\ 0 & 1 & - (1 - \frac{h_{r} Δ x}{k}) \end{matrix}] [\begin{matrix} T_{1} \\ T_{2} \\ T_{3} \end{matrix}] = [\begin{matrix} - T_{0} \\ 0 \\ \frac{h_{r} Δ x}{k} T_{s u r r} - α {q_{s}}^{'^{″}} \end{matrix}] \end{matrix}$

حل عددی در هندسه‌های متغییر دایمی یک بعدی با روش حجم محدود

$\begin{matrix} \frac{d (k A (x) \frac{d T}{d x})}{d x} + q_{g e n}^{'^{″}} A (x) = 0 \\ q_{l} = - k A {(x)}_{i - \frac{1}{2}} {\frac{d T}{d x}}_{i - \frac{1}{2}} = - k A (x_{i} - \frac{Δ x}{2}) \frac{T_{i} - {T_{i - 1}}_{}}{Δ x} \\ q_{l} = k A {(x)}_{i - \frac{1}{2}} \frac{T_{i - 1} - T_{i}}{Δ x} \\ q_{r} = k A {(x)}_{i + \frac{1}{2}} \frac{T_{i + 1} - T_{i}}{Δ x} \\ \Rightarrow q {_{g e n}^{'^{″}}}_{, i} A_{i} + k A_{i - \frac{1}{2}} \frac{T_{i - 1} - T_{i}}{Δ x^{2}} + k A_{i + \frac{1}{2}} \frac{T_{i + 1} - T_{i}}{Δ x} = 0 \end{matrix}$

حل عددی هدایت دایمی دو بعدی با روش حجم محدود

هدف: حل عددی معادله پخش گرما در حالت دایمی

راه حل: المانی را به صورت زیر در نظر می‌گیریم. و معادله پخش گرما در حالت دایمی را برای ان می‌نویسیم. از معادله پخش گرما برای حالت دایمی روی حجم انتگرال می‌گیریم.

پرونده:۱.jpg

معادله پخش گرما برای حالت دایمی:

$\vec{\nabla} . (k \vec{\nabla} T) + q_{g e n}^{'^{″}} = 0$

انتگرال گیری روی حجم از معادله پخش گرما:

$\int_{c . v .} [\vec{\nabla} . (k \vec{\nabla} T) + q_{g e n}^{'^{″}}] d V = 0 \to \int_{c . s .} k \vec{\nabla} T . d \vec{s} + \int_{V} q_{g e n}^{'^{″}} d V = 0$

$\sum (k \vec{\nabla} T)_{f} . \vec{s_{f}} + q_{g e n}^{'^{″}} V = 0$

که f همان W، E، S و N می‌باشد.

${(k \frac{\partial T}{\partial x})}_{E} S_{E} + {(k \frac{\partial T}{\partial y})}_{N} S_{N} - {(k \frac{\partial T}{\partial x})}_{W} S_{W} - {(k \frac{\partial T}{\partial y})}_{S} S_{S} + q_{g e n}^{'^{″}} V = 0$

با فرض ثابت بودن k، معادلات را می‌توان به صورت زیر نوشت.

$k (\frac{T_{i + 1, j} - T_{i, j}}{Δ x}) b Δ y - k (\frac{T_{i, j} - T_{i - 1, j}}{Δ x}) b Δ y + k (\frac{T_{i, j + 1} - T_{i, j}}{Δ y}) b Δ x - k (\frac{T_{i, j} - T_{i, j - 1}}{Δ y}) b Δ x + q_{g e n}^{'^{″}} Δ x Δ y b = 0$

روش دوم:

می‌توانیم برای یک المان همه انتقال حرارتی را که با اطراف دارد، به صورت ورودی در نظر بگیریم و با هم جمع کنیم:

$\begin{matrix} q_{e} + q_{w} + q_{n} + q_{s} + q_{g e n, i, j}^{'^{″}} Δ x Δ y b = 0 \\ q_{e} = k Δ y b \frac{T_{i + 1, j} - T_{i, j}}{Δ x} \\ q_{w} = k Δ y b \frac{T_{i - 1, j} - T_{i, j}}{Δ x} \\ q_{n} = k Δ x b \frac{T_{m, n + 1} - T_{m, n}}{Δ y} \\ q_{s} = k Δ x b \frac{T_{m, n - 1} - T_{m, n}}{Δ y} \\ \Rightarrow k Δ x Δ y [\frac{T_{m + 1, n} - 2 T_{m, n} + T_{m - 1, n}}{Δ x^{2}} + \frac{T_{m, n - 1} - 2 T_{m, n} + T_{m, n + 1}}{Δ y^{2}}] + q_{g e n_{m, n}}^{'^{″}} Δ x Δ y b = 0 \\ k (\frac{\partial^{2} T}{\partial^{2} x} + \frac{\partial^{2} T}{\partial^{2} y}) + q_{g e n}^{'^{″}} = 0 \\ i f : Δ x = Δ y \to k [T_{m + 1, n} + T_{m - 1, n} + T_{m, n - 1} + T_{m, n + 1} + - 4 T_{m, n}] + q_{g e n}^{'^{″}} Δ x^{2} = 0 \end{matrix}$

مثال3

در شکل زیر با توجه به اطلاعات داده شده:

الف: دمای نقاط را بدست آورید.

ب: شار ورود به سیستم را بدست آورید.

$T_{1} = T_{2} = T_{3} = T_{4} = T_{0}$

نقطه ۵:

{

$\begin{matrix} q_{w} = 0 \\ q_{e} = k Δ y b \frac{T_{6} - T_{5}}{Δ x} \\ q_{n} = k Δ x b \frac{T_{9} - T_{5}}{Δ y} \\ q_{s} = k Δ x b \frac{T_{1} - T_{5}}{Δ y} \\ Δ x = Δ y \\ q_{w} + q_{e} + q_{n} + q_{s} = 0 \\ k b [0 + 2 T_{6} - T_{1} + T_{9} - 4 T_{5}] = 0 \\ \to T_{6} + T_{9} - 4 T_{5} = T_{0} \end{matrix}$

نقطه ۶:

$T_{5} + T_{10} + T_{7} + T_{2} - 4 T_{6} = 0$

نقطه ۷:

$T_{6} + T_{11} + T_{8} + T_{3} - 4 T_{7} = 0$

نقطه ۸:

$2 T_{7} + T_{12} + T_{4} - 4 T_{8} = 0$

نقطه ۹:

$T_{13} + 2 T_{10} + T_{5} - 4 T_{9} = 0$

نقطه ۱۰:

$T_{9} + T_{14} + T_{11} + T_{6} - 4 T_{10} = 0$

نقطه ۱۱:

بقای انرژی:

$q_{w} + q_{e} + q_{n} + q_{s} + q = 0$
$\begin{matrix} k Δ y b \frac{T_{10} - T_{11}}{Δ x} + k Δ x b \frac{T_{7} - T_{11}}{Δ y} + k \frac{Δ y}{2} b \frac{T_{12} - T_{11}}{Δ x} + h (Δ x + Δ y) \frac{b}{2} (T_{\infty} - T_{11}) + k \frac{Δ x}{2} b \frac{T_{15} - T_{11}}{Δ y} = 0 \\ \to k (T_{10} + T_{7} + T_{12} + T_{15} - 3 T_{11}) + h Δ x (T_{\infty} - T_{11}) = 0 \\ \div \frac{1}{k} \to (T_{10} + T_{7} + T_{12} + T_{15} - 3 T_{11}) + \frac{h Δ x}{k} (T_{\infty} - T_{11}) = 0 \end{matrix}$

نقطه ۱۲:

بقای انرژی:

$\begin{matrix} q_{w} + q_{e} + q_{n} + q_{s} = 0 \\ \frac{h Δ x}{2} (T_{\infty} - T_{12}) + 0 + k \frac{Δ x}{2} b \frac{T_{8} - T_{12}}{Δ y} + k \frac{Δ y}{2} b \frac{T_{11} - T_{12}}{Δ x} = 0 \\ \to T_{8} + T_{11} - (2 + \frac{h Δ x}{k}) T_{12} = - \frac{h Δ x}{k} T_{\infty} \end{matrix}$

نقطه ۱۳:

$T_{14} + T_{9} - (2 + \frac{h Δ x}{k}) T_{13} = - \frac{h Δ x}{k} T_{\infty}$

نقطه ۱۴:

$\frac{1}{2} (T_{13} + T_{15}) + T_{10} - (2 + \frac{h Δ x}{k}) T_{14} = - \frac{h Δ x}{k} T_{\infty}$

نقطه ۱۵:

$\frac{1}{2} (T_{14} + T_{11}) - (1 + \frac{h Δ x}{k}) T_{15} = - \frac{h Δ x}{k} T_{\infty}$

حال با حل ماتریس زیر به روش ماتریس معکوس یا هر روشی که معتبر باشد می‌توان دما در نقاط مختلف را بدست آورد:

$A T = B$

A: ماتریس ضرایب
T: ماتریس دما
B: ماتریس معلومات

حال برای بدست آوردن شار ورودی به سیستم، المانی از پایین جسم جدا می‌کنیم.(از آنجا که انتقال حرارت پایا می‌باشد پس شار ورودی و شار خروجی یکسان می‌باشند.)

$\begin{matrix} q_{w} + q_{e} + q_{n} + q_{s} = 0 \\ 0 + k (\frac{Δ x}{2} b \frac{T_{5} - T_{0}}{Δ y} + Δ x b \frac{T_{6} - T_{0}}{Δ y} + Δ x b \frac{T_{7} - T_{0}}{Δ y} + \frac{Δ x}{2} b \frac{T_{8} - T_{0}}{Δ y}) + 0 + q_{s} = 0 \\ q_{s} = - k b (\frac{T_{5}}{2} + T_{6} + T_{7} + T_{8} + - 3 T_{5}) \end{matrix}$

مثال4

در شکل نشان داده شده اگر T₁=430 و T₃=394 و T₆=492 و ضریب رسانندگی برابر ۱w/m.k و ضریب انتقال جابجایی برابر ۵۰w/m² .k باشد:

الف) دمای سایر نقاط را بدست آورید.

ب) میزان کل انتقال حرارت را بدست آورید.

$\begin{matrix} Δ x = Δ y = 0.01 m \\ T_{1} = 430 k \\ T_{3} = 394 k \\ T_{6} = 492 k \\ K = 1 \frac{w}{m k} \\ h = 50 \frac{w}{m k^{2}} \\ T_{i n} = 600 k \\ T_{\infty} = 300 k \\ b=1 m \end{matrix}$

با توجه به شکل می‌توان دمای نقطهٔ ۸ و ۹ را به راحتی محاسبه کرد:

$\begin{matrix} T_{8} = T_{9} = T_{in} = 600 k \end{matrix}$

حال از روش تفاضل محدود معادلات را بدست می‌آوریم، یعنی برای هر نقطه تمام انتقال حرارتی را که به اطراف دارد به صورت وارد به نقطه مورد نظر، در نظر می‌گیریم.

پس داریم:

$\begin{matrix} point 2; K \frac{T_{1} - T_{2}}{Δ x} b \frac{Δ y}{2} + K \frac{T_{3} - T_{2}}{Δ x} b \frac{Δ y}{2} + K \frac{T_{6} - T_{2}}{Δ y} b Δ x + hb Δ x (T_{\infty} - T_{2}) = 0 \\ point 4; K \frac{T_{3} - T_{4}}{Δ x} b \frac{Δ y}{2} + hb \frac{Δ x}{2} (T_{\infty} - T_{4}) = 0 \\ point 5; K \frac{T_{1} - T_{5}}{Δ x} b Δ y + K \frac{T_{8} - T_{5}}{Δ y} b \frac{Δ x}{2} + K \frac{T_{1} - T_{5}}{Δ y} b \frac{Δ x}{2} = 0 \\ point 7; K \frac{T_{6} - T_{7}}{Δ x} b Δ y + K \frac{T_{3} - T_{7}}{Δ y} b Δ x = 0 \end{matrix}$

حال با حل دستگاه ۴ معادله ۴ مجهول از طریق حل ماتریس AX=B که

A ماتریس ضرایب

B ماتریس معلومات

X ماتریس دما

است به روش معکوس داریم:

$\begin{matrix} \to T_{2} = 422 k \\ T_{4} = 363 k \\ T_{5} = 503.5 k \\ T_{7} = 443 k \end{matrix}$

حال برای محاسبهٔ انتقال حرارت به دلیل پایا بودن انتقال حرارت، هم می‌توان انتقال حرارت از نقاط ۱و ۲و ۳و ۴ بدست آورد و هم از نقاط ۸و ۹. پس داریم:

$\begin{matrix} q_{8} + K \frac{T_{9} - T_{8}}{Δ x} b \frac{Δ y}{2} + K \frac{T_{5} - T_{8}}{Δ y} b \frac{Δ x}{2} +0=0 \\ q_{9} + K \frac{T_{8} - T_{9}}{Δ x} b \frac{Δ y}{2} + K \frac{T_{6} - T_{9}}{Δ y} b Δ x =0 \end{matrix}$

میدانیم که دمای نقطهٔ ۸و۹ با هم برابر پس:

$\begin{matrix} q_{8} = 48.25 w \\ q_{9} = 108 w \\ q = 8 (q_{8} + q_{9}) = 8 (156.25) = 1250 w \end{matrix}$

روش تفاضل محدود برای هدایت گذرا

روش صریح

گسسته‌سازی معادله گرما

در روش صریح تفاضل محدود دمای هر گره در زمان $t + Δ t$ را از اطلاعات مربوط به دمای هر گره و گره‌های مجاور در زمان قبلی t می‌توان محاسبه کرد (از اندیس P برای نشان دادن وابستگی T به زمان استفاده می‌شود)

$\frac{1}{α} \frac{T_{m, n}^{P + 1} - T_{m, n}^{P}}{Δ t} = \frac{T_{m + 1, n}^{P} + T_{m - 1, n}^{P} - 2 T_{m, n}^{P}}{(Δ x)^{2}} + \frac{T_{m, n + 1}^{P} + T_{m, n - 1}^{P} - 2 T_{m, n}^{P}}{(Δ y)^{2}}$

با حل گره در زمان جدید (P+1) و با فرض اینکه $Δ x = Δ y$ داریم:

$T_{m, n}^{P + 1} = (1 - 4 F o) T_{m, n}^{P} + F o (T_{m + 1, n}^{P} + T_{m - 1, n}^{P} + T_{m, n + 1}^{P} + T_{m, n - 1}^{P})$

که در آن Fo شکل تفاضل محدود عدد فوریه است.

$F o = \frac{a Δ t}{(Δ x)^{2}}$

اگر سیستم یک بعدی و در جهت x باشد شکل صریح معادله تفاضل برای گره داخلی m به صورت زیر است:

$T_{m}^{P + 1} = (1 - 2 F o) T_{m}^{P} + F o (T_{m - 1}^{P} + T_{m + 1}^{P})$

در معادله‌های صریح دماهای مجهول گره‌ها در زمان جدید فقط توسط دماهای معلوم گره‌ها در زمان قبل تعیین می‌شوند. لذا محاسبه دماهای مجهول سر راست است.

با کاهش مقادیر $Δ x$ و $Δ t$ می‌توان دقت حل تفاضل محدود را افزایش داد. البته با کاهش $Δ x$ تعداد گره‌هایی که باید بررسی شوند افزایش و با کاهش $Δ t$ تعداد بازه‌های زمانی برای نیل به زمان نهایی افزایش می‌یابد

یک عیب روش صریح این است که پایداری آن بی قید و شرط نیست. در مسایل گذرا با گذشت زمان دمای گره‌ها باید به طور پیوسته به مقادیر نهایی نزدیک تر شوند ولی با روش صریح این حل با نوسان‌های عددی (که بطور فیزیکی غیر ممکن هستند) مشخص می‌شوند. نوسان‌ها ممکن است ناپایدار باشند و در نتیجه جواب‌ها از شرایط پایای واقعی دور شوند. برای جلوگیری از این نتایج گول زننده $Δ t$ را باید از یک حد مشخص کمتر در نظر گرفت. این وابستگی زمان را معیار پایداری گویند.

معیار پایداری ایجاب می‌کند که ضریب مربوط به گره مورد نظر در زمان قبل بزرگ تر یا مساوی با صفر باشد.

در نتیجه برای یک گره داخلی یک بعدی Fo باید کوچکتر یا مساوی $\frac{1}{2}$ باشد

و برای یک گره داخلی دو بعدی Fo باید کوچکتر یا مساوی $\frac{1}{4}$ باشد.

مطابق شکل بالا یک دیوار با انتقال حرارت یک بعدی گذرا مفروض می باشد و مقادیر دمای نقطه صفر و ضریب جابه جایی معلوم می‌باشد بنا به حل عددی صریح می‌توان مسئله را به صورت زیر برای نقاظ نشان داده شده حل نمود:
گره صفر: $𝑇_{0}$

گره ۱:

$ρ c A Δ x \frac{\partial T_{1}}{\partial t} = ρ c A Δ x \frac{T_{1}^{n + 1} - T_{1}^{n}}{Δ t} = K A \frac{T_{0}^{n} - T_{1}^{n}}{Δ x} + K A \frac{T_{2}^{n} - T_{1}^{n}}{Δ x}$

$\Rightarrow T_{1}^{n + 1} - T_{1}^{n} = (1 - {2 f}_{o} {) T}_{1}^{n} + f_{o} T_{0}^{n} + f_{o} T_{2}^{n} \Rightarrow 1 - {2 f}_{o} > 0$

$f_{o} T_{0}^{n} + f_{o} T_{2}^{n} \Rightarrow 1 - {2 f}_{o} > 0$

گره ۲:

$T_{2}^{n + 1} = (1 - {2 f}_{o} {) T}_{2}^{n} + f_{o} {(T}_{1}^{n} + T_{3}^{n})$

گره ۳:

$\frac{ρ_{c}}{Δ t} \frac{Δ x A}{2} {(T}_{3}^{n + 1} - T_{3}^{n}) = K A \frac{T_{2}^{n} - T_{3}^{n}}{Δ x} + {h A (T}_{\infty} - T_{3}^{n})$

$\Rightarrow T_{3}^{n + 1} = (1 - {2 f}_{o} - {2 f}_{o} {B i) T}_{3}^{n} + {2 f}_{o} T_{2}^{n} + {2 f}_{o} {B i T}_{\infty}$

همانطور که قبلا ذکر شد روش حل صریح پایداری مشروط دارد و گام زمانی آن را مشابه زیر برای تمامی نقاط می‌توان بدست آورد. کمترین بازه زمانی شرط پایداری می‌باشد.

$\Rightarrow 1 - {2 f}_{o} (1 + B i) > 0 \Rightarrow f_{0} < \frac{1}{2 (B i + 1)} \Rightarrow Δ t < \frac{{Δ x}^{2}}{2 (B i + 1) α}$

مثال‌ها

مثال 5

نمای بالای دودکشی را ملاحظه می‌کنید که می‌تواند مصداقی از انتقال حرارت دو بعدی گذرا باشد و دمای گازهای خارج شده از دودکش و محیط و خواص در ادامه آمده است با فرض مقدار ضریب جابه جایی مطابق زیر حداقل گام زمانی برای پایداری را بدست آورید.
حل: با توجه به متقارن بودن دهانه خروجی دودکش می‌توان مسئله را برای $\frac{1}{8}$ حل نمود و تعمیم داد:

$Δ x = Δ y = 0.01$

$ρ = 2000 \frac{k g}{m^{3}} {, T}_{i} = 300 k$

$h = 50 \frac{w}{m^{2} k}$

$k = 1 \frac{w}{m k}$

$ρ c Δ V \frac{T^{n + 1} - T^{n}}{Δ t} = \sum q_{i n}$

گره ۱:

$T_{1}^{n + 1} - T_{1}^{n} = {2 f}_{o} {(T}_{5}^{n} - T_{1}^{n}) + {2 f}_{o} {(T}_{2}^{n} - T_{1}^{n}) + {2 B i f}_{o} {(T}_{\infty} - T_{1}^{n})$

گره ۲:

$T_{2}^{n + 1} - T_{2}^{n} = {2 f}_{o} {(T}_{6}^{n} - T_{2}^{n}) + f_{o} {(T}_{1}^{n} + T_{3}^{n} + {2 T}_{3}^{n}) + {2 B i f}_{o} {(T}_{\infty} - T_{2}^{n})$

گره ۳ همانند گره ۲

گره ۴:

$T_{4}^{n + 1} - T_{4}^{n} = {4 f}_{o} {(T}_{3}^{n} - T_{4}^{n}) + {4 B i f}_{o} {(T}_{\infty} - T_{4}^{n})$

گره ۵:

$T_{5}^{n + 1} - T_{5}^{n} = {4 f}_{o} {(T}_{6}^{n} - T_{5}^{n}) + f_{o} {(T}_{1}^{n} + T_{8}^{n} - {2 T}_{5}^{n})$

گره ۶:

$T_{6}^{n + 1} - T_{6}^{n} = f_{o} {(T}_{5}^{n} + T_{7}^{n} + T_{2}^{n} - T_{9}^{n} - {4 T}_{6}^{n})$

گره ۷:

$T_{7}^{n + 1} - T_{7}^{n} = {2 f}_{o} {(T}_{6}^{n} + T_{3}^{n} - {2 T}_{7}^{n})$

گره ۱:

$f_{o} < \frac{1}{2 (2 + B i)}$

گره ۲ و۳:

$f_{o} < \frac{1}{2 (2 + B i)}$

گره ۴:

$f_{o} < \frac{1}{4 (1 + B i)}$

گروه ۵ و ۶ و ۷:

$f_{o} < \frac{1}{4}$

$B i = \frac{50 \times 0.01}{1} = 0.5$

$\Rightarrow f_{o} < \frac{1}{4 (1 + B i)} \Rightarrow Δ t < \frac{{(0.01)}^{2}}{10^{- 7} (1 + 0.5)} = 166 s$

روش ضمنی

گسسته سازی معادله گرما

در روش صریح تفاضل محدود دمای هر گره در زمان $t + Δ t$ را از اطلاعات مربوط به دمای هر گره وگره‌های مجاور در زمان قبلی $t$ می‌توان محاسبه کرد. لذا تعیین دمای یک گره در هر لحظه مستقل از دمای سایر گره‌ها در همان لحظه است، گر چه محاسبه با این روش ساده است ولی محدودیت‌هایی برای انتخاب $Δ t$ وجود دارد.

برای یک نمو مکانی معین بازه زمانی باید با شرایط پایداری سازگار باشد. در نتیجه اغلب باید از مقادیر بسیار کوچک $Δ t$ استفاده شود و برای رسیدن به حل مسئله تعداد بسیار زیادی بازه زمانی لازم می‌شود. اغلب با استفاده از روش ضمنی تفاضل محدود به جای روش صریح می‌توان زمان محاسبه را کاهش داد.

با ارزیابی سایر دماها در زمان جدید $p + 1$ به جای زمان قبلی $p$ می‌توان شکل ضمنی معادله تفاضل محدود را بدست آورد.

شکل ضمنی برای یک گره داخلی یک سیستم دو بعدی عبارت است از:

$\frac{1}{α} \frac{T_{m, n}^{P + 1} - T_{m, n}^{P}}{Δ t} = \frac{T_{m + 1, n}^{P + 1} + T_{m - 1, n}^{P + 1} - 2 T_{m, n}^{P + 1}}{(Δ x)^{2}} + \frac{T_{m, n + 1}^{P + 1} + T_{m, n - 1}^{P + 1} - 2 T_{m, n}^{P + 1}}{(Δ y)^{2}}$

با مرتب کردن و به فرض $Δ x = Δ y$ نتیجه می‌شود:

$(1 + F o) T_{m, n}^{P + 1} - F o (T_{m + 1, n}^{P + 1} + T_{m - 1, n}^{P + 1} + T_{m, n + 1}^{P + 1} + T_{m, n - 1}^{P + 1}) = T_{m, n}^{P}$

از معادله فوق واضح است که دمای جدید $(m, n)$ به دماهای جدید گره‌های مجاورش که معمولا مجهولند بستگی دارد. لذا برای تعیین دماهای مجهول گره‌ها در $t + Δ t$ معادله‌های گره‌ای متناظر را باید همزمان حل کرد.

این حل با استفاده از روش گوس-سیدل یا با معکوس سازی ماتریس انجام می‌شود. لذا حل پیش رو شامل حل همزمان معادله‌های گره‌ای در لحظات $t = Δ t, 2 Δ t$ است که تا نیل به زمان نهایی دلخواه ادامه یابد.

فرمول بندی ضمنی نسبت به روش صریح این مزیت مهم را دارد که بدون قید و شرط همواره برای تمام بازه‌های مکانی و زمانی پایدار است. و هیچ قیدی برای $Δ x$ و $Δ t$ وجود ندارد.

چون در روش ضمنی می‌توان از مقادیر بزرگتر $Δ t$ استفاده کرد زمان محاسبه کاهش یافته و دقت نیز کمی کاهش می‌یابد.

با وجود این برای دستیابی به دقت زیاد $Δ t$ باید به اندازه کافی کوچک باشد تا نتایج از کاهش بیشتر $Δ t$ مستقل شوند.

شکل ضمنی را با روش موازنه انرزی نیز می‌توان بدست آورد. به عنوان مثال برای جابجایی و رسانش گذرای یک بعدی به صورت زیر است.

برای گره سطحی

$(1 + 2 F o + 2 F o B i) T_{0}^{P + 1} - 2 F o T_{1}^{P + 1} = 2 F o B i T_{\infty} + T_{0}^{P}$

برای هر گره داخلی

$(1 + 2 F o) T_{m}^{P + 1} - F o (T_{m - 1}^{P + 1} + T_{m + 1}^{P + 1}) = T_{m}^{P}$

انتقال حرارت یک بعدی

توضیحات داده شده در بالا را می‌توان مطابق شکل زیر برای دیواری یک بعدی با انتقال حرارت گذرا که با گره‌های ۰و۱و۲و۳ نشان داده شده است مطابقت داد.

گره ۱:

$ρ c \frac{\partial T}{\partial t} A Δ x = K A \frac{T_{0} - T_{1}}{Δ x} + K A \frac{T_{2} - T_{1}}{Δ x} = {ρ c}_{} A Δ x \frac{T_{1}^{n + 1} - T_{1}^{n}}{Δ t}$

$T_{1}^{n + 1} - T_{1}^{n} = \frac{α Δ t}{{Δ x}^{2}} {(T}_{0} + T_{2} - {2 T}_{1})^{n + 1} + q^{‴} \frac{Δ t}{{ρ c}_{}}$

$T_{1}^{n + 1} - T_{1}^{n} = f_{o} {(T}_{0} + T_{2} - {2 T}_{1})^{n + 1} \Rightarrow T_{1}^{n + 1} (1 + {2 f}_{o}) - f_{o} {(T}_{0}^{n + 1} + T_{2}^{n + 1}) = T_{1}^{n}$

گره ۲:

$T_{2}^{n + 1} (1 + {2 f}_{o}) - f_{o} {(T}_{1}^{n + 1} - T_{3}^{n + 1})$

گره ۳:

$ρ c \frac{Δ x}{2} A \frac{T_{3}^{n + 1} - T_{3}^{n}}{Δ t} = K A \frac{T_{2}^{n + 1} - T_{3}^{n + 1}}{Δ x} + {h A (T}_{\infty} - T_{3}^{n + 1}) \Rightarrow$

$T_{3}^{n + 1} (1 + {2 f}_{o} + {2 B i f}_{o}) - {2 f}_{o} T_{2}^{n + 1} = {2 B i f}_{o} T_{\infty}$

انتقال حرارت دو بعدی

مثال6

یک فین الومینیومی باk=180 و سطح مقطع مثلثی که طول آن L=5 cm و مقطع آن b=1 cm و عرض خیلی زیاد را در نظر بگیرید. مبنای فین با دمای $T_{0} = 200^{\circ} c$ درجهٔ سانتیگراد است. فین دمای خود را در معرض دمای محیط ۲۵ درجهٔ سانتیگراد با $h = 15 w / m^{2}^{\circ} k$ از دست می‌دهد با استفاده از روش تفاضل محدود دمای تمام نقاط را بیابید.

حل:
برای نوشتن معادلات finite diffrence برای گره‌های مختلف ازنوشتن معادلهٔ انرژی برای المان حجم استفاده می‌کنیم.

ابتدا $Δ x$ را محاسبه می‌کنیم:

$Δ x = \frac{L}{M - 1} = 0.01$

پس از ساده سازی معادله خواهیم داشت:

$\sum \overset{\cdot}{Q} = 0 \to K A_{l e f t} \frac{T_{m - 1} - T_{m}}{Δ x} + K A_{r i g h t} \frac{T_{m - 1} - T_{m}}{Δ x} + h A_{c o n v} (T_{\infty} - T_{m}) = 0$

$2 k w [L - (m - \frac{1}{2}) Δ x] \tan θ \frac{T_{m - 1} - T_{m}}{Δ x} + 2 k w [L - (m + \frac{1}{2}) Δ x] \tan θ \frac{T_{m - 1} - T_{m}}{Δ x} + h \frac{2 w Δ x}{\cos θ} (T_{\infty} - T_{m}) = 0$

$\to [1 - (m - \frac{1}{2}) \frac{Δ x}{L}] (T_{m - 1} - T_{m}) + [1 - (m + \frac{1}{2}) \frac{Δ x}{L}] (T_{m - 1} - T_{m}) + h \frac{Δ x^{2}}{k l \sin θ} (T_{\infty} - T_{m}) = 0$

$\tan θ = \frac{^{b}_{2}}{L} = 0.1 \to θ = 5.71$

$\begin{matrix} (5.5 - m) T_{m - 1} - (10.00838 - m) T_{m} + (4.5 - m) T_{m + 1} = - 0.209 \\ m = 1 \to - 8.00838 T_{1} + 3.5 T_{2} = - 900.209 \\ m = 2 \to 3.5 T_{1} - 6.00838 T_{2} + 2.5 T_{3} = - 0.209 \\ m = 3 \to 2.5 T_{2} - 4.00838 T_{3} + 1.5 T_{4} = - .209 \\ m = 4 \to 1.5 T_{3} - 2.00838 T_{4} + 0.5 T_{5} = - 0.209 \\ \to K A_{l e f t} \frac{T_{4} - T_{5}}{Δ x} + h A_{c o n v} (T_{\infty} - T_{5}) = 0 \to T_{4} - 1.00838 T_{5} = - .0209 \end{matrix}$

باحل چند معادله و چند مجهول دمای نقاط مختلف به دست می‌آید.

$T_{1} = 198.6 & T_{2} = 197.1 & T_{3} = 195.7 & T_{4} = 194.3 & T_{5} = 192.9$