شبيه سازي عددي شكست موج تنها بر روي ساحل شيب دار به روش
هيدروديناميك ذرات هموار نسبتاً تراكمپذير
محمد جواد كتابداري1* ، امين نعمتي روزبهاني2

دانشيار، دانشكده مهندسي دريا، دانشگاه صنعتي اميركبير
دانشجوي كارشناسي ارشد، دانشكده مهندسي دريا، دانشگاه صنعتي اميركبير
چكيده در اين مقاله از يك مدل عددي بدون شبكه بنام هيدروديناميك ذرات هموار نسبتاً تراكم پذير (WCSPH) براي شبيهسـازيروند شكست موج تنها بر روي ساحل شيبدار استفاده شده است. اين مدل دو بعدي بوده و سيال را به صورت كمـي تـراك مپـذي ر درنظر م يگيرد و علاوه بر حل معادلات حاكم بر سيال لزج براي به دست آوردن ميدان سرعت و چگالي، از حل معادله حالت براي به دست آوردن فشار استفاده ميكند. اين مسئله باعث كاهش حجم محاسبات نسبت بـه روش پايـ ه مـدل هيـدروديناميك ذراتهموار مي شود. براي شبيهسازي آشفتگي سيال در روند شكست موج از مدل آشفتگي SPS كـ ه بـه وسـيله تئـوري شـبيه سـازيگردابه هاي بزرگ (LES) به دست آمده، استفاده شده است. در تحقيق حاضر، براي بررسي دقت مدل در شبيهسـازي سـطح آزادابتدا مدلسازي شكست سد انجام شد و نتايج مدل با نتايج آزمايشگاهي Martin و Moyce(1952) مقايسه شد. سـپس بـرايمعتبرسازي مدل آشفتگي، نتايج پديده شكست سد توسط مدل با نتايج آزمايشگاهي Issa (2005) مورد مقايسه قرار گرفت. اين مقايسه ها نشان داد كه مدل تهيه شده ابزاري قوي جهت شبيهسازي رفتار جريان سيال آشفته است. در انتها تغييرات يـك مـوجتنها در حال شكست بر روي ساحل شيب دار مدلسازي شده و نتايج مورد بحث قرار گرفت.
كلمات كليدي: موج تنها، شكست موج، هيدروديناميك ذرات هموار نسبتاً تراكمپذير، ساحل شيب دار

Numerical Simulation of Breaking Solitary Wave on the Sloping Beach Using WCSPH Method

M. J. Ketabdari1, A. Nemati Roozbahani2

1-Associate Professor, Faculty of Marine Technology, Amirkabir University of Technology
2-M.Sc. student, Faculty of Marine Technology, Amirkabir University of Technology

Abstract
In this article, a numerical meshless method called Weakly Compressible Smoothed Particle Hydrodynamic (WCSPH) is used to simulate the solitary wave breaking process on the beach. The present model is a two dimensional model that considers the fluid weakly compressibility. This model solves the viscous fluid equations to obtain velocity field and density and solves the equation of state to obtain the pressure field. This leads to computational efficiency of the model in comparison with basic SPH model. To simulate the fluid turbulence in wave breaking process, the Sub Particle Scale (SPS) model that is obtained from Large Eddy Simulation (LES) theory of turbulence is used. In this research to consider the accuracy of free surface simulation, first a dam break test is performed to compare the model results with experimental data of
Martin and Moyce (1952). Afterwards to verify the turbulence model, a comparison between

ketabdar@aut.ac.ir : *
سال هشتم/ شماره 51/ بهار و تابستان 91 92

model results of dam break phenomenon and experimental results of Issa (2005) is carried out. These comparisons showed that the developed model is a powerful tool for simulation of the turbulent flow behavior. Finally, the transform of a breaking solitary wave on the sloping beach is modeled and the results were discussed.
Keywords: Solitary Wave, Wave Breaking, Weakly Compressible Smoothed Particle

1- مقدمه
شكست موج در نزديكي ساحل، يك مسئله تحقيقاتي مهم در مهندسي دريا و سواحل است.
خاصيت قوي غير خطي شكست موج باعث حركت پيچيده ذرات آب شده و از اين رو شبيهسازي عددي روند شكست موج بر روي ساحل را بسيار پيچيده م يكند. ميدان سرعت در هنگام شكست موج طبيعتي تصادفي و متغير داشته و از سويي ديگر با ورود حباب هاي هوا به همراه جت شيرجهاي، پيچيدگي هاي ديگري در اندازه گيري سرعت بوجود ميآيد .اندازه گيري موج در حال شكست در آزمايشگاه و طبيعت نياز به هزينه بالا و دستگاههاي خاص دارد. از طرفي پيشرفت تكنيكهاي عددي و ماشينهاي محاسبه گر در چند دهه اخير باعث شده كه مدل سازي عددي شكست موج مورد توجه زيادي قرار بگيرد. در گذشته به خاطر محدوديت در اندازه شبكه ها، استفاده از روش هاي برپايه شبكه باعث كاهش دقت حل اين مسئله مي شد. به همين دليل امروزه روش هاي بدون شبكه براي شبيهسازي رفتار سيال بيشتر مورد توجه قرار گرفته است .
30 سال هشتم/ شماره 51/ بهار و تابستان 91

روش هاي عددي كمي براي رهگيري سطح آزاد آب هنگامي كه موج به سمت ساحل حركت مي كند، وجود دارد. تعدادي از اين روش ها عبارتند ازMAC [1] و VOF [2]. اگرچه اين روش ها ب هطور موفقيت آميزي براي شبيه سازي امواج به كار رفتند ([3] و [4]) ولي استفاده از شبكه بندي اويلري، باعث ايجاد مشكلات زيادي در شبيه سازي شكست امواج و تغييرات ناگهاني سطح آزاد آب مانند پخش عددي1 ب ه دليل وجود ترم انتقال2 در معادلات ناوير- استوكس3 و يا استفاده از تعداد بسيار زياد المان ها براي شبيه سازي شكست موج م يشود .اخيرا روش هاي بدون شبكه بر پايه روش لاگرانژي4 ب همنظور غلبه بر اين مشكلات توسعه پيدا
Hydrodynamic, Slopped Beach
كرده اند. از ميان اين روش ها، هيدروديناميك ذرات هموار5، به عنوان يك روش قدرتمند جاي خود را در ميان ديگر روش هاي عددي باز نموده است. در اين روش ذرات در مختصات لاگرانژي حركت ميكنند و ترم انتقال در معادلات ناوير-استوكس ب هطور مستقيم به وسيله حركت ذرات محاسبه شده بدون اينكه مشكل پخش عددي وجود داشته باشد. هيدروديناميك ذرات هموار ابتدا براي مدلسازي هاي فضايي توسط لوسي6 [5] و همچنين توسط گينگولد و موناگان7 در سال
1977 ابداع گرديد [6]. سپس در زمينه هاي مختلف سيالاتي مثل مدل سازي هاي جريان سيال با فرض كردن سيال به صورت كمي تراكمپذير و استفاده از معادله حالت براي ب هدست آوردن فشار به نام “هيدروديناميك ذرات هموار نسبتاً تراكم پذير”8 گسترش يافت [7و8]. همچنين اين روش براي بسياري از شبيه سازي ها از جمله شكست سد9، پيشروي موج [9]، شكست موج بر روي سازه [10] و شكست موج تنها10 بر روي ساحل با شيب كم [11] گسترش داده شد بهطوري كه اكنون براي بسياري از مسائل پيچيده سطح آزاد آب به كار مي رود [12].
در اين مقاله به منظور درك رفتار سلسلة امواجي كه تك تك بر روي ساحل ميشكنند، شبيه سازي توليد يك موج تنها ب هوسيله موجساز پيستوني و سپس شكست آن بر روي ساحل ب هوسيله روش بدون شبكه هيدروديناميك ذرات هموار نسبتاً تراكم پذير صورت
گرفته است.

2- معادلات اساسي حاكم بر سيال لزج
معادلات حاكم بر حركت سيال، معادلات بقاء جـرمو ممنتوم11 بوده و به فرم لاگرانژي به صورت روابط زيـرنوشته مي شوند:

(1)
ρ Dt

163059-50441

Du1

=− ∇P+ g +Θ ()2
Dt ρ

rدر روابط بالا، ρ چگالي،u بردار شتاب ثقل ،t زمان
وΘ نيز ترم انتقال معادله حركت مي باشند.

2-1- گسسته سازي معادلات
در اين مقاله از دو نـوع تـرم انتقـال اسـتفاده شـدهاست. ابتدا از لزجت مصنوعي12 براي معتبر سازي اوليـهمدل و سپس از لزجت لايهاي به همراه يـك نـوع تـنشآشفتگي بـراي معتبـر سـاز ي مـدل آشـفتگي بـا انجـامآزمايش شكست سد و مقايسه آن با نتايج آزمايشـگاهياستفاده گرديد. با اعمال لزجت مصنوعي به عنـوان تـرمانتقال معادله حركت ،Θ=Πكه Π لزجت مصنوعي است و با توجه به روابط روش WCSPH، معادلات بقاء جرم و حركت به فرم روابط زير گسستهسازي ميشوند:

(Dρ)a =∑mb(ua −ub).∇raWab
217934-69196

Dtb

Du
()a =
Dt
643129129595

PaPbrr
−∑mb(2 + 2 +Πab).∇aWab + g
b ρa ρb

در روابط بالا،W تابع درونياب (كرنل13) و انديس b مربوط به ذرات مجاور ذره a است و سري هاي فوق بر روي تمامي ذرات همسايه ذره a جمع زده م يشوند.Πab يك نمونه لزجت مصنوعي است كه براي جلوگيري از نفوذ ذرات در يكديگر و همچنين كاهش پخش عددي در ترم انتقال معادله حركت به كار م يرود و ب هصورت روابط ذيل تعريف م يشود:

Πab µ
4296966399

−αcabab
=ρab

0
r vab.rab < 0  if (5)

otherwise

µab = huab.

2 ab2 ()6
rab +η

در روابط بالا،h طول هموار14 تابع درونياب است كه شعاع تاثير يك ذره نسبت به ذره مجاورش م يباشد[7].
13106511028

،cab = (ca +cb)/ 2 همچنين
1514863-1282

894586230367

2057407455917

/(rab =ra −rb ،ρab = (ρa +ρb و uab =ua −ub كه به ترتيب cab سرعت متوسط صوت ،ρab چگالي متوسط ،rab فاصله بين ذرات و
uab سرعت ذره a نسبت به ذره b است. η جهت
جلوگيري از صفر شدن مخرج به كار مي رود و داريم 2η2 =0/01h وα نيز ضريب ثابتي است كه در اين مقاله براي آن مقدار 3/0 انتخاب شده است .

2-2- تابع درونياب
اساس روش WCSPH بر پايه تئوري درونيابي بنا نهاده شده است. در اين تحقيق از يك نوع تابع درونياب به نام “مارپيچ درجه 3″15 كه توسط موناگان و لاتانزيو16 ارائه گرديد، استفاده شده است [13]:

 103 23 3
2 (1− q + q )q<1 
580649-105792

7π10h234
W(r,h) =2 (2−q)1<q< 2  (7)
28πh


0q> 2 

در رابطه بالا، q=r/h پارامتر بدون بعد طول است.

2-3- معادله حالت
در روش WCSPH، سيال تراكم پذير عمل مي كنـدكه به همين خاطر ميتوان از معادله حالت نوع تيـ ت17 براي ب هدست آوردن فشار استفاده كرد [14]:

58827198144

P=B ργ−1، B=ρ0cs2 /γ ()8

ρ0 

10134581532106

در رابطه (8)، B پـارامتر ثـابتي اسـت كـه بـه مـدولحجمي الاستيسيته مربوط است .γ عدد ثابتي بـي ن 1 و 7 است كه در اين مقالـه بـراي آن مقـدار ثابـت7 در نظر گرفته شده ،P فشار نسبي، 0ρ دانسيته در فشار اتمسفر ،v *(04 تا01) =cs كه cs سرعت صـوت وv سرعت بيشينه ذرات سيال در طول محاسبات اسـتكه در ايـن تح قيـق،cs = 10v درنظـر گرفتـه شـدهاســت. همچنــين (v= g(d+h كــه در آن h ارتفاع موج و d عمق آب است.

4-2- تكنيك حركت ذرات
حركت ذرات توسط تكنيك “هيدروديناميك ذرات هموار-ايكس”18 طبق رابطه( 9) تعريف م يشود:

99514585099

dra = ua +ε∑ mb uabWab ()9
dtb ρab

در رابطه بالا، ε مقداري ثابت است كه در اينجا برابر 5/0 در نظر گرفته م يشود. با استفاده از اين روش، به ذرات اجازه داده م يشود كه ب هخوبي سازماندهي شوند و در سرعتهاي زياد سيال، از نفوذ ذرات در يكديگر جلوگيري م يشود [15].

2-5- الگوريتم حل زماني معادلات
بــراي حــل معــا دلات حــاكم، از الگــوريتم زمــاني پيشبيني-تصحيح19 استفاده شده اسـت [15]. روابـط
(3) ، (4) و( 9) را م يتوان به فرم خلاصه زير بازنويسـيكرد:
 Dρ
(

)a = Da
 Dt
dx (

)a = Fa (10)
 Du
87478392755

Pan+1/ 2 = Bρaρn+01/ 2 γ−1 (12)


با استفاده از مقادير ب هدست آمده در اين مرحله، مقادير سرعت، چگالي و مكان ذره a ب هصورت رابطه زير اصلاح ميگردند:

 n+1/ 2n ∆tn+1/ 2
ρa=ρa + Da
865639-67870

2
 n+1/ 2n ∆tn+1/ 2
ua= ua +Fa
2
 n+1/ 2n ∆tn+1/ 2
ra= ra + 2 Ua

در نهايت پس از طي زمان t∆ و در پايان گام زماني 1n + ام و استفاده از مقادير اصلاحي، طبق روابط زير خواهيم داشت:

ρan+1 =2ρan+1/ 2 −ρan
uan+1 =2uan+1/ 2 −uan
 n+1n+1/ 2n
ra =2ra−ra

77572388999

Pan+1 = Bρρan0+1 γ−1


براي تعيين مقادير گام زماني بهينه از شـرط كورانـت20 طبق رابطه زير، استفاده م يشود:

 Dt  Dr
(

)a =Ua
 Dt

چنانچه روند حل از گام زماني n آغاز گردد، در گام زماني 5/0n+ خواهيم داشت:

 n+1/ 2n ∆tn
ρa=ρa +

Da
2
 n+1/ 2n ∆tn
ua= ua +

Fa (11)
2
 n+1/ 2n ∆tn
ra= ra +

2 Ua
32 سال هشتم/ شماره 51/ بهار و تابستان 91

∆t ≤ C

(16) Vmax

در رابطه بالا، C عدد كورانت بوده و معمولا 1/0 در نظر گرفته م يشود. dx فاصله اوليه بين ذرات و Vmax نيز سرعت حداكثر ذرات سيال در هر گام زماني است [16].

6-2- فيلتر چگالي
شبيه سازي ديناميكي در اين روش منجر به تغييرات زيادي در ميدان فشار ذرات سيال م يشود.
تلاش هايي به منظور غلبه بر اين مشكل صورت گرفته است كه يكي از آنها طراحي فيلتر مناسبي براي چگالي يا دوباره مقدا ردهي چگالي براي ذرات است. براي اين
منظور، از اصلاح درجه يك ام-ال-اس21، استفاده م يگردد. اين ترم در هر 03 گام زماني به فرم رابطه (17) يكبار بر روي ميدان چگالي ذرات اعمال م يشود
:[17]

WabMLS =
[β0(ra)+β1x(ra).(xa −xb)+β1z(ra).(za −zb)]Wab
β0 1
−1 
780286473897

β(ra)=β1x = A 0 β1z 0 A=∑Wb(ra)mbA~
bρb
~  1 (xa −xb)2 (za −zb)  A=(xa −xb) (xa −xb) (za −zb)(xa −xb)
(za −zb) (za −zb)(xa −xb)(za −zb)2
(17)

2-7- مدل آشفتگي
در اي ن تحقي ق از دو ن وع ت رم انتق ال در معادل ه حركت استفاده شده است. اولـين آن لزجـت مصـنوعيبود كه درروابط( 5) و( 6) توضيح داده شد. دومين آن اسـتفاده از لزجـت لايـه اي ب ههمـراه يـك نـوع ت نشآشفتگي است و با اين ترتيب رابطه 2() به رابطـه( 18) تبديل م يشود:

88542127036

كه 0ν لزجت سينماتيكي آب و τij نيز تنسور تنش اس- پي-اس22 است و طبق رابطه زير تعريف م يشود:

810769-499089296596233

124510996233

τijtijij2I2 (19) = 2νS −kδ − C ∆l δ S
ρ3

1600185-5447422
ij
ij

2

ij

ij



قیمت: تومان

دسته بندی : مهندسی دریا و بندر

پاسخ دهید