skip to Main Content

حل معادلات دیفرانسیل مرتبه اول  ODE  به کمک نرم افزار متلب

نرم افزار متلب یکی از پر کاربرد ترین نرم افزارهای محاسباتی برای دانشجویان رشته های مهندسی می باشد. سادگی محیط و وجود توابع آماده زیاد باعث شده است دانشجویان برای انجام تمرین ها و پروژه های درسی خود از این نرم افزار کمک بگیرند. در این پست حل معادلات دیفرانسیل مرتبه اول ODE به کمک نرم افزار متلب را در قالب یک مثال بررسی می کنیم. عموما فهم درس انتقال حرارت نسبت به سایر دروس مهندسی شیمی راحت تر می باشد و به همین دلیل مثال های این بخش را از درس حرارت انتخاب کرده‌ایم.

مثال اول:

توزیع دما را در حالت پایدار در یک دیواره ی آلومینیومی با ضخامت ۱ متر و ارتفاع زیاد که تحت تاثیر شار حرارتی ثابت به میزان \frac{W}{m^{2}} ۲۰۰۰ قرار دارد را بدست آورید. به دلیل اعمال شار ثابت، دمای سطح دیواره در ۴۰۰K ثابت می شود. همچنین ضریب هدایت حرارتی دیواره نیز  \frac{W}{mK} ۲۲۰ می باشد.

تصویر مثال اول حل معادله مرتبه اول عددی به کمک متلب

حل:

سوال در ساده ترین حالت ممکن قرار دارد. از قانون فوریه می دانیم که  q=-k\frac{\text{d}T}{\text{d}x} . خواسته ی مسئله بدست آوردن T(x)  می باشد پس باید معادله مرتبه اول \frac{q}{-k}=\frac{\text{d}T}{\text{d}x} را حل کنیم. با نگاه به معادله و یک انتگرال گیری ساده می دانیم جواب به صورت یک معادله خط خواهد بود اما می خواهیم روند حل را به کمک نرم افزار و به صورت عددی بررسی کنیم.

برای حل عددی معادلات دیفرانسیل روش های مختلفی وجود دارد. ما در اینجا از روش اولر برای حل معادلات استفاده می‌کنیم. دقت کنید پروژه ها و تمرین های درسی عموما مسائل پیچیده و سنگینی نیستند اما در شبیه سازی های واقعی و پروژه های مربوط به پایان نامه ها، بسته به نوع و شرایط خاص مسئله، ملاحظاتی پیرامون روش حل مطرح می شود و پارامتر هایی نظیر همگرایی و پیوستگی باید مورد بررسی قرار گیرند. این مباحث در دروسی مثل محاسبات عددی، CFD و ریاضیات پیشرفته در دوره ی کارشناسی ارشد مطرح می شوند و فعلا هدف ما تشریح این مباحث نمی باشد. بر گردیم به حل مسئله.

روش اولر برای حل معادلات دیفرانسیل:

در این روش پیشنهاد می شود به جای جمله مشتق مرتبه اولی یعنی  \frac{\text{d}T}{\text{d}x} جمله \frac{T_{2}-T_{1}}{\triangle x} را قرار دهیم. با این کار معادله دیفرانسیل مرتبه اول ODE به یک معادله جبری ساده تبدیل می شود. بعد از جا گذاری عبارت مشتق ، معادله به شکل زیر در می آید:

(T_2-T_1)/∆x=q/(-k)

و در نهایت با جا به جایی و مرتب کردن، معادله نهایی ما به صورت زیر در می آید:

T_2=T_1+∆x q/(-k)

اما T1 و T2وچیست؟ برای بهتر فهمیدن این پارامتر ها مجددا به شکل سوال دقت کنید. فرقی نمی کند در این مسئله شما دیواره را عمودی، افقی، مایل و … فرض کنید، مهم انتخاب دستگاه مختصات می باشد که مبنای حل ما است. در اینجا ما از مختصات کارتزین استفاده می کنیم. با توجه به شکل ما جهت مثبت محور ایکس را از سطح دیواره به داخل دیواره فرض می کنیم.

در حل عددی، ما مسئله را باید گسسته سازی کنیم. به زبان ساده تر ما طول یک متر را به صورت مجموعه تعدادی نقاط ریز با فاصله x∆ در نظر می گیریم و با این کار یک طول پیوسته را به مجموعه ای از نقاط گسسته تبدیل می کنیم. مثلا در اینجا طول یک متر را معادل ۱۰۰۰ نقطه با فاصله ۰٫۰۰۱ از یکدیگر درنظر می گیریم. همان طور که در شکل مشخص است T1 نقطه ی اول شروع ماست و بر روی سطح دیواره قرار دارد.پس نقطه T2 به فاصله x∆ از T1 قرار می گیرید و به همین ترتیب بقیه نقاط تعیین می شوند. الگوریتم حل در روش اولر به این صورت می باشد که برای T1مقدار اولیه تعیین می کنیم و با استفاده از معادله جبری که بدست آوردیم T2,T3,T4 ,…. را بدست می آوردیم.

روند حل معادله دیفرانسیل به روش اولر

حال به سراغ نرم افزار می رویم. از منوی HOME روی گزینه new script کلیک می کنیم. بسته به سلیقه هر کس روش خاصی برای برنامه نویسی دارد. در اینجا ما فعلا کد ها را به صورت تابع نمی نویسیم. برای حل عمق دیواره را متشکل از ۱۰۰۰ نقطه به فاصله ۰٫۰۰۱ در نظر می گیریم.

ابتدا شرایط اولیه و داده های مسئله را وارد می کنیم:

1
2
3
4
5
6
k=220; %w/mk
q=2000; %w/m^2
dx= 0.001; %meter
T0=400; %initial T ,kelvin
x=0:0.001:1; %delta X
T(1)=T0;

حال الگوریتم حل را وارد می کنیم.توجه کنید که شمارنده حلقه for باید برابر با تعداد نقاط دامنه ما یعنی ۱۰۰۰ نقطه باشد. برای بررسی نتایج پروفایل دما را بر حسب عمق دیواره رسم می کنیم:

1
2
3
4
5
6
7
for i=1:1000

T(i+1)= T(i) + dx*(-q/k); %kelvin

end

plot(x,T);

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
clear
clc
k=220; %w/mk
q=2000; %w/m^2
dx= 0.001; %meter
T0=400; %kelvin
x=0:0.001:1; %delta X
T(1)=T0;

for i=1:1000
T(i+1)= T(i) + dx*(-q/k); %kelvin
end

plot(x,T);

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

 

در این نمودار محور افقی عمق دیواره و محور عمودی دما را نشان می دهد که همان چیزی بود که انتظار داشتیم یعنی معادله یک خط صاف. حال کمی مسئله را واقعی تر کنیم.

مثال دوم:

مسئله قبل را با فرض متغیر بودن ضریب هدایت حرارتی حل کنید. در این حالت ضریب هدایت حرارتی از معادله  k=k_{0}(1+\beta T) بدست می آید.ثوابت این معادله نیز به میزان داده شده در زیر است.

حل:

کافی است در معادله ی   T_{2}=T_{1}-\triangle x \frac{q}{k} به جای k ،مقدار داده شده ی جدید را قرار دهیم. با جا گذاری معادله ی ما به صورت T_{2}=T_{1}-\triangle x \frac{q}{k_{0}(1+\beta T)} در میاید. سوالی که پیش می آید این است که مقدار T در جمله سمت راست را کدام دما بگیریم؟ T2 یا T1 ؟

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

در کد نوشته شده نیز کافی است جمله ضریب هدایتی و هم چنین معادله جبری را به صورت زیر به کد حل قسمت قبل اضافه کنیم و الگوریتم حل را نیز کمی تغییر دهیم:

1
2
3
4
5
6
7
8
9
10
11
beta=1.315e-4;

k0=228;

for i=1:1000

k(i)=k0*(1+beta*T(i)); %w/mk

T(i+1)= T(i) + dx*(-q/k(i)); %kelvin

end

حل دقیق این مسئله در کتاب روش های عددی در پدیده های انتقال دکتر نیک آذر، در فصل سوم مثال ۳٫۵آمده است و به صورت زیر می باشد.( T0دمای سطح صفحه می باشد)

برای بررسی نتایج حاصل از کد بالا، جواب دقیق و جواب عددی بدست آمده را در یک نمودار رسم می کنیم.

پس در نهایت کد ما به صورت زیر در می آید:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
clear
clc
beta=1.315e-4;
k0=228;
q=2000; %w/m^2
dx= 0.001; %meter
T0=400; %kelvin
x=0:0.001:1; %delta X
T(1)=T0;

for i=1:1000
k(i)=k0*(1+beta*T(i)); %w/mk
T(i+1)= T(i) + dx*(-q/k(i)); %kelvin
end

Texact=-1/beta + ( (T0+1./beta).^2 -q.*2.*x./(beta.*k0) ).^0.5;
plot(x,T,x,Texact);
error=(T-Texact)./Texact;

بحث و نتیجه گیری

تا اینجا حل عددی معادلات مرتبه اول ODE به کمک نرم افزار متلب را مشاهده کردیم. ولی همیشه باید بعد از حل عددی نسبت به صحت جواب های بدست آمده اطمینان حاصل کنیم. مشاهده می شود که در این مثال خاص هر دو نمودار یکی شده اند و روی هم افتاده اند و هم چنین خطای نسبی روش ما همواره کمتر از ۰٫۰۰۱ درصد می باشد. حالت پیش آمده به این دلیل می باشد که خطای روش اولر استفاده شده در اینجا از مرتبه یک می باشد و همچنین معادله نیز از مرتبه یک است، پس خطای ما بسیار بسیار ناچیز خواهد بود. نتیجه دیگری که از این نمودار می گیریم این است که برای آلومینیوم، در این شرایط دمایی، میزان تغییرات ضریب هدایت حرارتی تاثیر گذار نمی باشد، چرا که در حالت قبل نیز تقریبا به همین دما در فاصله ۱ متری رسیدیم.

دیدگاه بگذارید

avatar
  خبر دار شدن  
از
VimeoTelegramInstagramYoutube
Back To Top