Finite difference တွက်နည်းအကြောင်းရေးလာတာ အပိုင်း (၃) ကို ရောက်လာပါပြီ။ (အပိုင်း ၁ နဲ့ ၂ အပြီးမှာ မအားလပ်တာနဲ့ ဆက်မရေးဖြစ်တော့တာ တော်တော်တောင်ကြာသွား😁)
အပိုင်း (၁) မှာ ရူပဗေဒဆိုင်ရာ သဘောတရားအနည်းငယ် (အပူကူးပြောင်းပုံ ၃မျိုး အကြောင်း) ကို ယေဘုယျသဘောဆွေးနွေးဖြစ်ခဲ့ပါတယ်။ အပိုင်း (၂) မှာတော့ finite difference method ရဲ့ အရေးကြီး ညီမျှခြင်းတွေဖြစ်တဲ့ forward ၊ backward နဲ့ central difference ပုံသေနည်းတွေအကြောင်း တင်ပြခဲ့ပါတယ်။ သင်္ချာရှုထောင့်ကဲနေတော့ ကောင်းကောင်းနားလည်ဖို့ အတန်အသင့် အားစိုက်ထုတ်ဖို့ လိုကောင်းလိုပါမယ်။ အခု အပိုင်း (၃) မှာတော့ အပူလျှောက်ကူးခြင်း (heat conduction) ပုစ္ဆာကို finite difference နည်းလမ်းအသုံးပြုပြီး ဖြေရှင်းပါမယ်။ လောလောဆယ် excel ရဲ့ iterative function ကို သုံးပြီးတော့ပဲ ဖြေရှင်းပြပါမယ်။ အပိုင်း (၄) နဲ့ (၅) မှာတော့ matrix inversion တို့ Jacobi စတဲ့ နည်းလမ်းတွေကို အသုံးပြုပြီး ဖြေရှင်းပါမယ်။ ဒါတွေက numerical လေ့လာသူတို့အတွက် သိထားသင့်တဲ့ တွက်နည်းတွေမို့ပါ။ အဖြေက အနီးစပ်ဆုံး တူနေရမှာ ဖြစ်သလို ၊ မိမိအဖြေရဲ့ တိကျမှုအပေါ်မှာ ဝေဖန်ဆန်းစစ်ကြည့်နိုင်ဖို့လည်း အရေးကြီးပါတယ်။
Heat Equation
အရင်ဆုံး Heat equation (အပူညီမျှခြင်း) ခေါ် Heat diffusion equation လေးကို အရင်ဆုံး ကြည့်လိုက်ရအောင်။ wikipedia မှာ သွားရှာကြည့်မယ်ဆိုရင် အောက်ပါညီမျှခြင်းပုံစံမျိုးကို တွေ့ရမှာပါ။
where \alpha = diffusivity of a medium (ကြားခံနယ်ရဲ့ ပျံ့နှံ့ကိန်းသေ – သူက ကြားခံနယ်ရဲ့ ဒြပ်ပစ္စည်းအမျိုးအစားပေါ်မူတည်ပါတယ်)။
အထက်ပါ ညီမျှခြင်းထဲက T ကတော့ temperature ကို ကိုယ်စားပြုပါတယ်။ အခုပြထားတာက three dimension ညီမျှခြင်းပါ။ နောက်တစ်ချက်တွေ့ရမှာက temperature ဟာ နေရာပေါ်ရော အချိန်ပေါ်ရောပါလိုက်ပြီး ပြောင်းလဲပါတယ်။ သင်္ချာနည်းနဲ့ရေးရင် T = f(x, y, z, t) ပေါ့။ အဲ့ဒီထဲက x, y, z ဆိုတာ တည်နေရာ (spatial variables) ဖြစ်ပြီး t ကတော့ time ပါ။
Finite difference ကို စပြီးလေ့လာကြတဲ့အခါ ခုနက 3D ညီမျှခြင်းမျိုးကို ယူဆချက် (assumptions) တွေ ထပ်ပေးပြီးတော့ ရိုးရိုးရှင်းရှင်းပဲ လေ့လာလေ့ရှိကြတယ်။ ယူဆချက်တစ်ခုကတော့ steady state condition ပါ (တည်ငြိမ် သို့ တည်မြဲအခြေအနေလို့ ဘာသာပြန်ရပါမယ်)။ သဘောကတော့ steady state condition မှာ temperature ဟာ အချိန်နဲ့ လိုက်ပြီး မပြောင်းတော့ပါ။ ဒါ့ကြောင့် အထက်ပါ ညီမျှခြင်း (ဘယ်ဖက်ခြမ်း) ထဲက \frac{\partial T}{\partial t} = 0 ြဖစ်သွားပါမယ်။ Derivative of t ကို သုညပေးလိုက်တာပါ။
ညီမျှခြင်း ညာဖက်ခြမ်းက spatial ပြောင်းလဲခြင်း terms တွေက နေရာကို လိုက်ပြီး အပူချိန်တွေ ပြောင်းလဲတဲ့ သဘောသဘာဝကို ပြပါတယ်။ ဥပမာ – အမှတ်တစ်နေရာမှာ ရှိတဲ့ အပူချိန်ဟာ သူ့ဘေးပတ်ပတ်လည်က အမှတ်တွေရဲ့ အပူချိန်တွေပေါ်မူတည်ပြီး ပြောင်းလဲပါတယ်။ အဲ့ဒီမှာ အရေးကြီးတဲ့ assumption နောက်တစ်ခုကတော့ One dimension (တစ်ဖက်မြင်) အနေနဲ့ စဉ်းစားလိုက်တာပါပဲ။ တကယ်က 3D ညီမျှခြင်းကို သုံးရမှာ။ ဒါပေမယ့် cross-section area ကလည်း တသမတ်တည်းရှိနေမယ်၊ နောက်တခါ အလျားလိုမျိုး (အတိုင်းအတာ dimension တစ်ခုက) ကျန်တဲ့ကောင်တွေထက် သိသိသာသာကြီးနေမယ်ဆိုရင် 1D ညီမျှခြင်းကို အသုံးပြုလို့ရပါတယ်။ ဥပမာ မြင်သာအောင် rod (သို့မဟုတ်) bar တန်းလိုမျိုး geometry မျိုး (ပုံ ၁) ကို မြင်ယောင်ကြည့်ပါ။
ပုံထဲမှာ အပူချိန်က x direction မှာပဲ ပြောင်းလဲပါမယ်။ ကျန်တဲ့ y နဲ့ z direction မှာတော့ ပြောင်းလဲမှုမရှိပါ (တစ်နည်းပြောရရင် y နဲ့ z dimension က အပူချိန်တွေဟာ x က အပူချိန်အတိုင်းပဲ ရှိပါမယ်)။ ဒါကို one dimension (တစ်ဖက်မြင်) ပုစ္ဆာလို့ခေါ်ပါတယ်။ သင်္ချာလိုရေးရင် temperature T = f(x) ပါ။ y နဲ့ z အတွက် variables ေတွ ထည့်စဉ်းစားစရာ မလိုတော့ ခေါင်းရှုပ်စရာမလိုတော့ဘူးပေါ့။ အဲ့ဒီအခါ Heat equation ဟာလည်း အင်မတန်မှ ရိုးရှင်းလှတဲ့ ordinary differential equation လေးဖြစ်သွားပါတယ်။ အောက်ပါအတိုင်း ရေးလို့ရပါပြီ။
\frac{d^2T}{dx^2} = 0
———- Eq. (1)
မှတ်ချက်။ \alpha က အပေါင်းတန်ဖိုးရှိကိန်းသေဖြစ်တာကြောင့် ညာဖက်ကိုသွားစားတဲ့အခါ သုညပဲ ပြန်ရပါမယ် (အဆင့်ကျော်ပြီး ရေးထားတာပါ 😁)
ယူဆချက်တွေနဲ့တကွ Simplification တွေ လုပ်အပြီးမှာတော့ ပုစ္ဆာက တော်တော်လေးကို လွယ်ကူသွားတာကို တွေ့ရပါမယ်။ ပုံမှန်သင်ခန်းစာတွေမှာတော့ 2D plate အတွက် Laplace ညီမျှခြင်း ၊ ဒါမှမဟုတ်လည်း 1D steady state heat equation ကိုမှ heat generation (အပူထုတ်လွှတ်မှု) ကိုပါထည့်ပြီးစဉ်းစားလေ့ရှိပါတယ်။ အခုဆောင်းပါးမှာတော့ မိတ်ဆက်အဆင့်သာမို့ 1D (အလွယ်ဆုံး အရှင်းဆုံး) ညီမျှခြင်းကိုသာ ဖြေရှင်းသွားပါမယ်။
ဒါနဲ့ စကားမစပ် ဒီဆောင်းပါးထဲမှာ heat equation ဆိုပြီး wiki မှာတွေ့တဲ့ ညီမျှခြင်းကိုပဲ အလွယ်တကူကောက်ထည့်ထားတာပါ။ ဆောင်းပါးအရမ်းရှည်သွားမှာစိုးလို့ တွက်ထုတ်ပုံအစုံအလင်ကို ချမပြတော့ပါဘူး။ သဘောတရားအနည်းငယ်ကိုသာ ရှင်းပြပါ့မယ် (သိချင်ရင်တော့ စာရေးသူကိုသာ တိုက်ရိုက်ဆက်သွယ်မေးမြန်းနိုင်ပါတယ်) ။ အဆိုပါညီမျှခြင်းလေးဟာ control volume (differential volume သေးသေးလေးကို ဆိုလို) တစ်ခုရဲ့ energy balance ကို အခြေခံပြီးတော့မှ ပေါ်ပေါက်လာတာပါ။ ဒီနေရာမှာ energy balance ဆိုတာ control volume ထဲကိုဝင်လာတဲ့ စွမ်းအင်စုစုပေါင်းဟာ control volume ထဲကနေ ပြန်ထွက်သွားမယ့် စွမ်းအင်စုစုပေါင်းနဲ့ ညီပါတယ်။ (ပုံ ၂ မှာကြည့်ပါ)
ဒီနေရာမှာ စွမ်းအင်ဆိုတာက အပူစီးဝင်မှု (heat flow in) နဲ့ အပူစီးထွက်မှု (heat flow out) တို့ကို ခေါ်တာပါ။ သူတို့နှစ်ခုဆက်သွယ်ချက်ကို Taylor series အသုံးပြုပြီး ဖြန့်ရေးလို့ရပါတယ်။ တစ်ခါ အပူစီးဆင်းမှု heat flow ကို Fourier law အရ အပူချိန်ပြောင်းလဲခြင်း (temperature gradient) နဲ့ ပြန်ဆက်စပ်လို့ရပါတယ်။ ပြောရရင် အပူတစ်နေရာကနေ တစ်နေရာကို စီးတယ်ဆိုတာ တကယ်တော့ ကြားခံနယ်ရဲ့ အပူချိန်တွေ ပြောင်းလဲခြင်းပါ (ပုစ္ဆာတွက်ကြည့်တဲ့အခါ ပိုနားလည်သွားပါလိမ့်မယ်)။ နောက်ထပ် concept တွေဖြစ်တဲ့ တစ်ယူနစ်ဒြပ်ထုမှာ ရှိမယ့် internal energy ပြောင်းလဲခြင်း၊ တစ်ယူနစ်ထုထည်မှာရှိမယ့် အပူလွှတ်ခြင်း (heat generation per unit volume) စသဖြင့် ထည့်ပြီး စဉ်းစားလိုက်တဲ့အခါ လိုချင်တဲ့ heat equation ကို ရရှိမှာဖြစ်ပါတယ်။ စိတ်ဝင်စားသူတွေအတွက် ဆက်စပ်လေ့လာကြည့်လို့ရမယ့် လင့်ခ်လေးတွေကို ဆောင်းပါးအဆုံးမှာ ဖော်ပြပေးထားပါတယ်။ လေ့လာကြည့်ကြပါ။ နားမလည်တာ ဘာညာရှိရင်လည်း insight ကို (မြန်မာလို သို့ english လို) အချိန်မရွေးလာမေးလို့ ရပါတယ်။ ဒါကတော့ derivation အနှစ်ချုပ်ကို စာသားနဲ့ အတိုဆုံးချုပ်ပြီး ရှင်းပြလိုက်တာပဲ ဖြစ်ပါတယ်။
Boundary conditions
အထက်ပါ ညီမျှခြင်း (Eq. 1) ကို boundary value problem လို့ခေါ်တယ်။ ဆိုပါတော့ ပုံ ၁ ထဲက တုတ်ချောင်းရဲ့ ထိပ် ၂ ဖက်က အခြေအနေကို သိတယ် (အခြေအနေဆိုတာ temperature ကို ပြောတာပါ။)။ ကြားထဲက domain ရဲ့ temperature distribution တွေကို သိချင်တယ်။ ဘယ်လိုတွက်ကြမလဲပေါ့။ အောက်ဖော်ပြပါ ပုံမျိုး အကြမ်းလေးဆွဲပြီး စဉ်းစားပါမယ်။
ပုံ ၃ မှာပြထားတာက bar လေးတစ်ချောင်း။ သူ့ရဲ့ အလျားက L မီတာ ရှည်တယ်။ (နောက်မှ L ကို တန်ဖိုးတစ်ခုခုပေးတွက်ကြည့်တာပေါ့)။ အမှတ် A နဲ့ အမှတ် B မှာ ရှိတဲ့ temperature ကို သိတယ်။ အဲ့ဒီ A နဲ့ B အမှတ်က temperature condition ကို boundary condition ဆိုပြီးခေါ်တာပါ။ ဒါကို နောက်တစ်မျိုး Dirichlet boundary condition ဆိုပြီးတော့လည်း ခေါ်ကြပါတယ်။ အဲ့ဒီ Boundary condition တွေကတော့ အဖြေရဲ့ တန်ဖိုးကိုပါ ပေးပြီးသားမျိုးပေါ့။ (နောက်ထပ် Neumann condition ဆိုပြီးတော့ ရှိသေးတယ်။ ဒါကိုတော့ မပြောပြတော့ဘူးနော်။ သိချင်လို့ရှိရင် စာရေးသူကို ဆက်သွယ်မေးကြည့်လို့ရပါတယ်။ ခုထိတော့ ဘယ်သူမှ လာမမေးကြဖူးဘူး 😂 ။ )
ညီမျှခြင်းပုံစံနဲ့ ဖော်ပြရမယ်ဆိုရင်
အမှတ် A မှာ x တန်ဖိုးက သုညဆိုတော့ T(0) = T_0
အမှတ် B မှာ x = L ဆိုတော့ T(L) = T_e
ဒီနေရာမှာ T_0 နဲ့ T_e ဆိုတာက temperature တန်ဖိုးတစ်ခုခုဆိုပါတော့။
ဥပမာပုစ္ဆာတွက်ကြည့်ခြင်း
ကဲ ဥပမာအနေနဲ့ အခုပြောသွားတာတွေကို တန်ဖိုးတစ်ခုခု သတ်မှတ်ပြီး ဖြေရှင်းကြည့်ရအောင်။
L = length = 1 m,
K = thermal conductivity = 45 W/K/m
T0 = temperature at point A = 100 degC
Te = temperature at point B = 200 degC
ဖြေရှင်းရမယ့် ညီမျှခြင်းက \frac{d^2T}{dx^2} = 0 ပါ။ (K တန်ဖိုးတောင်မလိုတော့တာတွေ့ရပါမယ်)
ဒါကို finite difference ညီမျှခြင်း သုံးပြီးရှင်းပါမယ်။ အပိုင်း ၂ မှာရှင်းပြခဲ့တဲ့ second order derivative ပုံသေနည်းလေးကို ယူသုံးပါမယ်။ (မှတ်ချက်။ f နေရာမှာ T တွေ အစားထိုးပြီး ပြောင်းရေးထားပါတယ်။ temperature အတွက် function ပါ)
\frac{d^2T}{dx^2} \approx \frac{T(x +\Delta x) + T(x – \Delta x) – 2 T(x)}{\Delta x^2}
———- Eq. (2)
အထက်ပါ ညီမျှခြင်းကို differential equation ထဲကို အစားသွင်းလိုက်ရင်
\frac{T(x +\Delta x) + T(x – \Delta x) – 2 T(x)}{\Delta x^2} = 0
———- Eq. (3)
ရပါမယ် (ဒီ့ထက်ပိုရှင်းတဲ့ နည်းနဲ့ ရေးလို့ရပါသေးတယ်။ အောက်မှာ ထပ်ဖတ်ကြည့်ပါ)။
ပုံ ၃ ထဲက bar တန်းလေးကို အပိုင်း ၁၀ ပိုင်း ပိုင်းလိုက်မယ်။ တပိုင်းချင်းစီအတွက် nodal point လေးတွေရလာမယ်။ အဲ့ဒီ nodal point တွေပေါ်မှာပဲ တွက်ချက်မှုတွေ လုပ်ပါမယ်။ မိမိလိုအပ်ရင် လိုအပ်သလောက် nodal point များများလိုချင်ရင် အပိုင်းကို ၁၀ ပိုင်းထက်ပိုပြီး ပိုင်းပေးရပါမယ်။ အခုတော့ အပိုင်း ၁၀ ပိုင်းပိုင်းတွက်ရင် လုံလောက်ပါပြီ။ (ပုံမှာကြည့်ပါ)
ပုံထဲက \Delta x = \frac{L}{10} = 0.1 m ပါ။ (၁၀ပိုင်းပိုင်းဖို့ အလျားကို ၁၀ နဲ့စား)
T(x) ဆိုတာ မိမိစဉ်းစားချင်တဲ့အမှတ်။ ဒါကို လွယ်အောင် T(x) = T_i ဆိုပြီးရေးမယ်။ဘယ်ဖက်ကအမှတ် T(x - \Delta x) = T_{i-1} ၊ ညာဖက်ကအမှတ်ကို T(x + \Delta x) = T_{i+1}
အောက်ပါပုံမျိုးလေးနဲ့ ယေဘုယျပုံစံမျိုးပြလေ့ရှိတယ်
အပေါ်က ညီမျှခြင်း နံပါတ် (၃) ကို ဒီလိုပြန်ရေးလိုက်ပါ။
\frac{T_{i+1} + T_{i-1} – 2 T_i}{\Delta x^2} = 0
T_{i+1} + T_{i-1} – 2 T_i = 0
T_i = \frac{T_{i+1} + T_{i-1}}{2}
Domain ထဲက အမှတ်တိုင်း (nodal point) ရဲ့ temperature ကိုလိုချင်ရင် ဘယ်အမှတ်နဲ့ ညာအမှတ်က temperature ကို ပေါင်းပြီး ၂ နဲ့စားပေးရပါမယ်။ boundary အမှတ် (အမှတ် A) မှာရှိတဲ့ temperature က သိပြီးသား၊ သူ့ဘေးကအမှတ်ကို T2 ၊ T2 ရဲ့ ဘေးကအမှတ်ကို T3 စသဖြင့်ရေးလိုက်ရင်
T1 က သိပြီးသား (boundary condition ကနေ) ဟုတ်တယ်ဟုတ်။ ( T_1 = 100^{\circ}C )
T2 ကို လိုချင်ရင် သူ့ရဲ့ ဘယ်နဲ့ ညာကိုပေါင်း ၂ နဲ့စား ၊ ဆိုတော့ T_2 = \frac{T_1 + T_3}{2}
စသဖြင့် ညီမျှခြင်းတွေကို ရေးသွားရမှာပါ။ ဒါကို excel နဲ့ အလွယ်တကူ ဖြေရှင်းလို့ရပါတယ်။ မြန်အောင် excel ထဲမှာ တွက်ပြထားတာကို video အသေးစားလေးပြန်လုပ်ပေးလိုက်ပါတယ်။
ဗီဒီယိုထဲက excel file ကို download ဆွဲရန် နှိပ်ပါ
နောက်ဆုံးအဖြေက အောက်ပါအတိုင်း linearly varying ဖြစ်နေတဲ့ graph လေးကို ရပါမယ်။
သူ့ကို ပြန်ချိန်ကိုက်ကြည့်ချင်ရင် exact method နဲ့ integral ယူပြီး စစ်ကြည့်လို့ရပါတယ်။
\frac{d^2T}{dx^2} = 0
နှစ်ဖက်လုံးကို dx နဲ့လိုက်ပြီး integral ယူလိုက်ရင်
\frac{dT}{dx} = C
နှစ်ဖက်လုံးကို dx နဲ့လိုက်ပြီး integral ထပ်ယူလိုက်ရင်
T(x) = Cx + D
C နဲ့ D ကတော့ arbitrary constants တွေပါ။ သူတို့ကို boundary condition ကနေပြန်ရှာရပါမယ်။
At x = 0,
T(0) = C.(0) + D = 100
D = 100
At x = L = 1,
T(1) = C.(1) + D = 200
C = 100
ဒါဆိုရင်
T(x) = 100 x + 100
ပုံထဲမှာ ပြထားတဲ့ ညီမျှခြင်းလေးနဲ့ တူနေတာကို တွေ့နိုင်ပါတယ်။ ဒါကတော့ 1D steady state head conduction ညီမျှခြင်း (Dirichlet boundary condition) ကို ဖြေရှင်းပြထားတာပဲ ဖြစ်ပါတယ်။
အပိုင်း ၄ နဲ့ ၅ မှာလည်း ဒီပုစ္ဆာကိုပဲ နောက်ထပ်မတူတဲ့ တွက်နည်းတွေနဲ့ တွက်ကြည့်ပြီး ပိုမိုနားလည်အောင် လေ့လာသုံးသပ်ကြည့်ကြပါမယ်။ အားလုံး see you next time ပါ။
#yp
ဆက်စပ်လေ့လာကြည့်ရန်