(سایت پایان نامه ) –108- تحقیق علمی

1-2 دسته بندی کانال هادسته بندی کانال ها براساس قطر هیدرولیکی آنها به عنوان یک راه حل ساده و در عین حال کارا جهت بررسی رنج های مختلف مورد نظر می باشد. کاهش ابعاد کانال در فرایندهای مختلف، اثرات مختلفی داشته و ارائه یک شرط خاص منطبق بر متغیرهای فرایند ممکن است در نگاه اول گزینه ای بسیار نامناسب به نظر برسد، اما چنانچه تعداد فرایندها و متغیرهای حاکم بر انتقال از ابعاد معمولی به پدیده های مقیاس میکرو را در نظر بگیریم همان روش ساده دسته بندی کانال ها براساس ابعادشان را یک روش مناسب می بینیم که در کارهای علمی انجام شده در این زمینه نیز کاملا پذیرفته شده است. جدول (1-1) رنج های مختلف کانال را تحت جریان ها با نوع های متفاوت نشان می دهد:

شکل 1- SEQ شکل_1- * ARABIC 2 اثر قطر هیدرولیکی بر انتقال حرارت در میکرو کانال ها [3]وجود تعادل مابین نرخ انتقال حرارت و افت فشار موضوع مهمی در طراحی گذرگاه های جریان های خنک کننده برای دفع شارهای حرارتی بسیار بالایی است که در خنک کاری چیپ های ریز پردازنده ها با آن روبرو می شویم.
توسعه سیستم های میکر الکترو مکانیکی به طور عمومی نیارمند توسعه سیستم های دفع حرارتی است که متناسب با آن سیستم ها کوچک باشند. خنک کاری آینه هایی که در سیستم های لیزر به کار می روند شامل سیستم های خنک کننده ای است که باید حرارت بسیار زیادی را از سطح مقطعی بسیار کوچک دفع نمایند.
پیشرفت هایی که در بیولوژی و مهندسی ژنتیک رخ می دهد نیازمند انتقال تحت کنترل سیال و همچنین کنترل شرایط دمایی آن در گذرگاه هایی است که ابعاد آن از چندین میکرون بیشتر نمی باشد بنابراین داشتن درک درستی از جریان سیال و انتقال حرارت در این سیستم های با مقیاس میکرو جهت طراحی و بهره برداری از این سیستم ها ضروری می باشد.
حالت های متعددی از مسائل وجود دارد که نیاز به مطالعات بیشتر جهت اعلام نظر دارد. که برخی از آنها عبارتند از:
1. بررسی آزمایشگاهی اعتبار معادله های انتقال جریان لایه ای و جریان مغشوش، ضریب اصطکاک جریان لایه ای و معادلات انتقال حرارتی که براساس فرضیات کلاسیک برای کاربرد در میکرو کانال هایی که در آنها شاهد تغییری در پدیده های انتقال یا بروز پدیده فیزیکی جدیدی نیستیم، همچنین آزمایشاتی برای بررسی انتقال جرم و مومنتم و حرارت و فرایندهای انتقال جرم.
2. بررسی پدیده گذار از جریان لایه ای به جریان مغشوش، در مقیاس میکرو و با دقت به ارزیابی های آزمایشگاهی مربوط به این کار.
3. اثرات زبری نسبی سطح برای جریان با توجه به این که مقادیر مرتبط با زبری سطح که در میکرو کانال ها با آن سر و کار داریم بسیار بیشتر از مسائل با مقیاس معمولی می باشد همچنین اثر زبری بر گذار از جریان لایه ای به مغشوش، ضرایب اصطکاک و انتقال حرارت در این گذار.
4. بررسی اعتبار مقادیر ثابتی که به روش های آزمایشگاهی برای مقیاس معمولی به دست آمده اند، مانند افت های مربوط به تغییر سطح مقطع افت در خمدیگی لوله ها و نظایر آن که این مقادیر که از آزمایشات جریان سیال در مقیاس معمولی به دست آمده اند. مانند افت های مربوط به تغییر سطح مقطع افت در خمیدگی لوله ها و نظایر آن که این مقادیر که از آزمایشات جریان سیال در مقیاس معمولی حاصل شده اند باید با توجه به شرایط حاکم بر مقیاس میکرو نسبت به کار بردن آنها با توجه زیادی عمل کرد.
به چند دلیل ابعاد گذرگاه جریان در کاربردهای انتقال حرارت گرایش به سمت مقادیر کمتر و در مقیاس میکرو دارد:
1. نقش انتقال حرارت در مقیاس میکرو بسیار پر رنگ تر و موثرتر است.
2. نیاز تجهیزات میکرو الکترونیک برای انتقال بیشتر حرارت هم زمان با توسعه و تحول آن تجهیزات.
3. اضطراری که به واسطه کاربرد روز افزون و تنوع رو به گسترش تجهیزات مقیاس میکرو که نیاز به خنک کاری دارند، به وجود می آید.
با استفاده از کانال های با ابعاد کوچکتر به انتقال حرارت با کارآیی بالاتری دست می یابیم، هر چند که بر واحد طول افت فشار بیشتری را نیز شاهد خواهیم بود. چگالی حجمی بیشتر انتقال حرارتی- که لازمه تکنیک های تولید پیشرفته و طراحی های مسیرهای جریان پیچیده تر است- بر ضرورت توسعه میکرو کانال ها برای انتقال حرارت تاکید دارد، به طوری که بهینه سازی هر یک از کاربردهای متنوع میکرو کانال ها نتایج جدیدی برای ابعاد کانال به دست می دهد، به عنوان مثال در صنعت تبرید و سردخانه استفاده از تیوپ های به قطر 6 تا 8 میلیمتر با استفاده از میکرو پره ها دیگر جایگزین تیوپ های تخت با قطر زیاد شده است. در کاربردهای صنایع خودرو ابعاد رادیاتورها و اوپراتورها به حدود یک میلیمتر رسیده تا ما بین توان مورد نیاز جهت پمپ کردن، انتقال حرارت و تمیزی کل سیستم، توازن خوبی برقرار شود. همچنین در کاربردهای تهویه مطبوع ساختمان امکان اتصال سیستم های خنک کن تجهیزات الکترونیکی و میکرو الکترونیکی اتاق سرور به سیستم تهویه مطبوع ساختمان نیز در حال اجرایی شدن می باشد.

فصل دومجریان سیال در میکرو کانال
2-1 پیشگفتاربه عنوان یک تعریف میکرو کانال ها، کانال هایی هستند که ابعاد آنها از یک میلیمتر کمتر و از یک میکرومتر بیشتر می باشد. در ابعاد بالای یک میلیمتر، جریان خصوصیاتی مشابه بیشتر جریان های ماکروسکپیک را از خود بروز می دهد. در ابعاد پایین تر از یک میکرومتر، جریان با مشخصات جریان های مقیاس نانو بیشتر منطبق است. میکرو کانال ها از مواد گوناگونی چون شیشه، پلیمرها، سیلیکون و فلزات و با روش های گوناگون همچون میکرو ماشین کاری سطح، میکرو ماشین کاری حجمی، ماشین کاری، میکرو کاترها و روش هایی چون قالب سازی و ریختگی ساخته می شوند.
میکرو کانال ها به واسطه نسبت بالای سطح به حجمی که دارند و نیز حجم کوچک خود دارای مزایایی هم می باشند. نسبت بالای سطح به حجم به میکرو کانال ها نرخ بالایی از انتقال حرارت و جرم بخشیده و آنها را به ابزارهایی قدرتمند برای استفاده در مبدل های حرارتی کوچک تبدیل می سازد.
اثرات ساختار مولکولی در مایعات و گازها بسیار متفاوت است اگر عدد نودسن که به صورت Kn مساوی λ تقسیم بر Ls تعریف می شود که در آن λ مسیر توسط آزاد در گاز می باشد از عدد 3-10 بزرگتر باشد اثرات عدم تعادل شروع به نمایان شدن می کند. با افزایش عدد نودسن دیگر فرضیات محیط پیوسته و تئوری سیالات، غیر قابل استفاده خواهند بود و تجزیه و تحلیل این قبیل جریان ها نیازمند در نظر گرفتن پدیده های فیزیکی مختلف می باشد.
Kn=λLs (2.1)مایعات عموماً خاصیت تراکم ناپذیری دارند به این دلیل چگالی مایع در جریان داخل میکرو کانال ها، به عنوان تابعی از موقعیت نسبی به ابتدای کانال بسیار نزدیک به عددی ثابت باقی می ماند که این برخلاف تعداد زیاد گرادیان های فشاری است که مشخصات جریان با مقیاس میکرو را شکل می دهند و تعریف می کنند. این رفتار چگالی مایعات تحلیل جریان های مایع را به نسبت جریان های گازی بسیار ساده سازی می نماید. جایی که افت فشارهای زیاد در کانال ها منجر به انبساط زیاد و تغییرات زیاد در ظرفیت حرارتی می شوند.
2-2 قطر هیدرولیکیقطر هیدرولیکی یک عبارت مصطلح می باشد که برای جریان هایی که از درون کانال های با سطح مقطع غیر دایروی عبور می کنند بکار می رود. در عمل قطر هیدرولیکی به صورت زیر تعریف می شود:
Dh=4AP (2.2)که A سطح مقطع کانال و P محیط تر شده کانال می باشد.
برای کانال های با سطح مقطع دایروی داریم:
Dh=4AP=4πD24πD=D (2.3)همانطور که انتظار داشتیم قطر هیدرولیکی کانال با سطح مقطع دایروی برابر با قطر همان کانال می باشد.
برای محاسبه قطر هیدرولیکی در کانال های حلقوی داریم:
Dh=4AP=4(πro2-πri2)2πro+2πri=2(ro-ri) (2.4)که ro شعاع خارجی کانال و ri شعاع داخلی آن می باشد.
همچنین برای کانال های مستطیلی داریم:
Dh=4AP=4ab2(a+b)=2ab(a+b) (2.5)که a طول کانال و b عرض کانال می باشد.
2-3 طول توسعه یافتگی جریانجریان سیال پس از وارد شدن به لوله نیاز به طولی از لوله دارد تا منحنی سرعتش به شکل نهایی درآید و پس از آن تغییر نخواهد کرد. در این منطقه جریان بسیار شبیه به جریان لایه مرزی در حال رشد بوده که به سمت پایین دست جریان، ضخامت لایه مرزی در حال افزایش می باشد. در نهایت این لایه ها در مرکز کانال و در انتهای طول ورودی به یکدیگر ملحق می شوند. افت فشار از نقطه ورودی کانال تا موقعیت x از رابطه زیر به دست می آید:
P0-Px=fxDh+K(x)ρU22 (2.6)که در آن K(x) پارامتر افت فشار بوده که در شکل (2-1) برای یک مجرای دایره ای شکل و برای دو صفحه موازی آورده شده است [4].
همچنین طول توسعه یافتگی در جریان لایه ای تقریباً توسط رابطه زیر بدست می آید:
xDh=0.065Re (2.7)این امر در شکل (2-1) نیز بطور کامل قابل مشاهده است.

شکل 2- 1 پارامتر طول ورودی جریان K برای جریان لایه ای در ورودی مجرا2-4 حالت های انتقال حرارتمکانیزمی که حرارت به وسیله آن در یک سیستم تبدیل انرژی منتقل می شود بسیار پیچیده می باشد، با این حال سه حالت انتقال حرارت شناخته شده عبارتند از هدایت، جابجایی و تشعشع. جابجایی مکانیزمی از انتقال حرارت است که در اثر اختلاط یک بخش سیال با بخش دیگر تحت تاثیر حرکت جرم سیال اتفاق می افتد در حالی است که فرآیند واقعی، انتقال انرژی بین یک ذره سیال یا یک مولکول آن به دیگری همچون هدایت حرارتی می باشد انرژی ممکن است از نقطه ای از فضا به نقطه ای دیگر توسط جابجایی خود سیال نیز انجام شود به همین دلیل بررسی انتقال حرارت جابجایی به همراه حرکت سیال مورد بررسی قرار گرفته است. حرکت سیال ممکن است توسط عوامل خارجی مکانیکی مانند پمپ انجام پذیرد که در این حالت بدان جابجایی اجباری گفته می شود و اگر حرکت سیال به دلیل تفاوت در چگالی ناشی از اختلاف درجه حرارت موجود در جرم سیال باشد فرآیند را جابجایی آزاد یا جابجایی طبیعی می نامیم. انتقال حرارت های مهم در فاز مایع- بخار یعنی تبخیر و میعان نیز به عنوان مکانیزم های جابجایی در نظر گرفته می شوند چرا که در آنها نیز سهم حرکت سیال در انتقال حرارت بسیار است، البته در این دو حالت، پیچیدگی های دیگری مانند تبادل حرارت نهان نیز وجود خواهد داشت [5].
2-5 فرضیه پیوستگیدر تحلیل انتقال حرارت جابجایی در یک سیال، حرکت سیال بایستی به دقت همراه با فرآیند انتقال حرارت مورد مطالعه قرار گیرد. توصیف حرکت سیال در شکل بنیادی خود شامل مطالعه رفتار کلیه ذرات مجزای (مانند مولکول ها) تشکیل دهنده آن سیال می باشد. اساسی ترین رهیافت در تحلیل حرارت جابجایی نوشتن معادلات مکانیکی و ترمودینامیکی برای هر یک از ذرات تشکیل دهنده به طور مجزا و یا برای یک گروه آماری از ذرات تشکیل دهنده آن و به همراه در نظر گرفتن شرایط اولیه حاکم بر مسئله می باشد.
در بسیاری از کاربردها توجه اصلی به رفتار مولکولی سیال نیست، بلکه به اثرات متوسط یا اثرات ماکروسکپیک تعداد زیادی از مولکول ها توجه می شود، چرا که ما معمولاً همین مقادیر ماکروسکپیک را اندازه گرفته و مورد توجه قرار می دهیم. بنابراین در مطالعه انتقال حرارت جابجایی سیال به عنوان یک ماده با محیط پیوسته و قابل تقسیم تا بی نهایت در نظر گرفته شده و از ساختار مولکولی آن صرفنظر می شود، مدل محیط پیوسته تا زمانی که ابعاد و مسیر متوسط آزاد مولکول ها در قیاس با ابعاد مسئله و محیط به حد کافی کوچک باشند یعنی متوسط های آماری معنی دار باشند، دارای اعتبار می باشد.
با این وجود فرض محیط پیوسته زمانی که مسیر متوسط آزاد مولکول ها از لحاظ مرتبه بزرگی با کوچکترین طول موثر مسئله برابر می نماید دیگر دارای اعتبار نخواهد بود. در جریان های گازی انحراف حالت سیال از پیوستگی با عدد نودسن Kn نمایش داده می شود یعنی Kn=λL، مسیر متوسط آزاد λ همان متوسط فاصله ای است که مولکول بین دو برخورد موثر با مولکول های دیگری می پیماید و L طول مشخصه جریان می باشد. مدل های مناسب جهت جریان و انتقال حرارت برحسب عدد نودسن برای رژیم های مختلف گازی توسط Schaaf and Chamber که در جدول (2-1) ارائه گردیده اند [6]:
جدول 2-1
جریان پیوسته Kn<10-3جریان لغزشی 10-3<Kn<10-1جریان در مرحله گذار 10-1<Kn<10+1جریان آزاد مولکولی 10+1<Knدر رژیم جریان لغزشی مدل جریان پیوسته هم چنان برای محاسبه خواص جریان در نقاط دور از مرز جامد معتبر می باشد لیکن شرایط مرزی را می بایستی به منظور به حساب آوردن فعل و انفعال ناکامل بین مولکول های گاز و مرز جامد اصلاح نمود، در شرایط عادی 1/0 بیشتر جریان های گازی داخل میکرو کانال های چاه های حرارتی با طول مشخصه ای از مرتبه یک میکرومتر، عدد نودسن هائی کمتر از 1/0 دارند، بنابراین در اینگونه مسائل ما فقط رژیم جریان لغزشی داشته لیکن هنوز از مرحله جریان در حال گذار و جریان آزاد مولکولی فاصله داریم، مشخص است که فرض پیوستگی برای جریان های مایع داخل میکرو کانال های چاه های حرارتی معتبر می باشد.
2-6 اصول ترمودینامیکمناسب ترین چهارچوبی که مسائل انتقال حرارتی از این دست در آن می گنجد دیدگاه سیستم است که شامل مقادیری است که الزاماً ثابت نبوده و توسط یک مرز احاطه شده اند. مرزی می تواند کاملاً فیزیکی قسمتی فیزیکی قسمت دیگر مجازی یا به کل مجازی باشد. قوانین فیزیکی که مورد بررسی قرار می گیرد همواره برحسب مولفه های سیستم ارائه می شود.
یک حجم کنترل ناحیه مشخصی از فضا است که جرم، مومنتوم و انرژی می تواند از مرزهای آن عبور کرده همچنین جرم مومنتوم و انرژی در آن ذخیره می شود و نیروهای خارجی نیز می تواند بر آن اثر کند، تعریف دقیق یک تقسیم یا حجم کنترل می بایستی حداقل شامل تعریف سیستم مختصات باشد که این سیستم مختصات خود می تواند ایستا یا در حال حرکت باشد. مشخصه مورد توجه یک سیستم حالات آن است که شرطی از سیستم است که با خواص آن سیستم قابل توصیف است، یک خاصیت سیستم ممکن است براساس هر یک از مقادیری باشد که به حالت سیستم بر می گردد و منفصل از مسیری است که سیستم را به یک حالت مشخص می رساند اگر کلیه خواص یک سیستم بدون تغییر باقی بماند، چنین بیان می شود که سیستم در حالت تعادل است.
هر تغییری در یک یا بیشتر از خواص سیستم به معنی این است که تغییری در حالت سیستم روی داده است. مسیری از حالت های آتی که سیستم از آنها می گذرد فرآیند نامیده می شود. هنگامی که یک سیستم در حالت اولیه داده شده است تغییرات مختلفی را در حالت یا فرآیندهای خود می بینند و در نهایت به همان نقطه آغازین بر می گردد، چنین بیان می شود که سیستم در یک چرخه قرار گرفته است. خواص تنها زمانی قادر به توصیف حالت یک سیستم هستند که سیستم درحالت تعادل باشد، چنانچه انقال حرارتی بین دو سیستم که در تماس با یکدیگر قرار دارند رخ ندهد چنین گفته می شود که سیستم ها در تعادل حرارتی هستند. هر دو سیستمی که با هم در تعادل حرارتی باشند دارای درجه حرارت مساوی می باشند و چنانچه سیستم ها در تعادل حرارتی نباشند دماهای مختلفی خواهند داشت، در این حالت ممکن است انتقال انرژی از یک سیستم به سیستم دیگر صورت گیرد، بنابراین دما خاصیتی است که سطح حرارتی سیستم را اندازه می گیرد.
هنگامی که ماده ای وجود داشته باشد که قسمتی از آن مایع و قسمت دیگری و در حالت اشباع باشد کیفیت آن x توسط نسبت جرم بخار به جرم کل تعریف می شود کیفیت خاصیتی است که مقدار آن مابین صفر و یک می باشد. کیفیت تنها هنگامی معنی دارد که ماده در حالت اشباع یعنی دما و فشار اشباع قرار داشته باشد مقدار انرژی که می بایست به شکل حرارت به یک ماده در فشار ثابت انتقال یابد تا تغییر فاز صورت گیرد گرمای نهان نامیده می شود، این شامل تغییر در آنتالپی ماده به عنوان یک خاصیت ماده است که در شرایط اشباع برای هر یک از دوفاز مایع و بخار می باشد. گرمای تبخیر و جوش حرارت مورد نیاز جهت تبخیر کامل یک واحد جرم از مایع اشباع می باشد.
2-7 قوانین کلیقوانین کلی که به یک سیستم باز، مانند میکرو کانال های چاه های حرارتی اشاره دارند، می توانند به هر یک از فرم های انتگرالی یا دیفرانسیلی نوشته شوند. قانون بقای جرم به زبان ساده بیان می دارد زمانی که تبدیل جرم- انرژی نداشته باشیم جرم سیستم ثابت می ماند، لذا چنانچه چشمه یا چاه نداشته باشیم، Q=0 و نرخ تغییر جرم در حجم کنترل مساوی شار جرمی روی سطح کنترل خواهد بود. قانون دوم نیوتن درباره حرکت بیان می دارد که نیروی خالص F که در یک سیستم مختصات بر سیستم اثر می کند برابر با نرخ زمانی تغییرات مومنتوم خطی کل سیستم است. همچنین قانون بقای انرژی برای حجم کنترل بیان می دارد که نرخ تغییرات انرژی کل سیستم E برابر با مجموع نرخ زمانی تغییرات انرژی حجم کنترل و شار انرژی گذرنده از سطح کنترل می باشد.
قانون اول ترمودینامیک که بیان دیگری از بقای انرژی می باشد بیان می دارد که نرخ تغییرات انرژی برابر با اختلاف بین نرخ انتقال حرارت به سیستم و کار انجام شده توسط سیستم است. قانون دوم ترمودینامیک به معرفی آنتروپی، S به عنوان یک خاصیت سیستم می پردازد. این قانون بیان می دارد که نرخ تغییرات آنتروپی سیستم برابر یا بزرگتر از نرخ انتقال حرارت به سیستم تقسیم بر دمای سیستم در طی فرآیند انتقال حرارت می باشد، حتی در شرایطی که محاسبات آنتروپی مورد توجه نمی باشد، قانون دوم ترمودینامیک همچنان مهم بوده و معادل این مسئله است که نشان دهیم حرارت نمی تواند از سیستم با درجه حرارت کمتر به سیستمی با درجه حرارت بیشتر انتقال پیدا کند.
2-8 قوانین خاصقانون فوریه برای هدایت حرارتی براساس اصل پیوستگی بیان می گردد که شار حرارتی ناشی از هدایت در یک راستای معین (به عنوان مثال نرخ انتقال حرارت در واحد سطح) داخل یک محیط (خواه جامد، مایع یا گاز) به گرادیان دما در همان راستا بستگی دارد یعنی:
q''=-k∇T (2.8)که در آن q'' بردار شار حرارتی، K هدایت حرارتی و T دما می باشد. قانون نیوتون درباره سرد شدن چنین بیان می دارد که شار حرارتی از یک سطح جامد به وسیله جابجایی به سیال محیطی منتقل می شود، q'' به اختلاف دمای بین دمای سطح جامد Twو دمای جریان آزاد سیال T∞به شرح زیر وابسته می باشد:
q''=h( Tw- T∞) (2.9)که در این رابطه h ضریب انتقال حرارت است.
2-9 ساختار جریانرژیم های جریان ویسکوز بسته به ساختار جریان با لایه ای و مغشوش دسته بندی می شوند. در رژیم لایه ای ساختار جریان به وسیله حرکت ملایم درون لایه یا لایه ها مشخص می شود جریان در رژیم مغشوش با حرکت های سه بعدی اتفاقی ذرات سیال تحت تاثیر حرکت متوسط طبقه بندی می شود. این نوسانات مغشوش باعث تقویت بسیار زیاد انتقال حرارت جابجایی می گردند. با وجود که از لحاظ عملی جریان مغشوش هنگامی به وجود می آید که عدد رینولدز Re=ρUDhμ از یک مقدار بحرانی Recrit بیشتر شود. مقدار عدد رینولدز بحرانی به شرایط ورودی به کانال زبری سطح، لرزش هایی که به دیواره های کانال اعمال می شوند و هندسه سطح مقطع کانال بستگی دارد. Bhatti and Shah مقادیر Recrit را برای اشکال گوناگون از سطح مقطع کانال ارائه نمودند [7]. در کاربردهای عملی عدد رینولدز بحرانی را می توان به کمک رابطه زیر تخمین زد:
Re=ρUDhμ=2300 (2.10)که در آن Uسرعت متوسط جریان، Dh=4AS قطر هیدرولیکی کانال A و Sبه ترتیب مساحت سطح مقطع کانال ومحیط تر شده کانال می باشند. میکرو کانال ها معمولاً از لحاظ طول بیش از 1000 میکرومتر و از لحاظ قطر هیدرولیکی در حدود 10 میکرومت می باشند، سرعت متوسط جریان گاز با افت فشاری در حدود MPa 5/0 کمتر از ms 100 و عدد رینولدز مربوطه نیز کمتر از 100 می باشد. عدد رینولدز برای جریان های مایع به دلیل قوی تر بودن نیروهای ویسکوز در مایعات نسبت به گازها از این مقدار نیز کمتر باشد. به همین دلیل در بسیاری از کاربردها جریان داخل میکرو کانال ها را به صورت لایه ای در نظر گیرند. جریان مغشوش ممکن است در کانال هایی کوتاه و با قطر هیدرولیکی زیاد و تحت تاثیر افت فشارهای زیاد به وجود بیایند، بنابراین در بررسی میکرو کانال ها چندان مورد توجه نمی باشند.
2-10 طول ورودی حرارتیچنانچه دیواره های کانال گرم یا سرد شود لایه مرزی حرارتی نیز در طول سطوح داخلی کانال شکل خواهد گرفت. طول ورودی حرارتی Lt طول مورد نیاز کانال جهت توسعه یافتن جریان تا رسیدن به شرایط توسعه یافته خواهد بود. این طول به عدد رینولدز (Re) و عدد پرانتل (Pr) و نوع شرایط مرزی تحمیل شده به دیواره کانال بستگی دارد.

شکل 2- 2 جریان در حال توسعه هیدرودینامیکی و حرارتیرابطه تجربی برای طول ورودی کانال به شکل زیر می باشد:
LtDh=0.05 Re Pr (2.11)
فصل سومتشریح و حل پروژه
3-1پیش گفتار4381502651760X
Z
Y
00X
Z
Y
بررسی جریان های لغزشی امروزه بسیار مورد توجه قرار گرفته است [8-18]. همچنین مطالعات متعددی در مورد تولید انتروپی در میکرو کانال ها انجام شده است[19-23]. در این پروژه تلاش می شود تا جریان یک سیال گازی تراکم پذیر نیوتونی در یک میکروکانال مستطیل در حالت پایدار بررسی شود و سپس به بررسی تولید انتروپی در این کانال خواهیم پرداخت. در ابتدا معادلات حاکم بر جریان و معادلات حاکم بر انتقال حرارت ر بررسی می کنیم. از آنجایی که در نهایت می خواهیم از این معادلات در حل عددی استفاده کنیم لذا از صورتی از معادلات استفاده می کنیم که مناسب برای حل های عددی باشند و به راحتی گسسته سازی بر روی آنها انجام شود.
معادلات حاکم:
معادلات حاکم شامل معادلات حاکم بر جریان و معادلات حاکم بر انتقال حرارت و همچنین معادلات حاکم بر حالت سیال می باشند. معادلات حاکم بر جریان شامل معادلات پیوستگی و معادلات مومنتم می شوند که به قرار زیرند:
معادله پیوستگی:
divρu=0معادله مومنتم:
x-momentum divρuu=divμ g-- u-∂P∂xy-momentum divρvu=divμ g-- v-∂P∂yw-momentum divρwu=divμ g-- w-∂P∂zمعادلات حاکم بر انتقال حرارت شامل معادلات انرژی و معادله حالت میشوند. این معادلات نیز در فرم مناسب مورد استفاده به قرار زیرند:
divρh0u=divk g-- T+[∂uτxx∂x+∂uτyx∂y+∂uτzx∂z+∂vτxy∂x+∂vτyy∂y+∂vτzy∂z+∂wτxz∂x+∂wτyz∂y+∂(wτzz)∂z]P=ρRTو معادلات جانبی نیز برای بسته شدن دستگاه معادلات مورد نیاز هستند که به شرح زیرند:
h0=h+12(u2+v2+w2)h=cpTبا کمی دقت در معادلات مومنتم و انرژی در می یابیم که تمام این معادلات از دو بخش اساسی کانوکشن، عبارت سمت چپ تساوی و دیفیوژن، عبارت اول سمت راست تساوی تشکیل شده اند. به همین سبب می توان صورت کلی از معادلات شامل عبارت های کانوکشن و دیفیوژن را به صورت زیر نوشت:
divρϕu=divΠ g--ϕکه ضریب Π به سادگی بر اساس نوع معادله قابل تعیین است، لذا برای حل عددی ابتدا می بایست معادله بالا را که یک معادله کلی است برای سه معادله مومنتم و معادله انرژی است گسسته کنیم. به دلیل حجم بالای مطالب در این پروژه تنها گسسته شده معادلات آورده می شود و برای نحوه گسسته سازی و روش های آن خواننده به مراجع[30و31] ارجاع داده می شود.
3-2 گسسته سازی معادلات مومنتمبر اساس روش ارائه شده در مرجع [30] معادله گسسته شده به صورت زیر ارائه می شود و نکته قابل توجه در این گسسته سازی استفاده از فلاکس عبوری جریان است.
FeΦe-FwΦw+FnΦn-FsΦs+FtΦt-FbΦb=DeΦE-ΦP-DwΦP-ΦW+DnΦN-ΦP-DsΦP-ΦS+DtΦT-ΦP-DbΦP-ΦBمقادیر ϕ اگر با زیر نویس بزرگ باشد در آن نقاط مقادیر دقیق موجود می باشد ولی با زیر نویس حروف کوچک نیاز به تقریب داریم. از آنجایی که این تقریب تاثیر بسیار زیادی در همگرایی حل عددی دارد می بایست به دقت مورد توجه قرار گیرد. برای تقریب این نقاط در این پروژه از روش TVD استفاده کردیم. با این روش تقریب معادله گسسته شده به فرم زیر تبدیل می شود.
FeΦP+12ψre+ΦE-ΦPFeΦE+12ψre-ΦP-ΦE-FwΦW+12ψrw+ΦP-ΦWFwΦP+12ψrw-ΦW-ΦP+FnΦP+12ψrn+ΦN-ΦPFnΦN+12ψrn-ΦP-ΦN-FsΦS+12ψrs+ΦP-ΦSFsΦP+12ψrs-ΦS-ΦP+FtΦP+12ψrt+ΦT-ΦPFtΦT+12ψrt-ΦP-ΦT-FbΦB+12ψrb+ΦP-ΦBFbΦP+12ψrb-ΦB-ΦP=DeΦE-ΦP-DwΦP-ΦW+DnΦN-ΦP-DsΦP-ΦS+DtΦT-ΦP-DbΦP-ΦBدر اعمال تقریب TVD ردیف اول مربوط به فلاکس مثبت و ردیف دوم مربوط به فلاکس منفی می باشد. بعد از مرتب سازی متغییرها و ضرایب به معادله زیر خواهیم رسید:
aPϕP=aWϕW+aEϕE+aSϕS+aNϕN+aBϕB+aTϕT+SDCدر این معادله ضرایب به صورت زیر تعریف شده است:
aW=Dw+max⁡(Fw , 0)aE=De+max⁡(-Fe , 0)aS=Ds+max⁡(Fs , 0)aN=Dn+max⁡(-Fn , 0)aB=Db+max⁡(Fb , 0)aT=Dt+max⁡(-Ft , 0)و همچنین
aP=aW+aE+aS+aN+aB+aT+Fe-Fw+Fn-Fs+(Ft-Fb)و با توجه به تقریب استفاده شده برای نقاطی که در آنها مقدار قطعی وجود ندارد عبارت SDC نیز به این شکل تعریف می شود:
SDC=12Fe1-αeψ(re-)-αeψ(re+)(ϕE-ϕP)+12Fwαwψ(rw+)-1-αwψ(rw-)(ϕP-ϕW)+12Fn1-αnψ(rn-)-αnψ(rn+)(ϕN-ϕP)+12Fsαsψ(rs+)-1-αsψ(rs-)(ϕP-ϕS)+12Ft1-αtψ(rt-)-αtψ(rt+)(ϕT-ϕP)+12Fbαbψrb+-1-αbψrb-ϕP-ϕBقابل ذکر است که در عبارت بالا ψr تابع ون آلبادا در این تقریب است و به صورت زیر تعریف می شود:
ψr=r+r21+r2و همچنین داریم:
Fw>0 → αw=1 Fe>0 → αe=1Fw<0 → αw=0 Fe<0 → αe=0Fs>0 → αs=1 Fn>0 → αn=1Fs<0 → αs=0 Fn<0 → αn=0Fb>0 → αb=1 Ft>0 → αt=1Fb<0 → αb=0 Ft<0 → αt=0برای r داریم:
re+=ϕP-ϕWϕE-ϕP re-=ϕE-ϕEEϕP-ϕErw+=ϕW-ϕWWϕP-ϕW rw-=ϕP-ϕEϕW-ϕPrn+=ϕP-ϕSϕN-ϕP rn-=ϕN-ϕNNϕP-ϕNrs+=ϕS-ϕSSϕP-ϕS rs-=ϕP-ϕNϕS-ϕPrt+=ϕP-ϕBϕT-ϕP rt-=ϕT-ϕTTϕP-ϕTrb+=ϕB-ϕBBϕP-ϕB rb-=ϕP-ϕTϕB-ϕPبا توجه به اینکه گسسته سازی بر اساس حجم کنترل انتخاب شده برای مولفه های سرعت انجام می شود و در حل عددی به روش حجم محدود حجم های کنترل برای هر مولفه سرعت متفاوت است لذا با وجود اینکه معادله اصلی گسسته شده برای سه جهت مومنتم شکل یکسانی دارند ولی در تعاریف فلاکس جریان و در نتیجه ضرایب معادله بسیار با هم متفاوتند. به همین سبب برای روشن شدن روند محاسبات تمام ضرایب معادلات مومنتم به تفصیل آورده می شود. قابل ذکر است که از این به بعد از بالا نویس پریم برای مومنتم در جهت y و از بالا نویس زگوند برای مومنتم در جهت z استفاده می کنیم. همچنین از فرم کلی که تا کنون آن را تعریف کردیم برای مومنتم در جهت x استفاده می شود.
با توجه به توضیحات ارائه شده ابتدا فلاکس جربان را در معادله مومنتم در جهت x ارائه می کنیم.
Fe=12Fi+1,J,Z+Fi,J,Z=12(ρI+1,J,Z+ρI,J,Z2)ui+1,J,Z+(ρI,J,Z+ρI-1,J,Z2)ui,J,Z)Fw=12Fi,J,Z+Fi-1,J,Z=12(ρI,J,Z+ρI-1,J,Z2)ui,J,Z+(ρI-1,J,Z+ρI-2,J,Z2)ui-1,J,Z)Fn=12FI,j+1,Z+FI-1,j+1,Z=12(ρI,J,Z+ρI,J+1,Z2)vI,j+1,Z+(ρI-1,J,Z+ρI-1,J+1,Z2)vI-1,j+1,Z)Fs=12FI,j,Z+FI-1,j,Z=12(ρI,J,Z+ρI,J-1,Z2)vI,j,Z+(ρI-1,J-1,Z+ρI-1,J,Z2)vI-1,j,Z)Ft=12FI,J,z+1+FI-1,J,z+1=12(ρI,J,Z+ρI,J,Z+12)wI,J,z+1+(ρI-1,J,Z+ρI-1,J,Z+12)wI-1,J,z+1)Fb=12FI,J,z+FI-1,J,z=12(ρI,J,Z-1+ρI,J,Z2)wI,J,z+(ρI-1,J,Z-1+ρI-1,J,Z2)wI-1,J,z)دقت شود این فلاکس ها بر روی حجم کنترل u تعریف می شوند بنابر این زیر نویس های e,w,n,s,t,b به ترتیب بیانگر فلاکس های ورودی از پایین، بالا، جنوب، شمال، غرب و شرق حجم کنترل می باشند. فلاکس های ورودی بر روی حجم کنترل v به صورت زیر می باشند:
Fe'=12Fi+1,J,Z'+Fi+1,J-1,Z'=12(ρI+1,J,Z+ρI,J,Z2)ui+1,J,Z+(ρI,J-1,Z+ρI+1,J-1,Z2)ui+1,J-1,Z)Fs'=12FI,j-1,Z'+FI,j,Z'=12(ρI,J-1,Z+ρI,J-2,Z2)vI,j-1,Z+(ρI,J,Z+ρI,J-1,Z2)vI,j,Z)Fw'=12Fi,J,Z'+Fi,J-1,Z'=12(ρI,J,Z+ρI-1,J,Z2)ui,J,Z+(ρI-1,J-1,Z+ρI,J-1,Z2)ui,J-1,Z)Fn'=12FI,j,Z'+FI,j+1,Z'=12(ρI,J,Z+ρI,J-1,Z2)vI,j,Z+(ρI,J+1,Z+ρI,J,Z2)vI,j+1,Z)Ft'=12FI,J,z+1'+FI,J-1,z+1'=12(ρI,J,Z+1+ρI,J,Z2)wI,J,z+1+(ρI,J-1,Z+ρI,J-1,Z+12)wI,J-1,z+1)Fb'=12FI,J,z'+FI,J-1,z'=12(ρI,J,Z+ρI,J,Z-12)wI,J,z+(ρI,J-1,Z-1+ρI,J-1,Z2)wI,J-1,z)و همچنین فلاکس ورودی برای حجم کنترل w به شکل زیر خواهد بود:
Fe''=12Fi+1,J,Z''+Fi+1,J,Z-1''=12(ρI+1,J,Z+ρI,J,Z2)ui+1,J,Z+(ρI,J,Z-1+ρI+1,J,Z-12)ui+1,J,Z-1)Fs''=12FI,j,Z''+FI,j,Z-1''=12(ρI,J,Z+ρI,J-1,Z2)vI,j,Z+(ρI,J,Z-1+ρI,J-1,Z-12)vI,j,Z-1)Fw''=12Fi,J,Z''+Fi,J,Z-1''=12(ρI,J,Z+ρI-1,J,Z2)ui,J,Z+(ρI-1,J,Z-1+ρI,J,Z-12)ui,J,Z-1)Fn''=12FI,j+1,Z''+FI,j+1,Z-1''=12(ρI,J+1,Z+ρI,J,Z2)vI,j+1,Z+(ρI,J+1,Z-1+ρI,J,Z-12)vI,j+1,Z-1)Ft''=12FI,J,z+1''+FI,J,z''=12(ρI,J,Z+1+ρI,J,Z2)wI,J,z+1+(ρI,J,Z+ρI,J,Z-12)wI,J,z)Fb''=12FI,J,z''+FI,J,z-1''=12(ρI,J,Z+ρI,J,Z-12)wI,J,z+(ρI,J,Z-1+ρI,J,Z-22)wI,J,z-1)اکنون می توان معادلات مومنتم در سه جهت را با توجه به تعریف فلاکس های جریان در هر حجم کنترل به صورت گسسته نوشت. تعاریف ضرایب به صورت کلی پیش از این ارئه شده است و تنها با جایگذاری فلاکس جریان مربوطه می توان به تعاریف مناسب ضرایب برای هر جهت از حجم کنترل دست یافت.
با تعاریف بالا برای حجم کنترل در سه جهت داریم:
aPuP=aWuW+aEuE+aSuS+aNuN+aBuB+aTuT+SuDC+(Pw-Pe)aP'vP=aW'vW+aE'vE+aS'vS+aN'vN+aB'vB+aT'vT+SvDC+(Ps-Pn)aP''wP=aW''wW+aE''wE+aS''wS+aN''wN+aB''wB+aT''wT+SwDC+(Pb-Pt)می توان معادلات مومنتم را به شکل زیر نوشت که در آن nb معرف نقاط همسایگی نقطه ای است که بر روی حجم کنترل در آن نقطه معادله ممنتم به صورت گسسته در آمده است.
ai,J,Z ui,J,Z=anbunb+bi,J,Z+(PI-1,J,Z-PI,J,Z) aI,j,Z'vI,j,Z=anb'vnb+bI,j,Z'+(PI,J-1,Z-PI,J,Z)aI,J,z''wI,J,z=anb''wnb+bI,J,z''+(PI,J,Z-1-PI,J,Z)با کمی تغییر در معادلات بالا آن را به شکل مطلوب تری در می آوریم تا در مراحل بعدی محاسبات از آن استفاده کنیم.
ui,J,Z=ui,J,Z+di,J,Z(PI-1,J,Z-PI,J,Z)vI,j,Z=vI,j,Z+di,J,Z'(PI,J-1,Z-PI,J,Z)wI,J,z=wI,J,z+dI,J,z''(PI,J,Z-1-PI,J,Z)نکته قابل ذکر در معادلات بالا وجود دو دسته متغییر جدید است که آنها را به صورت زیر تعریف می کنیم:
ui,J,Z=anbunb+bi,J,Zai,J,Z di,J,Z=1ai,J,ZvI,j,Z=anb'vnb+bI,j,Z'aI,j,Z' di,J,Z'=1 aI,j,Z'wI,J,z=anb''wnb+bI,J,z''aI,J,z'' dI,J,z''=1aI,J,z''همچنین دو دسته از متغییرهای دیگر را به صورت زیر تعریف می کنیم:
ui,J,Z=u*i,J,Z+u'i,J,ZvI,j,Z=v*I,j,Z+v'I,j,ZwI,J,z=w*I,J,z+w'I,J,zکه در این تعریف مولفه های نشان داده شده با ستاره متغییر های کمکی هستند که در حل عددی در نهایت به مولفه های اصلی سرعت همگرا می شوند. همچنین مولفه های سرعت نشان داده شده با پرایم به عنوان تصحیح برای هر مرحله از حل عددی استفاده می شوند. از این رو در نهایت این متغییر ها به صفر همگرا خواهند شد. با توجه به این تعریف و با کم کردن معادلات مومنتم اصلی و معادلات مومنتم ستاره دار از هم و با یک سری از ساده سازی ها و تقریب ها که در مراجع به تفصیل آورده شده است به رابطه زیر خواهیم رسید که در حل های عددی بسیار به روند حل کمک خواهد کرد.
ui,J,Z=u*i,J,Z+di,J,Z(P'I-1,J,Z-P'I,J,Z)vI,j,Z=v*I,j,Z+di,J,Z'(P'I,J-1,Z-P'I,J,Z)wI,J,z=w*I,J,z+dI,J,z''(P'I,J,Z-1-P'I,J,Z)و به سادگی قابل درک است که:
u'i,J,Z=di,J,Z(P'I-1,J,Z-P'I,J,Z)v'I,j,Z=di,J,Z'(P'I,J-1,Z-P'I,J,Z)w'I,J,z=dI,J,z''(P'I,J,Z-1-P'I,J,Z)3-3 گسسته سازی معادله پیوستگیمعادله پیوستگی نقش بسیار پر رنگی در حل عددی معادله مومنتم ایفا می کند. در حقیقت از آنجاییکه در معادلات مومنتم سرعت و فشار با هم کوپل شده اند لذا نیاز به یک معادله کمکی برای حل فشار می باشد. از معادله پیوستگی به همین منظور استفاده می شود.
برای گسسته سازی معادله پیوستگی می بایست به این نکته دقت شود که حجم کنترل برای بررسی پیوستگی باید بر روی نود ها قرار بگیرد. در حقیقت حجم کنترل برای یک متغییر اسکالر بررسی می شود.
ρue-ρuw+ρvn-ρvs+ρwt-ρwb=0در نتیجه داریم:
ρI,J,Z+ρI+1,J,Z2ρi+1ui+1,J,Z-ρI-1,J,Z+ρI,J,Z2ρiui,J,Z+ρI,J+1,Z+ρI,J,Z2ρj+1vI,j+1,Z-ρI,J-1,Z+ρI,J,Z2ρjvI,j,Z+ρI,J,Z+1+ρI,J,Z2ρz+1wI,J,z+1-ρI,J,Z-1+ρI,J,Z2ρzwI,J,z=0با جایگذاری معادله () در معادله گسسته شده پیوستگی خواهیم داشت:
ρi+1ui+1,J,Z+di+1,J,Z(PI,J,Z-PI+1,J,Z)-ρiui,J,Z+di,J,Z(PI-1,J,Z-PI,J,Z)+ρj+1vI,j+1,Z+dI,j+1,Z'(PI,J,Z-PI,J+1,Z)-ρjvI,j,Z+dI,j,Z'(PI,J-1,Z-PI,J,Z)+ρz+1wI,J,z+1+dI,J,z+1''(PI,J,Z-PI,J,Z+1)-ρzwI,J,z+dI,J,z''(PI,J,Z-1-PI,J,Z)=0از آنجاییکه هدف ما از استفاده از معادله پیوستگی حل توزیع فشار است بنابراین می بایست معادله بالا را بر اساس فشار مرتب کنیم. معادله مورد نظر به این شکل خواهد بود:
ρi+1di+1,J,Z+ρidi,J,Z+ρj+1dI,j+1,Z'+ρjdI,j,Z'+ρz+1dI,J,z+1''+ρzdI,J,z''PI,J,Z=ρidi,J,ZPI-1,J,Z+ρi+1di+1,J,ZPI+1,J,Z+ρjdI,j,Z'PI,J-1,Z+ρj+1dI,j+1,Z'PI,J+1,Z+ρzdI,J,z''PI,J,Z-1+ρz+1dI,J,z+1''PI,J,Z+1+[ρiui,J,Z-ρi+1ui+1,J,Z+ρjvI,j,Z-ρj+1vI,j+1,Z+(ρzwI,J,z-ρz+1wI,J,z+1)]با تعریف متغییر های جدید برای ضرایب معادله بالا به صورت زیر
λI,J,Z=ρi+1di+1,J,Z+ρidi,J,Z+ρj+1dI,j+1,Z'+ρjdI,j,Z'+ρz+1dI,J,z+1''+ρzdI,J,z''λI-1,J,Z=ρidi,J,ZλI+1,J,Z=ρi+1di+1,J,ZλI,J-1,Z=ρjdI,j,Z'λI,J+1,Z=ρj+1dI,j+1,Z'λI,J,Z-1=ρzdI,J,z''λI,J,Z+1=ρz+1dI,J,z+1''bPI,J,Z=ρiui,J,Z-ρi+1ui+1,J,Z+ρjvI,j,Z-ρj+1vI,j+1,Z+ρzwI,J,z-ρz+1wI,J,z+1معادله مورد نیاز برای فشار به صورت زیر حاصل می شود.
λI,J,ZPI,J,Z=λI-1,J,ZPI-1,J,Z+λI+1,J,ZPI+1,J,Z+λI,J-1,ZPI,J-1,Z+λI,J+1,ZPI,J+1,Z+λI,J,Z-1PI,J,Z-1+λI,J,Z+1PI,J,Z+1+bPI,J,Zبه همین ترتیب با جایگذاری معادله () در معادله پیوستگی معادله مشابهی برای P' بدست می آید. معادله P' نیز همانند معادله P است با این تفاوت که در bP'I,J,Z بجای u ، v ، w از u* ، v* ، w* استفاده شده است. باقی ضرایب با معادله P مشابه است.
3-4 شرایط مرزیاعمال شرایط مرزی یکی از اساسی ترین قسمت های حل های عددی است. با اعمال شرط مرزی درست و صحیح است که مسئله فیزیکی تعریف صحیح خود را پیدا می کند و لذا در لحاظ کردن شرط های مزری دقت بسیار زیادی نیاز است.
در مسئله پیش رو از آنجا که با سه بعد مواجهیم نکات ظریفی در شرط های مرزی وجود دارد که می بایست به دقت مد نظر قرار گیرد که به تدریج به آنها اشاره می کنیم. نکته دیگر در اعمال شرط مرزی slip بر روی دیواره است. برای اعمال این شرط مرزی نیز از روش خاصی استفاده شده که هم این شرط مرزی را ارضاع کند و هم به نحوه مناسبی قابل استفاده در معادلات گسسته شده باشد.
3-4-1 اعمال شرط مرزی slip :شرط مرزی slip به صورت زیر تعریف می شود :
usurf-uwall=2-σvσvλ∂u∂nکه در آن usurf نشان دهنده سرعت سیال بر روی سطح و uwall نشان دهنده سرعت سطح می باشد که در مسئله مطرح شده برابر صفر است. σv در اکثر محاسبات مهندسی مشابه برابر واحد در نظر گرفته شده و همچنین λ بیانگر میانگین طول آزاد می باشد.
با تعریف
Kn=λDhداریم
usurf-uwall=Dh Kn ∂u∂nبا تقریب مناسبی برای ∂u∂n داریم:
∂u∂n=uP-usurfδ2که uP نمایان گر نقطه مرکزی حجم کنترل نزدیک به دیوار برای مومنتم در جهت x است. با جایگذاری معادله () در معادله () داریم:
usurf-uwall=2 Dh Kn δ×(uP-usurf)با مرتب سازی برای usurf و uP-usurf داریم :
usurf=2 Dh Kn δ1+2 Dh Kn δ uP+11+2 Dh Kn δ uwalluP-usurf=11+2 Dh Kn δ(uP-uwall)با تعریف
βv=2 Dh Kn δخواهیم داشت
usurf=βv1+βv uP+11+βv uwalluP-usurf=11+βv(uP-uwall)توضیح و تعاریف فوق برای مولفه های دیگر سرعت نیز کاربرد دارد، به همین سبب به شرح نحوه اعمال شرط مرزی برای حجم کنترل های سرعت می پردازیم.
3-4-2 شرط های مرزی برای u
با توجه به نوع شبکه بندی تعریف شده برای مسئله بعد طول به m قسمت تقسیم شده است. همچنین به علت وجود تقارن در سطح مقطع کانال تنها یک چهارم کانال حل می شود و لذا نیم بعد عرضی کانال به n قسمت و نیم بعد ارتفاع کانال به q قسمت تقسیم شده است. با این تعاریف برای مولفه سرعت در جهت x شروط مرزی را ارائه می کنیم.
شرایط مرزی در دو قسمت تقارن به صورت زیر تعیین می شود:
uJ=0=usouth=u1uZ=0=ubott=uZ=1برای اعمال شرط مرزی در دیواره با توجه به توضیحات ارائه شده برای شرط مرزی slip داریم:
uJ=n+1+uJ=n2=usurf⇒uJ=n+1+uJ=n=2βv1+βvuJ=nuJ=n+1=βv-1βv+1uJ=nبه همین ترتیب این رابطه برای دیواره دیگر نیز صادق است، بنابراین داریم:
uJ=n+1=unorth=-uJ=n slip uJ=n+1=βv-1βv+1uJ=nuZ=q+1=utop=-uZ=q slip uZ=q+1=βv-1βv+1uZ=qبرای شرط مرزی ورودی داریم:
ui=1=uwest=uinletui=0=uwestwest=2uinlet-ui=2برای شرط مرزی خروجی باید به این نکته دقت کرد که سرعت محاسبه شده در خروجی می بایست پیوستگی کل را ارضا کند به همین دلیل شرط مرزی خروجی به صورتی تعریف شده که شار ورودی جرم به کانال در ورودی به اصلاح سرعت خروجی بپردازد. با این توضیح داریم:
ui=m+1=ueast=ui=m×n×q×uinletJ=1nZ=1qui=m,J,Zui=m+2=ueasteast=2ueast-ui=m3-4-3 شرط های مرزی برای v
شرایط مرزی برای دو مرز تقارن به صورت زیر است:
vj=1=vsouth=0vj=0=vsouthsouth=-vj=2vZ=0=vbott=vZ=1با کمی دقت در تعریف شرط مرزی تقارن در مولفه سرعت v و مقایسه ی آن با شرط مرزی در u به این نکته پی می بریم که یک شرط مرزی اضافی برای v تعریف شده است. این تعریف به آن دلیل است که با توجه به استفاده از روش TVD و نوع تعریف توابع در این روش در یکی از مرز ها نیازمند تعریف شرط مرزی دیگری می باشیم.
شرط مرزی بعدی که باید مورد بررسی قرار گیرد شرط مرزی برای دیواره است. نکته مهم در شرط مرزی دیواره آنست که مولفه سرعت v بر یکی از دیوار ها عمود و با دیواره دیگر موازی است. به همین دلیل بر خلاف u که بر روی هر دو دیوار شرط مرزی slip برقرار است برای v و همچنین w که در ادامه به شرط های مرزی آن می پردازیم بر روی دیواری که سرعت موازی با آن دیوار است شرط مرزی slip برقرار است.
vj=n+1=vnorth=vwall=0vj=n+2=vnorthnorth=-vj=nvZ=q+1=vtop=-vZ=q slip vZ=q+1=βv-1βv+1vZ=qبرای شرط مرزی ورودی داریم:
vI=0=vwest=0vI=-1=vwestwest=-vI=1و برای شرط مرزی خروجی نیز از شرط گرادیان برابر صفر استفاده می کنیم. خواهیم داشت :
vI=m+1=veast=vI=mvI=m+2=veasteast=vI=m3-4-4 شرط های مرزی برای w
برای w نیز همانند v برای شرط مرزی تقارن داریم :
wJ=0=wsouth=wJ=1wz=1=wbott=wsym=0wz=0=wbottbott=-wz=2همچنین برای مرز با دیواره خواهیم داشت :
wJ=n+1=wnorth=-wJ=n slip wJ=n+1=βv-1βv+1wJ=nwz=q+1=wtop=wwall=0wz=q+2=wtoptop=-wz=qو همانند v برای مرز ورودی و خروجی داریم :
wI=0=wwest=0wI=-1=wwestwest=-wI=1wI=m+1=weast=wI=mwI=m+2=weasteast=wI=m3-5 اعمال شرط های مرزی در معادلات مومنتم
با توجه به وجود شرط های مرز ی مختلف در مسئله و معادلات گسسته شده مومنتم در هر جهت لازم است که تاثیر هر یک از مرز ها بر معادلات به دقت بررسی شود و ضرایب معادلات گسسته شده مجددا محاسبه گردد. به همین سبب بررسی کامل معادلات مومنتم در نزدیکی مرز ها و تغییرات اعمال شده برای هر ضریب به صورت کامل در ضمیمه ارائه شده است.

3-6 گسسته سازی معادله انرژیبا توجه به مطالب ارائه شده در بخش قبل معادله انرژی حاکم بر جریان بدین شکل خواهد بود.
∂∂xkρh+12ρujujuk=∂∂xiujτij-∂qj∂xjکه در آن داریم:
τij=λδij∂uk∂xk+μ∂ui∂xj+∂uj∂xiλ=-23μh=cPTهرچند می توان با ساده سازی شکل های ساده تری از معادله انرژی را برای این مسئله یافت اما باید به این نکته توجه شود که این صورت از معادله انرژی برای استفاده در حل عددی به روش حجم محدود حیاتی است. با توجه به معادله انرژی درمی یابیم این معادله شامل عبارت کانوکشن، دیفیوژن و اتلاف ویسکوز است. به همین دلیل به صورت مجزا به بررسی و گسسته سازی این عبارت ها می پردازیم.
در ابتدا برای عبارت کانوکشن داریم:
∂∂xkρhuk=∂∂xρcPTu+∂∂yρcPTv+∂∂zρcPTwلذا
∂∂xρcPTudV=∂∂xρcPTuAdx=[ρcPTue-(ρcPTu)w]Aو به همین ترتیب خواهیم داشت
∂∂yρcPTvdV=[ρcPTvn-(ρcPTv)s]A∂∂zρcPTwdV=[ρcPTwt-(ρcPTw)b]Aبرای عبارت دیفیوژن در معادله انرژی داریم
-∂qj∂xj=-∂∂x-k∂T∂x-∂∂y-k∂T∂y-∂∂z-k∂T∂zلذا
∂∂xk∂T∂xdV=∂∂xk∂T∂xAdx=k∂T∂xe-k∂T∂xwAو به همین ترتیب هواهیم داشت
∂∂yk∂T∂ydV=k∂T∂yn-k∂T∂ysA∂∂zk∂T∂zdV=k∂T∂zt-k∂T∂zbAبا گسسته سازی معادلات بالا داریم
ρecPTeue-ρwcPTwuw+ρncPTnun-ρscPTsus+ρtcPTtut-ρbcPTbub=kδTE-TP-kδTP-TW+kδTN-TP-kδTP-TS+kδTT-TP-kδTP-TB+Sدر معادله بالا S نمایان گر تمام عباراتی است مانند عبارت اتلاف ویسکوز که تا کنون در معادله انرژی به آنها نپرداخته ایم. با تغییر کوچکی در معادله بالا خواهیم داشت:
cPFeTTe-cPFwTTw+cPFnTTn-cPFsTTs+cPFtTTt-cPFbTTb=kδTE-TP-kδTP-TW+kδTN-TP-kδTP-TS+kδTT-TP-kδTP-TB+Sکه با دقت در شبکه بندی داریم
FTe=(ρI,J,Z+ρI+1,J,Z2)ui+1,J,ZFTw=(ρI-1,J,Z+ρI,J,Z2)ui,J,ZFTn=(ρI,J,Z+ρI,J+1,Z2)vI,j+1,ZFTs=(ρI,J-1,Z+ρI,J,Z2)vI,j,ZFTt=(ρI,J,Z+ρI+1,J,Z+12)wI,J,z+1FTb=(ρI,J,Z-1+ρI,J,Z2)wI,J,zهمچنین با توجه به اینکه شبکه بندی انتخاب شده یک شبکه بندی یکنواخت است می توان تعریف زیر را ارائه کرد:
Γe=Γw=Γn=Γs=Γt=Γb=Γ=kδcPکه معادله انرژی را به فرم زیر خواهیم داشت:
FeTTe-FwTTw+FnTTn-FsTTs+FtTTt-FbTTb=ΓeTE-TP-ΓwTP-TW+ΓnTN-TP-ΓsTP-TS+ΓtTT-TP-ΓbTP-TB+1cPSبا توجه به شبکه بندی برای مسئله بعد از گسسته سازی ترم کانوکشن با عبارت هایی مواجه می شویم که در آن نقاط مقدار قطعی برای دما موجود نمی باشد لذا نیاز مند استفاده از تقریب های مناسبی می باشیم که با شرایط مسئله سازگار باشد و به همگرایی حل کمک کند. مانند روش استفاده شده در معادلات ممنتم از روش TVD استفاده می کنیم.
با استفاده از تقریب TVD داریم:
aTPTP=aTWTW+aTETE+aTSTS+aTNTN+aTBTB+aTTTT+STDCکه در این معادله انرژی گسسته شده ضرایب عبارتند از:
aTW=Γw+max⁡(FTw , 0)aTE=Γe+max⁡(FTe , 0)aTS=Γs+max⁡(FTs , 0)aTN=Γn+max⁡(FTn , 0)aTB=Γb+max⁡(FTb , 0)aTT=Γt+max⁡(FTt , 0)و همچنین داریم :
aTP=aTW+aTE+aTS+aTN+aTB+aTT+FTe-FTw+FTn-FTs+(FTt-FTb)با توجه به استفاده از تقریب TVD داریم:
STDC=12FTe1-αTeψ(rTe-)-αTeψ(rTe+)(TE-TP)+
12FTwαTwψ(rTw+)-1-αTwψ(rTw-)(TP-TW)+
12FTn1-αTnψ(rTn-)-αTnψ(rTn+)(TN-TP)+
12FTsαTsψ(rTs+)-1-αTsψ(rTs-)(TP-TS)+
12FTt1-αTtψ(rTt-)-αTtψ(rTt+)(TT-TP)+
12FTbαTbψrTb+-1-αTbψrTb-TP-TB+1cpSدر رابطه بالا عبارت 1cpS نماینده اتلاف ویسکوز در معادله انرژی می باشد. در ادامه به تفصیل در مورد گسسته سازی عبارت های اتلاف ویسکوز بحث خواهد شد.
می بایست دقت شود که در معادله بالا داریم :
FTw>0 → αTw=1 FTe>0 → αTe=1FTw<0 → αTw=0 FTe<0 → αTe=0FTs>0 → αTs=1 FTn>0 → αTn=1FTs<0 → αTs=0 FTn<0 → αTn=0FTb>0 → αTb=1 FTt>0 → αTt=1FTb<0 → αTb=0 FTt<0 → αTt=0همچنین داریم:
rTe+=TP-TWTE-TP rTe-=TE-TEETP-TErTw+=TW-TWWTP-TW rTw-=TP-TETW-TPrTn+=TP-TSTN-TP rTn-=TN-TNNTP-TNrTs+=TS-TSSTP-TS rTs-=TP-TNTS-TPrTt+=TP-TBTT-TP rTt-=TT-TTTTP-TTrTb+=TB-TBBTP-TB rTb-=TP-TTTB-TPحال به بررسی عبارت اتلاف ویسکوز در معادله انرژی می پردازیم. برای S داریم:
S=∂∂xuτxx+∂∂xvτxy+∂∂xwτxz+∂∂yuτyx+∂∂yvτyy+∂∂ywτyz+∂∂zuτzx+∂∂zvτzy+∂∂zwτzz-∂∂xu12ρu2+12ρv2+12ρw2-∂∂yv12ρu2+12ρv2+12ρw2-∂∂zw12ρu2+12ρv2+12ρw2به دلیل تعداد زیاد عبارات به بررسی آنها به صورت مجزا می پردازیم. در عبارت اول در ترم اتلاف ویسکوز داریم:
τxx=2μ∂u∂x-23μ∂u∂x+∂v∂y+∂w∂zبا مرتب سازی رابطه بالا می توان آن را به شکل زیر نوشت :
τxx=23μ2∂u∂x-∂v∂y+∂w∂zبا انتگرال گیری از عبارت اول در اتلاف ویسکوز بر روی حجم کنترل دما در شبکه بندی ارائه شده داریم:
∂∂xuτxxAdx=uτxxe-uτxxw=23μu2∂u∂x-∂v∂y+∂w∂ze-23μu2∂u∂x-∂v∂y+∂w∂zwبرای محاسبه عبارت بالا می بایست تک تک اجزا آن را با بررسی شبکه بندی به دقت محاسبه کرد. با توجه به شبکه بندی و اینکه در حجم کنترل دما این عبارت ها می بایست تعیین شوند داریم:
ue=ui+1,J,Z∂u∂xe=ui+2,J,Z-ui,J,Z2δ∂v∂ye=vI,j+1,Z+vI+1,j+1,Z-(vI,j,Z+vI+1,j,Z)2δ∂w∂ze=(wI,J,z+1+wI+1,J,z+1)-(wI,J,z+wI+1,J,z)2δو همچنین برای کروشه دوم داریم:
uw=ui,J,Z∂u∂xw=ui+1,J,Z-ui-1,J,Z2δ∂v∂yw=vI-1,j+1,Z+vI,j+1,Z-(vI-1,j,Z+vI,j,Z)2δ∂w∂zw=(wI-1,J,z+1+wI,J,z+1)-(wI-1,J,z+wI,J,z)2δدر عبارت دوم از اتلاف ویسکوز داریم :
τxy=μ∂u∂y+∂v∂xبا انتگرال گیری بر حجم کنترل دما بدست می آید :
∂∂xvτxyAdx=vτxye-vτxyw=μv∂u∂y+∂v∂xe-μv∂u∂y+∂v∂xwو خواهیم داشت:
ve=vI,j,Z+vI,j+1,Z+vI+1,j,Z+vI+1,j+1,Z4∂u∂ye=ui+1,J+1,Z-ui+1,J-1,Z2δ∂v∂xe=vI+1,j+1,Z+vI+1,j,Z-(vI,j+1,Z+vI,j,Z)2δهمچنین برای کروشه دوم داریم:
vw=vI-1,j,Z+vI-1,j+1,Z+vI,j,Z+vI,j+1,Z4∂u∂yw=ui,J+1,Z-ui,J-1,Z2δ∂v∂xw=vI,j+1,Z+vI,j,Z-vI-1,j+1,Z+vI-1,j,Z2δدر عبارت سوم از اتلاف ویسکوز دازیم:
τxz=μ∂u∂z+∂w∂xو با انتگرال گیری بر روی حجم کنترل خواهیم داشت:
∂∂xwτxzAdx=wτxze-wτxzw=μw∂u∂z+∂w∂xe-μw∂u∂z+∂w∂xwبرای عبارات در کروشه اول داریم:
we=wI,J,z+wI,J,z+1+wI+1,J,z+wI+1,J,z+14(∂u∂z)e=ui+1,J,Z+1-ui+1,J,Z-12δ(∂w∂x)e=(wI+1,J,z+1+wI+1,J,z)-(wI,J,z+1+wI,J,z)2δو در کروشه دوم:
ww=wI-1,J,z+wI-1,J,z+1+wI,J,z+wI,J,z+14(∂u∂z)w=ui,J,Z+1-ui,J,Z-12δ(∂w∂x)w=(wI,J,z+1+wI,J,z)-(wI-1,J,z+1+wI-1,J,z)2δعبارت چهارم در اتلاف ویسکوز در بردارنده عبارت زیر می باشد:
τyx=μ∂u∂y+∂v∂xکه با انتگرال گیری داریم:
∂∂yuτyxAdy=uτyxn-uτyxs=μu∂u∂y+∂v∂xn-μu∂u∂y+∂v∂xsحال برای عبارت های داخل کروشه خواهیم داشت:
un=ui,J,Z+ui+1,J,Z+ui,J+1,Z+ui+1,J+1,Z4(∂u∂y)n=(ui,J+1,Z+ui+1,J+1,Z)-(ui,J,Z+ui+1,J,Z)2δ(∂v∂x)n=vI+1,j+1,Z-vI-1,j+1,Z2δو
us=ui,J-1,Z+ui+1,J-1,Z+ui,J,Z+ui+1,J,Z4(∂u∂y)s=(ui,J,Z+ui+1,J,Z)-(ui,J-1,Z-ui+1,J-1,Z)2δ(∂v∂x)s=vI+1,j,Z+vI-1,j,Z2δبرای عبارت پنجم داریم:
τyy=23μ2∂v∂y-(∂u∂x+∂w∂z)که با انتگرال گیری داریم:
∂∂yvτyyAdy=vτyyn-vτyys=23μv2∂v∂y-∂u∂x+∂w∂zn-23μv2∂v∂y-∂u∂x+∂w∂zsو عبارت های درون کروشه ها اینچنینند:
vn=vI,j+1,Z(∂v∂y)n=vI,j+2,Z-vI,j,Z2δ∂u∂xn=ui+1,J,Z+ui+1,J+1,Z-ui,J,Z+ui,J+1,Z2δ(∂w∂z)n=wI,J,z+1+wI,J+1,z+1-(wI,J,z+wI,J+1,z)2δو همپنین
vs=vI,j,Z(∂v∂y)s=vI,j+1,Z-vI,j-1,Z2δ∂u∂xs=ui+1,J-1,Z+ui+1,J,Z-ui,J-1,Z+ui,J,Z2δ(∂w∂z)s=wI,J-1,z+1+wI,J,z+1-(wI,J-1,z+wI,J,z)2δدر عبارت ششم از اتلاف ویسکوز داریم:
τyz=μ∂v∂z+∂w∂yو در تنیجه خواهیم داشت
∂∂ywτyzAdy=wτyzn-wτyzs=μw∂v∂z+∂w∂yn-μw∂v∂z+∂w∂ysلذا داریم:
wn=wI,J,z+wI,J,z+1+wI,J+1,z+wI,J+1,z+14(∂v∂z)n=vI,j+1,Z+1-vI,j+1,Z-12δ(∂w∂y)n=wI,J+1,z+1+wI,J+1,z-(wI,J,z+1+wI,J,z)2δو
ws=wI,J-1,z+wI,J-1,z+1+wI,J,z+wI,J,z+14(∂v∂z)s=vI,j,Z+1-vI,j,Z-12δ(∂w∂y)s=wI,J,z+1+wI,J,z-(wI,J-1,z+1+wI,J-1,z)2δدر عبارت هفتم از اتلاف ویسکوز داریم :
τzx=μ∂u∂z+∂w∂xبا توجه به جهت انتگرال گیری داریم:
∂∂zuτzxAdz=uτzxt-uτzxb=μu∂u∂z+∂w∂xt-μu∂u∂z+∂w∂xbو در کروشه ها خواهیم داشت:
ut=ui,J,Z+ui+1,J,Z+ui,J,Z+1+ui+1,J,Z+14∂u∂zt=ui,J,Z+1+ui+1,J,Z+1-(ui,J,Z+ui+1,J,Z)2δ∂w∂xt=(wI+1,J,z+1-wI-1,J,z+1)2δو
ub=ui,J,Z-1+ui+1,J,Z-1+ui,J,Z+ui+1,J,Z4∂u∂zb=ui,J,Z+ui+1,J,Z-(ui,J,Z-1+ui+1,J,Z-1)2δ∂w∂xb=(wI+1,J,z-wI-1,J,z)2δعبارت هشتم و انتگرال گیری آن به قرار زیر است:
τzy=μ∂v∂z+∂w∂y∂∂zvτzyAdz=vτzyt-vτzyb=μv∂v∂z+∂w∂yt-μv∂v∂z+∂w∂ybو در نیجه برای کروشه ها داریم:
vt=vI,j,Z+vI,j+1,Z+vI,j,Z+1+vI,j+1,Z+14∂v∂zt=vI,j+1,Z+1+vI,j,Z+1-(vI,j+1,Z+vI,j,Z)2δ∂w∂yt=(wI,J+1,z+1-wI,J-1,z+1)2δو همچنین
vb=vI,j,Z-1+vI,j+1,Z-1+vI,j,Z+vI,j+1,Z4∂v∂zb=vI,j+1,Z+vI,j,Z-(vI,j+1,Z-1+vI,j,Z-1)2δ∂w∂yb=(wI,J+1,z-wI,J-1,z)2δو در نهایت برای عبارت نهم داریم:
τzz=23μ2∂w∂z-∂u∂x+∂v∂y∂∂zwτzzAdx=wτzzt-wτzzb=23μw2∂w∂z-∂u∂x+∂v∂yt-23μw2∂w∂z-∂u∂x+∂v∂ybکه در آن عبارت های درون کروشه ها به قرار زیرند:
wt=wI,J,z+1∂w∂zt=(wI,J,z+2-wI,J,z)2δ∂u∂xt=ui+1,J,Z+ui+1,J,Z+1-ui,J,Z+ui,J,Z+12δ∂v∂yt=vI,j+1,Z+vI,j+1,Z+1-(vI,j,Z+vI,j,Z+1)2δو همچنین
wb=wI,J,z∂w∂zb=(wI,J,z+1-wI,J,z-1)2δ∂u∂xb=ui+1,J,Z-1+ui+1,J,Z-ui,J,Z-1+ui,J,Z2δ∂v∂yb=vI,j+1,Z-1+vI,j+1,Z-(vI,j,Z-1+vI,j,Z)2δهمچنین برای سه عبارت پایانی داریم:
-∂∂xu12ρu2+12ρv2+12ρw2Adx=-u12ρu2+12ρv2+12ρw2e--u12ρu2+12ρv2+12ρw2w-∂∂yv12ρu2+12ρv2+12ρw2Ady=-v12ρu2+12ρv2+12ρw2n--v12ρu2+12ρv2+12ρw2s-∂∂zw12ρu2+12ρv2+12ρw2Adz=-w12ρu2+12ρv2+12ρw2t--w12ρu2+12ρv2+12ρw2bکه در رابطه بالا تمامی عبارت های درون کروشه ها در صفحات پیشین آورده شده است.
3-6-1 نحوه اعمال شرط مرزی پرش دما بر روی دیواره هابرای اعمال شرط مرزی پرش دما بر روی دیواره ها همانند شرط مرزی لغزش عمل می کنیم. در این مسئله برای شرط مرزی پرش از رابطه زیر استفاده شده است:
Tsurf-Twall=2-σtσtλPr2γ1+γ∂T∂nکه در بسیاری از محاسبات مهندسی σt=1 در نظر گرفته شده است . برای∂T∂n از تقریب زیر استفاده می کنیم و در ادامه خواهیم داشت:
∂T∂n=TP-Tsurfδ2Tsurf-Twall=2-σtσtDh Kn Pr2γ1+γ∂T∂nو در نتیجه
Tsurf-Twall=2δDh Kn Pr2γ1+γ×(TP-Tsurf)با در نظر گرفتن
βt=2δDh Kn Pr2γ1+γداریم
Tsurf=βt1+βt TP+11+βt TwallTP-Tsurf=11+βt(TP-Twall)3-6-2 شرایط مرزی برای دمابا توجه به رابطه بدست آمده در صفحه قبل و این نکته که برای دما بر روی دیوار داریم:
TP+Tnorth2=Tsurfمی توان به رابطه زیر برای نود خارج از محدوده مسئله برای دما دست یافت:
TP+Tnorth=2βt1+βt TP+21+βt TwallnorthTJ=n+1=Tnorth=βt-11+βt TP+21+βt Twallnorthهمچنین برای مرز تقارن داریم
TJ=0=Tsouth=TJ=1همانند روابط بالا برای مرز دیگر دیواره داریم:
TZ=q+1=Ttop=βt-11+βt TP+21+βt Twalltopو برای مرز تقارن هم ارز آن داریم:
TZ=0=Tbott=TZ=1برای ورودی و خروجی نیاز به دو شرط خارج از منطقه حل خواهیم بود که به قرار زیرند :
TI=0=Twest=TinletTI=-1=Twestwest=2Tinlet-TI=1TI=m+1=Teast=TI=mTI=m+2=Teasteast=TI=mبا توجه به توضیحات صفحات پیشین بدیهی است که این شرایط مرزی تاثیر زیادی بر معادله انرژی در نود های نزدیک مرز خواهد گذاشت لذا به بررسی این نود ها به صورت کاملا مشروح در ضمائم می پردازیم.

3-7 الگوریتم حلدر این پروژه از الگوریتم SIMPLER برای حل معادلات مومنتم و انرژی استفاده شده است. با توجه به توضیحات کامل موجود در منابع مختلف شرح مختصری از مراحل مختلف این الگوریتم ارائه می شود.
3-7-1مراحل SIMPLER
مرحله اول :
ui,J,Z=anbunb+bi,J,Zai,J,ZvI,j,Z=anb'vnb+bI,j,Z'aI,j,Z'wI,J,z=anb''wnb+bI,J,z''aI,J,z''مقادیر unb ، vnb وwnb از مقادیر u ، v و wبدست آمده از دور قبل محاسبات یا از مقادیر اولیه برای مولفه های سرعت استفاده می کنند.
ضرایب a'',b'',a',b',a,b نیز به همین ترتیب عمل میکنند.
مرحله دوم: حل معادله Pضرایب معادله P نیز از w,v,u های دور قبل استفاده می کنند. تنها bp از vوu محاسبه شده در مرحله اول بدست می آید.
مرحله سوم: حل معادلات مومنتم
ai,J,Zui,J,Z*=anbunb*+PI-1,J,Z-PI,J,Z+bi,J,ZaI,j,Z'vI,j,Z*=anb'vnb*+PI,J-1,Z-PI,J,Z+bI,j,Z'aI,j,Z''wI,j,Z*=anb''wnb*+PI,J,Z-1-PI,J,Z+bI,j,Z''در مرحله سوم با حل معادلات بالا v*,u* بدست می آیند.
در این معادلات از P های محاسبه شده در Step 2 استفاده می شود. تمام ضرایب از u های دور قبل بدست می آیند.

Related posts: