بررسي شک لگيري موج خودنگهدار شکافان هسته اي در يک راکتور سريع بحراني

*
محسن حسامي و مهدي نصري نصرآبادي
گروه مهندسي هست هاي، دانشکده علوم و فناوري هاي نوين، دانشگاه اصفهان

(دريافت مقاله: ١٥/٠٥/١٣٩١ – دريافت نسخه نهايي: ٢٧/٠٧/١٣٩٢)

چكيده – اين مقاله به بررسي پديده اي به نام موج پاياي شکافان هست هاي که با برقراري سيکلU-Pu در يک راکتـور هسـت هاي سـريع ايجـادمي شود، مي پردازد. از موارد جالب توجه در مورد اين نوع راکتور، ايمني حاصل از حالت خودبخودي موج و همچنـين بـازد ة حـدود ۵۰ درصـديمصرف سوخت آن است که در راکتورهاي حرارتي امروزي در حدود ۱ تا ۲ درصد مي باشد. جهت بررسي امکان تشکيل اين موج، سيستم حاصل از معادلة پخش نوترون که با معادلات مصرف سوخت و سينتيک کوپل شده است درنظر گرفته شد. روش خطي سازي به کار گرفته شده در حل ايـنسيستم نه تنها زيرسيستم هاي پخش و واکنش را از يکديگر جدا مي کند بلکه مسئله روبرو شدن با چند صد روز زمان واقعي شبي هسازي را نيـزبسيار ساده تر مي سازد. گسسته سازي فضائي و زماني زيرسيستم خطي شدة پخش به ترتيب با استفاده از روش المان محدود سـه بعـدي و روشکرانک – نيکلسون و نهايتاﹰ همگرائي آن به کمک تکرارهاي مبتني بر روش بهينة BiCGStab(L) انجام شد. روش المان به المان به کار گرفته شده در کنار BiCGStab(L) حافظة مورد نياز را به اندازة يک دهم مقدار حافظة استفاده شده در غياب اين روش کاهش داد. با توجه به تغييرات انـدکغلظت عناصر موجود در زيرسيستم واکنش، روش دقيق رانگ – کوتا (۵) براي حل معادلات ديفرانسيل مرتبه اول ايـن زيـر سيسـتم اسـتفادهگرديد. با توجه به زمان طولاني شبيه سازي، سه بعدي بودن مسئله و محاسبات سنگين حاصل از آن، از روش OpenMp جهت موازي کردن برنامة نوشته شده به زبان فرترن ۹۵ استفاده شد.

واژگان كليدي : موج خودنگهدار شکافان هسته اي، سيستم نوترونيک، خطي سازي، BiCGStab(L)، روش المان محدود سه بعدي المان به المان، .OpenMp.

Analysis of Self-Sustained Nuclear Fission Wave Evolution in
a Fast Critical Reactor

M. Hesami and M. N. Nasrabadi

Department of Nuclear Engineering, Faculty of Advanced Sciences & Technologies,
University of Isfahan, Isfahan, Iran

Abstract: Present paper deals with the phenomenon known as self-sustained nuclear fission wave that is set up with establishment of the U-Pu cycle in a fast nuclear reactor. The safety arised from self-sustained state of the wave and also 50% the possibility of initiation and evolution of the wave, the system comprised of neutron diffusion equation coupled with burn-up and kinetic equations have been considered. The linearization method used to solve this system not only separates diffusion subsystem
* : مسئول مكاتبات، پست الكترونيكي: mnnasrabadi@ast.ui.ac.ir
from reaction one but also simplifies greately dealing with several real hundred days of simulation. Space and time discretizations of the linearized diffusion subsystem were performed by 3D finite element and crank-nicolson methods, respectively. The convergence of equation was carried out using the iterations based on the optimized BiCGStab(L) method. Rung-Kutta (5) method was used to solve the first order differential equations of the reaction subsystem since the concentrations of elements in the reaction subsystem changed slowly. Using the element by element FEM and BiCGStab(L) used for solving the system dropped the required memory down to 1/10. Considering the long period of the simulation, 3D case of the problem and heavy computations, the OpenMp method was used to parallelize the code written in FORTRAN 95.

Keywords: Nuclear reactor, self-sustained nuclear fission wave, neutronic system, linearization, BiCGStab(L), 3D element by element FEM, OpenMp.

شار نوترون (n/cm2.s) ϕ ضريب پخش نوترون(1-cm) D
ثابت واپاشي دختر هسته هاي شش گروهي (1-s) λ
i چگالي اورانيوم ٢٣٨ (3n/cm) N۸
سرعت نوترون (m/s) ν چگالي اورانيوم ٢٣٩ (3n/cm) N۹
سطح مقطع جذب ميکروسکوپيک اورانيوم ٢٣٨ (2cm) σ8 a چگالي پلوتونيوم ٢٣٩ (3n/cm) NPu
سطح مقطع جذب ميکروسکوپيک پلوتونيم ٢٣٩ (2cm) σPu
a چگالي تک گروه مؤثر پاره هاي شکافت (3n/cm)

N
سطح مقطع شکافت ميکروسکوپيک پلوتونيم ٢٣٩ (2cm) σfPu چگالي منبع نوترون (n/cm3.s) q
نيمه عمر مؤثر واپاشي اورانيوم ٢٣٩ از طريق دو کاهش بتا (s) τ زمان (s) t
فهرست علائم

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

Terra Power به پشتيباني مالي بيل گيتس روي آن انجام شـود .
اين گروه در سال ۲۰۱۰ طي مقاله اي نتايج حاصل از تحقيقـاتخود را که با استفاده از روش مونت کارلو انجام شده بود منتشر کردند. آنها در اين مقاله حداقل ميزان سـوزش سـوخت جهـتشکل گيري موج را ۱۴% به دست آورده اند [۴].
خيلي قبل تر از عنوان شدن چنين رژيـم کـاري بـراي يـکراکتور هسته اي، افرادي نظيـر فينبـرگ۷ و دريسـکول ۸ مسـائليمانند استفاده از سوخت طبيعي (غني نشده) و هـم چنـين نقـشزايش و تبديل در قلب راکتـور را بررسـي نمـوده بودنـد ولـياحتمالا فوکتيستوف۹ اولين کسي است که موج پايـاي شـکافانهسته اي بر مبناي زايش، تبـديل و سـوزش را مطـرح نمـوده وشرايط لازم جهت شـکل گيـري آن را مـورد بررسـي قـرار دادهاست [۹ -۸]. افرادي مانند سکيموتو۱۰، فومين۱۱ و ونـدام ۱۲ نيـزپژوهش هايي در اين زمينه انجام داده انـد . سـکيموتو بـا فـرضثابت بـودن مـوج در طـول قلـب و بـا درنظـر گـرفتن مسـائلترموهيدروليک، راکتوري بـر مبنـاي ايـن مـوج طراحـي کـرد ه است [۲]. فومين با درنظـر گـرفتن معـادلات چنـد گروهـي در حالت دو بعدي شکل گيري موج را بررسـي کـرده اسـت [۱۰].
وندام نيز بـا تکيـه بـر اصـول رياضـي و حضـور يـک فـاکتورفيدبک، جوابي تحليلي براي موج به دسـت آورده و بـه بررسـيخصوصياتي نظير سرعت موج و عـرض فعـال مـوج پرداختـهاست [۱۲ -۱۱].
براي بررسي چنـين مـوجي لازم اسـت کـه معادلـة ترابـردنوترون، معادلات سوزش و سينتيک۱۳ به صورت کوپل شـده بـايکديگر حل شوند. در ابتداي کار اين راکتور، جهت راه انـدازيموج از يک چشمة خارجي استفاده مـي شـود . چشـم ة نـوترونخارجي، اورانيوم ۲۳۸ موجود درون قلب راکتور را طبق رابطـ ة زير به پلوتونيوم ۲۳۹ تبديل مي کند:
238 U(n, )γ fi 23992 U23.5β-minfi
92
23993 Np fi2.35β-day23994 Pu(n, fission) (۱)
با توجه به نحوة توزيع اولية شار حاصل از چشمة نـوترونخارجي، به دليل انـرژي نسـبتﹰا بـالاي نـوترون هـاي ورودي، آن نواحي از قلب که سرشار از اورانيوم ۲۳۸ و نزديک تـر بـه مـرزچشمه مي باشند، بيشتر به پلوتونيوم ۲۳۹ تبديل خواهند شد. در اين ناحيه عمل شکافت نيز انجام مي گيـرد ولـي فرآينـد غالـبتوليد پلوتونيوم ۲۳۹ خواهد بود. با گذشت زمان و انباشته شدن پلوتونيوم ۲۳۹ طبق شرايط لازم و کافي کـه فوکتيسـتوف بيـانکرده است، يک ناحيـ ة بحرانـي موضـعي در سـمت چشـمه ونزديک آن ايجاد خواهـد شـد. مقـداري جلـوتر از ايـن ناحيـ ة بحراني به طرف داخل قلب شرايط زيـر بحرانـي بـوده و دارايضريب تکثير۱۴ کوچک تـر از يـک اسـت . ايـن ناحيـة بحرانـيموضعي طبق شرايط گفته شده توسط فوکتيستوف مـي توانـد در حکم چشمه اي براي ناحية زير بحراني روبروي خود عمل کرده و به اين صورت موج شکافان هسته اي موردنظر ايجـاد خواهـدشد. با توجه به استراتژي گفته شده، اين موج به صورت خود به خود به سمت جلو حرکت خواهد کرد و به هيچ گونـه دخالـتخارجي نياز نخواهد داشت. امتياز تئوري فوکتيستوف در مـورداين پديده در همين خودکار بودن موج و عـدم نيـاز بـه کنتـرلخارجي جهت ادامة حرکت آن است.
۲ – مدل نوترونيک
همواره در بررسي يک پديدة نوتروني جديد سعي مـي شـود ازمدل هاي سادة نوتروني ولي به اندازة کافي دقيق و مؤثر استفاده شود. از اين جهت در اين تحقيق از معادلة پخش تـک گروهـيبراي بررسي پديده اي به نام موج پاياي شکافان هسته اي استفاده شده است. صورت کلي اين معادله به صورت زير است:

.Dq (۲)
که در آن q چگالي منبع نوترون، D ضريب پخش و v سـرعتنوترون هاي تک گروهي اسـ ت. شـکلq بـا توجـه بـه سـيکلاســتفاده شــده، عناصــر شــرکت کننــده در ســوزش، دختــرهسته ها۱۵، نحوة برخورد با نوترون هاي تـ أخيري و… مـي توانـدشکل هاي متفاوتي به خود بگيـر د. در تحقيقـات گذشـته [۴، ۶، ۱۰] دو ســيکل معــروفU-Pu و Th-U اســتفاده شــده انــد.
به طورکلي شکل گيري موج با به کار گيري سيکل U-Pu آسان تـرانجام مي شود. اين به خاطر سطح مقطع هاي مناسب تر اورانيـوم۲۳۸ در اين رژيم کـاري، کوتـاه تـر بـودن دورة زمـاني تبـديلاورانيوم ۲۳۸ به پلوتونيوم ۲۳۹ به نسبت يک دهـم دورة زمـانيتبديل توريم ۲۳۲ به اورانيوم ۲۳۳ و امکان رسيدن بـه چگـاليسوخت عملي بالاتر با استفاده از اورانيوم ۲۳۸ است. هم چنـين پسماند موجود راکتورهاي هسته اي فعلي عمـدتاﹰ اورانيـوم ۲۳۸ مي باشد که با توجـ ه بـه قابـل اسـتفاده بـودن آن در ايـن نـوعراکتور، مناسب تر است که تحقيقات روي سـيکل کـاريU-Pu انجام شود [۴].
اگر به اين سيکل و پارامترهـاي زمـاني آن توجـه کنـيم در خواهيم يافت که با توجه به کوتاه بودن زمان حضـور اورانيـوم۲۳۹ و هم چنين نزديک بودن خصوصيات اين عنصر با نپتونيوم
۲۳۹، جابه جا کردن اين دو با يک عنصر، نه تنهـا خطـاي قابـلتوجهي در محاسبات به وجود نخواهد آورد، بلکه طول محاسبات نيز کمتر خواهد شد. تعداد دختر هسته هاي ايجاد شده از شکافت پلوتونيوم ۲۳۹ نسبتاﹰ زياد مي باشد؛ درنظر گرفتن تمام آنها، هم به جهت افزايش ميزان محاسبات و هم به جهت اينکه بعضي از آنها احتمال شکل گيري کمي داشته و نقش ناچيزي در مسئله دارند، معقول و اقتصادي به نظر نمي رسد. به اين ترتيب بـراي بررسـينقش پاره هاي شکافت و هم چنين نقش نوترون هاي تـ أخيري از شش گروه پاره هاي شکافت استفاده شده است و مـابقي دختـرهسته ها به صورت يک گروه مؤثر در معادلات مربوط به سوزش وارد شده اند. با توجه به موارد بيان شـده معـادلات مربـوط بـهسوزش به صورت زير خواهند بود:
∂∂N8t =−ϕσ 8a N8 (۳)

97536-259439

∂∂N9t =ϕσ 8a N8− 1

τ N9 (۴)
∂NPu∂t = 1τ N9−ϕ σ +σ ( PuaPuf ) NPu (۵)
با فرض اينکه شکافت به دو پارة شکافت، حداکثر احتمال را داشته باشد، براي يک گروه مؤثر پاره هاي شکافت خواهيم داشت:

ifNPui N– i (۶)
و در نهايت معادلات مربوط به سينتيک نيز بهصورت زير در خواهند آمد:
9753696189

∂∂N–ti =β ϕσi fPu NPu −λi N– i i =1,…,6 (۷)

در معادلات فوق ۸N ،NPu ، N۹ ، N و N– i به ترتيب غلظـتاورانـيم ۲۳۸، اورانـيم ۲۳۹، پلوتونيـوم ۲۳۹، يـک گـروه مـؤثر پاره هاي شکافت و دختر هسته هاي شش گروهي مـي باشـند . τ زمان مؤثر تبديل اورانيوم ۲۳۹ از طريق کاهش بتا به پلوتونيـوم۲۳۹ است. سطح مقطع هاي σaPu ، σ8a وσfPu به ترتيـب سـطحمقطع هاي جذب اورانيوم ۲۳۸، پلوتونيـوم ۲۳۹ و سـطح مقطـعشکافت پلوتونيوم ۲۳۹ مي باشند. βi کسر نوترون هاي تأخيري و λi ثابت واپاشي دختر هسته هاي شـش گروهـي هسـتند . در نهايت شکل qي ظاهر شده در معادلـ ة پخـش بـه صـورت زيـرخواهد بود:
6 q = υ −β ϕσ[ (1)] fPu NPu +∑λi N– i −ϕ×
(۸)
283464-297465

[a Nia N N]Sext
2270760-16787

در اين رابطه σ سطح مقطع مؤثر درنظر گرفته شده براي مابقي پاره هاي شکافت، Sext منبع نوترون خارجي وυ تعداد متوسـطنـوترون هـاي آزاد شـده بـه ازاي هـر شـکافت پلوتونيـوم ۲۳۹ مي باشد.

۳ – حل سيستم معادلات
۳ -۱ – استراتژي کلي
سيستم تشکيل شده از معادلات به دست آمده در قسمت قبل در واقع يک سيستم پخش – واکنش۱۶ اسـت . بـا توجـه بـهشـکل q و ميـزان پيچيـدگي آن، راه ح ل هـاي متف اوتي را مي توان براي حل آن در پـيش گرفـت [۱۳]. معادلـ ة پخـشنوترون از نوشتن قانون پيوسـتگي بـراي نـوترون بـه دسـتمي آيد. اگر به معادلة پخش دقت شود (رابطة ۲)، براي اينکه معادله فقط برحسب شار نوترون باشد، طبق رابطـ ة شـار بـانوترون (در حالت تک انرژي) سـمت چـپ ايـن معادلـه درمعکوس سرعت نوترون ها ضرب شده اسـت [۱۴]. سـرعتنوترون ها حتي در سطح انرژي حرارتـي آنهـا عـدد بزرگـياست. اين ضريب حساسيت سيسـتم را بـالا بـرده و هنگـامبررسي همگرايي سيستم بهناچـار بايسـتي پلـه هـاي زمـانيکوچک انتخاب شوند. در غير اين صـورت آشـفتگي۱۷ وارد شده به سيستم زياد بوده و تکرارها جهت همگرايي سيسـتمبسيار زياد خواهند شد و يا سيستم اصلاﹰ همگرا نخواهد شد
.[۱۵]
نکتة مهم اين است که شار نوترون در حالت پاياي سيسـتميعني وقتي که تغييرات شار فقط در اثـر تغييـر غلظـت عناصـرباشد، تغييرات اندک و آهسته اي را به خود خواهد ديد. مقيـاسزماني اين نوع تغييرات بسيار بزرگ تر از مقياس زماني تغييراتي است که بر اثر آشفتگي سريع ايجـاد مـ يشـود . اگـر بخـواهيمسيستم را در حالت کاملاﹰ غيرخطي حل کنيم، با توجه بـه نکتـة گفته شده، شبيه سازي بسيار زمان گير خواهد بود. مي توان گفـتبه شرط اينکه تغييرات شديدي در شار وجـود نداشـته باشـد وتغييرات شار فقط در اثر تغييرات غلظـت عناصـر باشـد، دقـت
استراتژي کام ﹰلا غيرخطي نسبت به دقت اسـتراتژي هـاي خطـ ي خيلي متفاوت نخواهد بود. بـا توجـه بـه نـوع و ذات مسـئله و
ج

شکل۱ – شکل قلب و شرايط مرزي آن

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

با انرژي ۱ که به صورت عمود وارد قلب مي شوند) فـرضشده است کـه در آن 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 کـه در ايـن روش بـه تکرارهاي داخلي معروف است) ايـن روش در حـدود ۴ تـا ۵ باقي ماند. روش هاي مبتني بر زيرفضاي کرايلـف نيـاز بـه يـکحدس اوليه براي جواب دارند و هرچه اين حدس بهتر انتخاب شود تکرارها در جهت رسيدن به محدودة خطـاي مجـاز کمتـرخواهد شد. در اين تحقيق از پيش حلگر قطري۳۷ استفاده شـدهاست.

۳ -۳ – حل زير سيستم واکنش
زيرسيستم واکنش شامل مابقي معادلات مرتبـ ة اول مربـوط بـهتغيير غلظـت عناصـر و نـوترون هـاي تـأخيري کـه در نهايـتتغييرات مربوط به q را رقم مـي زننـد، مـي شـود . بـا توجـه بـهتکرارهاي فوق العاده زياد موجود در مسير حل مسـ ئله و امکـانگسترش خطا، تغييرات اندک غلظت عناصر و هم چنـين نکـاتاشاره شده در قسمت (۳ -۱)، يکی از بهترين راهکارها اسـتفادهاز روش رانگ – کوتا مي باشد. دقت اين روش کمـک خواهـدکرد که تغييرات غلظت عناصر را که دليل اصلي تغييـرات شـاراست بسيار بهتر ثبت کنيم. هزينة محاسباتي اين روش نسبت به روشي مانند تتا زيادتر است ولـي بـا انتخـاب تعـداد گـام هـايمناسب، مي توان به دقت به مراتب بهتر دست پيدا کرد. در ابتـدااز روش صريح اولر ( 0θ= ) براي اين کار استفاده شد که نتيجة مناسبي در بر نداشت و به واگرائي جواب ها منجر مـي شـد . در اين تحقيق به طور خاص از روش ۵ گامي رانگ – کوتا استفاده شده است. ضرائب مربوط بـه ايـن روش در کتـب مربـوط بـهمعادلات ديفرانسيل معمولي آمده است که جهـت جلـوگيری ازاطناب کلام از آوردن آنها صرف نظر مي شود.
هر يک از



قیمت: تومان


پاسخ دهید