د و در این فاصله زمـانی ، تغییـراتعناصر و نوترو نهای تأخیری با یک شار ثابت به دست آمد.

با انرژی ۱ که به صورت عمود وارد قلب می شوند) فـرضشده است کـه در آن 0n .D∇ϕ=J مـی باشـد . در ناحیـ ه دوم شرط مرزی به صورت ALBEDO (رابطه ۱۵) درنظر گرفته شده است [۱۷].
شرایط مرزی و اولیه درنظر گرفته شده به صورت روابط زیر می باشند: پخش – واکنش مربوطه شرح داده شد. در این قسمت بـه زیـرسیستم پخش پرداخته خواهد شد. این زیر سیستم یـک معادلـه دیفرانسیل با مشتقات جزئی از مرتبه دو است. حل ایـن معادلـهرا با گسسته سازی فضایی۱۸ شـرو ع کـرده و بـرای ایـن کـار ازروش شناخته شده المان محدود۱۹ اسـتفاده مـی کنـیم . بـر طبـق
: ( )
3129281660104

27675841768210

در قسمت قبـل روش اتخـاذ شـده بـرای حـل زمـانی سیسـتم ۱). در ناحیه MeVاول شرط مرزی به صورت شار ثابت (نوترون هایی
رابطه۲داریم
 n .D

.Dq (۹)
ϕ(r,t)y 0,t 0= = = J , 0ϕ(r,t)y 0,t 0≠ = = 0 (۱۶) برای راحتی کار q را به صورت زیر بازنویسی می کنیم: (۱۰) q = ϕ+a b
1106424110957

که در آن a و b به راحتی از رابطه(۸) بهدست می آیند. بـه ایـن (۱۷) N8(r,t)

ρ NA
۳ -۲ – حل زیر سیستم پخش
ترتیب معادله پخش به صورت زیر در خواهد آمد:

.Da b (۱۱)
با توجه به روند کلی روش المان محدود، دو طرف رابطه (۱۱) را در یک تابع تست ضرب کرده و روی فضای مربوط به مسئله انتگرال گیری می کنیم. بدین ترتیب معادله پخش بـه رابطـه زیـرتبدیل خواهد شد:

(۱۲)
∫ v a dXϕ+∫ v b dX
ΩΩ
که در آن X = (x,y,z) و Ω فضای موردنظر است. با توجـهبه توابع گرین۲۰ داریم:

و هم چنین داریم [۱۶]:
D ∂ϕ∂n = n .D∇ϕ (۱۴)
که در این دو رابطه ∂Ω مرز سیسـتم وn۲۱
بـردار نرمـال آن می باشد. با توجه به هندسه راکتور و شرایط مرزی خاص درنظر گرفته شده برای مسئله، ∂Ω به دو ناحیه تقسیم می شود. ناحیه اول (1∂Ω) صفحه دایره ای شکل منتهی الیه سمت چپ قلب و ناحیه دوم (2∂Ω) مابقی مرزهای خارجی سیستم است (شـکل

شکل ۲ – نمای یک چهارم قلب

N9(r,t)N(r,t) =t 0 ==t 0= 00 , N (r,t , NPu(r,ti) =t 0)= =t 0= 0 0 (۱۸)
که در آن 22α= JJ−+∂Ω∂Ω اسـت . 2J−∂Ω جریـان جزئـی خـارج
شونده از یـک سـطح کوچـک مـرز و 2J+∂Ω جریـان جزئـیمنعکس شده (داخل شونده) از یک سطح کوچک مرز خـارجیبه طرف داخل قلب است.
با توجـه بـه تقـارن موجـود در هندسـه قلـب و هـم چنـین تغییرات مکانی شار، یک چهارم از قلب در شبیه سـازی درنظـر گرفته شده است (شکل ۲). این عمل به نوبه خود باعث ایجـادنوعی دیگر از شرط مرزی می شود که در محاسبات عـددی بـهشرط مرزی آزاد۲۲ معروف است [۱۸]. خوشبختانه برای اعمـالاین نوع شرط مرزی نیازی به انجام تغییرات در ماتریس ضرائب نمی باشد. هنگام استفاده از این نوع شرط مرزی بایـد از وجـودتقارن اطمینان حاصل شود، در غیر این صورت جواب مناسـببه دست نخواهد آمد. با ترکیب روابـط (۱۳، ۱۴، ۱۵) و معادلـ ه
پخش و هم چنین اضافه کردن شرط مرزی ناحیه اول، به رابطـ ه زیر خواهیم رسید:
∂ϕ−α
∫ v J dX0+ ∫ v a ϕdX + ∫ v b dX
∂Ω1ΩΩ
(۱۹)
در ایـن رابطـه بـه دنبـال جـوابی از ϕ هسـتیم کـه در فض ای سوبولف۲۳ از مرتبه یک H (10 Ω) صـدق کنـد. یعنـیϕ بایـدروی این فضا پیوسته بوده و مشتق آن پیوسته تکه ای باشد[۱۶].
در ادامهϕ را با مقـدار تقریبـی آن (فـرم ضـعیف روش المـانمحدود) جایگزین می کنیم:
N
(۲۰)
با یکسان فرض کـردن شـکل فضـاییϕ (تـابع وزن در روشباقیمانده وزن دار شده۲۴) و تابع تست v طبق روش گالرکین۲۵
[۱۶]، معادله پخش به شکل زیر در خواهد آمد:
∑N ∫ϕ ϕ dX dΦj =−∑N ∫∇ϕ ∇ϕ DdXΦ −

N
∑∫ϕi a i ϕj dXΦ + ϕj ∫ i b dXi i =1,….,N
j 1= ΩΩ
(۲۱)
که در آن N تعداد گـره هـا ی۲۶ درنظـر گرفتـه شـده (یـا تعـدادمعادلات موجود در دسـتگاه کلـی) در فضـای قلـب در جهـتمش بندی هندسه موردنظر است. با فرض:
∫ϕ ϕi j dX = Mij
Ω (۲۲)
∫∇ϕ ∇ϕi Dj dX =Kij
Ω (۲۳)
∫ϕi b dXi= b

i
Ω (۲۴)

∫ϕi a i ϕj dX = aij (۲۵)

سیستم گسسته شده فضایی زیر به دست خواهد آمد:
1210056428842

∑j 1N= M ij

ddtΦj =−∑j 1N= K ij Φ −j1 121+α−α∑j 1N= Mij∂Ω2 Φ +j N
∫ ϕi J0i dX +∑a ij Φ +jbi i =1,…,N
∂Ω1j 1=
(۲۶)
نقطه اثر قسمت های مربوط به شرایط مرزی کـه در اثـر اعمـالتوابع گرین ایجاد شده اند، بستگی به اندیس های اختصاص داده شده به گره ها و المان های مرزی خواهد داشـت . شـرط مـرزیناحیه اول با استفاده از روش پنالتی۲۷ و شرط مرزی ناحیـ ه دوم با محاسبه یک سری انتگرال های بـرداری سـه بعـدی و انجـامتصحیحات لازم در ماتریس ضرائب اعمال خواهـد شـد [۱۸]. فعلاﹰ برای ساده بودن شکل معادله و اینکه بتوانیم مسئله زمان و نحوه برخورد با آن را ساده تر بیان کنـیم، ایـن دو قسـمت را ازمعادله حذف می کنیم. بـا بسـط رابطـه (۲۶) روی انـدیسi در نهایت سیستم ماتریس گونه زیر را خواهیم داشت:
[M].{ddtΦ}+[K].{ }Φ ={a} .{ } {b

T Φ +

} (۲۷)
که در آن [ ]، { } و T به ترتیب بیانگر ماتریس، بردار و ترانهاده می باشند. با بهکارگیری روش تتا۲۸ می توان شار در لحظه جدید را با استفاده از متوسط وزن دار شده دو گرادیـان زمـانی متـوالیشار به دست آورد [۱۹]. طبق ایـن روش بـرای شـار در لحظـه جدید خواهیم داشت:
{ }Φ = Φ +∆1{ }0t((1−θ) {ddtΦ}0 +θ{ddtΦ}1) (۲۸)
گرادیان های زمانی موجود در این رابطه با استفاده از رابطه (۲۷) قابل محاسبه هستند:
1408176-501

2036064-501

[M].{ddtΦ}0 +[K].{ }Φ =0{a} .{ }T Φ +0{b} (۲۹)
[M].{ddtΦ}1 +[K].{ }Φ =1{a} .{ }T Φ +1{b} (۳۰)
در نهایت سیستم گسسته شده کامل به دست خواهد آمد:
([M]+θ∆t[K]).{ }Φ =1([M]−∆ −θt(1) [K]).{ }Φ +0
643128-1147

1240536-1147

2057400-1147

2673096-1147

θ∆t({a} .{ }T Φ +1{b})+∆t (1−θ)({a} .{ }T Φ +0{b})
(۳۱)
تغییرات θ از صفر تا یـک مـا را از روش صـریح کامـل۲۹ بـهروش ضمنی کامل۳۰ هدایت خواهد کرد [۱۹]. بهترین انتخـاببا قـرار دادن .05θ= و رسـیدن بـه روش معـروف کرانـک – نیکلسون۳۱ حاصل خواهد شد. این روش به لحـاظ پا یـداری ومسائل مربوط به همگرائی، شرایط خوبی را برای مسـ ئله ایجـادمی کند. انتخاب این مقدار برای θ ، ظاهراﹰ میـزان محاسـبات راافزایش خواهد داد ولی با این انتخاب مرتبه دقت بالاتر رفتـه وتعداد تکرارها در جهت رسیدن به محدوده خطای۳۲ موردنظر و تأمین شدن معیار توقف۳۳ تکرارکننده کمتر خواهد شد [۱۵].
در حالت پایای سیستم، استفاده از روش صریح کامل یعنـی0θ= نیز با توجه به t∆ ی کوچک درنظر گرفتـه شـده، دقـتمناسبی را در پی خواهد داشت. این موضـوع حالـت غیرخطـیحاصل شده از روش تتا برای /0 5θ= را از بین خواهـد بـرد. موارد مربوط به تحلیل خطا و تعیین پله زمانی در مراجـع [۱۵] و [۱۹] آمده است که برای طولانی نشدن مطلب از طـرح آنهـاخودداری می شود. با مشخص شدن θ و پلـ ه زمـانی، سیسـتمخطی آشنای زیر به دست خواهد آمد:
A . X = B (۳۲)
این سیستم را می توان بـا اسـتفاده از روش هـایی کـه براسـاسزیرفضای کرایلف۳۴ شکل گرفته اند، به مراتب سری عتر، بهینه تر و دقیق تر از رو شهای سنتی حل کرد. از میان آنهـا مـی تـوان بـهروش های Lanczos ،Arnoldi ،GMRES و CG اشاره کرد. در این روش ها به جای تشکیل سیستم کلی ماتریس ضرائب و حل کردن مستقیم آن، تنها از یک سری ضرب های بردار در ماتریس استفاده می شود. این روش ها خصوصاﹰ وقتی که با مسائل بزرگ روبرو هستیم مؤثرتر و دقیق تر عمل خواهند کرد و میزان حافظ ه و محاسبات مورد نیاز آنها کمتر خواهد بود.
در این تحقیق بهطـور خـاص از روشBiCGStab(L) کـهنمونه توسعه یافته روش CG می باشـد و مشـکلات اولیـه ایـنروش را به همراه ندارد [۲۰] استفاده شده است. از نکـات مهـمدر مورد این روش، ساده تر بودن اجـرای مـوازی آن نسـبت بـهروش های سنتی مانند گوس – جردن جهـت انجـام محاسـباتموازی می باشـد . وقتـی از روش BiCGStab(L) در کنـار روشالمان محدود استفاده شود، حل سیسـتم عنـوان شـده در رابطـه (۳۲) به انجام یـک سـری حلقـه هـایfor و ضـرب بـردار درماتریس مربوط به تک المان ها تبـدیل خواهـد شـد [۱۸]. ایـنحلقه ها چه در سـطح حافظـه تسـهیم شـده۳۵ و چـه در سـطححافظه توزیع شده۳۶ به راحتی قابل موازی سازی هستند. تعیـینمقدار L در روش BiCGStab(L) به نـوع مسـئله بسـتگی دارد .
به طورکلی با افزایش L میزان حافظه مورد نیاز افزایش و میـزانمحاسبات برای رسیدن به یک دقت مشـخص کـاهش خواهـدیافت [۲۰]. در این تحقیق مقدار L برابر ۴ در نظر گرفتـه شـدهاست. برای برآورده شـدن تلـورانس در محاسـب ه خطـا، تعـدادتکرارهـای خـارجی (در مقابـل عـدد L کـه در ایـن روش بـه تکرارهای داخلی معروف است) ایـن روش در حـدود ۴ تـا ۵ باقی ماند. روش های مبتنی بر زیرفضای کرایلـف نیـاز بـه یـکحدس اولیه برای جواب دارند و هرچه این حدس بهتر انتخاب شود تکرارها در جهت رسیدن به محدوده خطـای مجـاز کمتـرخواهد شد. در این تحقیق از پیش حلگر قطری۳۷ استفاده شـدهاست.

  • 2