نوع مقاله : مقاله پژوهشی
نویسندگان
1 دانشجوی دکتری، دانشکده مهندسی عمران و محیطزیست، دانشگاه صنعتی امیرکبیر، تهران
2 ....
چکیده
در این مقاله، مدل لاگرانژی دوبعدی قائم با روش نیمه ضمنی ذره متحرک (MPS) جهت مدلسازی مسائل موج و برخورد با سازههای دریایی توسعه داده شده است. بهمنظور اعمال مدل به جریان ویسکوز، لازم است اثر تلاطم در نظر گرفته شود. با این هدف، از ویسکوزیته گردابی ثابت استفاده شده است. مدل عددی حاضر در چند مثال روند یابی موج در کانال و برخورد موج به دیواره به کار رفته است. مقایسه نتایج حاصل از مدل عددی حاضر با اطلاعات آزمایشگاهی و عددی موجود در ادبیات فنی نشاندهنده افزایش توانایی مدل در برآورد سطح آب بهخصوص در مقایسه با سایر مدلهای غیر ویسکوز است. همچنین میدان فشار برخلاف اکثر مدلهای لاگرانژی بهخوبی محاسبه شده است.
کلیدواژهها
1- مقدمه
جریانهای با سطح آزاد، بهخصوص در مناطق دریایی یا کانالهای موجدار به سبب نقش تجاری و اقتصادی خود، اهمیت هیدرولیکی خاصی دارند. شبیهسازی هیدرودینامیکی این مناطق به خاطر پیچیدگیهایی که در شرایط مرزی و سطح آزاد آنها وجود دارد، تا حدودی پیچیده میباشد. با وجود پیشرفتهای زیادی که در مدلسازی عددی سطوح آزاد آب به دست آمده است، همچنان دشواریهایی در شبیهسازی سطوح آزاد پیچیده یا اندرکنش سازه و آب وجود دارد. اخیرا روشهای لاگرانژی درزمینهٔ مدلسازی سطح آزاد بسیار موردتوجه قرار گرفتهاند. یکی از جدیدترین این روشها، روش نیمه ضمنی ذره متحرک (MPS) است که اولین بار توسط Koshizuka & Oka (1996) معرفی شد [1]. در گذشته روش MPS برای مدلسازی تعدادی از پدیدههای هیدرولیکی چون شکست سد، شکست موج، خروج جت و جریان روی سرریز بکار گرفته شده است [2، 3].
Gotoh and Sakai (2006) یک مدل چند فازه MPS را برای شبیهسازی مسائل با فازهای مایع و گاز یا مایع و جامد، انتقال رسوب و اجسام شناور توسعه دادند [4]. Ataie-Ashtiani and Farhadi (2006) توایع کرنل مختلف را با هم مقایسه کرده و رابطهای جهت افزایش پایداری مدل لاگرانژی ارائه دادند [5]. Shibata and Koshizuka (2007) مدل MPS سهبعدی را جهت شبیهسازی برخورد موج به عرشه کشتی و پیشبینی فشار ناشی از برخورد بکار بردند [6]. Suzuki et al. (2007) روش MPS را با روش هامیلتونی ترکیب کرده و روش جدیدی پیشنهاد نمود [7]. فیاض (1386) متغیرهای دیگری مانند تنش سطحی به مدل چند فازه MPS افزود و پایداری و دقت مدل را ارتقا بخشید [8]. Khayyer and Gotoh (2009) روی بخش بقای ممنتوم مدل کار کردند و رابطه جدیدی برای تغییرات فشار پیشنهاد دادند [3]. آنها همچنین برای غلبه بر نوسانات فشار کمی تراکمپذیری را برای سیال قائل شدند. سپس Khayyer and Gotoh (2010) مدل مرتبه بالاتری برای پایدارسازی و ارتقای محاسبات فشار در مدل MPS معرفی کردند [9]. Shakibaeinia and Jin (2010) تراکمپذیری کمی را برای مدل MPS قائل شده و در شرایط مرزی، استراتژی جایگزینی ذرات را پیشنهاد نمودند [10]. Tanaka and Masunaga (2010) برای هموارسازی نوسانات فشار در زمان و مکان تلاشهایی انجام دادند [11]. Lee et al. (2010) مسائل برخورد بین دو فاز مایع یا مایع و جامد را با روش MPS مدلسازی کردند [12]. بهمنظور غلبه بر نوسانات غیرفیزیکی فشار، Kondo and Koshizuka (2010) رابطه جدیدی برای جمله منبع در معادله پواسن فشار پیشنهاد دادند [13]. بهطورکلی روش MPS تاکنون بیشتر برای مدلسازی جریانهای غیر ویسکوز بکار رفته و اثر آشفتگی در مدلهای چند فازه MPS دیده نشده است [5، 7]. فیاض (1386) تلاشهایی برای افزودن اثر آشفتگی به مدل چند فازه MPS با استفاده از مفهوم آشفتگی گردابی ثابت انجام داد [8]. همچنین شیرازپور (1388) اثر مدل آشفتگی طول اختلاط پرانتل و مدل دو معادلهای k-e را بر روی تغییرات سطح آزاد در مسئله شکست سد بررسی کرد [14]. در پژوهش حاضر، روش MPS برای شبیهسازی سطح آزاد در جریان ویسکوز توسعه داده شده است. به این صورت که از ویسکوزیته گردابی ثابت استفاده شده و اثر آشفتگی با استفاده از این مدل در مسائل موج دیده شده است.
2- روش تحقیق
2-1- معادلات حاکم
معادلات حاکم در جریان ویسکوز شامل معادلات پیوستگی و ممنتوم بهصورت زیر میباشد:
که در آن u بردار سرعت، t زمان، ρچگالی سیال، P فشار، νt ویسکوزیته گردابی سیال و f بردار شتاب گرانشی و تنش سطحی است.
در این روش برخلاف روشهای اولری سطح آزاد جریان بهطور خودکار برای ذرات به دست میآید. ازاینرو میتوان ترم فشار را بهصورت متغیر مستقل شامل مقادیر فشار استاتیکی و دینامیکی در نظر گرفت و با حل معادله پواسون فشار، مقدار آن را برای تمامی ذرات محاسبه کرد. معادله پواسون فشار بهکاررفته بهصورت رابطه (3) میباشد [1]:
2-2- منقطع سازی مدل
در روش MPS معادلات پیوستگی و ممنتوم به معادلات اندرکنش ذرات با اپراتورهای مختلف تبدیل میشوند. همه اندرکنشها در فاصله مشخصی به نام شعاع تأثیر محدود میگردند. وزن دهی به ذرات همسایه موجود در شعاع تأثیر یک ذره مشخص، با استفاده از تابع کرنل انجام میشود. در تحقیق حاضر از رابطهQuartic Spline بهعنوان تابع کرنل استفاده شده است.
که در آن r فاصله بین دو ذره، re شعاع تأثیر و w تابع کرنل است.
تابع گرادیان بین دو ذره همسایه با استفاده از تابع وزنی بهصورت زیر بیان میگردد [1]:
که در آن کمیت فیزیکی، d ابعاد مسئله (2 یا 3 برای مسائل دو یا سهبعدی)، n0 چگالی استاندارد که برای سیال غیرقابل تراکم ثابت میماند، ri بردار مکان ذره i و حداقل مقدار مربوط به ذرات همسایه در شعاع تأثیر است.
تابع لاپلاس نیز بهصورت زیر بیان میگردد [1]:
در معادلات بالا،ni چگالی ذره i در موقعیت ri بهصورت زیر به دست میآید:
و ضریبλ بهصورت زیر تعریف میگردد [1]:
2-3- روش حل عددی
روش MPS مورداستفاده در پژوهش حاضر بهگونهای است که معادلات حاکم بهصورت نیمه صریح حل میشوند. برای افزایش سرعت همگرایی مدل از روش پروجکشن استفاده شده است. در این روش معادلات ناویر استوکس در دو نیم گام زمانی حل میشوند. در نیم گام زمانی اول، معادلات اول با حضور جملات ویسکوزیته و ثقل حل میشوند درحالیکه از سایر جملات صرفنظر میشود. در این حالت موقعیت و سرعت هر ذره محاسبه میشود. سپس معادله پواسن فشار حل شده و مقادیر فشار برای کلیه ذرات به دست میآید. در نیم گام زمانی دوم مقادیر بهدستآمده از نیم گام زمانی اول با حضور سایر جملات معادله ناویر استوکس مانند فشار و کشش سطحی، تصحیح میگردد [13].
در ابتدای حل لازم است ویسکوزیته گردابی محاسبه گردد. پس از انتخاب ویسکوزیته آشفتگی برای هر ذره با تحلیل حساسیت، حل پروجکشن آغاز میشود. معادله نیم گام زمانی اول به بیان ریاضی بهصورت زیر خواهد بود:
پس از حل معادله بالا بهصورت صریح، مقادیر تغییرات سرعت برای هر ذره به دست میآید. با استفاده از آن میتوان موقعیت و سرعت تصحیحشده ذرات را نیز محاسبه کرد.
که در آن ut، rt، ut+1/2 و rt+1/2 به ترتیب سرعت و موقعیت هر ذره در گام زمانی فعلی t و نیم گام زمانی آینده t+1/2 است.
برای یافتن مقادیر فشار معادله پواسن بهصورت ضمنی حل میشود:
در نهایت با دانستن سرعت و موقعیت در نیم گام زمانی اول و مقادیر سرعت، گام دوم پروجکشن با حل معادله ناویر استوکس در حضور سایر جملات اجرا میشود.
این معادله بهطور ضمنی حل شده و مقادیر جدید سرعت و موقعیت برای هر ذره در گام زمانی آینده به دست میآید.
3- تجزیهوتحلیلدادههاوبیاننتایج
3-1- صحت سنجی و ارزیابی مدل
جهت روندیابی موج در یک کانال مستقیم، از مدل دوبعدی قائم MPS برای پیشبینی تغییرات پروفیل موج بهخصوص مقدار و موقعیت قله آن استفاده شده است. در این مثال کانالی به طول 7.5 متر و با عمق آب اولیه 0.5 متر در نظر گرفته شده است. پروفیل اولیه موج با معادله زیر تعریف شده است [17]:
که در آن h تراز
سطح آزاد موج، H ارتفاع موج، h عمق آب و x فاصله جهت طولی کانال است.
برای تولید موجی با ارتفاع 0.2 متر و طول 3 متر، معادله پروفیل اولیه موج با سرعت قائم صفر و سرعت افقی U با رابطه زیر بهعنوان ورودی مدل عددی در نظر گرفته شده است [17]:
که در آن g شتاب ثقل است.
در هر مسئله پارامترهای بهینه بر اساس عدد کورانت طبق رابطه زیر انتخاب شدهاند:
که در آن r0 فاصله اولیه بین ذرات، ∆t گام زمانی و Vmax حداکثر سرعت جریان است.
تغییرات تراز سطح آب در حالتهای بدون ویسکوزیته و با ویسکوزیته ثابت توسط مدل عددی MPS به دست آمده است. مقایسه نتایج با رابطه تئوری انتقال موج بدون استهلاک در شکل 1 نشان داده شده است. در این مثال قطر ذرات 0.045 متر در نظر گرفته شده و جمعا از 2138 ذره در محاسبات استفاده شده است. عدد کورانت برابر 0.059 به دست آمده است.
الف)
ب)
شکل(1): مقایسه پروفیل موج مدل عددی با پروفیل تئوری بدون استهلاک، از چپ به راست پس از زمان 0.24، 0.51، 0.75 و 1 ثانیه در حالت الف) بدون ویسکوزیته ب) ویسکوزیته گردابی ثابت برابر 0.002
همانگونه که در شکل 1 قابلمشاهده است شبیهسازی در حالتهای الف ( بدون ویسکوزیته) سطح آب ناهموار بوده و تعدادی از ذرات در هوا پراکنده شدهاند. ولی در حالت ب که مدلسازی آشفتگی انجام شده است، سطح آب بهصورت هموارتری به دست آمده است. مقایسه موقعیت و مقدار قله موج پیشبینی شده توسط مدل عددی MPS در زمانهای مختلف با مقدار تئوری آن نشان میدهد موقعیت قله موج در هر دو حالت در کلیه زمانها با پروفیل تئوری تطابق دارد. ارتفاع موج در حالت الف مدل عددی (بدون ویسکوزیته) با پروفیل تئوری منطبق است ولی در حالت ب به علت آنکه اثر ویسکوزیته جریان و استهلاک ناشی از آشفتگی توسط مدل آشفتگی دیده شده است، ارتفاع کمتری برای موج به دست آمده است.
در شکل 2 تغییرات میدان فشار در حالتهای مختلف در نظرگیری ویسکوزیته در زمان0.81 ثانیه نمایش داده شده است. مشاهده میشود به غیر از حالت بدون ویسکوزیته، حالت با ویسکوزیته گردابی ثابت بهخوبی توانسته است لایههای افزایش فشار از کف کانال تا سطح آزاد را پیشبینی نماید. همچنین یک منطقه پرفشار زیر قله موج به علت فشارهای دینامیکی ناشی از تغییرات شدید سطح آزاد آبدیده میشود. این موضوع با انتظار ما برای تغییرات فشار تطابق داشته و میتواند نحوه پیشبینی فشار را تأیید نماید.
الف) ب)
شکل (2): تغییرات فشار زیر موج پس از گذشت زمان 0.81 ثانیه در حالتهای الف) بدون ویسکوزیته ب) با ویسکوزیته ثابت برابر 0.002
3-2- حرکت تک موج در محفظه ساده
مدل لاگرانژی برای مطالعه حرکت تک موج در یک محفظه ساده با دیوارههای قائم بکار گرفته شده است. موج در حالت اولیه با معادله پروفیل h و با سرعت عمودی صفر و سرعت افقی U توسط روابط (17) و (18) تولید شده است [17].
طول محفظه 3 متر و عمق آب 0.5 متر در تظر گرفته شده است. موج با ارتفاعهای مختلف از 0.15 تا 0.3 متر در سمت چپ محفظه تشکیل شده و به سمت راست پیشروی میکند. پس از برخورد موج به دیواره عمودی محفظه بهصورت قائم روی آن بالا میرود. میزان بالاروی روی دیواره قائم به کمک عمق آب بیبعد شده است. مقادیر محاسبهشده توسط مدل عددی لاگرانژی در حالت بدون ویسکوزیته با نتایج آزمایشگاهی Camfield and Street (1968) و مدل عددی SPH Monaghan and Kos (1999) در شکل 3 مقایسه شده است [17]. برای نسبت ارتفاع موج به عمق آبهای 0.3 تا 0.6، از 900 تا 1800 ذره محاسباتی با قطر 0.033 تا 0.05 در مدل عددی استفاده شده است.
شکل (3): مقایسه نتایج مدل عددی MPS حاضر با دادههای آزمایشگاهی و شبیهسازی عددی SPH
مقادیر خطای محاسباتی نسبت به مقادیر آزمایشگاهی با استفاده از رابطه RMSE زیر برای روش SPH برابر 1% و برای مدل MPS حاضر برابر 3.7% به دست آمده است.
4- نتیجهگیری
مدل MPS دوبعدی قائم برای شبیهسازی مسائل موج توسعه داده شده است. به منظور اعمال مدل به جریان ویسکوز، مدلسازی آشفتگی به روش ویسکوزیته گردابی ثابت با مدل هیدرودینامیک ترکیب شده است. کاربرد مدل حاضر در مثالهای موجود در ادبیات فنی و مقایسه نتایج بهدستآمده با دادههای تئوری، تجربی و عددی، بیانگر ارتقای توانایی مدل در پیشبینی تغییرات سطح آزاد آب میباشد. همچنین افزودن مدلسازی آشفتگی به مدل عددی سبب افزایش پایداری و هموارسازی سطح آزاد آب شده است. بهخصوص محاسبه میدان فشار برخلاف اکثر مدلهای لاگرانژی بهخوبی انجام شده است.