پروژه شماره 1 درس دینامیک سیالات محاسباتی 1
مقدمه :
روشهای عددی از کم هزینه ترین روشهای حل مسایل مهندسی بوده و با صرف زمان اندک و همچنین حل هندسه های گوناگون مزایا و برتری آن بر سایر روشها به اثبات رسیده است .دینامیک سیالات محاسباتی عبارت است از شبیه سازی عددی مسایل مهندسی به همراه روشهای عددی فوق ، در ادامه حل عددی معادله موج که از اولین و ساده ترین مسایل آغازین CFD بوده و با روشهای گسسته سازی متفاوتی بررسی گردیده است ارائه خواهد گردید. پیش از ارائه نرم افزار های مختلف برای انجام محاسبات موارد فوق(که امروزه به وفور یافت می شوند) ، کد نویس در محیط های فرترن ، مطلب ، سی و... معمول بوده و در این مقوله با استفاده از کدنویسی فرترن 90 به حل مساله اقدام می نماییم.
صورت مساله :
معادله موج زیر را در نظر بگیرید :
+a =0 a : constant
الف : معادله فوق را با استفاده از روش های زیر :
1 – روش اویلر صریح FTBS
2 – روش اویلر ضمنی FTCS
3 – روش لکس – فردریک
4 – روش لکس – وندروف
برای حالتی که :
a=1 ∆x = 0.5 ∆t = 0.01
و شرایط مرزی و اولیه :
U(0,t) = 0
U (x,0) = 1 for 10∆x≤x≤15∆x
<10 ∆x & x>15∆x for U(x,0) = 0
و با تعداد نقاط 61 حل کنید.
نتایج را بر حسب برای t=0 و t=100∆t و t=200∆t رسم کنید.
بررسی مساله :
بر اساس مطالعات انجام شده و بررسی معیار پایداری ، معادله موج برای سه روش 1 و 3 و4 در حالت صریح پایدار است . اما برای روش 2 در حالت صریح ناپایدار است .لذا از روش 2 که برای حالت ضمنی پایدار است استفاده می کنیم.
الف. حل روشهای صریح :
ابتدا معادله را بر اساس روابط استاندارد هر روش که در مراجع حل عددی موجود می باشد گسسته نموده و پس از چک کردن مجدد پایداری مساله(با استفاده از معادله اصلاح شده) بر اساس داده های صورت مساله به ادامه حل می پردازیم. (معادله پیوسته ای که به صورت معادله تفاضل محدود در آمده و تنها یک مجهول را بر حسب مقادیر معلوم بیان می کند روش صریح می گویند) .
ب. حل روشهای ضمنی : در این حالت هرگاه معادله را گسسته کنیم ،بیش از یک مجهول در معادله تفاضل محدود ظاهر می شود. در نتیجه مجموعه ای ازطریق حل دستگاه معادلات را باید حل کرد که این حل مستلزم مدت زمانی طولانی محاسبات در هر گام زمانی می شود.(در این روش از الگوریتم توماس جهت حل دستگاه معادلات استفاده شده است).
مزیت مهم روش های ضمنی پایداری معادلات تفاضل محدود می باشد.
بنابراین گام های زمانی بزرگتر در این روش مجاز است.
فیزیک و هندسه مساله :
مساله یک بعدی در راستای x بوده و دراین مسئله می خواهیم مقدار را برحسب یعنی ( ( را برای 61 نقطه در سه زمان t=0 و t=100∆t و t=200∆t را بدست آوریم سپس رسم کنیم .(مقدار تابع در ترسیم نمودارها دما فرض شده است).
چنانکه در صورت مساله مشخص است معادله خطی می باشد.
گسسته سازی :
1- FTBS :
+a =0 → + a =0 →
→
می توان نشان داد که خطای برشی عبارت از O(∆t , ∆x ) می باشد . همانطور که مشاهده می شود در معادله اصلاح شده این روش
(1 – ) مشتق زوج بالاترین مرتبه از نقطه نظر خطا را در خطای برشی دارا بوده پس باعث تحمیل اتلاف مصنوعی به جواب وتغییر شکل موج می گردد.جهت بررسی پایداری روش با استفاده از روش فون نیومن خواهیم داشت :
=
= 1 - (1- )
پس :
= 1 –
= ( 1 - 2 )2 + ??????
پس با توجه به اینکه همواره می باشد برای اینکه 1 باشد باید گردد یعنی پایداری روش در حالت 0≤ امکان پذیر است و روش عددی FTBS مورد استفاده برای حل عددی معادله موج در این حالت پایدار می باشد.
2- روش اویلر ضمنیFTCS :
معادله موج را از روش اویلر ضمنی گسسته می کنیم.
می توان نشان داد که خطای برشی عبارت از O( ) می باشد. رابطه فوق تشکیل یک دستگاه از معادلات جبری سه قطری می دهد. جهت تحلیل خطا از معادله اصلاح شده استفاده می شود.
در کل، روش های ضمنی نیازمند زمان محاسباتی بیشتری در هر گام زمانی هستند ولی مطمینا" دارای این خاصیت هستند که می توان از گام زمانی بزرگتری برای حل استفاده نمود زیرا معمولا" بدون هیچگونه قید و شرطی پایدارهستند .پایداری روش را می توان با استفاده از روش فون نیومن مورد بررسی قرار داد. داشتیم: g=
پس به ازای جمیع مقادیر است و روش عددی فوق برای حل معادله موج همواره پایدار است.
3-روش لکس فردریک :
+ =
می توان نشان داد که خطای برشی عبارت از می باشد. این روش همواره سازگار نیست به علت اینکه اگر به سمت صفر میل کند ممکن است صفر نشود.این روش با خطای اتلافی زیادش برای مشهور می باشد. جهت بررسی پایداری روش، با استفاده از روش فون نیومن خواهیم داشت:
بنابراین پایداری این روش در حالت امکان پذیر است و روش عددی لکس مورد نظر برای حل عددی معادله موج در این حالت پایدار می باشد.
4-روش لکس – وندروف( LAX-Wendroff):
با بسط سری تیلور برای داریم:
می توان نشان داد که خطای برشی عبارت از O( ) می باشد. همانطور که مشاهده می شود در معادله اصلاح شده این روش:
مشتق فرد بالاترین مرتبه از نقطه نظر خطا را در خطای برشی دارا بوده،پس این روش دارای خاصیت پراکندگی است.ولی با توجه به داشتن مشتق زوج در جملات ابتدایی خطای برشی،در کل دارای خاصیت پخش می باشد(خاصیت پراکندگی و اتلافی ).جهت بررسی پایداری روش،با استفاده از روش فون نیومن خواهیم داشت:
g =1-
=1-2
پس این روش به ازاء ،آنگاه می شود و روش پایدار می گردد.روش دو مرحله ای لکس – وندروف نیز وجود دارد که از نظر خطای برشی و پایداری همانند روش فوق می باشد.
شروط مرزی :
شروط مرزی در این مساله شرط دریکله می باشد (
شرط مرزی دیریکله یعنی مقدار خود تابع روی مرز مشخص شده باشد) و نیازی به گسسته کردن ندارد.
فلوچارت و شرح کد :
1-ارائه توضیحات نویسنده کد و شرح کارایی کد
2-تعریف پارامترها (عموما در برنامه ها این پارامترها به صورت نامحدود تعریف می گردد ولی به واسطه مشخص بودن تعداد
GRID این پارامترها به صورت ثابت در این برنامه تعریف گردیده است (به خصوص در خصوص ماتریس ها).
3-تعریف پارامترهای اولیه نظیر طول مساله ، ضرائب ثابت معادله ، گام زمانی و... .
4- تولید شبکه مساله که در ماتریس X(i) ذخیره می گردد.
5-تعریف شرایط اولیه که در این قسمت بر اساس دستور مساله بالا به هر مقدار مکانی یک مقدار
U تخصیص می یابد.
6- چاپ مقادیر اولیه(
t=0) که از خواسته های مساله می باشد
.
7-تعریف 2 حلقه زمانی برای تکرار محاسبات در 2 زمان خواسته شده
8- محاسبه زمان
n+1 برای هریک از مقادیر
U در هر گره در حلقه های فوق برای روشهای 1و3و4 (این مقدار بدست آمده در طی تکرار زمانی جهت زمانی آتی مقدار زمان گذشته فرض می گردد).
9-محاسبه ضرایب ماتریس 3 قطری و ماتریس مقادیر ثابت و اجرای الگوریتم توماس در طی حلقه های زمانی براسی روش 2
10-چاپ نتایج به همراه اطلاعات مورد نیاز نرم افزار
TECPLOT
کد نوشته شده:
به دلیل تفاوت فرمت کد و نرم افزار ورد جهت استفاده به فایلهای پیوست (شامل کلیه فایلهای مورد نیازجهت اجرا) مراجعه شود.
نتایج و نمودارها :
با تمهیدات در نظر گرفته شده در کد خروجی کلیه محاسبات در پایان هر حلقه زمانی بر حسب متد محاسبه در یک فایل جداگانه ذخیره گردیده است و به دلیل تفاوت فرمت از ذکر مقادیر خوداری می گردد و جهت مشاهده کافی است فایلهای زیر در یک محیط ویرایشگر متن گشوده شوند.
1-فایل با نام
Output.DAT حاوی مقادیر نتایج در زمان
t=0
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
2-فایل با نام
Output2.DAT حاوی مقادیر نتایج روش
FTBS در زمان
t=100dt
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
3-فایل با نام
Output3.DAT حاوی مقادیر نتایج روش
FTBS در زمان
t=200dt
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
4-فایل با نام
Output4.DAT حاوی مقادیر نتایج روش
FTCS در زمان
t=100dt
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
5-فایل با نام
Output5.DAT حاوی مقادیر نتایج روش
FTCS در زمان
t=200dt
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
6-فایل با نام
Output7.DAT حاوی مقادیر نتایج روش
Lax–
Friedrich در زمان
t=100dt
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
7-فایل با نام
Output8.DAT حاوی مقادیر نتایج روش
Lax–
Friedrich در زمان
t=200dt
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
8-فایل با نام
Output9.DAT حاوی مقادیر نتایج روش
Lax–
Wendroff در زمان
t=100dt
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
9-فایل با نام
Output10.DAT حاوی مقادیر نتایج روش
Lax–
Wendroff در زمان
t=200dt
این مقادیر در نمودار زیر توسط نرم افزار
TECPLOT ترسیم گردیده است.
بحث و بررسی :
چنانکه در نمایش اولیه (حالت زمان صفر ) مشخص گردیده است مساله مذکور حل گذرای یک معادله موج در بازه مکانی 61 گره در امتداد محور
X در 2 بازه زمانی به روشهای متفاوت می باشد و نظر به اینکه حالت تعادلی مساله از خواسته های مساله نمی باشد نمودارهای هر حالت برای زمانهای درخواست شده ارائه گردیده است. همچنین با توجه به اینکه کلیه روشهای فوق پایدار بوده و یا با رعایت شروط پایداری حل و بررسی شده اند تفاوت نمودارها تنها ناشی از الگوریتم های حل می باشد که قطعا در ادامه و تا زمان تعادل تقریبا همه به یک شکل تغییر خواهند یافت و فقط زمان رسیدن به زمان تعادل متفاوت خواهد بود.در شکل زیر گراف ترکیبی حالتهای فوق برای زمان نهایی
200dt ارائه گردیده است.
