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

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

شرایط CFL

همانطور که اشاره شد بررسی پایداری روش‌های مختلف عددی، یکی از مباحث بسیار مهم است که کاربرد زیادی در علوم توربوماشین و آیرودینامیک دارد و پایداری حل‌های عددی انجام شده در این علوم اهمیت بسیار زیادی دارد. بنابراین هدف اصلی این مطلب بررسی پایداری روش تفاضل محدود است. بر این اساس در قدم اول، شرایط Courant-Friedrichs-Levy را به منظور بررسی پایداری روش‌های تفاضل محدود برای معادلات هذلولوی مورد مطالعه قرار می‌دهیم. توجه کنید که این شرایط یعنی Courant-Friedrichs-Levy را با نماد اختصاری CFL نیز نمایش می‌دهند.

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

روش FTBS، گسسته سازی معادلات را به صورت پیشرو و گسسته سازی در مکان را به صورت پسرو انجام می‌دهد. برای آشنایی بیشتر و دقیق‌تر با شیوه گسسته سازی معادلات و آشنایی با روش FTBS به مطلب «روش تفاضل محدود (Finite Difference Method) — از صفر تا صد» مراجعه کنید.

البته حالت دیگری نیز برای ناپایدار شدن حل وجود دارد. در این حالت با کاهش سایز شبکه و ثابت ماندن گام زمانی، روش FTBS در نهایت ناپایدار می‌شود. بنابراین واضح است که انتخاب گام زمانی مستقل از سایز شبکه نیست و بنابراین به منظور حفظ پایداری حل، اندازه گام زمانی و سایز شبکه باید با یکدیگر به صورت هماهنگ تغییر کنند.

نکته دیگری که باید به آن اشاره کرد این است که در یک حل عددی می‌توان ناحیه‌ای را تعیین کرد که حل در یک نقطه مانند $$ (x_ i , t ^ n )$$ به مقادیر موجود در آن ناحیه بستگی دارد. شکل زیر «ناحیه وابستگی حل عددی» (Numerical Domain of Dependence) در روش FTBS نشان می‌دهد.

ناحیه وابستگی حل عددی در روش FTBS
شکل ۱: ناحیه وابستگی حل عددی در روش FTBS

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

یکی از شرط‌های اساسی برای همگرایی یک روش تفاضل محدود که برای معادله دیفرانسیل با مشتقات جزئی (PDE) هذلولوی نوشته شده، این است که ناحیه وابستگی حل عددی شامل «ناحیه وابستگی ریاضی» (Mathematical Domain of Dependence) باشد. این موضوع را با استفاده از شرط Courant-Friedrichs-Levy یا شرط کوارنت مورد بررسی قرار می‌دهند. این شرط را به صورت اختصاری با نماد CFL نیز نمایش می‌دهند. این شرط به افتخار نویسنده آن که برای اولین بار این شرایط را مورد بررسی قرار داد، به این نام شناخته می‌شود.

به عنوان مثال، برای معادله یک بعدی جابه‌جایی که با استفاده از «الگوریتم بادسو» (Upwind Scheme) یا آپویند گسسته شده، شرط کورانت یا CFL رابطه زیر را برای پایداری حل الزامی می‌داند.

عدد کورانت CFL
رابطه ۱

عددی که در بالا نشان داده شده است یعنی $${|x| \Delta t} \over {\Delta x}$$ را عدد کورانت یا عدد CFL می‌نامند. همچنین همانطور که اشاره شد، این عدد را به صورت خلاصه با نماد CFL نیز نمایش می‌دهند.

به صورت کلی می‌توان بیان کرد که شرط پایداری برای روش‌های تفاضل محدود «صریح» (Explicit) این است که عدد کورانت یا عدد CFL با استفاده از یک عدد ثابت محدود شده باشد. این عدد ثابت به روش عددی خاصی که استفاده شده، وابستگی شدیدی دارد.

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

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

معادله نیمه گسسته
رابطه ۲

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

فرض کنید که (U(t پاسخ دقیق ریاضی معادله نیمه گسسته فوق باشد. در این حالت، آشفتگی و اختلالی را به فرم (e(t برای این معادله در نظر بگیرید که با پاسخ دقیق ((U(t) جمع می‌شود. با استفاده از این فرضیات، فرم نهایی پاسخ به شکل زیر در می‌آید.

آشفتگی و اختلال موجود در جریان
رابطه ۳

در ادامه به دنبال یافتن پاسخ برای دو سوال اساسی در روش تفاضل محدود هستیم. سوال اول این است که آشفتگی و اختلال (e(t برای یک مسئله نیمه گسسته، در طول زمان رشد می‌کند؟ سوال بعدی این است که برای ارضای پایداری یک روش انتگرالی زمانی چه محدودیت‌هایی وجود دارد و گام زمانی چگونه باید انتخاب شود؟

بنابراین در گام اول و به منظور یافتن پاسخی مناسب برای سوال‌های بالا، (V(t را در رابطه ۲ قرار می‌دهیم. بنابراین رابطه ۲ به شکل زیر در می‌آید.

معادله نیمه گسسته روش تفاضل محدود
رابطه ۴

حال از رابطه ۳ برای نوشتن (V(t استفاده می‌کنیم و مقادیر موجود در این رابطه را در معادله فوق (رابطه 4) جایگذاری می‌کنیم.

معادله نیمه گسسته روش تفاضل محدود
رابطه ۵

با مقایسه معادله فوق (رابطه ۵) با معادله ابتدای این بخش (رابطه ۲)، ساده سازی‌ها را انجام می‌دهیم و در نهایت رابطه زیر به دست می‌آید.

معادله نیمه گسسته روش تفاضل محدود
رابطه ۶

بنابراین اغتشاش و اختلال موجود در جریان باید معادله همگن بالا را ارضا کند. مطالعه معادلات خطی مانند معادله بالا، نشان می‌دهد که مقدار (e(t با افزایش زمان بدون محدودیت افزایش می‌یابد. البته توجه کنید که شرط این موضوع این است که تمام قسمت‌های حقیقی «مقادیر ویژه» (Eigenvalue) ماتریس A، مثبت باشند.

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

همانطور که در تحلیل مقدار ویژه مربوط به روش انتگرالی «معادلات دیفرانسیل معمولی» (Ordinary Differential Equation) مشاهده شد، روش انتگرالی باید برای تمام مقادیر ویژه موجود در مسئله، پایدار باشد. یکی از روش‌هایی که با استفاده از آن می‌توانیم پایداری را مورد مطالعه قرار دهیم، این است که مقادیر ویژه مقیاس شده به وسیله گام زمانی را در صفحه مختلط λ Δt رسم کنیم و در نهایت با ناحیه مطلوب برای پایداری انتگرال معادله دیفرانسیل معمولی مقایسه کنیم. در این شرایط مقدار گام زمانی برای پایدار شدن تمام مقادیر ویژه انتگرال معادله دیفرانسیل معمولی مشخص می‌شوند.

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

مثال 1

در این مثال به بررسی ماتریس پایداری روش FTCS برای «جابه‌جایی یک بعدی» ($$1-D \enspace Convection$$) پرداخته می‌شود.

به صورت مختصر می‌توان بیان کرد که روش FTCS خلاصه شده عبارت Forward Time, Central Space است که با اعمال آن، مشتق زمانی به صورت پیش‌رو و مشتق مکانی به صورت مرکزی در سیستم معادلات اعمال می‌شود.

روش FTCS در مطلب «روش تفاضل محدود (Finite Difference Method) — از صفر تا صد» به صورت کامل مورد مطالعه قرار گرفته است و شیوه نمایش آن برای معادله جابه‌جایی یک بعدی نیز بررسی شد. بنابراین با اعمال روش FTCS به معادله جابه‌جایی یک بعدی، فرم نهایی معادله به شکل زیر در می‌آید.

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

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

همانطور که مشاهده می‌شود با استفاده از معادله 7، اطلاعات جدید برای هر U مربوط به نقطه i قابل محاسبه است. این موضوع با استفاده از کد متلب مربوط به آن به خوبی نشان داده شده است. در این شرایط، تنها در شرایطی نیاز به محاسبه ماتریس A است که بررسی پایداری و محاسبه مقادیر ویژه موضوع مهمی برای ما باشد. در ادامه دقیقا همین عمل انجام خواهد شد.

به عبارت دیگر، روندی که در ادامه طی می‌شود به ترتیب این است که در ابتدا ماتریس A و مقادیر ویژه ماتریس A محاسبه می‌شوند و در ادامه مقادیر ویژه مقیاس شده با استفاده از گام زمانی (Δt)، توسط ناحیه پایداری اویلر پیشرو پوشش داده می‌شود.

توجه کنید که این کد می‌تواند شرایط مرزی مانند «پریودیک» (Periodic) و «جریان ورودی» (Inflow) – «جریان خروجی» (Outflow) را مورد مطالعه قرار دهد. این شرایط در مثال موجود در مطلب «روش تفاضل محدود» بررسی شد. بنابراین در ادامه به بررسی پایداری و مقادیر ویژه در هر دو حالت شرایط مرزی پرداخته می‌شود.

در ادامه منحنی λ Δt برای حالتی رسم می‌شود که در آن مقدار CFL برابر با یک قرار داده شده است. این موضوع در شکل‌های ادامه این مطلب به خوبی نمایش داده شده است. توجه کنید که برای یک مسئله یک بعدی عدد CFL به شکل زیر تعریف می‌شود.

عدد کورانت CFL
رابطه ۸

همانطور که در شکل‌های بالا مشاهده می‌شود، شرط مرزی پریودیک، مقادیر ویژه موهومی تولید می‌کند که این مقادیر با فاصله گرفتن از مرکز به سمت $$\pm i$$ حرکت می‌کند. توجه کنید که در حقیقت شرط مرزی پریودیک مقدار ویژه برابر با صفر تولید می‌کند، بنابراین ماتریس A منفرد و تکین است.

نکته دیگر این است که در شرایط مرزی ورودی و خروجی، مقادیر ویژه درون ناحیه با قسمت حقیقی منفی حرکت می‌کند. اما در این حالت نیز با فاصله گرفتن از مرکز، به مقادیر موهومی $$\pm i$$ میل می‌کند.

توجه کنید که برخی از مقادیر λ Δt وجود دارند که در خارج از ناحیه پایداری اویلر قرار گرفته‌اند. در نتیجه ما می‌توانیم با استفاده از کاهش مقدار گام زمانی تمام مقادیر λ Δt را به درون ناحیه پایداری منتقل کنیم. اما اثبات می‌شود که این کار به صورت عملی غیر ممکن است؛ زیرا در نهایت مقادیر ویژه به سمت $$\pm i$$ میل می‌کنند. بنابراین می‌توان نتیجه گرفت که هیچ مقدار محدودی برای گام زمانی Δt وجود ندارد که به کمک آن بتوان مقادیر ویژه را به داخل ناحیه پایداری دایروی مربوط به روش اویلر پیشرو منتقل کرد.

توضیحات بالا نشان می‌دهد که روش FTCS برای جابه‌جایی ناپایدار است.

بررسی پایداری روش تفاضل محدود در شرایط مرزی پریودیک
شکل ۲: شرایط مرزی پریودیک
بررسی پایداری روش تفاضل محدود در شرایط مرزی جریان ورودی و خروجی
شکل ۳: شرایط مرزی دریکله جریان ورودی و بادسو جریان خروجی

نکته مهم دیگری که در ادامه به صورت دقیق مورد مطالعه قرار می‌گیرد، بررسی شیوه تغییرات طیف مقادیر ویژه ماتریس A با تغییر اندازه شبکه Δx است.

در بخش قبل مشاهده شد که برای معادله «جابه‌جایی» (Convection) با تغییر و اصلاح مقدار Δx باید مقدار Δt نیز اصلاح شود تا عدد CFL ثابت باقی بماند و پایداری حفظ شود. بنابراین تحت شرایط بیان شده و برای یک CFL ثابت، انتظار داریم که طیف مقدار ویژه با تغییر تعداد نقاط تغییر خاصی نکند.

شکل‌های زیر مقادیر ویژه هنگام حضور شرایط مرزی پریودیک و دریکله را در حالتی رسم کرده که پارامتر Δx با ضریب 10 اصلاح شده و به عبارت دیگر تعداد نقاط شبکه حل ۱۰ برابر شده است.

بررسی پایداری روش تفاضل محدود در شرایط مرزی پریودیک
شکل 4: شرایط مرزی پریودیک
بررسی پایداری روش تفاضل محدود در شرایط مرزی جریان ورودی و خروجی
شکل 5: شرایط مرزی دریکله جریان ورودی و بادسو جریان خروجی

همانطور که مشاهده می‌شود، در حالتی که شرایط مرزی به صورت پریودیک است، ما همچنان مقادیر ویژه موهومی که به $$\pm i$$ میل می‌کند را مشاهده می‌کنیم و زمانی که شرایط مرزی دریکله داریم، همچنان مقادیر به صورت اندک در ناحیه منفی قرار می‌گیرند ولی نسبت به مش درشت تر یعنی شکل 3 مقادیر به سمت محور موهومی نزدیک‌تر شده‌اند.

توجه کنید که در شکل‌های بالا نیز Δt طوری در نظر گرفته شده است که همانند حالت قبل، پارامتر CFL برابر با یک باشد.

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

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

«ماتریس دوری یا گردشی» (Circulant Matrix)
رابطه ۹

برای آنکه به صورت دقیق ویژگی‌های ماتریس بالا را متوجه شوید، کافی است که به ستون‌های آن، عناصر موجود در قطر اصلی، شروع و پایان هریک از ستون‌ها دقت کنید. این نوع از ماتریس را «ماتریس دوری یا گردشی» (Circulant Matrix) نیز می‌نامند و مقادیر ویژه آن با استفاده از رابطه زیر قابل محاسبه هستند.

مقادیر ویژه روش تفاضل محدود
رابطه ۱۰

مثال 2

همانطور که در مثال قبل مشاهده شد، زمانی که شرط مرزی پریودیک اعمال می‌شود، گسسته سازی مرکزی معادله مربوط به پدیده یک بعدی جابه‌جایی، نتایجی را به صورت موهومی ارائه می‌کند و زمانی که مقادیر ویژه با استفاده از گام زمانی برای حالتی که CFL برابر با یک است، مقیاس می‌شود، مقادیر ویژه در طول محور بین مقادیر $$\pm i$$ دیده می شوند. این موضوع در شکل‌های بالا به صورت دقیق نشان داده شده است.

توجه کنید که برای یک مسئله «جابه‌جایی» (Convection) که سرعت در آن مقداری ثابت دارد و شرط مرزی پریودیک اعمال شده است، ماتریس به صورت دوری یا گردشی خواهد بود. بنابراین ما می‌توانیم از معادله 10 استفاده کنیم و مقادیر ویژه را به صورت تحلیلی محاسبه کنیم. بنابراین در ابتدا ضرایب aj را محاسبه می‌کنیم.

برای گسسته سازی مرکزی داریم:

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

توجه کنید که مقادیر ضرایب aبرای سایر jها برابر با صفر است. با قرار دادن مقادیر ضرایب aدر معادله ۱۰، مقادیر ویژه به فرم زیر محاسبه می‌شوند.

مقدار ویژه ماتریس ضرایب روش تفاضل محدود
رابطه ۱2

برای ساده سازی معادله بالا از این نکته استفاده می‌شود که رابطه $$e^{i 2 \pi n} = 1$$ برقرار است. بنابراین رابطه فوق به شکل زیر در می‌آید.

مقدار ویژه ماتریس ضرایب روش تفاضل محدود
رابطه ۱۳

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

بررسی پایداری روش تفاضل محدود
رابطه ۱۴

بنابراین همانطور که در مثال 1 نیز به صورت دقیق مورد بررسی قرار گرفت، مقادیر ویژه مربوط به حل یک بعدی جابه‌جایی به صورت موهومی هستند. این موضوع در اینجا نیز نشان داده شد. بنابراین در  شرایط خاص این مسئله نیز زمانی که شرط مرزی به صورت پریودیک است، مقادیر ویژه به صورت موهومی هستند. نکته دیگری که باید به آن اشاره کرد این است که وقتی مقدار CFL برابر با یک است ($$CFL = |u| \Delta t / \Delta x = 1$$)، مقادیر ویژه در محدوده $$\pm i$$ قرار می‌گیرند.

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

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

در صورتی که به مباحث ارائه شده، علاقه‌مند هستید و قصد یادگیری در زمینه‌های مطرح شده در مکانیک سیالات را دارید، آموز‌ش‌های زیر به شما پیشنهاد می‌شود:

^^

بر اساس رای 1 نفر

آیا این مطلب برای شما مفید بود؟

نظر شما چیست؟

نشانی ایمیل شما منتشر نخواهد شد. بخش‌های موردنیاز علامت‌گذاری شده‌اند *

برچسب‌ها