KUCA Modeling problems

User Support
1

Hello everybody, recently I have carried out the KUCA modeling(Case Ⅱ-1).

image
image1066×820 62.5 KB

Here are my Python script


import openmc
import os
os.chdir('kuka')
# ==================================================
fuel=openmc.Material(name='93 fuel')
fuel.add_nuclide('U235',1.50694E-03)
fuel.add_nuclide('U238',1.08560E-04)
fuel.add_nuclide('Al27',5.56436E-02)
fuel.set_density('sum')#
# ==================================================
air=openmc.Material(name='air')
air.add_element('N' ,0.000039)
air.add_element('O' ,0.000011)
air.set_density('atom/b-cm', 5.000E-05)
# ==================================================
al=openmc.Material(name='Al')
al.add_nuclide('Al27',1)
al.set_density('atom/b-cm',6.00385E-02)
# ==================================================
# PE1/8
pe1_8=openmc.Material(name='PET1/8')
pe1_8.add_element('C',4.01084E-02)
pe1_8.add_element('H',8.02167E-02)
pe1_8.set_density('atom/cm3',(4.01084E-02+8.02167E-02)*1e24)
# ==================================================
# PE1/4
pe1_4=openmc.Material(name='PET1/4')
pe1_4.add_element('C',4.04356E-02)
pe1_4.add_element('H',8.08711E-02)
pe1_4.set_density('sum')
# ==================================================
# PE1/2
pe1_2=openmc.Material(name='PET1/2')
pe1_2.add_element('C',4.03280E-02)
pe1_2.add_element('H',8.06560E-02)
pe1_2.set_density('sum')
# ==================================================
# PE19
pe19=openmc.Material(name='PET19')
pe19.add_element('C',4.00042E-02)
pe19.add_element('H',8.00083E-02)
pe19.set_density('sum')
# ==================================================
# Pb-Bi
pb_bi=openmc.Material(name='pb_bi')
pb_bi.add_nuclide('Pb204',1.87461E-04)
pb_bi.add_nuclide('Pb206',3.25860E-03)
pb_bi.add_nuclide('Pb207',3.00266E-03)
pb_bi.add_nuclide('Pb208',7.15378E-03)
pb_bi.add_nuclide('Bi209',1.67670E-02)
pb_bi.set_density('sum')
# ==================================================
# Pb-Bi coating 1
pb_bi_c1=openmc.Material(name='pb_bi coating 1st')
pb_bi_c1.add_nuclide('H1',2.83301E-03)
pb_bi_c1.add_nuclide('H2',4.25015E-07)
pb_bi_c1.add_nuclide('O16',2.06885E-03+7.88039E-07+4.14757E-06)
pb_bi_c1.add_element('C',2.27058E-03)
pb_bi_c1.add_nuclide('Ti46',2.50941E-05)
pb_bi_c1.add_nuclide('Ti47',2.28983E-05)
pb_bi_c1.add_nuclide('Ti48',2.31493E-04)
pb_bi_c1.add_nuclide('Ti49',1.72522E-05)
pb_bi_c1.add_nuclide('Ti50',1.69385E-05)
pb_bi_c1.add_nuclide('S32',4.16818E-04+8.77326E-08)
pb_bi_c1.add_nuclide('S33',3.28997E-06)
pb_bi_c1.add_nuclide('S34',1.84677E-05)

pb_bi_c1.add_nuclide('Ba132',4.43050E-07)
pb_bi_c1.add_nuclide('Ba134',1.06025E-05)
pb_bi_c1.add_nuclide('Ba135',2.89167E-05)
pb_bi_c1.add_nuclide('Ba136',3.44526E-05)
pb_bi_c1.add_nuclide('Ba137',4.92619E-05)
pb_bi_c1.add_nuclide('Ba138',3.14521E-04)
pb_bi_c1.set_density('sum')
# ==================================================
# Pb-Bi coating 1
pb_bi_c2=openmc.Material(name='pb_bi coating 2nd')
pb_bi_c2.add_nuclide('H1',3.78991E-03)
pb_bi_c2.add_nuclide('H2',5.64072E-07)
pb_bi_c2.add_nuclide('O16',4.58782E-04+1.74753E-07+9.19753E-07)
pb_bi_c2.add_element('C',4.03671E-03)
pb_bi_c2.add_nuclide('Ti46',5.36896E-05)
pb_bi_c2.add_nuclide('Ti47',4.89918E-05)
pb_bi_c2.add_nuclide('Ti48',4.95287E-04)
pb_bi_c2.add_nuclide('Ti49',3.69116E-05)
pb_bi_c2.add_nuclide('Ti50',3.62405E-05)
pb_bi_c2.add_nuclide('Si28',1.86243E-05)
pb_bi_c2.add_nuclide('Si29',9.43026E-07)
pb_bi_c2.add_nuclide('Si30',6.25992E-07)
pb_bi_c2.set_density('sum')
# ==================================================
# control rod
rod=openmc.Material(name='B2O3')
rod.add_nuclide('B10',3.87448E-03)
rod.add_nuclide('B11',1.68447E-02)
rod.add_nuclide('O16',3.10787E-02)
rod.set_density('sum')
# ==================================================
# moderator PE1_8
m_pe1_8=openmc.Material(name='moderator PE1/8')
m_pe1_8.add_element('C',3.95860E-02)
m_pe1_8.add_element('H',7.77938E-02)
m_pe1_8.set_density('sum')
# ==================================================
# moderator PE10
m_pe10=openmc.Material(name='moderator PE10')
m_pe10.add_element('C',4.08960E-02)
m_pe10.add_element('H',7.97990E-02)
m_pe10.set_density('sum')
# ==================================================
# geometry
# region 
r_up1_l1_in=openmc.rectangular_prism(5.08,5.08,'z',    (-0.05,0.05))
r_up1_l1_gap=openmc.rectangular_prism(5.13,5.13,'z',  (-0.05,0.05) )
r_up1_l1_clad=openmc.rectangular_prism(5.43,5.43,'z',  (-0.05,0.05))

# region pb_bi
r_up1_l1_in_pi=openmc.rectangular_prism(5.00,5.00,'z',    (-0.05,0.05))
r_up1_l1_gap_pi=openmc.rectangular_prism(5.08,5.08,'z',  (-0.05,0.05) )
r_up1_l1_clad_pi=openmc.rectangular_prism(5.16,5.16,'z',  (-0.05,0.05))

#=======================================================

# control rod region
r=openmc.ZCylinder(r=1.69)
r1=openmc.ZCylinder(r=2)
r2=openmc.ZCylinder(r=2.25)
r3=openmc.ZCylinder(r=2.5)
r_0=openmc.ZPlane(z0=0)
r_1=openmc.ZPlane(z0=147)

rod_in=openmc.Cell(fill=rod,region=-r ) #& +r_0 & -r_1
rod_clad_in=openmc.Cell(fill=al,region=+r & -r1)
rod_gap=openmc.Cell(fill=None,region=+r1 & -r2)
rod_clad_out=openmc.Cell(fill=al,region=+r2 & -r3)
rod_out=openmc.Cell(fill=None,region=+r3)
u_rod=openmc.Universe(cells=[rod_in,rod_clad_in,rod_gap,rod_clad_out,rod_out])

#=======================================================
# Air
c_air=openmc.Cell(fill=air)
u1=openmc.Universe(name='air',cells=[c_air])

#  Al 
c_al=openmc.Cell(fill=al,region=r_up1_l1_in)
c_al_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_al_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_al_out=openmc.Cell(fill=air)
u_al=openmc.Universe(cells=[c_al,c_al_in,c_al_clad,c_al_out])

#  m_pe10
c_m_pe10=openmc.Cell(fill=m_pe10,region=r_up1_l1_in)
c_m_pe10_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_m_pe10_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_m_pe10_out=openmc.Cell(fill=air)
u_m_pe10=openmc.Universe(cells=[c_m_pe10,c_m_pe10_in,c_m_pe10_clad,c_m_pe10_out])

#  m_pe1_8
c_m_pe1_8=openmc.Cell(fill=m_pe1_8,region=r_up1_l1_in)
c_m_pe1_8_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_m_pe1_8_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_m_pe1_8_out=openmc.Cell(fill=air)
u_m_pe1_8=openmc.Universe(cells=[c_m_pe1_8,c_m_pe1_8_in,c_m_pe1_8_clad,c_m_pe1_8_out])

#  pe19
c_pe19=openmc.Cell(fill=pe19,region=r_up1_l1_in)
c_pe19_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_pe19_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_pe19_out=openmc.Cell(fill=air)
u_pe19=openmc.Universe(cells=[c_pe19,c_pe19_in,c_pe19_clad,c_pe19_out])


#  fuel
c_fuel=openmc.Cell(fill=fuel,region=r_up1_l1_in)
c_fuel_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_fuel_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_fuel_out=openmc.Cell(fill=air)
u_fuel=openmc.Universe(cells=[c_fuel,c_fuel_in,c_fuel_clad,c_fuel_out])

#  pb_bi
c_pi=openmc.Cell(fill=pb_bi,region=r_up1_l1_in_pi)
c_pi_in=openmc.Cell(fill=pb_bi_c1,region=~r_up1_l1_in_pi & r_up1_l1_gap_pi)
c_pi_clad=openmc.Cell(fill=pb_bi_c2,region=~ r_up1_l1_gap_pi & r_up1_l1_clad_pi)
c_pi_out=openmc.Cell(fill=air)
u_pi=openmc.Universe(cells=[c_pi,c_pi_in,c_pi_clad,c_pi_out])

#  F assembly
z_s_F=[0, 2, 24.5425 ,49.9425 ,89.9475 , 98.52 ,149.32, 152.40102]
z_list_F=[openmc.ZPlane(z0=i)for i in z_s_F]

lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_F=+z_list_F[0] & -z_list_F[1] & lat_box
r_low_F18=+z_list_F[1] & -z_list_F[2] & lat_box
r_low_F10=+z_list_F[2] & -z_list_F[3] & lat_box
r_fuel_F=+z_list_F[3] & -z_list_F[4] & lat_box
r_up_F18=+z_list_F[4] & -z_list_F[5] & lat_box
r_up_F10=+z_list_F[5] & -z_list_F[6] & lat_box
r_viod_F=+z_list_F[6] & -z_list_F[7] & lat_box

c_al_F=openmc.Cell(fill=al,region=r_al_F)

a=openmc.RectLattice()
a.pitch=(5.53,5.53,0.3175)
a.lower_left=(-5.53/2,-5.53/2,2)
a.universes=[[[u_m_pe1_8]]]*71

c_low_F18=openmc.Cell(fill=a,region=r_low_F18)

b=openmc.RectLattice()
b.pitch=(5.53,5.53,25.4)
b.universes=[[[u_m_pe10]]]
b.lower_left=(-5.53/2,-5.53/2,24.5425)

c_low_F10=openmc.Cell(fill=b,region=r_low_F10)

c=openmc.RectLattice()
c.pitch=(5.53,5.53,0.15875)
c.lower_left=(-5.53/2,-5.53/2,49.9425)
c.universes=[[[u_fuel]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]
            ,[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*36

c_fuel_F=openmc.Cell(fill=c,region=r_fuel_F)

d=openmc.RectLattice()
d.pitch=(5.53,5.53,0.3175)
d.lower_left=(-5.53/2,-5.53/2,89.9475)
d.universes=[[[u_m_pe1_8]]]*27

c_up_F18=openmc.Cell(fill=d,region=r_up_F18)

e=openmc.RectLattice()
e.pitch=(5.53,5.53,25.4)
e.lower_left=(-5.53/2,-5.53/2,98.52)
e.universes=[[[u_m_pe10]]]*2

c_up_F10=openmc.Cell(fill=e,region=r_up_F10)

c_void_F=openmc.Cell(fill=None,region=r_viod_F)
u_F=openmc.Universe(cells=[c_al_F,c_low_F18,c_low_F10,c_fuel_F,c_up_F18,c_up_F10,c_void_F])

# ================================================================
#  f assembly
z_s_f=[0, 2, 25.1775 ,50.5775, 60.1025 , 79.1525 , 88.6775 , 97.25 , 148.05 , 151.13102]
z_list_f=[openmc.ZPlane(z0=i)for i in z_s_f]
lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_f=+z_list_f[0] & -z_list_f[1] & lat_box
r_low_f18=+z_list_f[1] & -z_list_f[2] & lat_box
r_low_f10=+z_list_f[2] & -z_list_f[3] & lat_box
r_low_p  =+z_list_f[3] & -z_list_f[4] & lat_box
r_fuel_f=+z_list_f[4] & -z_list_f[5] & lat_box
r_up_p  =+z_list_f[5] & -z_list_f[6] & lat_box
r_up_f18=+z_list_f[6] & -z_list_f[7] & lat_box
r_up_f10=+z_list_f[7] & -z_list_f[8] & lat_box
r_viod_f=+z_list_f[8] & -z_list_f[9] & lat_box

c_al_f=openmc.Cell(fill=al,region=r_al_f)

aa=openmc.RectLattice()
aa.pitch=(5.53,5.53,0.3175)
aa.lower_left=(-5.53/2,-5.53/2,2)
aa.universes=[[[u_m_pe1_8]]]*73

c_low_f18=openmc.Cell(fill=aa,region=r_low_f18)

bb=openmc.RectLattice()
bb.pitch=(5.53,5.53,25.4)
bb.universes=[[[u_m_pe10]]]
bb.lower_left=(-5.53/2,-5.53/2,z_s_f[2])

c_low_f10=openmc.Cell(fill=bb,region=r_low_f10)

cc=openmc.RectLattice()
cc.pitch=(5.53,5.53,9.525)
cc.universes=[[[u_m_pe10]]]
cc.lower_left=(-5.53/2,-5.53/2,z_s_f[3])

c_low_p=openmc.Cell(fill=cc,region=r_low_p)


dd=openmc.RectLattice()
dd.pitch=(5.53,5.53,19.05)
dd.lower_left=(-5.53/2,-5.53/2,z_s_f[4])
dd.universes=[[[u_pi]]]

c_fuel_f=openmc.Cell(fill=dd,region=r_fuel_f)

ee=openmc.RectLattice()
ee.pitch=(5.53,5.53,9.525)
ee.lower_left=(-5.53/2,-5.53/2,z_s_f[5])
ee.universes=[[[u_m_pe10]]]

c_up_p=openmc.Cell(fill=ee,region=r_up_p)

ff=openmc.RectLattice()
ff.pitch=(5.53,5.53,0.3175)
ff.universes=[[[u_m_pe1_8]]]*27
ff.lower_left=(-5.53/2,-5.53/2,z_s_f[6])

c_up_f18=openmc.Cell(fill=ff,region=r_up_f18)

gg=openmc.RectLattice()
gg.pitch=(5.53,5.53,25.4)
gg.universes=[[[u_m_pe10]]]*2
gg.lower_left=(-5.53/2,-5.53/2,z_s_f[7])

c_up_f10=openmc.Cell(fill=gg,region=r_up_f10)

c_void_f=openmc.Cell(fill=None,region=r_viod_f)

u_f=openmc.Universe(cells=[c_al_f,c_low_f18,c_low_f10,c_low_p,c_fuel_f,c_up_p,c_up_f18,c_up_f10,c_void_f])


# ================================================================
#  16 assembly
z_s_16=[0, 2, 13.1125 , 35.655 , 61.055 , 78.835 , 87.4075 , 138.2075 , 149.32 , 152.40102]
z_list_16=[openmc.ZPlane(z0=i)for i in z_s_16]
lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_16    =+z_list_16[0] & -z_list_16[1] & lat_box
r_low_al_16=+z_list_16[1] & -z_list_16[2] & lat_box
r_low_1618 =+z_list_16[2] & -z_list_16[3] & lat_box
r_low_1610 =+z_list_16[3] & -z_list_16[4] & lat_box
r_fuel_16  =+z_list_16[4] & -z_list_16[5] & lat_box
r_up_1618  =+z_list_16[5] & -z_list_16[6] & lat_box
r_up_1610  =+z_list_16[6] & -z_list_16[7] & lat_box
r_up_al_16 =+z_list_16[7] & -z_list_16[8] & lat_box
r_viod_16  =+z_list_16[8] & -z_list_16[9] & lat_box

c_al_16=openmc.Cell(fill=al,region=r_al_16)

aaa=openmc.RectLattice()
aaa.pitch=(5.53,5.53,0.15875)
aaa.lower_left=(-5.53/2,-5.53/2,2)
aaa.universes=[[[u_al]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*10

c_low_al_16=openmc.Cell(fill=aaa,region=r_low_al_16)

bbb=openmc.RectLattice()
bbb.pitch=(5.53,5.53,0.3175)
bbb.universes=[[[u_m_pe1_8]]]*71
bbb.lower_left=(-5.53/2,-5.53/2,z_s_16[2])

c_low_1618=openmc.Cell(fill=bbb,region=r_low_1618)

ccc=openmc.RectLattice()
ccc.pitch=(5.53,5.53,25.4)
ccc.universes=[[[u_m_pe10]]]
ccc.lower_left=(-5.53/2,-5.53/2,z_s_16[3])

c_low_1610=openmc.Cell(fill=ccc,region=r_low_1610)

ddd=openmc.RectLattice()
ddd.pitch=(5.53,5.53,0.15875)
ddd.lower_left=(-5.53/2,-5.53/2,z_s_16[4])
ddd.universes=[[[u_fuel]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*16

c_fuel_16=openmc.Cell(fill=ddd,region=r_fuel_16)

eee=openmc.RectLattice()
eee.pitch=(5.53,5.53,0.3175)
eee.lower_left=(-5.53/2,-5.53/2,z_s_16[5])
eee.universes=[[[u_m_pe1_8]]]*27

c_up_1618=openmc.Cell(fill=eee,region=r_up_1618)

fff=openmc.RectLattice()
fff.pitch=(5.53,5.53,25.4)
fff.universes=[[[u_m_pe10]]]*2
fff.lower_left=(-5.53/2,-5.53/2,z_s_16[6])

c_up_1610=openmc.Cell(fill=fff,region=r_up_1610)

ggg=openmc.RectLattice()
ggg.pitch=(5.53,5.53,0.15875)
ggg.universes=[[[u_al]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*10
ggg.lower_left=(-5.53/2,-5.53/2,z_s_16[7])

c_up_al_16=openmc.Cell(fill=ggg,region=r_up_al_16)

c_void_16=openmc.Cell(fill=None,region=r_viod_16)

u_16=openmc.Universe(cells=[c_al_16,c_low_al_16,c_low_1618,c_low_1610,c_fuel_16,c_up_1618,c_up_1610,c_up_al_16,c_void_16])
# ============================================================================================================================
# A’-I left
l_ai_lat=openmc.RectLattice()
l_ai_lat.pitch=(5.53,5.53)
l_ai_lat.lower_left=(0,0)
l_ai_lat.universes=[[u_al]*4]*7

l_ai_box=openmc.rectangular_prism(5.53*4,5.53*7,'z',(5.53*4/2,5.53*7/2))
l_ai_cell=openmc.Cell(fill=l_ai_lat,region=l_ai_box)
# ============================================================================================================================
# J-P left
l_jp_lat=openmc.RectLattice()
l_jp_lat.pitch=(5.53,5.53)
l_jp_lat.lower_left=(0,39.01)
l_jp_lat.universes=[[u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_al,u_pe19,u_al],
                    [u_al,u_al,u_pe19,u_pe19],
                    [u_al]*4,
                    [u_al]*4,]

l_jp_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2,5.53*6/2+39.01))
l_jp_cell=openmc.Cell(fill=l_jp_lat,region=l_jp_box)
# ============================================================================================================================
# Q-T left
l_qt_lat=openmc.RectLattice()
l_qt_lat.pitch=(5.53,5.53)
l_qt_lat.lower_left=(0,72.49)
l_qt_lat.universes=[[u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_pe19,u_al  ,u_al  ],
                    [u_al,u_pe19,u_al  ,u_al  ]]

l_qt_box=openmc.rectangular_prism(5.53*4,5.53*3,'z',(5.53*4/2,5.53*3/2+72.49))
l_qt_cell=openmc.Cell(fill=l_qt_lat,region=l_qt_box)
# ============================================================================================================================
# U-Z left
l_uz_lat=openmc.RectLattice()
l_uz_lat.pitch=(5.53,5.53)
l_uz_lat.lower_left=(0,89.38)
l_uz_lat.universes=[[u_al]*4]*6

l_uz_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2,5.53*6/2+89.38))
l_uz_cell=openmc.Cell(fill=l_uz_lat,region=l_uz_box)
# ============================================================================================================================
# U-Z left1
l_uz1_lat=openmc.RectLattice()
l_uz1_lat.pitch=(5.53,5.53)
l_uz1_lat.lower_left=(22.42,89.38)
l_uz1_lat.universes=[[u_al]*6,
                     [u_pe19]*6,
                     [u_pe19,u_pe19,u_al,u_al,u_pe19,u_pe19],
                     [u_pe19,u_pe19,u_al,u_al,u_pe19,u_pe19],
                     [u_pe19]*6,
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10]]

l_uz1_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+22.42,5.53*6/2+89.38))
l_uz1_cell=openmc.Cell(fill=l_uz1_lat,region=l_uz1_box)
# ============================================================================================================================
# Q-T left1
l_qt1_lat=openmc.RectLattice()
l_qt1_lat.pitch=(5.53,5.53)
l_qt1_lat.lower_left=(22.42,72.49)
l_qt1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_rod   ,u_F]]

l_qt1_box=openmc.rectangular_prism(5.53*6,5.53*3,'z',(5.53*6/2+22.42,5.53*3/2+72.49))
l_qt1_cell=openmc.Cell(fill=l_qt1_lat,region=l_qt1_box)
# ============================================================================================================================
# J-P left1
l_jp1_lat=openmc.RectLattice()
l_jp1_lat.pitch=(5.53,5.53)
l_jp1_lat.lower_left=(22.42,39.01)
l_jp1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_rod   ,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_rod   ,u_F]]

l_jp1_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+22.42,5.53*6/2+39.01))
l_jp1_cell=openmc.Cell(fill=l_jp1_lat,region=l_jp1_box)
# ============================================================================================================================
# A’-I left1
l_ai1_lat=openmc.RectLattice()
l_ai1_lat.pitch=(5.53,5.53)
l_ai1_lat.lower_left=(22.42,0)
l_ai1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_al]*6]

l_ai1_box=openmc.rectangular_prism(5.53*6,5.53*7,'z',(5.53*6/2+22.42,5.53*7/2))
l_ai1_cell=openmc.Cell(fill=l_ai1_lat,region=l_ai1_box)
# ============================================================================================================================
# A’-I left2
u2=openmc.Universe(name='empty uni',cells=[openmc.Cell()])
l_ai2_lat=openmc.RectLattice()
l_ai2_lat.pitch=(5.53,5.53)
l_ai2_lat.lower_left=(55.9,0)
l_ai2_lat.universes=[[u_F,u_m_pe10,u_F],
                     [u_F,u_m_pe10,u_F],
                     [u_F,u_m_pe10,u_F],
                     [u_16,u_m_pe10,u_16],
                     [u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_al,u_m_pe10,u_al]]

l_ai2_box=openmc.rectangular_prism(5.53*3,5.53*7,'z',(5.53*3/2+55.9,5.53*7/2))
l_ai2_cell=openmc.Cell(fill=l_ai2_lat,region=l_ai2_box)
# ============================================================================================================================
# J-P left2
l_jp2_lat=openmc.RectLattice()
l_jp2_lat.pitch=(5.53,5.53)
l_jp2_lat.lower_left=(55.9,39.01)
l_jp2_lat.universes=[[u_F,u_F,u_F],
                     [u_f,u_f,u_f],
                     [u_f,u_f,u_f],
                     [u_f,u_m_pe10 ,u_f],
                     [u_F,u_m_pe10 ,u_F],
                     [u_F,u_m_pe10 ,u_F]]

l_jp2_box=openmc.rectangular_prism(5.53*3,5.53*6,'z',(5.53*3/2+55.9,5.53*6/2+39.01))
l_jp2_cell=openmc.Cell(fill=l_jp2_lat,region=l_jp2_box)
# ============================================================================================================================
# Q-T left2
l_qt2_lat=openmc.RectLattice()
l_qt2_lat.pitch=(5.53,5.53)
l_qt2_lat.lower_left=(55.9,72.49)
l_qt2_lat.universes=[[u_pe19,u_pe19,u_pe19],
                     [u_pe19,u_pe19,u_pe19],
                     [u_F   ,u_F   ,u_F   ]]

l_qt2_box=openmc.rectangular_prism(5.53*3,5.53*3,'z',(5.53*3/2+55.9,5.53*3/2+72.49))
l_qt2_cell=openmc.Cell(fill=l_qt2_lat,region=l_qt2_box)
# ===========================================================================================================================
# U-Z left2
l_uz2_lat=openmc.RectLattice()
l_uz2_lat.pitch=(5.53,5.53)
l_uz2_lat.lower_left=(55.9,89.38)
l_uz2_lat.universes=[[u_al]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_m_pe10]*3]

l_uz2_box=openmc.rectangular_prism(5.53*3,5.53*6,'z',(5.53*3/2+55.9,5.53*6/2+89.38))
l_uz2_cell=openmc.Cell(fill=l_uz2_lat,region=l_uz2_box)
# ===========================================================================================================================
# U-Z left3
l_uz3_lat=openmc.RectLattice()
l_uz3_lat.pitch=(5.53,5.53)
l_uz3_lat.lower_left=(72.79,89.38)
l_uz3_lat.universes=[[u_al]*6,
                     [u_pe19,u_pe19,u_pe19,u_pe19,u_pe19,u_al],
                     [u_pe19,u_pe19,u_al  , u_al ,u_pe19,u_al],
                     [u_pe19,u_pe19,u_al  , u_al ,u_pe19,u_al],
                     [u_pe19,u_pe19,u_pe19,u_pe19,u_pe19,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_al]]

l_uz3_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+72.79,5.53*6/2+89.38))
l_uz3_cell=openmc.Cell(fill=l_uz3_lat,region=l_uz3_box)
# ============================================================================================================================
# Q-T left3
l_qt3_lat=openmc.RectLattice()
l_qt3_lat.pitch=(5.53,5.53)
l_qt3_lat.lower_left=(72.79,72.49)
l_qt3_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F     ,u_rod   ,u_m_pe10,u_m_pe10,u_pe19,u_pe19]]

l_qt3_box=openmc.rectangular_prism(5.53*6,5.53*3,'z',(5.53*6/2+72.79,5.53*3/2+72.49))
l_qt3_cell=openmc.Cell(fill=l_qt3_lat,region=l_qt3_box)
# ============================================================================================================================
# J-P left3
l_jp3_lat=openmc.RectLattice()
l_jp3_lat.pitch=(5.53,5.53)
l_jp3_lat.lower_left=(72.79,39.01)
l_jp3_lat.universes=[[u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_rod,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_rod,u_m_pe10,u_m_pe10,u_pe19,u_pe19]]

l_jp3_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+72.79,5.53*6/2+39.01))
l_jp3_cell=openmc.Cell(fill=l_jp3_lat,region=l_jp3_box)
# ============================================================================================================================
# A’-I left3
l_ai3_lat=openmc.RectLattice()
l_ai3_lat.pitch=(5.53,5.53)
l_ai3_lat.lower_left=(72.79,0)
l_ai3_lat.universes=[[u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_al,u_al,u_al,u_al,u_pe19,u_pe19]]

l_ai3_box=openmc.rectangular_prism(5.53*6,5.53*7,'z',(5.53*6/2+72.79,5.53*7/2))
l_ai3_cell=openmc.Cell(fill=l_ai3_lat,region=l_ai3_box)
# ============================================================================================================================
# A’-I left3
l_ai4_lat=openmc.RectLattice()
l_ai4_lat.pitch=(5.53,5.53)
l_ai4_lat.lower_left=(106.27,0)
l_ai4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al]]

l_ai4_box=openmc.rectangular_prism(5.53*4,5.53*7,'z',(5.53*4/2+106.27,5.53*7/2))
l_ai4_cell=openmc.Cell(fill=l_ai4_lat,region=l_ai4_box)
# ============================================================================================================================
# J-P left4
l_jp4_lat=openmc.RectLattice()
l_jp4_lat.pitch=(5.53,5.53)
l_jp4_lat.lower_left=(106.27,39.01)
l_jp4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al]]

l_jp4_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2+106.27,5.53*6/2+39.01))
l_jp4_cell=openmc.Cell(fill=l_jp4_lat,region=l_jp4_box)
# ============================================================================================================================
# Q-T left4
l_qt4_lat=openmc.RectLattice()
l_qt4_lat.pitch=(5.53,5.53)
l_qt4_lat.lower_left=(106.27,72.49)
l_qt4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al],
                     [u_al,u_al,u_m_pe10,u_al]]

l_qt4_box=openmc.rectangular_prism(5.53*4,5.53*3,'z',(5.53*4/2+106.27,5.53*3/2+72.49))
l_qt4_cell=openmc.Cell(fill=l_qt4_lat,region=l_qt4_box)
# ===========================================================================================================================
# U-Z left4
l_uz4_lat=openmc.RectLattice()
l_uz4_lat.pitch=(5.53,5.53)
l_uz4_lat.lower_left=(106.27,89.38)
l_uz4_lat.universes=[[u_al]*4]*6

l_uz4_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2+106.27,5.53*6/2+89.38))
l_uz4_cell=openmc.Cell(fill=l_uz4_lat,region=l_uz4_box)

# ============================================================================================================================
box=openmc.rectangular_prism(150,150,'z',(150/2,150/2),'vacuum')
z1=openmc.ZPlane(z0=0)
z1.boundary_type='vacuum'
z2=openmc.ZPlane(z0=150)
z2.boundary_type='vacuum'
box_region= +z1& -z2 & box
all_cell=openmc.Cell(fill=air)
cells=[l_ai_cell,l_jp_cell,l_qt_cell,l_uz_cell,l_uz1_cell,l_qt1_cell,l_jp1_cell,l_ai1_cell,l_ai2_cell,l_jp2_cell,l_qt2_cell,l_uz2_cell
        ,l_uz3_cell,l_qt3_cell,l_jp3_cell,l_ai3_cell,l_ai4_cell,l_jp4_cell,l_qt4_cell,l_uz4_cell,all_cell]
u_all=openmc.Universe(cells=cells)
cell_all=openmc.Cell(fill=u_all,region=box_region)
geom = openmc.Geometry([cell_all])
geom.export_to_xml()

mats = list(geom.get_all_materials().values())
openmc.Materials(mats).export_to_xml()
point = openmc.stats.Point((61, 27, 51))
src = openmc.Source(space=point)
settings = openmc.Settings()
settings.source = src
settings.batches = 200
settings.inactive = 10
settings.particles = 1000
# batches = 100
# inactive = 50
# particles = 1000
# Instantiate a Settings object
# settings_file = openmc.Settings()
# settings_file.batches = batches
# settings_file.inactive = inactive
# settings_file.particles = particles
# # Create an initial uniform spatial source distribution over fissionable zones
# bounds = [55.3, 22.12, 50, 72, 77, 88]
# uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
# settings_file.source = openmc.Source(space=uniform_dist)
settings.export_to_xml()

openmc.run(geometry_debug=1)
p = openmc.Plot()
p.filename = 'pinplot'
p.width = (10, 10)
p.pixels = (2000, 2000)
p.color_by = 'cell'
# p.colors={fuel:'blue',m_pe10:'green',m_pe10:'red',air:'black',al:'yellow'}
# p.colors={fuel:'blue'}
p.basis='xy'
p.origin=(55.9,0,100)
# p.origin=(10*5.53+0.3,4*5.53+0.3,50.0425)
plots = openmc.Plots([p])
plots.export_to_xml()
openmc.plot_geometry()

And the k result only have 0.24544,leakage fraction have 0.83594,which is far from the benckmark.

image
image796×151 6.27 KB

I tried to find out if there was a geometric error, compared the Plot with the benchmark I have not find any differences.But when I use the debug mode ,the system alerts me that there have detected overlapping cells.
image
image696×109 2.84 KB

I’ve been checking for a long time and can’t find my mistake :frowning: , so if anyone would like to help me find my mistake or suggest where I might have made it I’d really appreciate it!Thanks in advance!

2

Hi @world3,

Neat benchmark model! Unfortunately, I don’t have time to look through the Python code anytime soon, but I wanted to point you to the interactive plotter (GitHub - openmc-dev/plotter: Native plotting GUI for model design and verification), which might help you to debug some of the problems you’re having with the geometry. It should make it easier to determine the regions in the model with overlaps by turning on the “overlaps” option in the color dialog box or using the Ctrl + P keyboard shortcut.

I hope this helps you out!

Best,

Patrick

1 Like
3

Thank you very much！ It’s help me alot :smiley:

2 Likes
4

By using the openmc-plotter GUI ,I have fixed my overlapping cells problem.But here is a new issue. The result doesn’t match the benchmark with the error about 1400pcm(Case Ⅱ-6).Even though my script can work fine, I found that there are sometimes getting the warnings like that:

image
image794×42 2.69 KB

and
image
image702×27 1.8 KB
.
From the coordinates of the generated particle file, the geometry is located but no strange shapes appear.
I found that @paulromano had dealt with similar problems, but I didn’t seem to solve mine by following his ideas :frowning: .Would you mind give me some advice? Thanks in advance!
Here are my latest model:

import openmc
import os
os.chdir('kuca')
# os.environ["OPENMC_CROSS_SECTIONS"]='/home/swagger/software/openmc/xs/endfb71_hdf5/cross_sections.xml'
# ==================================================
fuel=openmc.Material(name='93 fuel')
fuel.add_nuclide('U235',1.50682E-03)
fuel.add_nuclide('U238',9.25879E-05)
fuel.add_nuclide('U236',4.82971E-06)
fuel.add_nuclide('U234',1.13659E-05)
fuel.add_nuclide('Al27',5.56436E-02)
fuel.set_density('sum')
# ==================================================
air=openmc.Material(name='air')
air.add_element('N' ,0.000039) 
air.add_element('O' ,0.000011)
air.set_density('atom/b-cm', 5.000E-05)
# ==================================================
al=openmc.Material(name='Al')
al.add_nuclide('Al27',1)
al.set_density('atom/b-cm',6.00385E-02)
# ==================================================
# PE1/8
pe1_8=openmc.Material(name='PET1/8')
pe1_8.add_element('C',4.01084E-02)
pe1_8.add_element('H',8.02167E-02)
pe1_8.set_density('sum')
# ==================================================
# PE1/4
pe1_4=openmc.Material(name='PET1/4')
pe1_4.add_element('C',4.04356E-02)
pe1_4.add_element('H',8.08711E-02)
pe1_4.set_density('sum')
# ==================================================
# PE1/2
pe1_2=openmc.Material(name='PET1/2')
pe1_2.add_element('C',4.03280E-02)
pe1_2.add_element('H',8.06560E-02)
pe1_2.set_density('sum')
# ==================================================
# PE19
pe19=openmc.Material(name='PET19')
pe19.add_element('C',4.00042E-02)
pe19.add_element('H',8.00083E-02)
pe19.set_density('sum')
# ==================================================
# Pb-Bi
pb_bi=openmc.Material(name='pb_bi')
pb_bi.add_nuclide('Pb204',1.87461E-04)
pb_bi.add_nuclide('Pb206',3.25860E-03)
pb_bi.add_nuclide('Pb207',3.00266E-03)
pb_bi.add_nuclide('Pb208',7.15378E-03)
pb_bi.add_nuclide('Bi209',1.67670E-02)
pb_bi.set_density('sum')
# ==================================================
# Pb-Bi coating 1
pb_bi_c1=openmc.Material(name='pb_bi coating 1st')
pb_bi_c1.add_nuclide('H1',2.83301E-03)
pb_bi_c1.add_nuclide('H2',4.25015E-07)
pb_bi_c1.add_nuclide('O16',2.06885E-03+7.88039E-07+4.14757E-06)
pb_bi_c1.add_element('C',2.27058E-03)
pb_bi_c1.add_nuclide('Ti46',2.50941E-05)
pb_bi_c1.add_nuclide('Ti47',2.28983E-05)
pb_bi_c1.add_nuclide('Ti48',2.31493E-04)
pb_bi_c1.add_nuclide('Ti49',1.72522E-05)
pb_bi_c1.add_nuclide('Ti50',1.69385E-05)
pb_bi_c1.add_nuclide('S32',4.16818E-04+8.77326E-08)
pb_bi_c1.add_nuclide('S33',3.28997E-06)
pb_bi_c1.add_nuclide('S34',1.84677E-05)

pb_bi_c1.add_nuclide('Ba132',4.43050E-07)
pb_bi_c1.add_nuclide('Ba134',1.06025E-05)
pb_bi_c1.add_nuclide('Ba135',2.89167E-05)
pb_bi_c1.add_nuclide('Ba136',3.44526E-05)
pb_bi_c1.add_nuclide('Ba137',4.92619E-05)
pb_bi_c1.add_nuclide('Ba138',3.14521E-04)
pb_bi_c1.set_density('sum')
# ==================================================
# Pb-Bi coating 1
pb_bi_c2=openmc.Material(name='pb_bi coating 2nd')
pb_bi_c2.add_nuclide('H1',3.78991E-03)
pb_bi_c2.add_nuclide('H2',5.64072E-07)
pb_bi_c2.add_nuclide('O16',4.58782E-04+1.74753E-07+9.19753E-07)
pb_bi_c2.add_element('C',4.03671E-03)
pb_bi_c2.add_nuclide('Ti46',5.36896E-05)
pb_bi_c2.add_nuclide('Ti47',4.89918E-05)
pb_bi_c2.add_nuclide('Ti48',4.95287E-04)
pb_bi_c2.add_nuclide('Ti49',3.69116E-05)
pb_bi_c2.add_nuclide('Ti50',3.62405E-05)
pb_bi_c2.add_nuclide('Si28',1.86243E-05)
pb_bi_c2.add_nuclide('Si29',9.43026E-07)
pb_bi_c2.add_nuclide('Si30',6.25992E-07)
pb_bi_c2.set_density('sum')
# ==================================================
# control rod
rod=openmc.Material(name='B2O3')
rod.add_nuclide('B10',3.87448E-03)
rod.add_nuclide('B11',1.68447E-02)
rod.add_nuclide('O16',3.10787E-02)
rod.set_density('sum')
# ==================================================
# moderator PE1_8
m_pe1_8=openmc.Material(name='moderator PE1/8')
m_pe1_8.add_element('C',3.95860E-02)
m_pe1_8.add_element('H',7.77938E-02)
m_pe1_8.set_density('sum')
# ==================================================
# moderator PE10
m_pe10=openmc.Material(name='moderator PE10')
m_pe10.add_element('C',4.08960E-02)
m_pe10.add_element('H',7.97990E-02)
m_pe10.set_density('sum')
# ==================================================
# geometry
# region 
r_up1_l1_in=openmc.rectangular_prism(5.08,5.08,'z',    (-0.05,0.05))
r_up1_l1_gap=openmc.rectangular_prism(5.13,5.13,'z',  (-0.05,0.05) )
r_up1_l1_clad=openmc.rectangular_prism(5.43,5.43,'z',  (-0.05,0.05))

# region pb_bi
r_up1_l1_in_pi=openmc.rectangular_prism(5.00,5.00,'z',    (-0.05,0.05))
r_up1_l1_gap_pi=openmc.rectangular_prism(5.08,5.08,'z',  (-0.05,0.05) )
r_up1_l1_clad_pi=openmc.rectangular_prism(5.16,5.16,'z',  (-0.05,0.05))

#=======================================================

# control rod region
r=openmc.ZCylinder(r=1.69)
r1=openmc.ZCylinder(r=2)
r2=openmc.ZCylinder(r=2.25)
r3=openmc.ZCylinder(r=2.5)
r_0=openmc.ZPlane(z0=0)
r_1=openmc.ZPlane(z0=147)

rod_in=openmc.Cell(fill=rod,region=-r ) #& +r_0 & -r_1
rod_clad_in=openmc.Cell(fill=al,region=+r & -r1)
rod_gap=openmc.Cell(fill=None,region=+r1 & -r2)
rod_clad_out=openmc.Cell(fill=al,region=+r2 & -r3)
rod_out=openmc.Cell(fill=None,region=+r3)
u_rod=openmc.Universe(cells=[rod_in,rod_clad_in,rod_gap,rod_clad_out,rod_out])

#=======================================================
# Air
c_void=openmc.Cell(fill=air)
u2=openmc.Universe(name='void',cells=[c_void])

#  Al 
c_al=openmc.Cell(fill=al,region=r_up1_l1_in)
c_al_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_al_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_al_out=openmc.Cell(fill=air,region=~r_up1_l1_clad)
u_al=openmc.Universe(cells=[c_al,c_al_in,c_al_clad,c_al_out])

#  m_pe10
c_m_pe10=openmc.Cell(fill=m_pe10,region=r_up1_l1_in)
c_m_pe10_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_m_pe10_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_m_pe10_out=openmc.Cell(fill=air,region=~ r_up1_l1_clad)
u_m_pe10=openmc.Universe(cells=[c_m_pe10,c_m_pe10_in,c_m_pe10_clad,c_m_pe10_out])

#  m_pe1_8
c_m_pe1_8=openmc.Cell(fill=m_pe1_8,region=r_up1_l1_in)
c_m_pe1_8_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_m_pe1_8_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_m_pe1_8_out=openmc.Cell(fill=air,region=~r_up1_l1_clad)
u_m_pe1_8=openmc.Universe(cells=[c_m_pe1_8,c_m_pe1_8_in,c_m_pe1_8_clad,c_m_pe1_8_out])

#  pe19
c_pe19=openmc.Cell(fill=pe19,region=r_up1_l1_in)
c_pe19_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_pe19_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_pe19_out=openmc.Cell(fill=air,region=~r_up1_l1_clad)
u_pe19=openmc.Universe(cells=[c_pe19,c_pe19_in,c_pe19_clad,c_pe19_out])


#  fuel
c_fuel=openmc.Cell(fill=fuel,region=r_up1_l1_in)
c_fuel_in=openmc.Cell(fill=air,region=~r_up1_l1_in & r_up1_l1_gap)
c_fuel_clad=openmc.Cell(fill=al,region=~ r_up1_l1_gap & r_up1_l1_clad)
c_fuel_out=openmc.Cell(fill=air,region=~ r_up1_l1_clad)
u_fuel=openmc.Universe(cells=[c_fuel,c_fuel_in,c_fuel_clad,c_fuel_out])

#  pb_bi
c_pi=openmc.Cell(fill=pb_bi,region=r_up1_l1_in_pi)
c_pi_in=openmc.Cell(fill=pb_bi_c1,region=~r_up1_l1_in_pi & r_up1_l1_gap_pi)
c_pi_clad=openmc.Cell(fill=pb_bi_c2,region=~ r_up1_l1_gap_pi & r_up1_l1_clad_pi)
c_pi_out=openmc.Cell(fill=air,region=~r_up1_l1_clad_pi)
u_pi=openmc.Universe(cells=[c_pi,c_pi_in,c_pi_clad,c_pi_out])

#  F assembly
z_s_F=[0, 2, 24.5425 ,49.9425 ,89.9475 , 98.52 ,149.32, 152.40102]
z_list_F=[openmc.ZPlane(z0=i)for i in z_s_F]
# z_list_F[0].boundary_type='reflective'
# z_list_F[7].boundary_type='reflective'
lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_F=+z_list_F[0] & -z_list_F[1] & lat_box
r_low_F18=+z_list_F[1] & -z_list_F[2] & lat_box
r_low_F10=+z_list_F[2] & -z_list_F[3] & lat_box
r_fuel_F=+z_list_F[3] & -z_list_F[4] & lat_box
r_up_F18=+z_list_F[4] & -z_list_F[5] & lat_box
r_up_F10=+z_list_F[5] & -z_list_F[6] & lat_box
r_viod_F=+z_list_F[6]  & lat_box

c_al_F=openmc.Cell(fill=al,region=r_al_F)

a=openmc.RectLattice()
a.pitch=(5.53,5.53,0.3175)
a.lower_left=(-5.53/2,-5.53/2,z_s_F[1])
a.universes=[[[u_m_pe1_8]]]*71
a.outer=u2
c_low_F18=openmc.Cell(fill=a,region=r_low_F18)

b=openmc.RectLattice()
b.pitch=(5.53,5.53,25.4)
b.universes=[[[u_m_pe10]]]
b.lower_left=(-5.53/2,-5.53/2,z_s_F[2])
b.outer=u2
c_low_F10=openmc.Cell(fill=b,region=r_low_F10)

c=openmc.RectLattice()
c.pitch=(5.53,5.53,0.15875)
c.lower_left=(-5.53/2,-5.53/2,z_s_F[3])
c.universes=[[[u_fuel]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]
            ,[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*36
c.outer=u2
c_fuel_F=openmc.Cell(fill=c,region=r_fuel_F)

d=openmc.RectLattice()
d.pitch=(5.53,5.53,0.3175)
d.lower_left=(-5.53/2,-5.53/2,z_s_F[4])
d.universes=[[[u_m_pe1_8]]]*27
d.outer=u2
c_up_F18=openmc.Cell(fill=d,region=r_up_F18)

e=openmc.RectLattice()
e.pitch=(5.53,5.53,25.4)
e.lower_left=(-5.53/2,-5.53/2,z_s_F[5])
e.universes=[[[u_m_pe10]]]*2
e.outer=u2
c_up_F10=openmc.Cell(fill=e,region=r_up_F10)

c_void_F=openmc.Cell(fill=None,region=r_viod_F)
u_F=openmc.Universe(cells=[c_al_F,c_low_F18,c_low_F10,c_fuel_F,c_up_F18,c_up_F10,c_void_F])

# ================================================================
#  f assembly
z_s_f=[0, 2, 25.1775 ,50.5775, 60.1025 , 79.1525 , 88.6775 , 97.25 , 148.05 , 151.13102]
z_list_f=[openmc.ZPlane(z0=i)for i in z_s_f]
lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_f=+z_list_f[0] & -z_list_f[1] & lat_box
r_low_f18=+z_list_f[1] & -z_list_f[2] & lat_box
r_low_f10=+z_list_f[2] & -z_list_f[3] & lat_box
r_low_p  =+z_list_f[3] & -z_list_f[4] & lat_box
r_fuel_f=+z_list_f[4] & -z_list_f[5] & lat_box
r_up_p  =+z_list_f[5] & -z_list_f[6] & lat_box
r_up_f18=+z_list_f[6] & -z_list_f[7] & lat_box
r_up_f10=+z_list_f[7] & -z_list_f[8] & lat_box
r_viod_f=+z_list_f[8]  & lat_box

c_al_f=openmc.Cell(fill=al,region=r_al_f)

aa=openmc.RectLattice()
aa.pitch=(5.53,5.53,0.3175)
aa.lower_left=(-5.53/2,-5.53/2,2)
aa.universes=[[[u_m_pe1_8]]]*73
aa.outer=u2
c_low_f18=openmc.Cell(fill=aa,region=r_low_f18)

bb=openmc.RectLattice()
bb.pitch=(5.53,5.53,25.4)
bb.universes=[[[u_m_pe10]]]
bb.lower_left=(-5.53/2,-5.53/2,z_s_f[2])
bb.outer=u2
c_low_f10=openmc.Cell(fill=bb,region=r_low_f10)

cc=openmc.RectLattice()
cc.pitch=(5.53,5.53,0.15875*2)
cc.universes=[[[u_pe19]],[[u_fuel]]]*15
cc.lower_left=(-5.53/2,-5.53/2,z_s_f[3])
cc.outer=u2
c_low_p=openmc.Cell(fill=cc,region=r_low_p)


dd=openmc.RectLattice()
dd.pitch=(5.53,5.53,0.3175)
dd.lower_left=(-5.53/2,-5.53/2,z_s_f[4])
dd.universes=[[[u_pi]],[[u_fuel]]]*30
dd.outer=u2
c_fuel_f=openmc.Cell(fill=dd,region=r_fuel_f)

ee=openmc.RectLattice()
ee.pitch=(5.53,5.53,0.3175)
ee.lower_left=(-5.53/2,-5.53/2,z_s_f[5])
ee.universes=[[[u_pe19]],[[u_fuel]]]*15
ee.outer=u2
c_up_p=openmc.Cell(fill=ee,region=r_up_p)

ff=openmc.RectLattice()
ff.pitch=(5.53,5.53,0.3175)
ff.universes=[[[u_m_pe1_8]]]*27
ff.lower_left=(-5.53/2,-5.53/2,z_s_f[6])
ff.outer=u2
c_up_f18=openmc.Cell(fill=ff,region=r_up_f18)

gg=openmc.RectLattice()
gg.pitch=(5.53,5.53,25.4)
gg.universes=[[[u_m_pe10]]]*2
gg.lower_left=(-5.53/2,-5.53/2,z_s_f[7])
gg.outer=u2
c_up_f10=openmc.Cell(fill=gg,region=r_up_f10)

c_void_f=openmc.Cell(fill=None,region=r_viod_f)

u_f=openmc.Universe(cells=[c_al_f,c_low_f18,c_low_f10,c_low_p,c_fuel_f,c_up_p,c_up_f18,c_up_f10,c_void_f])


# ================================================================
#  16 assembly
z_s_16=[0, 2, 13.1125 , 35.655 , 61.055 , 78.835 , 87.4075 , 138.2075 , 149.32 , 152.40102]
z_list_16=[openmc.ZPlane(z0=i)for i in z_s_16]
lat_box=openmc.rectangular_prism(5.53,5.53,origin=[0]*2)
r_al_16    =+z_list_16[0] & -z_list_16[1] & lat_box
r_low_al_16=+z_list_16[1] & -z_list_16[2] & lat_box
r_low_1618 =+z_list_16[2] & -z_list_16[3] & lat_box
r_low_1610 =+z_list_16[3] & -z_list_16[4] & lat_box
r_fuel_16  =+z_list_16[4] & -z_list_16[5] & lat_box
r_up_1618  =+z_list_16[5] & -z_list_16[6] & lat_box
r_up_1610  =+z_list_16[6] & -z_list_16[7] & lat_box
r_up_al_16 =+z_list_16[7] & -z_list_16[8] & lat_box
r_viod_16  =+z_list_16[8] &   lat_box

c_al_16=openmc.Cell(fill=al,region=r_al_16)

aaa=openmc.RectLattice()
aaa.pitch=(5.53,5.53,0.15875)
aaa.lower_left=(-5.53/2,-5.53/2,2)
aaa.universes=[[[u_al]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*10
aaa.outer=u2
c_low_al_16=openmc.Cell(fill=aaa,region=r_low_al_16)

bbb=openmc.RectLattice()
bbb.pitch=(5.53,5.53,0.3175)
bbb.universes=[[[u_m_pe1_8]]]*71
bbb.lower_left=(-5.53/2,-5.53/2,z_s_16[2])
bbb.outer=u2
c_low_1618=openmc.Cell(fill=bbb,region=r_low_1618)

ccc=openmc.RectLattice()
ccc.pitch=(5.53,5.53,25.4)
ccc.universes=[[[u_m_pe10]]]
ccc.lower_left=(-5.53/2,-5.53/2,z_s_16[3])
ccc.outer=u2
c_low_1610=openmc.Cell(fill=ccc,region=r_low_1610)

ddd=openmc.RectLattice()
ddd.pitch=(5.53,5.53,0.15875)
ddd.lower_left=(-5.53/2,-5.53/2,z_s_16[4])
ddd.universes=[[[u_fuel]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*16
ddd.outer=u2
c_fuel_16=openmc.Cell(fill=ddd,region=r_fuel_16)

eee=openmc.RectLattice()
eee.pitch=(5.53,5.53,0.3175)
eee.lower_left=(-5.53/2,-5.53/2,z_s_16[5])
eee.universes=[[[u_m_pe1_8]]]*27
eee.outer=u2
c_up_1618=openmc.Cell(fill=eee,region=r_up_1618)

fff=openmc.RectLattice()
fff.pitch=(5.53,5.53,25.4)
fff.universes=[[[u_m_pe10]]]*2
fff.lower_left=(-5.53/2,-5.53/2,z_s_16[6])
fff.outer=u2
c_up_1610=openmc.Cell(fill=fff,region=r_up_1610)

ggg=openmc.RectLattice()
ggg.pitch=(5.53,5.53,0.15875)
ggg.universes=[[[u_al]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]],[[u_m_pe1_8]]]*10
ggg.lower_left=(-5.53/2,-5.53/2,z_s_16[7])
ggg.outer=u2
c_up_al_16=openmc.Cell(fill=ggg,region=r_up_al_16)

c_void_16=openmc.Cell(fill=None,region=r_viod_16)

u_16=openmc.Universe(cells=[c_al_16,c_low_al_16,c_low_1618,c_low_1610,c_fuel_16,c_up_1618,c_up_1610,c_up_al_16,c_void_16])
# ============================================================================================================================
# u2=openmc.Universe(name='empty uni',cells=[openmc.Cell()])
# A’-I left
l_ai_lat=openmc.RectLattice()
l_ai_lat.pitch=(5.53,5.53)
l_ai_lat.lower_left=(0,0)
l_ai_lat.universes=[[u_al]*4]*7

l_ai_box=openmc.rectangular_prism(5.53*4,5.53*7,'z',(5.53*4/2,5.53*7/2))
l_ai_cell=openmc.Cell(fill=l_ai_lat,region=l_ai_box)
# ============================================================================================================================
# J-P left
l_jp_lat=openmc.RectLattice()
l_jp_lat.pitch=(5.53,5.53)
l_jp_lat.lower_left=(0,39.01-0.1)
l_jp_lat.universes=[[u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_al,u_pe19,u_al],
                    [u_al,u_al,u_pe19,u_pe19],
                    [u_al]*4,
                    [u_al]*4,]

l_jp_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2,5.53*6/2+39.01-0.1))
l_jp_cell=openmc.Cell(fill=l_jp_lat,region=l_jp_box)
# ============================================================================================================================
# Q-T left
l_qt_lat=openmc.RectLattice()
l_qt_lat.pitch=(5.53,5.53)
l_qt_lat.lower_left=(0,72.49-0.1)
l_qt_lat.universes=[[u_al,u_pe19,u_pe19,u_pe19],
                    [u_al,u_pe19,u2  ,u2  ],
                    [u_al,u_pe19,u2  ,u2  ]]

l_qt_box=openmc.rectangular_prism(5.53*4,5.53*3,'z',(5.53*4/2,5.53*3/2+72.49-0.1))
l_qt_cell=openmc.Cell(fill=l_qt_lat,region=l_qt_box)
# ============================================================================================================================
# U-Z left
l_uz_lat=openmc.RectLattice()
l_uz_lat.pitch=(5.53,5.53)
l_uz_lat.lower_left=(0,89.38-0.1)
l_uz_lat.universes=[[u_al]*4]*6

l_uz_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2,5.53*6/2+89.38-0.1))
l_uz_cell=openmc.Cell(fill=l_uz_lat,region=l_uz_box)
# ============================================================================================================================
# U-Z left1
l_uz1_lat=openmc.RectLattice()
l_uz1_lat.pitch=(5.53,5.53)
l_uz1_lat.lower_left=(22.42-0.1,89.38-0.1)
l_uz1_lat.universes=[[u_al]*6,
                     [u_pe19]*6,
                     [u_pe19,u_pe19,u2,u2,u_pe19,u_pe19],
                     [u_pe19,u_pe19,u2,u2,u_pe19,u_pe19],
                     [u_pe19]*6,
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10]]

l_uz1_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+22.42-0.1,5.53*6/2+89.38-0.1))
l_uz1_cell=openmc.Cell(fill=l_uz1_lat,region=l_uz1_box)
# ============================================================================================================================
# Q-T left1
l_qt1_lat=openmc.RectLattice()
l_qt1_lat.pitch=(5.53,5.53)
l_qt1_lat.lower_left=(22.42-0.1,72.49-0.1)
l_qt1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u2   ,u_m_pe10]]

l_qt1_box=openmc.rectangular_prism(5.53*6,5.53*3,'z',(5.53*6/2+22.42-0.1,5.53*3/2+72.49-0.1))
l_qt1_cell=openmc.Cell(fill=l_qt1_lat,region=l_qt1_box)
# ============================================================================================================================
# J-P left1
l_jp1_lat=openmc.RectLattice()
l_jp1_lat.pitch=(5.53,5.53)
l_jp1_lat.lower_left=(22.42-0.1,39.01-0.1)
l_jp1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u2   ,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_F],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u2   ,u_m_pe10]]

l_jp1_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+22.42-0.1,5.53*6/2+39.01-0.1))
l_jp1_cell=openmc.Cell(fill=l_jp1_lat,region=l_jp1_box)
# ============================================================================================================================
# A’-I left1
l_ai1_lat=openmc.RectLattice()
l_ai1_lat.pitch=(5.53,5.53)
l_ai1_lat.lower_left=(22.42-0.1,0)
l_ai1_lat.universes=[[u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_pe19,u_pe19,u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_al]*6]

l_ai1_box=openmc.rectangular_prism(5.53*6,5.53*7,'z',(5.53*6/2+22.42-0.1,5.53*7/2))
l_ai1_cell=openmc.Cell(fill=l_ai1_lat,region=l_ai1_box)
# ============================================================================================================================
# A’-I left2

l_ai2_lat=openmc.RectLattice()
l_ai2_lat.pitch=(5.53,5.53)
l_ai2_lat.lower_left=(55.9-0.1,0)
l_ai2_lat.universes=[[u_m_pe10,u2,u_m_pe10],
                     [u_m_pe10,u2,u_m_pe10],
                     [u_m_pe10,u2,u_m_pe10],
                     [u_m_pe10,u2,u_m_pe10],
                     [u_m_pe10,u2,u_m_pe10],
                     [u_m_pe10,u2,u_m_pe10],
                     [u_al,u2,u_al]]

l_ai2_box=openmc.rectangular_prism(5.53*3,5.53*7,'z',(5.53*3/2+55.9-0.1,5.53*7/2))
l_ai2_cell=openmc.Cell(fill=l_ai2_lat,region=l_ai2_box)
# ============================================================================================================================
# J-P left2
l_jp2_lat=openmc.RectLattice()
l_jp2_lat.pitch=(5.53,5.53)
l_jp2_lat.lower_left=(55.9-0.1,39.01-0.1)
l_jp2_lat.universes=[[u_F,u_F,u_F],
                     [u_f,u_f,u_f],
                     [u_f,u_f,u_f],
                     [u_f,u2 ,u_f],
                     [u_F,u2 ,u_F],
                     [u_m_pe10,u2 ,u_m_pe10]]

l_jp2_box=openmc.rectangular_prism(5.53*3,5.53*6,'z',(5.53*3/2+55.9-0.1,5.53*6/2+39.01-0.1))
l_jp2_cell=openmc.Cell(fill=l_jp2_lat,region=l_jp2_box)
# ============================================================================================================================
# Q-T left2
l_qt2_lat=openmc.RectLattice()
l_qt2_lat.pitch=(5.53,5.53)
l_qt2_lat.lower_left=(55.9-0.1,72.49-0.1)
l_qt2_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_m_pe10,u_m_pe10,u_m_pe10],
                     [u_F   ,u_F   ,u_F   ]]

l_qt2_box=openmc.rectangular_prism(5.53*3,5.53*3,'z',(5.53*3/2+55.9-0.1,5.53*3/2+72.49-0.1))
l_qt2_cell=openmc.Cell(fill=l_qt2_lat,region=l_qt2_box)
# ===========================================================================================================================
# U-Z left2
l_uz2_lat=openmc.RectLattice()
l_uz2_lat.pitch=(5.53,5.53)
l_uz2_lat.lower_left=(55.9-0.1,89.38-0.1)
l_uz2_lat.universes=[[u_al]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_pe19]*3,
                     [u_m_pe10]*3]

l_uz2_box=openmc.rectangular_prism(5.53*3,5.53*6,'z',(5.53*3/2+55.9-0.1,5.53*6/2+89.38-0.1))
l_uz2_cell=openmc.Cell(fill=l_uz2_lat,region=l_uz2_box)
# ===========================================================================================================================
# U-Z left3
l_uz3_lat=openmc.RectLattice()
l_uz3_lat.pitch=(5.53,5.53)
l_uz3_lat.lower_left=(72.79-0.1,89.38-0.1)
l_uz3_lat.universes=[[u_al]*6,
                     [u_pe19,u_pe19,u_pe19,u_pe19,u_pe19,u_al],
                     [u_pe19,u_pe19,u2  , u2 ,u_pe19,u_al],
                     [u_pe19,u_pe19,u2  , u2 ,u_pe19,u_al],
                     [u_pe19,u_pe19,u_pe19,u_pe19,u_pe19,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_al]]

l_uz3_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+72.79-0.1,5.53*6/2+89.38-0.1))
l_uz3_cell=openmc.Cell(fill=l_uz3_lat,region=l_uz3_box)
# ============================================================================================================================
# Q-T left3
l_qt3_lat=openmc.RectLattice()
l_qt3_lat.pitch=(5.53,5.53)
l_qt3_lat.lower_left=(72.79-0.1,72.49-0.1)
l_qt3_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u2   ,u_m_pe10,u_m_pe10,u_pe19,u_pe19]]

l_qt3_box=openmc.rectangular_prism(5.53*6,5.53*3,'z',(5.53*6/2+72.79-0.1,5.53*3/2+72.49-0.1))
l_qt3_cell=openmc.Cell(fill=l_qt3_lat,region=l_qt3_box)
# ============================================================================================================================
# J-P left3
l_jp3_lat=openmc.RectLattice()
l_jp3_lat.pitch=(5.53,5.53)
l_jp3_lat.lower_left=(72.79-0.1,39.01-0.1)
l_jp3_lat.universes=[[u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u2,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_F,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u2,u_m_pe10,u_m_pe10,u_pe19,u_pe19]]

l_jp3_box=openmc.rectangular_prism(5.53*6,5.53*6,'z',(5.53*6/2+72.79-0.1,5.53*6/2+39.01-0.1))
l_jp3_cell=openmc.Cell(fill=l_jp3_lat,region=l_jp3_box)
# ============================================================================================================================
# A’-I left3
l_ai3_lat=openmc.RectLattice()
l_ai3_lat.pitch=(5.53,5.53)
l_ai3_lat.lower_left=(72.79-0.1,0)
l_ai3_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_m_pe10,u_pe19,u_pe19],
                     [u_al,u_al,u_al,u_al,u_pe19,u_pe19]]

l_ai3_box=openmc.rectangular_prism(5.53*6,5.53*7,'z',(5.53*6/2+72.79-0.1,5.53*7/2))
l_ai3_cell=openmc.Cell(fill=l_ai3_lat,region=l_ai3_box)
# ============================================================================================================================
# A’-I left4
l_ai4_lat=openmc.RectLattice()
l_ai4_lat.pitch=(5.53,5.53)
l_ai4_lat.lower_left=(106.27-0.1,0)
l_ai4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u2,u2,u_m_pe10,u_al],
                     [u2,u2,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al]]

l_ai4_box=openmc.rectangular_prism(5.53*4,5.53*7,'z',(5.53*4/2+106.27-0.1,5.53*7/2))
l_ai4_cell=openmc.Cell(fill=l_ai4_lat,region=l_ai4_box)
# ============================================================================================================================
# J-P left4
l_jp4_lat=openmc.RectLattice()
l_jp4_lat.pitch=(5.53,5.53)
l_jp4_lat.lower_left=(106.27-0.1,39.01-0.1)
l_jp4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u2,u2,u_m_pe10,u_al],
                     [u2,u2,u_m_pe10,u_al],
                     [u_m_pe10,u_m_pe10,u_m_pe10,u_al]]

l_jp4_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2+106.27-0.1,5.53*6/2+39.01-0.1))
l_jp4_cell=openmc.Cell(fill=l_jp4_lat,region=l_jp4_box)
# ============================================================================================================================
# Q-T left4
l_qt4_lat=openmc.RectLattice()
l_qt4_lat.pitch=(5.53,5.53)
l_qt4_lat.lower_left=(106.27-0.1,72.49-0.1)
l_qt4_lat.universes=[[u_m_pe10,u_m_pe10,u_m_pe10,u_al],
                     [u2,u2,u_m_pe10,u_al],
                     [u2,u2,u_m_pe10,u_al]]

l_qt4_box=openmc.rectangular_prism(5.53*4,5.53*3,'z',(5.53*4/2+106.27-0.1,5.53*3/2+72.49-0.1))
l_qt4_cell=openmc.Cell(fill=l_qt4_lat,region=l_qt4_box)
# ===========================================================================================================================
# U-Z left4
l_uz4_lat=openmc.RectLattice()
l_uz4_lat.pitch=(5.53,5.53)
l_uz4_lat.lower_left=(106.27-0.1,89.38-0.1)
l_uz4_lat.universes=[[u_al]*4]*6

l_uz4_box=openmc.rectangular_prism(5.53*4,5.53*6,'z',(5.53*4/2+106.27-0.1,5.53*6/2+89.38-0.1))
l_uz4_cell=openmc.Cell(fill=l_uz4_lat,region=l_uz4_box)

# ============================================================================================================================
box=openmc.rectangular_prism(160,160,'z',(50,50),'vacuum')
z1=openmc.ZPlane(z0=0)
z1.boundary_type='vacuum'
z2=openmc.ZPlane(z0=152.40102)
z2.boundary_type='vacuum'

ass_region= l_ai_box | l_ai1_box | l_ai2_box | l_ai3_box | l_ai4_box | l_jp_box | l_jp1_box | l_jp2_box | l_jp3_box | l_jp4_box | l_qt_box | l_qt1_box | l_qt2_box | l_qt3_box | l_qt4_box | l_uz_box  | l_uz1_box | l_uz2_box | l_uz3_box | l_uz4_box

box_region=    +z1& -z2 & box
all_cell=openmc.Cell(fill=air,region= ~ass_region & box_region)
cells=[l_ai_cell,l_jp_cell,l_qt_cell,l_uz_cell,l_uz1_cell,l_qt1_cell,l_jp1_cell,l_ai1_cell,l_ai2_cell,l_jp2_cell,l_qt2_cell,l_uz2_cell
        ,l_uz3_cell,l_qt3_cell,l_jp3_cell,l_ai3_cell,l_ai4_cell,l_jp4_cell,l_qt4_cell,l_uz4_cell,all_cell]
u_all=openmc.Universe(cells=cells)
cell_all=openmc.Cell(fill=u_all,region=box_region)
geom = openmc.Geometry([cell_all])
geom.export_to_xml()

mats = list(geom.get_all_materials().values())
openmc.Materials(mats).export_to_xml()

batches = 300
inactive = 50
particles = 10000
settings_file = openmc.Settings()
settings_file.batches = batches
settings_file.inactive = inactive
settings_file.particles = particles
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [0, 0, 0, 160, 160, 150]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings_file.source = openmc.Source(space=uniform_dist)
settings_file.export_to_xml()

openmc.run(geometry_debug=1)
p = openmc.Plot()
p.filename = 'pinplot'
p.width = (5, 5)
p.pixels = (5000, 5000)
p.color_by = 'material'
# p.colors={fuel:'blue',m_pe10:'green',m_pe10:'red',air:'black',al:'yellow'}
# p.colors={fuel:'blue'}
p.basis='xz'
p.origin=(56.373749634942726,57.04400653836615,79.76028052979754)
# p.origin=(10*5.53+0.3,4*5.53+0.3,50.0425)
plots = openmc.Plots([p])
plots.export_to_xml()
openmc.plot_geometry()
5

Regarding the lattice errors, one thing you may want to try is to ensure that every lattice in the problem has a defined outer universe (i.e., lattice.outer = ...). Other than that, I don’t see any obvious problems in your setup. Are you seeing good agreement with the benchmark keff value on the other cases for this problem?

6

No,I’m not.In the Case Ⅰ-1, these is also about 1400pcm error. :frowning:

7

Do you know if other codes do get good agreement on this benchmark? For some benchmarks, it’s not surprising to see large differences for a variety of reasons (large experimental uncertainties, unknown material compositions, etc.). However, if you know that other codes get good agreement on this particular benchmark and you’re not seeing good agreement with OpenMC, then something is clearly wrong.

8

Hello, Pual. The banchmark book have MCNP code with case Ⅰ-1,which do get good agreement with the experience result. Comparing with the OpenMC’s input, the MCNP code is so complex that I can’t figure out.Here is the MCNP code.
kuca3.py (43.3 KB)
I have the similar modeling method to the MCNP code, but I’m not get the good agreement to the benchmark :frowning: . In the terminal, the OpenMC sometimes output this:
image
Whether there is a way I can get which lattice couldn’t locate the pariticle?

9

One thing I notice now is that you don’t have any S(a,b) tables listed for your materials. For the PET materials, you need to specify a S(a,b) table so that thermal scattering is treated correctly. I believe the correct table to use is:

material.add_s_alpha_beta('c_H_in_CH2')

Also, I’ll mention that I recently released a tool for converting MCNP inputs, so that could serve as a good check on the OpenMC model you are generating manually.

1 Like
10

Hello, Pual , thanks for your help , the tool you released is awesome. By using the attribute add_s_alpha_beta(‘c_H_in_CH2’), the keff become to
image
which is about 1200pcm lower that the benchmark.
I used the tool to convert MCNP(case Ⅰ-1) code,and got good agreement with the benchmark.And I noticed that the Geometry XML file converted by your tool using the 2D lattice to define each assembly rather than 3D lattice(my building method) didn’t have any warning.I tried to open the particle 7555 HDF5 file, locating the xyz coordinate :
image
and I didn’t find anything strange :

image
image646×638 1.67 KB

I wonder if I have some non-standard definition in 3D lattice or something I have omitted. Would you tell me what is causing the warning and how to solve this warning? By the way, when I replace the ‘air’ material with ‘None’, there will be more warning.

1 Like