در کــار حاضــر جریــان ســیال آب، درون همزنــی بــا ابعــادآزمایشگاهی (با حجم ۱۰ لیتر) که دارای چهار پره ثابت و یـکپروانه شش پره راشتون می باشد، بـه روش عـددی شـبیه سـازیشده است. این هندسه، یک هندسه مرسوم در کارهـای عـددیاست و برای آن تعدادی نتایج عـددی و آزمایشـگاهی موجـوداسـت (شـکل ۱). عـدد رینولـدز در ایـن همـزن بـه صـورت Re = ND2 ν تعریـف شـده و نـوع جریـان داخـل همـزن رامشخص می کند. در این رابطه، N سرعت زاویـه ای پروانـه (در واحد D ،( revs قطر پروانه وn لزجت جنبشی سیال کـاریاست. به دلیل وجود هندسه پیچیده و دوران پروانه، عدد رینولدز جریان ۲۹۰۰۰ بوده و جریان ایجاد شده مغشوش می باشد. دلیل
جدول ۱ – مقایسه نتایج حاصل از روش میدان نیرو با روش کمانه کردن اصلاح شده
Re =100 روش میدان نیرو [۲۲] روش کمانه کردن اصلاح شده [۳۲] حد مجاز بالا و پایین [۳۱]
CD ۳/۳۳۳ ۳/۲۲۷۵ ۳/۲۲- ۳/۲۴
CL ۱/۰۵۱۱ ۱/۰۰۴۰ ۰/۹۹-۱/۰۱

انتخاب این عدد رینولدز، نتایج آزمایشگاهی موجود است [۲۳].

۳ – روش شبکه بولتزمن
روش شبکه بولتزمن، روشی با مبنای ذره ای – جنبشی است که برای شبیه سازی جریان های سیال مورد استفاده قرار می گیرد. در این روش، معادله جنبشی به منظور بـ ه دسـت آوردن تـابع توزیـعسرعت ذره f(x,u,t)r r حل می شود که در آن ur بـردار سـرعت،xr بردار مکان ذره و t زمان است. کمیت هـای م اکروسـکوپیکمانند چگالی ρ و ممنتوم ρur با استفاده از مقادیر تـابع توزیـعاحتمال حضور ذره f که از حل معادله جنبشی بهدست می آینـدمحاسبه می شوند. یک مدل جنبشی پر کاربرد، تقریـب گسسـتهشده معادله بولتزمن با یـک زمـان آرامـش اسـت کـه بـه مـدل۱۰BGK نیز موسوم است. شکل BGK معادله بولتزمن بهصورت زیر می باشد:

∂∂ft + u. fr∇r =−

λ1 (f −f (eq)) (١)
که در آن f(eq) ، تابع توزیـع تعـادلی (توزیـع تعـادلی ماکسـول – بولتزمن) و λ زمان آرامش می باشد. در حالت طبیعـی یـک ذرهسیال در بی نهایت جهت (ur ) مجاز به حرکت میباشـد . اولـینگــام بــرای حــل عــددی معادلــه (۱)، بــه منظــور محاســبه f، گسسته سازی سرعت ur (سرعت حرکت ذره) می باشـد . بـرایاین منظور بدون آن که قوانین بقاء نقض شوند، ذره به حرکـتبا سرعت های خاصی (urα ) محدود می گردد [۲۴]:

(٢)
در معادله بالاfα ≡ f(x,u ,t)r rα ، تابع توزیع مربوط بـهa امـین سرعت گسسته شده urα ذره است. fα بیانگر احتمـال وجـودذره ای در زمــان t در مکــان xr اســت کــه دارای ســرعتی

شکل ۱ – هندسه جریان [۷]، همزن(شکل سمت چپ) دارای چهار پره ثابت است که از دوران صلب سیال جلوگیری میکند، پروانه (شکل
سمت راست) یک پروانه شش پره راشتون است. ضخامت دیسک و ضخامت پر ههای ثابت و متحرکD 017.0 است

(
شکل ٢ – یک مدل سه بعدی، شبکه نوزده سرعتی ( 19D3Q

برابر urα است و fα(eq) تابع توزیع تعادلی متناظر با آن اسـت . می شود [٢۵]:
مـمدـدلل 19نــQو3Dزده معرســورفعت یاس تســ (ه بعـشکـدل۲ی) ، شــبیکیک ه از بــومدلتل زمهانی، کــسه ه بعدبــهی 013 , α=
408432-336669

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

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

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

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

است که با موفقیت در حل بسیاری از مسائل سه بعـدی جریـان(۴) 61wα =181 , α= −
 1
سیال مورد استفاده قرار گرفته است. در این مـدل سـرعت هـای18736 , α= −
گسسته شده ذرات، reα می باشند. با گسسته سازی فضای سرعت ذره در معادله جنبشی، چگالی و تابع توزیع تعادلی بـرای شـبکه 19D3Q طبـق رابطـه زیـر ممنتوم به شکل زیر قابل محاسبه می باشند:

4044696105725

5263896105725

ρ =ur ∑18 e frα α = ∑18 e frα α (eq) ((۶۵))f (eq) =ρw [1+ 3 e .ur r + 9 (e .urα r)2 − 32:شودu.ur r تع٣)ری ف [می)
ααc2 α2c42c
609600118594

w، ضریب وزنی است و طبق رابطه (۴) جایگذاری 1α=1 r α=و در آن a سرعت صوت در این مـدل برابـر اسـت بـاcs = c کـه در
r
3
2676145-23407

آنcr =δxrδt و δxr و δt ، به ترتیب ثابت شـبکه و انـدازه گـامزمانی است. شکل کاملاﹰ گسسته شده معادله (١) با گام زمانی δt و گام مکانی reαδt به صورت زیر خواهد بود:
fα (xri + erαδt,t +δ −t) fα (x ,tri ) =

1τfα (x ,tri )−fαeq (x ,tri ) (٧)
2365248-36219

که در آنτ=λδt زمان آرامش بدون بعـد و xri مختصـات یـکنقطه در فضای فیزیکی می باشد. به این معادله، معادلـه گسسـتهشده بولتزمن با تقریب BGK گفته می شود. این معادله معمولاﹰ در دو مرحله، مطابق رابطه (٨) حل می شود. در این رابطه، علامت ~ بیانگر مقدار تابع توزیع بعد از مرحله برخورد می باشد.
مرحله برخورد۱۱:
f (x ,t%α ri +δ =t)f (x ,t)α ri−

1τ[f (x ,tα ri)−fα(eq)(x ,tri)]
مرحله جاری شدن۱۲:
f (xα ri + δ +δ =erα t,tt)f (x ,t%α ri +δt)
(٨)
در روش شبکه بولتزمن، پارامتر زمان آرامش t طبق رابطه (٩)، از لزجت جنبشی سیال محاسبه می گردد: (٩) ν= τ−( 0.5)cs2δt
این روابط و فرضیات مرتبط با آنها باعـث مـی گـردد تـا مـدلBGK (رابطه (۷)) به یک روش با مرتبه دوم دقت، بـرای حـلجریان های تراکم ناپذیر سیال تبدیل شود [۲۶]. باید به این نکته توجه کرد که مرحله برخورد در معادله (۸) کاملاﹰ موضعی است و به داده های نقاط دیگر نیازی ندارد. مرحله جـاری شـدن نیـزفقط با نقاط همسایه تبادل دارد و لذا به محاسبه نیاز نـدارد . بـهدلایل فوق، این روش از قابلیت بالایی، بـرای پـردازش مـوازیبرخوردار است.

۴ – شیوه اعمال مرز جامد در روش شبکه بولتزمن
۴ -۱ – روش کمانه کردن
روش کمانه کردن برای اعمال شرط مرزی عدم لغـزش بـه کـارمی رود. در این روش، در کنـار دیـوار، مرحلـه برخـورد انجـاممی شود و در مرحله جار ی شدن به جای اینکه تابع توزیع در یک جهت را به گره بعد منتقل کننـد، در روی همـین نقطـه، جهـتورودی را برعکس می کنند؛ به عنوان مثال اگر در شکل ٢، تـابعتوزیع 1f قبل از انجام گام جاری شدن به سمت دیوار برود، در گام جاری شدن روی همین گره سیال، 2f برابـر 1f قـرار دادهمی شود. به عبارت دیگر، تابع توزیع به جهت دورشدن از دیوار، تغییر جهت میدهد و هیچ تغییری در نقطه جامد نیاز نیست. لازم به ذکر است که اگر دیوار متحرک باشد، علاوه بر تغییـر جهـ ت دادن تابع توزیع ورودی، مقدار خروجی نیز متناسب با سـرعتدیوار متفاوت خواهد بود و طبق رابطه (١٠) مقدار تابع توزیع در جهت خروجی معین می شود:
f (x ,t)α rf= f (x ,t)%α rf−6wα αρe .r urw (١٠)
در چگاایلـین راسیبالطـ، هreαα جهـجهـت تبـ رورعوکدیس،wαr ،wα ضـریب وزنـی، ρ u سـرعت حرکـت مـرزجامد و اندیس f مربوط به سیال است. لازم به ذکر اسـت کـهرابطه فوق، وقتی دیوار دقیقاﹰ وسط فاصله دو گـره قـرار داشـتهباشد، از مرتبه دوم دقت برخوردار است، در غیـر ایـ نصـورتمرتبه اول دقت را داراست. برای رفع این مشـکل و هـم چنـینبرای اعمال شرط مرزی کمانه کردن بـر روی مرزهـای منحنـی،روش کمانه کردن اصلاح شده در سال ١٩٩٨ توسـط فلیپـووا وهانل پیشنهاد شد [٢٧]. این روش در سال ١٩٩٩ توسـط مـی وهمکارانش با روابـط سـاده تـری بیـان گردیـد و تـاکنون مـورداستفاده قرار می گیرد [٢٨، ٢٩]. در روش کمانـه کـردن اصـلاحشده، برای یک مرز منحنی که داخل سـیال قـرار گرفتـه اسـت
(شکل ٣)، مقدار fa طبق رابطه (١١) محاسبه میشود [٢٨]:
f (x ,tα rf +∆ = −χt)(1)f (x ,t%α rf)
+χf (x ,t)α* rb − ρ6 w e .α αr urw (١١)
در رابطه فوق، پارامترهـای ذکرشـده بـ هصـورت ذیـل تعریـفمی شوند:
f (x ,t)α* rb=ρw [1α +3e .urα rbf + 4.5(e .urα rf )2 −1.5u .ur rff ]
(١٢)

شکل ۳ – تصویر یک شبکه یکنواخت دوبعدی و مرز منحنی داخل آن [۲۸]

414528114855

urbf = (∆−∆1)urf + ur∆w , χ=

2∆−τ 1 , ∆≥ 0.5 (١٣)
urbf = urff = u (xrf rff ,t) , χ=

, ∆<
r
310896-77498

∆ = xrf −−xrxrwb (١۵)
xf
پس از انجام گام جاری شدن مطابق روابط فوق، دیگر نیـاز بـهاعمال نیرو در نقاط شبکه نیست و شرط عدم لغزش و بـه تبـعآن اثر مرز جامد متحرک اعمال شده است.
مزایای روش کمانه کردن:
این روش، هم در مسائل جریان دائـم و هـم در مسـائلجریان گذرا قابل استفاده است.
به دلیـل اینکـه در ایـن روش، فقـط بـرای نقـاط سـیال،محاسبه و مبادله اطلاعـات وجـود دارد، تمـام نقـاط جامـد ازمحاسبات حذف می شوند و توابع توزیـع هـم در نقـاط جامـدمحاسبه نمی شوند، لذا زمان و هزینه محاسباتی به انـدازه سـهمنقاط جامد در بین کل نقاط کاهش می یابد.
عیوب روش کمانه کردن:
وجود میان یابی در این روش، موضعی بودن روش شبکه بولتزمن را تا حدودی از بین می برد.
بهدلیل اینکه برای محاسبه fa به روابـط (١١) تـا (١۵) نیاز است، استفاده از این روش در هندسه های پیچیـده راحـتنیست. با این وجود در صورتی که روابط مذکور به طور صحیح در یک هندسه اعمال شود، جـواب مسـئله سـریع تـر از روشمیدان نیرو به دست می آید زیرا به تکرار داخلی و ضـریب زیـر تخفیف١٣ نیازی نیست.

۴ -۲ – روش میدان نیرو یا روش مرز شناور
در این روش که هم برای مرز ساکن و هم بـرای مـرز متحـرکبه کار می رود، جسم جامد از درون دامنه حل حذف میگـردد واثر آن توسط یک سری نیـرو، در نقـاط ارتبـاط جامـد و سـیالاعمال می شود. ایـده اصـلی ایـن روش در سـال ۱۹۷۲ توسـطفردی به نام پسکین۱۴ برای شبیه سازی دریچه هـای قلـب مـورداستفاده قرار گرفت، اما این روش رسماﹰ در سـال ۱۹۹۳ توسـطگلداستین و همکارانش بـرای شـبیه سـازی شـرط مـرزی عـدملغزش مطـرح شـد [۳۰]. لازم بـه ذکـر اسـت کـه گلداسـتین وهمکارانش این ایده را در یـک کـارCFD مـورد اسـتفاده قـراردادند. در سال ۱۹۹۶ ایـن روش توسـط اگلـز در روش شـبکهبولتزمن مورد اسـتفاده قـرار گرفـت [۴]، و پـس از او در سـال۱۹۹۹، بـه بیـانی سـاده تـر، توسـط درکسـن و همکـارانش در شبیه سازی جریان داخل یک همـزن، بـه روش شـبکه بـولتزمنمورد استفاده قرار گرفت [۷].
روش کار نسبتاﹰ ساده است، امـا بـ هدلیـل اینکـه در مراجـعمذکور روش و فرمول بندی مورد نیاز خیلی فشرده عنوان شـدهاست، در اینجا به بیانی ساده برای یک دامنه دوبعدی، این روش توضیح داده می شود. همان طور که در اول این بخش عنوان شد، در روش میدان نیرو، جسـم جامـد از داخـل دامنـه محاسـباتیحذف می گردد و فقط اثر آن توسـط یـک سـری نیـرو اعمـالمی شود. برای این منظور ابتدا دامنه حل با شبکه مورد نیاز برای روش شــبکه بــولتزمن کــه یکنواخــت و مربعــی مــی باشــد، شبکه بندی می شود؛ سپس مرزهای جسم جامـد بـا یـک سـرینقاط مرجع، معین می شود و مختصـات ایـن نقـاط و شـمارنده

شکل ۴ – تصویر یک نقطه مرجع جامد بههمراه نقاط اطراف آن در
شبکه بولتزمن

هر نقطه ذخیره می گردد. مکان این نقاط مرجع بایـد روی مـرزجسـم جامـد بـا یـک سـری نقـاط مرجـع، معـین مـی شـود ومختصات این نقاط و شمارنده هر نقطه ذخیره می گـردد . مکـانایـن نقـاط مرجـع بایـد روی مـرز جسـم جامـد باشـند، امــامحدودیتی از لحاظ انطباق با نقاط شـبکه، و یـا یکسـان بـودنفاصله ندارند و کاملاﹰ آزادانه انتخاب مـی شـوند . مـثلاﹰ مـ یتـوانبرای نقاط حساس جسم جامـد، تعـداد نقـاط مرجـع بیشـتر ونزدیک تری را انتخاب نمود. فرض کنید نقطـهA در شـکل (۴) یک نقطه مرجع از جسم جامد باشد و نقـاط اطـراف آن، نقـاطمربوط به شبکه بندی شبکه بولتزمن باشـند . در ایـن روش ابتـدااز روی تابع توزیـع هـر نقطـه، سـرعت نقـاط شـبکه محاسـبهمی شود و با توجه به فاصـله نقطـهA از گـره هـای اطـراف آن،مقدار سرعت در نقطه A توسط میان یابی پیدا می شود. با توجـهبه شکل (۴) و نامگذاری های آن، مقدار سـرعت در نقطـهA از رابطه (۱۶) به دست می آید:
r
uA = −α(1)(1−β)urws + −α β(1) urwn
+αβuren +α −β(1)ures
حال با توجه به سرعت واقعی که در نقطه A مد نظر می باشـد،انحراف سرعت نقطه A از سـرعت واقعـی، طبـق رابطـه (۱۷) به دست می آید:
∆urA = urA,real − urA
اگر به نقاط اطراف نقطه A نیروهای مناسبی وارد شود، می توان انحراف سرعت در ایـن نقطـه را برطـرف نمـود. لـذا انحـرافسرعت در نقطه A به گره های اطراف برون یـابی مـی شـود و ازروی انحراف سرعت در هر نقطه شبکه، نیروی مورد نیاز بـرایاصلاح این انحراف به دست می آید. بنابراین مقدار نیروی مـوردنیاز در هر نقطه شبکه، طبق رابطه (۱۸) به دست میآید:
r Fws = (1-α)(1- )β ρ∆urA r
Fwn = (1-α βρ∆)urA
Fren =αβρ∆urA (١٨)
rFes =α(1- )β ρ∆urAلازم به ذکر است که این مقادیر فقـط ناشـی از نقطـه مرجـعA می باشند و برای بقیه نقاط مرجع هم، باید اثر انحراف سرعت ها محاسبه شده و به مقادیر فوق افزوده شود. لـذا رابطـه (۱۸) بـهرابطه (۱۹) تبدیل می شود که در آن مقـادیر نیـروی مربـوط بـهنقاط مرجع قبلی نیز حذف نمی شوند:
rr
Fws = Fws + −α(1)(1−β ρ∆)urA
rr
Fwn = Fwn + −α βρ∆(1)urA
Fren = Fren +αβρ∆urA
Fres = Fres +α −β ρ∆(1 ) urAبا توجه به اینکه تغییر شـدید در نیـرو باعـث واگرایـی برنامـهمی شود، در اینجا از ترکیب نیروی حاصل با نیروی گـام زمـانیقبلی، به همراه یک ضریب زیرتخفیف (h ) استفاده میشود [٧].
یعنی داریم:
r
Ft = +ηFr(Frt−∆t −F)r
که در آن Fr مقدار محاسبه شـده از رابطـه (١٩) اسـت . بـرایاعمال نیروی به دست آمده، این نیرو به عنوان یک نیروی حجمی متغیر درنظر گرفته شـده و مقـدار آن در نقـاط مختلـف شـبکهاعمال می شود. بدیهی است که وجود ضریب زیرتخفیـف، اثـردیواره جامد را با تأخیر در داخل دامنه حل وارد می کند، لذا این روش برای مسائل جریان دائم مفید میباشد. لازم بهذکـر اسـتکه برای افزایش دقت روش میدان نیرو، مـی تـوان در هـر گـامزمانی، چند بار تکرار داخلی انجام داد و نیروها را اصلاح نمود. با این کار، اثر مرز جامد نیز سریع تر در دامنه حل اعمال می شود [٧]. به عنوان مثال در مرجع [١۵] برای افزایش دقت شبیه سازی، در هر گام زمانی ده تکرار داخلی انجام شده است.
مزایای روش میدان نیرو:
این روش به سادگی در برنامـه کـامپیوتری قابـل اعمـالاست.
برنامه را برای هندسه های مختلف می تـوان بـه کـار بـرد. برای این کار فقط نقاط مرجع جامد را باید عوض کرد و نیازی به تغییر در برنامه نیست.
عیوب روش میدان نیرو:
وجود میان یابی در این روش، موضعی بودن روش شـبکهبولتزمن را تا حدودی از بین می برد.
وجود ضریب زیرتخفیف، کاربرد روش میدان نیرو را به حل جریان های دائم محدود می کند یا در صورت نیاز بـه حـلجریان های غیردائم به برنامه تکرار داخلی اضافه می شود.
وجود تکرار داخلی، زمان مـورد نیـاز بـرای رسـیدن بـهجواب را افزایش می دهد.
به دلیل اینکه در این روش نقاط جامد، با سـیال جـایگزینمی شود و روابط مربـوط بـه نقـاط سـیال در ایـن نقـاط اعمـالمی گردد، میزان محاسبات افـزایش مـی یابـد . بـرای اینکـه میـزانافزایش حجم محاسبات مربوط به این مشـکل بهتـر درک شـود،یک همزن استوانه ای به قطر و ارتفاع D درنظر بگیرید. با توجـهبه اینکـه در روش شـبکه بـولتزمن از شـبکه یکنواخـت مربعـیاستفاده می شود، برای شبیه سازی این همزن، دامنـه ای بـه شـکلمکعب بر آن محیط می شود. این مکعـب اضـلاعی بـه انـدازهD دارد و کل این حجم در محاسبات وارد می شود. از آنجا که نقاط بیرون استوانه در مسئله اصلی وجود نداشته، بهاندازه حجـم ایـننقاط اضافی به محاسبات افزوده شده است. بـا توجـه بـه اینکـهحجم مکعب 3D ، حجم استوانه 4

3πD و اختلاف حجم آنهـا4 3(4−π)D / می باشد، میزان محاسبات، بیست و هفت درصـداز یک استوانه تنها بیشتر شده است. لازم به ذکـر اسـت کـه ایـنمقدار غیر از حجم محور، پره های متحرک و ثابت می باشد، زیـرابهجای این اجسام جامد نیز سیال قرار می گیـرد و در محاسـباتبه عنوان سـیال وارد مـی شـون د. البتـه مـی تـوان ترکیبـی از روشکمانه کردن و روش میدان نیرو را بـ ه کـار بـرد و ایـن مشـکل راکاهش داد. به این ترتیب که در دیواره های همـزن، روش کمانـه کردن و برای پره هـای ثابـت، محـور و پـره هـای متحـرک روشمیدان نیرو را اعمال نمود. با این کار، حجم اضافی بیرون ظـرفحذف می شود و فقط نقاط حاوی پره های ثابت، پره های متحرک و محور به شکل سیال در محاسبات وارد می شوند.
در کار حاضر بـا توجـه بـه نکـاتی کـه در مـورد دو روشکمانه کردن اصلاح شده و میـدان نیـرو عنـوان شـد، و داد ههـایجدول ۱ که دقیق تر بـودن روش کمانـه کـردن اصـلاح شـده رانسبت به روش میدان نیرو تأییـد مـی کنـد، تـا زمـانی کـه روشکمانه کردن اصلاح شده قابـل اعمـال باشـد، ایـن روش نسـبتبه روش میدان نیرو در اعمال مرز جامد ترجیح داده میشود. لذا در کار حاضر در کنار روش شبکه بولتزمن از روش کمانه کردن اصلاح شده برای شبیه سازی مرزهای جامد بهره گرفته میشـود .

  • 1

پاسخ دهید