Subcritical multiplication problem with k-effective in 'fixed source' mode

Hello all, I am now simulating the ‘fixed source’ mode with several cylindrical sources,just study the influence of different components on the detector.
But when I ran the simulation, I got an error:

ERROR:
The secondary particle bank appears to be growing without bound. You are likely running a subcritical multiplication problem with k-effective close to or greater than one

I tried to change the energy spectrum of the source, reduce the number of sources and even switch to one point source, but nothing changed until I reduced the enrichment of the material.

But this does not fit my situation. Is there any other way to solve this error?

Here are my part of code:

# Create materials for the problem
fuel=[]
fuel.append(openmc.Material(name='UO2_1'))
fuel[-1].set_density(units='g/cm3', density=10.42)
fuel[-1].add_element('U', 1., enrichment=3.0)
fuel[-1].add_element('O',  2.)
fuel[-1].temperature = 900.0

fuel.append(openmc.Material(name='UO2_2'))
fuel[-1].set_density(units='g/cm3', density=10.42)
fuel[-1].add_element('U', 1., enrichment=4.5)
fuel[-1].add_element('O',  2.)
fuel[-1].temperature = 900.0

UO2 = openmc.Material(name='UO2')
UO2.set_density(units='g/cm3', density=10.42)
UO2.add_element('U', 1., enrichment=4.5)
UO2.add_element('O',  2.)
UO2.temperature = 900.0

Gd2O3 = openmc.Material(name='Gd2O3')
Gd2O3.add_element("Gd",2)
Gd2O3.add_element("O",3)
Gd2O3.set_density('g/cm3',7.41)
Gd2O3.temperature = 900.0

fuel.append(openmc.Material.mix_materials([UO2,Gd2O3],[0.92,0.08],'wo',name="fuel+Gd"))

# Define problem settings
r=openmc.stats.PowerLaw(0.0,0.25,1.0)
phi=openmc.stats.Uniform(0.0,2*pi)
z=openmc.stats.Uniform(-fuelLength/2,fuelLength/2)
x=[]
y=[]
m=12
n=4
for i in range(5):
    for j in range(9-i):
        if i==0:
            x.append((3 ** 0.5) * m * L2 / 2 + (j - (8 - i) / 2) * L1)
            y.append((n - 8 + 0.5 * m) * L2 + (3 ** 0.5) * i * L1 / 2)
        else:
            x.append((3 ** 0.5) * m * L2 / 2 + (j - (8 - i) / 2) * L1)
            y.append((n - 8 + 0.5 * m) * L2 + (3 ** 0.5) * i * L1 / 2)
            x.append((3 ** 0.5) * m * L2 / 2 + (j - (8 - i) / 2) * L1)
            y.append((n - 8 + 0.5 * m) * L2 - (3 ** 0.5) * i * L1 / 2)

energy_bins = [1.00E-05, 0.625, 100, 1000, 5.00E+03, 10000, 50000, 1.00E+05, 200000, 300000, 4.00E+05,
               500000, 600000, 700000, 800000, 900000, 1000000, 2000000, 3000000, 4000000, 5000000,
               6000000, 7000000, 8000000, 9000000, 10000000, 20000000]
probability_density = [0.114804354, 0.130632628, 0.073037242, 0.056694686, 0.025809602, 0.06903753,
                       0.038666972, 0.051744062, 0.041062851, 0.034657177, 0.02416415, 0.027668522,
                       0.027214132, 0.025025168, 0.021317823, 0.015144402, 0.108801774, 0.059257734,
                       0.026894671, 0.014246912, 0.007213418, 0.003657553, 0.001703832, 0.00081529,
                       0.00038822, 0.000339297,0]

sources=[]
for i in range(len(x)):
    source = openmc.Source()
    source.strength = 1/len(x)
    source.space = openmc.stats.CylindricalIndependent(r, phi, z, origin=(x[i], y[i], 0.0))
    source.angle = openmc.stats.Isotropic()
    source.energy = openmc.stats.Tabular(energy_bins, probability_density,'histogram')
    sources.append(source)

settings = openmc.Settings()
settings.run_mode = 'fixed source'
settings.source = sources

With thanks, best regards
Richard