مدلسازی انتقال حرارت پایا با تکنیک تفاضلات محدود در متلب + ویدیو

مقدمه

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

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

پیش‌نیازهای آموزش بررسی مدلسازی انتقال حرارت پایا

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

دقیق شدن در بررسی مدلسازی انتقال حرارت پایا: حل دقیق یا حل تقریبی؟

همان طور که می‌دانید برای بررسی انتقال حرارت پایا در اجسام می‌توان به دو صورت کلی عمل کرد:

• حل دقیق
• حل تقریبی

حل دقیق انتقال حرارت پایا در متلب

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

حل تقریبی انتقال حرارت پایا در متلب

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

یکی از راه‌های حل عددی معادلات در متلب استفاده از دستورات حل عددی مانند ode45 است. اما روشی که در این بخش قصد بررسی آن را داریم حل عددی به وسیله‌ی تفاضلات محدود است. برای فهم بهتر مسئله و کاربردی بودن آن، به بررسی فین‌ها یا پره‌ها می‌پردازیم.

استخراج معادلات و گسسته‌سازی در مدلسازی انتقال حرارت پایا با متلب

مدلسازی انتقال حرارت پایا را با نوشتن موازنه‌ی انرژی در یک بعد آغاز می‌کنیم. این‌کار در واقع برابری مجموع شارهای گرمایی وارده به جسم است. معادله‌ٔ موازنهٔ انرژی را برای یک پره‌ٔ یک بعدی با دمای پایه‌ی ثابت و انتقال حرارت جابه‌جایی با محیط، به شکل زیر در می‌آید:

چون در اینجا قصد داریم به کدنویسی انتقال حرارت پایا در متلب بپردازیم، وارد جزئیات نحوه‌ی نوشتن موازنه برای هر المان نمی‌شویم. برای خواندن مراحل به دست آوردن این معادله، می‌توانید به کتاب انتقال حرارت هولمن (Holman JP. Heat Transfer. Tata McGraw-Hill Education; 2002.) مراجعه کنید.

حالا باید از تکنیک تفاضلات محدود کمک گرفته و معادله مورد نظر را گسسته‌سازی کنیم. کتاب «کاربرد روش‌های عددی در مهندسی شیمی» نوشته‌ی دکتر شهره فاطمی مرجعی است که به خوبی شیوه‌ی گسسته‌سازی معادلات به کمک تفاضلات محدود را برای دانشجویان مهندسی شیمی شرح داده است. معادله گسسته شده در نهایت به این شکل خواهد بود.

در این جا تفاوت معادلهٔ اول که از حل تحلیلی به دست آمده با معادله دوم روشن می‌شود. در حل تحلیلی المان‌ها ضخامت‌های خیلی کم نزدیک به صفر دارند که آن ها را با dx نشان می‌دهند. ولی در حل عددی این المان‌ها دارای ضخامت‌های مشخصی هستند که معمولا برای سادگی مسئله ضخامت‌های برابر در نظر گرفته می­‌شود.

در این معادله، \( \Large T_{n} \) دمای المان \( \Large n \)اُم است. فرض می‌شود که دمای المان برابر با دمای نقطه‌ای در مرکز المان است. نقطه گذاری‌های ما معمولاً از مرزها شروع می‌­شوند. در ادامه با دو نوع نقطه مواجه می‌شویم:

• یکی نقاط داخلی که معادلات یکسانی دارند و فقط زیروندهای (اندیس) آن‌ها تغییر می‌­کند.
• دیگری نقاط مرزی که معادله‌ حاکم بر آن از شرایط مرزی مسئله حاصل می‌شود.

بنابراین تعداد معادلات برابر با تعداد نقاط مسئله است:

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

الگوریتم حل معادلات انتقال حرارت پایا

در بخش قبل، معادلات حاکم بر انتقال حرارت پایا را استخراج کرده و به حالت گسسته تبدیل کردیم. حال باید سعی کنیم تا این معادلات را حل کنیم. برای فهمیدن الگوریتم حل، از اولین نقطه مجهول یعنی نقطه دوم شروع می‌کنیم. در این معادله قصد ما به دست آوردن \( \Large T_2 \) یا دمای نقطه‌ی دوم است. پس آن را جدا کرده و به سمت دیگر معادله می‌بریم. (اگر نخواهیم یا نتوانیم چنین کاری کنیم، باید الگوریتم دیگری را پیش بگیریم که در قسمت آخر مجموعهٔ انتقال حرارت پایا به این دست موارد خواهیم پرداخت.)

در معادله‌ی بالا، مشخصات محیط و قطعه یعنی \( \Large P , h \) و \( \Large T_{\infty} \) مشخص و ثابت است. دمای نقطه اول \( \Large T_1 \) هم که از شرایط مرزی داریم. پس تنها مجهول باقی مانده دمای نقطه بعدی یعنی نقطه سوم \( \Large T_3 \) است. شاید تصور کنید که بتوان دمای نقطه‌ی سوم نیز از معادله‌ی بعدی به دست آورد. اما در معادله‌ی بعدی نیز دمای نقطه‌ی چهارم مجهول است. بنابراین این معادلات مستقلاً قابل حل نیستند و باید آن‌ها را به صورت هم‌زمان حل کرد.

حل معادلات دیفرانسیل حاکم بر انتقال حرارت پایا

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

1. از معادله‌ی اول شروع کرده و مقادیر حدس اولیه را در سمت راست معادله جایگذاری می‌کنیم. با این کار مقداری برای \( \Large T_2 \) به دست می‌آید.
2. سپس سراغ معادله‌ی دوم می‌رویم. در آن معادله، برای دمای نقطه‌ی چهارم از حدس اولیه و برای دمای نقطه‌ی دوم از مقدار جدید آن که در معادله‌ی قبل به دست آمده استفاده می‌کنیم. زیرا مقدار جدید \( \Large T_2 \) نسبت به مقدار قبلی دقیق‌تر است.
3. این روند را تا معادله‌ی آخر ادامه می‌دهیم.
4. حال از اعداد جدید به دست آمده، به عنوان حدس اولیه برای حل دوباره این معادلات از اول تا آخر استفاده می‌کنیم و این کار را تا جایی ادامه می دهیم که مقادیر به دست آمده با مقادیر قبلی خود اختلاف کم و قابل قبولی داشته باشند.

این الگوریتم که برای حل این دستگاه معادلات جبری چند معادله چند مجهول گفته شد، به روش حل عددی گاوس – سایدل مشهور است. در مجموعه‌ٔ مدلسازی انتقال حرارت پایا در متلب نحوه‌ٔ حل دستگاه معادلات انتقال حرارت پایا به کمک روش گاوس – سایدل، به همراه حل مثال بررسی شده است.

انتقال حرارت پایا در دو بُعد و مدلسازی آن در متلب

بعد از بررسی مسئله‌ی یک بعدی، به این می‌پردازیم که چه طور معادلات انتقال حرارت پایا را به هندسه‌ی دوبعدی تعمیم دهیم. با انتخاب المان شروع می‌کنیم. ساده‌ترین المان برای صفحه‌ی دو بعدی، مربع است. قدم بعدی نام‌گذاری نقطه‌ها است که می‌توانیم مثل کتاب‌های انتقال حرارت از یک تا n  نام‌گذاری کنیم. این کار ما را هم در کد نویسی و هم در بررسی جواب ما را دچار مشکل می‌کند. پس بهتر است مثل یک ماتریس نام‌گذاری کنیم که از \( \Large T_{1,1} \) شروع می شود و تا \( \Large T_{m,n} \) ادامه دارد.

سپس سراغ به دست آوردن معادله می‌رویم. یک المان میانی و المان‌هایی که روی آن تاثیر می‌گذارد را به عنوان نمونه برداشته و موازنه انرژی را با تکنیک تفاضلات محدود برای آن المان می‌نویسیم. با نوشتن موازنه، معادله روبرو حاصل می‌شود:

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

که R در معادله بالا نشان‌دهنده‌ی مقاومت انتقال حرارت بین دو المان مورد نظر است که برابر \( \Large \frac{\Delta x}{kA} \) یا \( \Large \frac{\Delta y}{kA} \) است. معادله بالا به شکل خلاصه‌تر به صورت زیر خواهد بود:

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

شرایط مرزی برای مدلسازی در متلب

در مسائل مدلسازی انتقال حرارت پایا معادلات دیفرانسیل از مرتبه‌ی دو هستند. یعنی مشتق دوم دما در معادلات حضور دارد. بنابراین برای حل در هر بُعد به دو شرط مرزی نیاز داریم. یک شرط مرزی در ابتدای بازه (مثلاً پایه‌ی پره) و شرط دوم در انتهای بازه (مثلاً انتهای پره). در بخش‌های قبلی فرض بر این بود که شرط مرزی در ابتدا و انتهای پره، از نوع اول است. به این معنی که مقدار دما در ابتدا و انتهای پره، مشخص است.

حال می‌خواهیم شرایط مرزی ممکن دیگر مانند عایق بودن یا وجود تشعشع را بررسی کنیم. در حالتی که انتهای میله عایق باشد، در معادله‌ٔ موازنه‌ٔ انرژی برای نقطه‌ٔ آخر تغییر کوچکی ایجاد می‌گردد. به این شکل که عبارت مربوط به تبادل حرارت با محیط از معادله حذف می‌شود و معادله حاصل به شکل زیر درمی‌آید:

در صورتی که انتهای میله تحت تشعشع حرارتی قرار داشته باشد، باید شار وارد شده به المان (نقطه) انتهایی را برابر تشعشع وارد شده قرار داد. رابطه‌ٔ مربوط به تشعشع به شکل زیر است:

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

دقت حل مدلسازی

برای تعیین دقت حل دو عامل موثر هستند:

• اولی تعداد المان‌ها است که بیشتر شدن آن ها باعث کاهش اندازه‌ی هر المان شده و در نتیجه به نمایش بهتر توزیع دما کمک می­‌کند.
• دیگری خطای محاسبه سعی و خطا که هر چه کمتر شدن آن اعداد ما را به حالت دقیق نزدیک‌تر می­‌کند.

در عکس زیر تأثیر تعداد المان‌ها بر دقت حل معادلات، نمایش داده شده است. همان طور که در تصویر مشخص است، افزایش تعداد المان‌ها به پیوستگی توزیع دما کمک می‌کند.

تعداد نقاط برابر ۱۰

تعداد نقاط برابر ۵۰

در شکل زیر نیز تأثیر کاهش مقدار خطا را مشاهده می‌کنید که باعث نزدیک شدن توزیع دما با حالت دقیق می‌­شود.

توزیع دما در میله با مقدار خطای ۱

توزیع دما در میله با مقدار خطای ۰.۱