بردار یکه نرمال eˆ n1 Psat فشار اشباع
بردار حالت آب یا گاز در p Up g چگالی بخار اشباع
شرایط چپ و راست H L ,H R l چگالی مایع اشباع

١- مقدمه
بسياري از پديده هـاي فيزيکـي در سـرعت هـاي بسـيار بـالا درجريان ايجاد می شـود و انتشـار امـواج تراکمـی و انبسـاطی درفرآيندهای فيزيکی ناشی از سرعت بالا رخ خواهد داد. جريـانسيال در اين پديده ها داراي ماهيتي تراکم پذير بوده و بسياري از مسائلي که در حوزه ي علوم مکانيک سيالات و يا هوا فضا قـرارمـی گيرنـد داراي ايـن شـرايط هسـتند. مسـائلي نظيـر حرکـت موشک ها، جريان مافوق صوت در کانال ها و يـا مسـألة انفجـارجزء اين دسته از مسائل قرار خواهند گرفت.
انفجار پديدة پيچيده اي است کـه در طـي آن در کسـري ازثانيه مادة منفجره تبديل به گاز با فشار و دماي بـالا مـي شـود وانرژي خود را به صـورت ناگهـاني در محـيط آزاد مـي کنـد . در صورت مهـار انـرژي عظـيم آزاد شـده ناشـي از يـک واکـنشانفجاري و استفاده از خصايص آن و هدايت آن مـي تـوان ايـنانرژي را در بسياري از فرآيند ها و صنايع مورد استفاده قرار داد.
از جمله کاربردهاي انفجار می توان به اسـتفاده از ايـن مـواد درکارهـای عمرانـی و معـدنی، شـکل دهـي فلـزات توسـط مـوجضربه اي در محيط هاي مانند آب و هوا، جـوش دادن قطعـات وپوشش دهي مواد نام برد.
مواردي که در بالا ذکر شد اثر انفجار برروی اجسـام جامـددر مجاورت مادة منفجره را نشان می دهند. طبيعي است که ايـناثر شديد و سريع برروی اين اجسام در اثر عبور امواج تراکمـياز محيط سيال حاصل مي شود و با توجه به اينکه انفجار در چه سيالي صورت پذيرد نحوة انتشار اين امـواج و سـرعت انتشـارامواج در محيط و در نهايت رسيدن موج فشاري به جسم مـوردنظر اهميت پيدا خواهد کرد. محيط آب، يکي از پرکـاربرد تـرينمحيط هايي است که انفجار در آن صورت مي گيرد. آب با توجه به خصوصيات ترموديناميکي ويـژ ة خـود، در شـرايط مختلـفرفتارهاي متفاوتي را نسبت به هوا بروز مي دهـد . يکـي از ايـنپديده ها که در بسـياري از مـوارد حـائز اهميـت اسـت، پديـدة کاويتاسيون مي باشد. شـناخت و مطالعـة ايـن پديـده بـه منظـوراستفاده مناسـب از انـرژي حاصـل از انفجـار در محـيط آب ازاهميت به سزايي برخوردار است. يکي از زمينه ها و کاربردهـاييکه فيزيک انفجار زير آب و کاويتاسيون بـه وضـوح در آن ديـدهمي شود صنايع دريايي و کشتي سازي و اژدرهاي مورد اسـتفادهاست. انفجار مواد منفجره از قبيل موشک ها، مين هـاي دريـايي،چاه هاي نفتي، در نزديکي سازه هـاي دريـايي ماننـد اسـکله هـا،خطوط انتقال زير آب و کشتي هـا نـوعي بارگـذاري را بـر ايـنسازه ها تحميل مي کند که به دليل فيزيک پيچيـد ة آن و همچنـين صدمات وارده از لحاظ علمی و عملـی دارای اهميـت بسـزايياست.
به طور کلي پديدة کاويتاسيون در جريآن هاي تراکم پـذير بـهمعناي تشکيل فاز بخار در سيال است که خـود معلـول کـاهشفش ار س يال ت ا ح د فش ار اش باع س يال م يباش د [١و ٢]. درطبقه بندي کاويتاسيون، اين پديده به دو گروه کاويتاسيون دائم يا متصل١ و کاويتاسيون غيردائم يا گذرا٢ تقسيم بندي مـي شـود.
کاويتاسيون غيردائم شـامل پديـده هـاي پيچيـدة فيزيکـي ماننـدتشکيل ناحية کاويتاسـيوني و سـطح تمـاس ديناميـک، توسـعة کاويتاسيون و درنهايت اضمحلال اين ناحيه مـي شـود . در دهـة اخير کاويتاسيون به طور گسترده اي توسط محققـين مختلـف بـاروش هاي آزمايشگاهي و عددي مورد مطالعه قرار گرفته اسـت [٣- ٥]. همچنين لازم بهذکر است که عموم مطالعـات صـورت گرفته برروی کاويتاسيون متصل تمرکز يافتـه اسـت. ايـن نـوعکاويتاسيون به صورت جريان دائمي بوده و تغييـر شـکل ناحيـة کاويتاسيون بسيار آهسته و يا به صورت تنـاوبي خواهـد بـود. از طرف ديگــر شکل و ســطح تماس ناحيــة کاويتاســيون در کاويتاسيون غيردائـم بسـيار متحـرک و پويـا بـوده و در پديـدة انفجار زير آب و خطوط لولة صنعتي قابل مشاهده است.
به منظور شبيه سازي جريان دوفازي در روش تکسيالي، با دو فاز موجود به صورت يک مخلوط برخورد خواهد شـد. از اينـرو بر خلاف روش هاي ديگر- نظير روش دوسيالي- تنها يک دسـتهمعادلة ديفرانسيل براي بيان حرکت سـيال کفايـت مـي نمايـد . بـاتوجه به آنکه يک دسته معادله ديفرانسيل در مدلسـازي اسـتفادهمي شود، تنها يک معادله حالـت نيـز بـراي کـل مخلـوط بـه کـار مي رود. يکي از مهم تـرين دغدغـه هـا در ايـن روش، اسـتفاده ازمعادله حالتی است که بتواند تغييرات سيال از حالت مايع اشـباعتا فاز بخار را پوشش دهد [6-٨]. همان طور که مشـخص اسـتدر اين روش با توجه به ماهيت آن، جزئيات گذار از يک فاز بـهفاز ديگر و نيز ميزان تبادل جرم و انرژي بين يک فـاز بـا ناحيـة کاويتاسيون محاسبه نمي شود. اما از جهـت ديگـر، در ايـن روشم دل س ازي ايج اد ناحي ة کاويتاس يون و توس عه و در نهاي ت اضمحلال آن به سادگي انجام مي پذيرد [6 و ٧] و اين مدل به طور گسترده اي در شبيه سازي پديدة انفجار زير آب مورد استفاده قرار گرفته است. در بين مدل هاي تک سـيالي، مـدل خـلاء٣ يکـي از اولين مدل ها ي به کار گرفته شده است که توسط تنگ و هانگ در سال ١٩٩6 ارائه شده است [٩]. اين مدل يک مـدل بـا يـک فـازخالص است و داراي مشکلاتي براي بهکارگيري آن در چند بعـدمي باشد. يکي ديگر مدل هاي محبـوب و پرکـاربرد، مـدل قطـع٤ مي باشد. اين مدل براساس يک ايدة ساده شکل گرفته است و در اين مدل هنگامي که ميزان فشار به مقداري کمتر از فشـار اشـباعمي رسد، فشار ناحية کاويتاسيون برابر بـا فشـار اشـباع قـرار داده مي شود [١٠]. مدل قطع نيز مانند مدل خلاء يک مدل با يک فـازخالص بوده و قوانين بقا را نقض مي کند.
يکي ديگر از مدل هاي کاويتاسيون، مدل اشـميت اسـت کـهنخستين بار توسـط اشـميت و همکـاران در سـال ١٩٩٠ بـرايجريان فشار بـالا و سـرعت بـالا در يـک نـازل کوچـک مـورداسـتفاده قـرار گرفـت [٧]. در ايـن مـدل، جريـان کاويتاسـيون به صورت مخلوطي همگن و باروتروپيک از گاز و مـايع در نظـرگرفته شده است. در اين روش رابطة مـورد اسـتفاده در معادلـة حالت٥ با استفاده از انتگرال گيـري از رابطـة سـرعت صـوت درجريان دوفازي به دست مي آيد. در حقيقت، تنها در حالت هـاييمي توان مدل اشميت را به کار بسـت کـه نسـبت چگـالي بسـيارکوچک باشـد و يـا اينکـه ميـزان پـرش فشـار در عـرض مـرزکاويتاسيون از يک حد خاصي بـالاتر باشـد. ايـن شـرايط تنهـازماني ارضاء مي شود کـه فشـار محـيط بـالا باشـد و يـا انـدازهکاويتاسيون کوچک باشد. به همين دليل است که اين مدل براي جريان با فشار بالا داخل يک نازل کوچک مدل مناسبی است.
پس از مدل اشميت يکي از مدل هـايي کـه بـراي شـبيه سـازي
جريان کاويتاسيون معرفي شده است، مدل ايزنتروپيک مي باشد کـهدر سال ٢٠٠٤ ميلادي توسط ليو و همکاران ارائه گرديـد [١١]. در اين مدل، مولفه هاي بخار مخلوط کاويتاسيون بـه صـورت همگـن،تراکم پذير و ايزنتروپيک درنظر گرفته مي شـود . ايـن مـدل قابليـتخوبي در پيش بيني اندازة پيک فشاري و سرعت انتشـار امـواج رادارد و در بسياري از مسائل انفجـار زيـر آب مـورد اسـتفاده قـرارگرفته است [١٢].
ژي و همکاران [١٣] در سال ٢٠٠٦ مدل ديگري بـه نـام مـدلاصلاح شدة اشميت را ارائه کردند. ايـن مـدل تـا حـدي توانسـتضعف هاي مدل اشميت را پوشش دهد و براي دستة وسيع تـري ازمسائل قابل استفاده باشد. همچنين اين مدل براي شبيه سازي انفجار يک بعدي و دو بعدي مورد استفاده قرار گرفت.
در مقاله حاضر مسـألة انفجـار در زيـر آب در نزديکـي جسـمجامد با استفاده از روش اويلري – لاگرانـژي ALE بـرروی شـبکة تطبيقي شبيه سـازي شـده اسـت و انتشـار امـواج در داخـل آب وتشکيل و اضمحلال ناحية کاويتاسيون مـورد بررسـي قـرار گرفتـهاست. به منظور حل جريان از يک معادلة حالت متناسب بـا فيزيـکحاکم بر انفجار زير آب به همراه حلگر دقيق ريمـان اسـتفاده شـدهاست. بخش هاي مختلفي که در ادامه آورده شده اسـت بـه ترتيـبعبارتند از: معادلة حاکم، مسأله ريمان و معادلة حالت، روش عددي و نتايج.

۲- معادلات حاکم
در جريآن هاي تراکم پذير با سرعت بالا بـهعلـت غالـب بـودن اثـرعبارت هاي همرفت6 به عبارت هاي پخش٧ از اثر عبارت هاي پخش در معادلات ناوير- استوکس صرف نظـر شـده و جريـان سـيال بـااستفاده از معادلات اويلر توصيف مي گردد. ايـن معـادلات در فـرمبقايي عبارتند از:
121920100928

Ut FxU GyU 0
 uv
U   u , F U    u2 P , G U   uv
  vuv v2 P
که U بردار کميت هـای بقـايي،F و G بردارهـای شـار مـی باشـند . همچنين چگالی، u و v سرعت هـای جريـان در راسـتایx و y بوده و P بيانگر فشار میباشد.
در مسألة انفجار در زير آب در حلALE ۸ با دو محـيط مجـزاشامل: محيط گازي داخل حبـاب انفجـاري و محـيط آب پيرامـونحباب مواجه هستيم. از اينرو به منظور بسته شدن دستگاه معـادلاتديفرانسيل حاکم بر جريان از معادلة حالت مناسب با هر يک از دو محيط بايد استفاده شود. بدين منظور از رابطة معادلـة گـاز ايـده آل براي محيط گاز استفاده شده است. براي محـيط آب نيـز از معادلـةحالت مناسب استفاده شده است که در قسـمت هـاي بعـد توضـيحداده خواهد شد.

٣- حل دقيق مسأله ريمان
در بيشتر روش هاي عددي در جريـآن هـاي تـراکم پـذير بـه منظـوربه دست آوردن شار عبوري از سطوح سلولي نيـاز بـه حـل مسـالهريمان در هر مرز سـلولي مـي باشـد [١٤]. در واقـع مسـاله ريمـانعبارت است از بررسي جريان حاصل در معادلة اويلر با يک شـرطاولية ناپيوسته در مرز. در حالت اوليه سمت راست سيال در حالـتثابت Hr و طرف چپ در حالت ثابتHl مي باشد:
Ut0 HHlRxx00 (٣)
لازم به ذکر است که حل دقيق مسألة ريمان تنها در موارد خـاص وبراي گاز ايده آل و آب با معادلة حالت مشخص موجود بـوده و دربسياري از موارد از حل تقريبي ريمـان بـه منظـور بـه دسـت آوردنکميت هاي برداري استفاده مـيشـود . در روش هـاي تـک سـيالي وبه منظور تعيين ناحية کاويتاسـيون در آب، حلگرهـاي دقيـق مسـألةريمان بدين صورت عمل مي کنند که از حل دقيق ريمان تک فازي استفاده مي شود و بردار کميت هاي بقايي به دست مي آيد [١٥]. پس از آن با استفاده از معادلة حالت و دانستن مقدار فشار اشباع سـيال،ناحية کاويتاسيون تعيين مي گردد. در اين پـژوهش از حلگـر دقيـقريمان که در مرجع [٩] توسط تانـگ و هوآنـگ ارائـه شـده اسـتاستفاده شده است.
در اين حلگر از معادلة حالت تيت٩ مطابق با يک فرآيند تغيير فاز هم دما در فاز مايع استفاده شده است. همان طور که مي دانـيمپديدة کاويتاسيون در جريان تراکم پذير با کاهش فشار و در يـکفرآيند دما ثابت رخ مي دهد. شکل (١) نشـان دهنـدة ايـن پديـده برروی دياگرام فازي مي باشد.
همان طورکه در شکل (٢) نشان داده شده است هنگام ايجـادکاويتاسيون در جريان تراکم پـذير بـرروی خـط دمـا ثابـت روينمودار p-v حرکـت کـرده و از حالـت مـايع اشـباع وارد ناحيـةدوفازي و سپس ناحية بخار خواهيم شد. بنابراين معادلـة حالـت مورد نظر با توجه به فيزيک حاکم در يـک فرآينـد هـم دمـا يـکمعادلة باروتروپيک بوده و با استفاده از رابطة (٤) تعيين ميگردد:
(Psat  B)

   Bl
l 
P  Psat gl (٤)
C2vapg

٤- روش عددي حل جريان دوفازي

3577590-2159423

شکل ۱- تغيير فاز در فرآيندهاي مختلفشکل ٢- تغيير فاز در يک فرآيند همدما

عمومًاً جريان هايي كه شامل انتشـار امـواج در دو محـيط غيرقابـلاختلاط با معادلات حالت متفاوت مـي شـود، دو ناپيوسـتگي مـوجضربه و مرز بين دو محيط وجود خواهـد داشـت. مـثًلاً در مسـالهانفجار زير آب، به وجود آمدن موج ضربه ناشـي از اخـتلاف فشـاربسيار زياد در دو سوي مرز، باعث حضـور ايـن دو ناپيوسـتگي درنزديكي يكديگر مي شود. برهم كنش اين دو و اثرات آنها بر يكديگر كه در هنگام انعكاس موج از مرزهاي محيط روي مي دهد نيز بايـدمورد توجه قرار گيرد. لازم به ذکر است در اين پژوهش با توجه بـهحضور ناحية کاويتاسيون عمًلاً با سه محيط مواجه خواهيم بود کـهشامل ناحية گازي، ناحية آب و ناحية دوفـازي کـه در اثـر کـاهشفشار آب رخ مي دهد خواهد بود و به منظور مدل سازي فيزيک ايـنجريان پيچيده بايد محل اين نواحي و سطوح تماس بين اين نواحي به کمک روش عددي شبيهسازي شود.

٤-١- روش حل عددي ALE
در كار حاضر از روش مرتبـه دوم ALE بـراي حـل معـادلاتحاكم و به دست آوردن خواص جريان استفاده مي شود. با توجه به وجود مرز متمايز بين ناحية حباب گازي و آب، مرز مشـترکبين اين دو ناحيه به صورت کاملا لاگرانژي شبيه سـازي شـده وساير نقاط شبکه نيز به صورت دلخواه حرکت داده خواهند شـد .
به منظور ايجاد کمترين پيچيدگي در شبکة محاسـباتي در هنگـامحرکت لاگرانژي نقاط شـبکة محاسـباتي در دو محـيط حبـابگازي و آب با استفاده از يک تابع هموار حرکت داده مي شـوند .
در واقع با به کار بردن روش حاضـر از مزيـت هـاي روش هـاياويلري و لاگرانژي به طور همزمان استفاده شده است.
ترتيب كار بدين صورت است كه ابتدا خواص ميـدان حـلدر هر المان با توجه به اين که اين المان در کداميک از دو محيط قرارگرفته است، مشخص مي شود. سپس با استفاده از حل دقيـقريمان براي دو محيط آب و گاز، سرعت در نقاطي که مرز بـيندو سيال را مشخص مي کنند تعيين مي گـردد . پـس از مشـخصشدن سرعت در نقاط مرزي سرعت حرکت بقيه نقاط شبکه بـه گونـه اي تعيـين مـ يشـود کـه کمتـرين درهـم پيچيـدگي را در المان هاي آن ايجاد کند. در شبيه سازي حاضر بـه منظـور تعيـينسرعت نقاط از حل معادله لاپلاس استفاده شده اسـت . شـرايطمرزي برروی ديواره ها سرعت صفر و بـرروی سـطح حبـاب ازسرعت لاگرانژي حرکت حباب که از حل دقيـق مسـأله ريمـانبه دست مي آيد استفاده مي شود. پس از آن محـل نقـاط داخلـيشبکه در گام زماني بعدي محاسبه مي شود.
استفاده از حل معادله لاپلاس به دليل همـوار بـودن آن و عـدموجود اکسترمم هاي نسبي در دامنه حل اين معادلـه اسـت. پـس ازتعيين سرعت و موقعيـت نقـاط شـبکه در گـام زمـاني مـورد نظـرمعادلات جريان براي هر المان به روش حجم محدود انتگرال گيـريمي شوند. براي بهدست آوردن شار جريان روي هر ضلع المان ها نيز مجددًا از حل دقيق معادله ريمان به همراه معادله حالت مناسب براي هر محيط استفاده مي شود. پس از تعيين متغيرهاي حالت در تمامي المان هاي شبکه لازم است که کيفيت شبکه بـراي اسـتفاده در گـامزماني بعد مورد آزمـون قرارگيـرد و درصـورت نيـاز، ميـدان حـلمجددًا شبکه بندي شده و اطلاعات از شبکه قبلي بـه شـبکه جديـدمنتقل شوند و سپس عمليات حل با استفاده از شـبکه جديـد ادامـهمي يابد.
براي راحتي كار، ابتدا روش عددي براي حالت يـك بعـديتوضيح داده مي شود. در اين حالت معادلات حاکم بـراي حجـم کنترل (t) با سطح کنترل (t) را مـي تـوان بـهصـورت زيـرنوشت

t  U d  Fd  0 (٥)
(t)(t)
u  x F(U)   u  x u  P  (6)
E(u  x) Puكه در آن U   [ u E]T تابع شار مطابق با رابطـة (6) و x سرعت نقاط شبکه مي باشد. درصورت درنظـر گـرفتن فـرمنيمه منفصل معادلات حاکم به صورت زير نوشته ميشود:
25289-10968

dtd U i  i  Fi 0

262599540842

114085740842

که Ui مقدار متوسط سلولي است وFi شاری است که از سطح کنترل i عبور ميکند:

U i Ud
721757298867

1843421298867

i Fi F U ,xi12 ,xi12 F U ,xi12 ,xi12  (٩) خـواهيم [tn, tn+1] با انتگرال گيري از اين معادلات در بازه زماني :داشت
tn1
83072313620

in1U

in1inUin   Fdti 0
tn
که به دليل تغييـر شـکل حجـم کنتـرل در بـازه زمـاني يادشـده،انتگرال زماني تابع شار به سرعت حركـت شـبكه بسـتگي پيـدامي كند و بهصورت زير خواهد بود:
1130189-65881

Fdti tFi Uk,xn12 ,xn12  براي تعميم روش به حالت دوبعدي، ابتدا به هر نقطه مرزي يـکمحور مختصات نسبت داده مي شود که جهات محورهـاي آن درمختصات کلي باتوجه بـه شـکل (٣) و از روابـط (١٢) محاسـبهمي شوند:
ˆt  eˆt12eˆt2 nˆ  eˆn12 eˆn2
که ˆn و tˆ بردار نرمال عمود و مماس بر ضلع المـان مـي باشـد . سپس خواص سيال آب و گاز در دوسوي اين نقطه با استفاده از خواص المان هاي هـم نـوعي کـه در ايـن نقطـه اشـتراک دارنـد،به روش متوسط گيري وزني براساس فاصله مراکز آن المـان هـا از نقطه مزبور به ترتيب زير به دست آورده ميشود (شكل (٤)):

UigU
12894471640

378095-6740

i dig i dilil Upg i d1ig Upl  i d1il (١٣)
کهUp بردار حالت آب يا گاز در اين نقطه و Ui بـردار حالـت

شکل ٤- انتقال خواص از مراکز المان ها به رئوس آنها
:
شکل ٣- ايجاد محورهای مختصات محلی

هريک از المان هاي هم نوع پيرامون آن است. همچنـينdi فاصـلهمرکز المان تا ايـن نقطـه مـي باشـد . بـا اسـتفاده از ايـن خـواص وبه کارگيري حل دقيق مساله ريمان، سرعت اين نقطـه در جهـت ˆn به دست ميآيد. مقدار سرعت در جهت tˆ نيز بـ ا توجـه بـه جهـتسرعت در امتداد ˆn به صورت زير محاسبه ميشود:
if *n   0 *t tl if *n   0 *t tg که v*n سرعتي است كه در جهت ˆn از حل دقيـق مسـاله ريمـان به دست مي آيـد و vtl و vtg مولــفه هـاي سـرعت در جهـتt ˆ هستند که از تصوير ميانگين سرعت هاي المان هاي هم نوع پيراموني در جهت tˆ در دوسوي نقطـه مزبـور بـه دسـت مـي آينـد . پـس ازمحاسبه v*t و v*n وتبديل آنها از مختصات محلـي بـه مختصـاتکلي x,y شرط مرزي در مرز بين دو محـيط بـه دسـت مـي آيـد .
علاوه بر آن شرط مرزي برروی ديواره ها نيز صـفر در نظـر گرفتـهمي شود. در نهايت معادلة لاپلاس برروی نقاط شبکه با استفاده از دو شرط مرزي به دست آمده حل شده و سرعت نقاط شبکه بـه دسـتمي آيد. سپس هندسه جديد شبکه با استفاده از روابط زير بـه دسـتآورده ميشود:
2594610500895

ypn1  ypn typn xnp1  xnp txpn همچنين با استفاده از روابط زير موقعيت و سرعت نقاط شـبکه درزمان 12tn محاسبه مي شود:

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

شما می توانید تکه های دیگری از این مطلب را با جستجو در همین سایت بخوانید

ولی برای دانلود فایل اصلی با فرمت ورد حاوی تمامی قسمت ها با منابع کامل

اینجا کلیک کنید

313944-42683

xpn12  xn1t xn ,ynp12  yn1t yn
218799-29629

xnp12  xn12 xn ,ynp12  yn12 yn در روابط فوق زيرنويس p مربوط به نقاط مي باشد. پس از به دست آوردن سرعت هر نقطه از شبکه، سرعت هـر ضـلع از المـان هـا ي شبکه از روابط (١٧) محاسبه ميشود:
1187301122306

288141122306

2ye  yp12 yp2 xe xp12xp که زيرنويس e مربوط به اضلاع شبکه است و نيـز 1p و 2p معـرفدو رأس طرفين هر ضلع هستند. اين سرعت ها از مختصات کلي به مختصات محلي که در جهات عمود و مماس بر هر ضـلع تعريـفمي شوند، منتقل مي گردند (vn edge,vt edge در شکل (٥).
از حل دقيق مساله ريمان برروي هر ضلع شبکه در زمان tn و در جهت عمود بر آن، بردار *U براي اين ضلع در اين گـام زمـانيبه دست ميآيد:
*


دیدگاهتان را بنویسید