روش تفاضل محدود (Finite Difference Method) – از صفر تا صد

تقریب‌ تفاضل محدود

در مطالب قبلی وبلاگ فرادرس، معادلات دیفرانسیل معمولی (Ordinary Differential Equation) به صورت دقیق و با استفاده از تقریب سری تیلور مورد بررسی قرار گرفتند. این تقریب برای مشتق در زمان $$\partial U \over \partial t$$ نوشته شد. با استفاده از این روش، تقریبی برای مشتق U بر حسب زمان و در زمان n به دست آمد که این تقریب با نماد زیر نشان داده می‌شود.

بنابراین همانطور که در ابتدای مطلب اشاره شد، با استفاده از روش تفاضل محدود می‌توان یک محیط پیوسته را به محیط گسسته تبدیل کرد و بازنویسی معادلات و محاسبه مشتق‌های جزئی را در آن انجام داد. فواصل نقاط مختلف گسسته شده در روش تفاضل محدود برابر با $$\triangle x$$ در نظر گرفته می‌شود (این فواصل می‌توانند با یکدیگر برابر یا متفاوت باشند) و شکل زیر به خوبی یک «شبکه» (Mesh) و ساختار یک بعدی برای تقریب تفاضل محدود را به تصویر کشیده است.

در ادامه هدف ما محاسبه یک تقریب تفاضل محدود برای $$\partial U \over \partial x$$ در نقطه $$x _ i$$ است. برای یک تابع مشتق‌پذیر U، مشتق جزئی در نقطه $$x _ i$$ را می‌توان با استفاده از رابطه زیر مورد محاسبه قرار داد.

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

این روش به نام روش «تفاضل مرکزی» (Central Difference) معروف است و به عملگر δ_2x، «عملگر تفاضل مرکزی» (Central Difference Operator) گفته می‌شود.

مشکل اصلی روش تفاضل محدود مرکزی این است که وقتی مقدار U در $$x _ i + \Delta x$$ و $$x _ i - \Delta x$$ با یکدیگر برابر باشند، مشتق در نقطه $$x _ i$$ را برابر با صفر محاسبه می‌کند. این در حالی است که در بسیاری از توابع نوسانی این گونه نیست و مشتق مقداری مثبت یا منفی دارد. بنابراین از روش‌هایی با مراتب بالاتر و یا روش‌های «یک طرفه» (One-Side) نیز برای محاسبه مشتق جزئی استفاده می‌شود. برای مثال یک «تقریب تفاضل پس‌رو» (Backward Difference Approximation) به صورت زیر قابل بیان است.

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

خطای برشی تقریب تفاضل محدود

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

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

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

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

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

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

مثال

تقریب تفاضل محدودی برای محاسبه مشتق جزئی مرتبه اول با استفاده از چهار نقطه بنویسید.

در این قسمت برای محاسبه تقریب تفاضل محدود مشتق مرتبه اول از چهار نقطه $$U _ {i}$$ ، $$U _ {i-1}$$ ، $$U _ {i+1}$$ و  $$U _ {i+2}$$ استفاده می‌شود. بنابراین در ابتدا بسط تیلور برای این چهار نقطه را به صورت زیر می‌نویسیم.

در ادامه و با استفاده از رابطه‌ای که در ابتدای بخش «تقریب تفاضل محدود برای مشتق‌های مراتب بالاتر» بیان شد، تقریب تفاضل محدود را می‌نویسیم. بنابراین داریم:

بنابراین با توجه به روابطی که در بالا ارائه شد، 4 ضریب مجهول در تفاضل محدود این عبارت حضور دارند. این چهار  ضریب مجهول به ترتیب  $$\alpha_ {i-1}$$ ، $$\alpha_ {i}$$،  $$\alpha_ {i+1}$$ و $$\alpha_ {i+2}$$ هستند.

در ادامه، ضرایب را طوری انتخاب می‌کنیم که ضریب U_xi برابر با ۱ و ضرایب سایر موارد موجود یعنی U_i، U_xxi و U_xxxi برابر با صفر باشد. بنابراین چهار معادله و چهار مجهول موجود است که می‌توان این دستگاه معادلات را به فرم ماتریس زیر بیان کرد.

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

بر این اساس با قرار دادن ضرایب بالا در معادله اولیه، تقریب تفاضل محدود با دقت مرتبه سه برای نقاط $$U _ {i}$$ ، $$U _ {i-1}$$ ، $$U _ {i+1}$$ و  $$U _ {i+2}$$ به شکل زیر در می‌آید.

پیشنهاد ما این است که یک دور مسئله بالا را خودتان حل کنید و نتیجه آن را با حل ارائه شده در این قسمت مقایسه کنید و بعد از آن و برای تمرین بیشتر، تقریب تفاضل محدود را برای مشتق سوم در نقطه i یعنی U_xxxi با استفاده از 5 نقطه  $$U _ {i}$$ ، $$U_ {i-2}$$، $$U _ {i-1}$$ ، $$U _ {i+1}$$ و  $$U _ {i+2}$$ بنویسید. در نهایت خطای برشی و دقت رابطه به دست آمده را محاسبه کنید.

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

تقریب تفاضل محدود در حالت دو بعدی

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

پارامتری مانند U را در نظر بگیرید که مقدار آن در مختصات (i, j) با U_i,j نشان داده می‌شود. در این شرایط، تقریب تفاضل مرکزی برای مشتق مرتبه اول در راستای x و y را می‌توان با استفاده از رابطه زیر نمایش داد.

نکته مهمی که باید به آن اشاره کرد این است که مجموعه نقاطی که در تفاضل محدود دو بعدی و سه بعدی در نظر گرفته می‌شود شامل نقاطی است که در اطراف U_i,j قرار دارند.

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

روش تفاضل محدود

معادله یک بعدی «نفوذ و جابه‌جایی» (Convection-Diffusion) را در نظر بگیرید. این معادله را می‌توان به شکل زیر بیان کرد.

فیلم آموزش دینامیک‌ سیالات‌ محاسباتی CFD 1 در فرادرس

کلیک کنید

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

رابطه فوق یک معادله دیفرانسیل معمولی (توجه کنید که علامت ∂ به d تبدیل شده) برای U_i را نشان می‌دهد و شامل مقادیر U در نقاط $$U _ {i-1}$$، $$U _ {i}$$ و $$U _ {i+1}$$ می‌شود.

در ادامه مقادیر پارامتر U در تمامی نقاط شبکه را با استفاده از یک بردار به شکل زیر نمایش می‌دهیم.

در نهایت این معادلات به صورت مجموعه‌ای از معادلات دیفرانسیل معمولی در می‌آیند که می‌توان آن‌ها را به شکل خلاصه شده زیر نمایش داد.

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

در رابطه بالا، ماتریس A را می‌توان به شکل زیر نمایش داد.

توجه کنید که سطر i ام ماتریس بالا، ضرایب مربوط به معادله دیفرانسیل معمولی پارامتر U در نقطه i را نشان می‌دهند. نکته مهم دیگری که باید به آن اشاره کرد این است که به غیر از حالاتی که در آن‌ها شرایط مرزی حضور دارد، مقدار A_i,j ضریب α_j موجود در تقریب تفاضل محدود را نمایش می‌دهد.

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

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

رابطه زیر که در بالا نیز مورد مطالعه قرار گرفت را دوباره در نظر بگیرید.

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

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

مثال

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

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

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

در این مثال، ناحیه حل را برای تولید شبکه در محدود 0 تا ۱ و شرایط مرزی را پریودیک در نظر می‌گیریم. کد متلب این مثال در ادامه آورده شده است.

نتایج این حل با استفاده از روش تفاضل محدود در شکل‌های زیر برای زمان t=0.5 ، t=0.25 و t=1.0 به تصویر کشیده شده است.

در مباحث بعدی وبلاگ فرادرس نشان داده می‌شود که الگوریتم FTCS برای هر مقداری از Δt برای معادله «جابه‌جایی خالص» (Pure Convection) ناپایدار است. این موضوع در این مثال نیز به خوبی نشان داده شده است.

شرایط مرزی

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

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

رابطه بالا، سرعت را در نقطه و «نود» (Node) یا گره اول یعنی i=0 نشان می‌دهد. این نود، مرز ورودی را نشان می‌دهد. توجه کنید که در حل‌های عددی نقاط موجود در شبکه و مش حل را نود یا گره می‌نامند.

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

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

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

یکی دیگر از موارد بسیار مهم در مسیر اجرای یک حل عددی تعیین پارامتر b در رابطه ۲۴ است. در واقع در مسائلی که شرط مرزی دریکله در آن به کار رفته، پارامتر b مقدار شرط مرزی را نیز در بر می‌گیرد. برای مثال رابطه بازنویسی شده به وسیله شرایط مرزی در مسئله نفوذ و جابه‌جایی (رابطه ۳۷) در نظر بگیرید. در این رابطه مقدار b را می‌توان به شکل زیر محاسبه کرد.

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

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

• دینامیک سیالات محاسباتی (CFD) — از صفر تا صد
• معادلات ناویر استوکس (Navier Stokes) — از صفر تا صد
• پایداری روش تفاضل محدود — به زبان ساده
• آیرودینامیک (Aerodynamics) چیست؟ — از صفر تا صد
• پیوستگی و بقای جرم در سیالات — از صفر تا صد
• مومنتوم خطی (Linear Momentum) در سیالات — از صفر تا صد
• توربوماشین (Turbomachinery) — به زبان ساده

^^