Variance reduction method based on adjoint discrete ordinate
-
摘要: 对于深穿透类型的屏蔽计算，为了得到较为可信的统计结果，蒙特卡罗方法(MC方法)需要模拟大量的粒子，巨大的计算时间是其存在的主要问题。源偏倚和权窗技巧能够有效降低深穿透问题的计数误差。开展了基于共轭离散纵标(SN)的MC减方差方法研究，根据SN方法的共轭注量率计算并生成了源偏倚和权窗参数，编写了JMCT程序的源抽样子程序，并且在秦山一期测量值基础上进行了验证，成功应用到CAP1400压力容器快中子注量率和堆腔中子和光子剂量率计算中。数值结果表明，对于深穿透屏蔽计算问题，和无偏的MC方法相比，基于共轭SN的MC减方差方法能够在保证结果精度的前提下，提高计算效率1~2个量级。Abstract: For deep-penetration shielding calculation, Monte Carlo method (MC method) requires modeling a great number of particles to obtain reliable results, thus huge computation time is the main problem of the MC method. Source biasing and weight window technique effectively decrease the tally error of deep penetration problem. This paper studies the variance reduction (VR) method based on adjoint Discrete Ordinate (SN), generates source biasing factors and weight window parameters for the MC method by using the adjoint fluence rates of the SN method, and develops source sampling subroutine for JMCT. The VR method was verified at phaseⅠ of Qinshan Nuclear Power Plant measurements. It was applied to CAP1400 pressure vessel fast neutron fluence rate and cavity neutron and photon dose rate calculations. Numerical results show that the VR method based on SN increases calculation efficiency by 1~2 orders for deep-penetration shielding calculation with high precision compared with unbiased MC method.
-
Key words:
- deep-penetration /
- discrete ordinate method /
- Monte Carlo method /
- source biasing /
- weight window
-
静止无功补偿器(SVC)相对于机械投切电容器组式的无功补偿装备，实现了无功功率精确、平滑补偿，是目前较为先进实用的无功补偿装置，已经得到了广泛应用。另外，SVC装置具有良好的分相补偿能力，对于电弧炉工作中造成的三相电压、电流不平衡有着独特的治理效果; 且晶闸管比全控型半导体开关器件造价低，有着更高的额定电压电流参数，适合拓展于高压、大容量场合。静止无功发生器(SVG)有着更快的响应时间，且可以不借助于无源支路，灵活的产生容性或感性无功功率进行调节。虽受限于目前全控型开关器件的电压容量及经济成本，然而在已经装设了SVC装置的工业现场，不改动已有的无源滤波支路，使用SVG配合SVC进行剩余无功功率补偿，具有非常大的应用拓展空间。本文针对广西某钢厂电弧炉现场，通过PSCAD/EMTDC对SVC和SVG联合方案进行了系统仿真验证[1-4]。
1. SVC与SVG无功补偿的基本原理
工业现场SVC有多种类型，主要有五种型式：磁阀控制电抗器型(MCR型)、可控硅控制空芯电抗器型(TCR型)、可控硅阀控制高阻抗变压器型(TCT型)、可控硅开关控制电容器型(TSC型)及自饱和电抗器型(SSR型)。基本原理及仿真主要围绕TCR型SVC。而本文所述的某钢厂电弧炉现场的负荷及SVC+SVG补偿装置的电压等级均为35 kV，并联于同一条母线下。
1.1 基于斯坦门茨理论的SVC补偿原理
TCR型SVC的基本原理是以母线三相电压矢量为基准，SVC装置发出与负荷基波无功电流幅值相等、相位相反的基波电流抵消负荷无功，而无源支路的容性无功容量与负荷最大感性无功基本持平，晶闸管阀控制单元根据负荷功率的变化调控电抗器感值灵活控制SVC需要补偿的无功功率[5]，其基本结构如图 1所示。
根据斯坦门茨(Steinmetz, Charles Proteus)平衡化补偿原理[6]，控制单元根据负荷无功计算出相应的电纳值，然后由电纳值推导出相应的触发角发送到晶闸管阀组。电纳值B与触发角α的关系式为
B(α)=1ωL2π+sin2α−2απ (1) 式中：ω为角频率; L为电抗值。
而斯坦门茨平衡化原理是由三相参考电压与负荷电流分别求出各相的电纳值。其基本关系式为
Bab=(Re(˙Ubc˜˙Ial)+Re(˙Uca˜˙Ibl)−Re(˙Uab˜˙Icl))√3U2ab (2) Bbc=(Re(˙Uca˜˙Ibl)+Re(˙Uab˜˙Icl)−Re(˙Ubc˜˙Ial))√3U2bc (3) Bca=(Re(˙Uab˜˙Icl)+Re(˙Ubc˜˙Ial)−Re(˙Uca˜˙Ibl))√3U2ca (4) 式中：˙Uab,˙Ubc,˙Uca为母线三相参考电压矢量; ˜˙Ial,˜˙Ib1,˜˙Icl是负荷三相相电流矢量。
根据斯坦门茨原理公式(2)，(3)和(4)求出的电纳B值，由式(1)查表，分别求出三相触发角，实现无功开环算法[7]，其结构如图 2所示。
1.2 基于瞬时功率理论的SVG补偿原理
SVG，相较于基于半控型晶闸管控制的SVC装置，有着更快的响应时间。不同于根据斯坦门茨原理来计算电路的等值电纳，SVG装置使用脉冲调制技术来驱动开关管器件。而目前常用的控制技术主要有基于快速傅里叶变换的电流控制技术和基于瞬时无功功率理论的控制方法。由于基于傅里叶变换的控制技术在原理上至少存在一个周期的延时，因此在快速响应的补偿场合，大多使用基于无功功率理论的控制方法。图 3所示典型SVG逆变器的电路拓扑。
求解上述电路的数学模型，并将态电路函数先进行PARK变换到d-q域，然后进行拉普拉斯变换，可有
{(sL+R)id=usd−urd+wLiq(sL+R)iq=usq−urq+wLidudc=32sC(Sdid+Sqiq) (5) 式中：id，iq为ia，ib，ic经过PARK变换在旋转p-q坐标系下的分量; Sa，Sb，Sc为三相桥臂的开关状态，S＝1表示该相桥臂上管开通，下管关断; S＝0表示该相下管导通，上管关断，经过p-q变换后为Sd，Sq。
针对参考电流，根据瞬时无功功率理论，如图 4所示，先经p-q变换求出其在旋转坐标系下的直流量，然后经过低通滤波器、p-q反变换得出其基波无功电流的参考量[8]。
根据瞬时无功功率理论得出参考电流量后，逆变器的控制器和实际控制框图如图 5所示。
在工业现场，通常将SVG装置直接并联于负荷母线下进行无功补偿[9-10]。本文提出的方案将SVG与SVC并联然后并联于同一母线的负荷下。
2. 广西某钢厂电弧炉的无功补偿现状
广西某钢厂110 kV高压变电站母线接入最小短路容量576 MV·A，安装50 MV·A主变压器一台，主变压器中压侧负荷主要为电弧炉和精炼炉负荷，低压侧主要为高压电动机负荷，目前场内35 kV母线侧安装TCR+FC型SVC，设置2次C型、3次、4次及5次单调谐滤波支路，总基波补偿容量25 Mvar(兆乏)，系统主接线如下图 6所示。
而该钢厂电弧炉经过改造后，发出的无功功率远超原始无功量，因此需重新计算，并于实际比对。无功需量计算应综合考虑初炼交流电弧炉(EAF)和精炼电弧炉(LF)的最大无功发生量和同时率，工程推荐值按照下式确定
QD=k1QEAF,max+k2QL,F,max (6) 式中：k1为EAF无功发生量的计算系数，工程推荐值为0.9~1.2;k2为LF无功发生量的计算系数，工程推荐值为0.4~0.6;QEAF, max为EAF最大无功发生量，单位为兆乏(Mvar); QLF, max为LF最大无功发生量，单位为兆乏(Mvar); QD为无功需量，单位为兆乏(Mvar)。
根据计算，EAF，LF同时运行时35 kV所需最大无功发生量为35.25 Mvar。
35 kV母线所带轧机总装机容量为11 MV·A，取功率因数为0.7，功率因数提高至0.92，按公式可知
Q=P(√1cosφ21−1−√1cosφ22−1)=4.57Mvar (7) 式中：φ1为EAF的功率因数角; φ2为LF的功率因数角。因此，35 kV母线总无功补偿容量为39.82 Mvar，考虑一定的补偿裕度并结合电能质量测试报告，最终确定主变35 kV母线SVC装置的基波补偿容量为40 Mvar。针对原先25 Mvar的SVC系统，出现将近15 Mvar的无功缺口。
3. 基于PSCAD/EMTDC的SVC/SVG联合应用仿真
SVC具备良好的分相补偿能力，能更好适应三相不平衡工况，但响应时间较长; 而SVG有着更快响应时间。结合二者优点，针对该钢厂电弧炉现场无功补偿实际，仍可利用原先的SVC装置，由SVG填补电弧炉改造产生的增多的无功功率。将SVG单元前置与SVC单元之前，两装置独立检测、独立工作。先投入SVC及SVG装置，稳定之后一定时间内投入负荷无功源，使用PSCAD/EMTDC进行系统仿真，电弧炉负荷利用无功波动源模拟，仿真步长设为156 μs。系统的仿真结构图如图 7所示，SVG，SVC和无功波动源都接于35 kV母线。
SVC控制算法是根据式(1)，(2)，(3)，(4)编写而成的C语言模块，然后由FORTRAN文件编写接口[11-12]。而SVG单元是由基本模块串联而成的三相桥[13]，控制模块根据式(5)及图 4、图 5所示的控制算法搭建。SVC及SVG控制模块如图 8和图 9所示。
无功波动源是最大0~40 Mvar可变的三相不平衡感性无功波动源，充分模拟该钢厂的电弧炉无功功率变化。图 10为投入无功源之后TCR环内电流，可以看出三相TCR环内电流并不同步，响应了相应的不平衡。
图 11为无功源投入之后，SVG发出的无功功率。SVG发出最大近16 Mvar的容性无功功率，且感性无功值基本均小于0.3 Mvar，从另一个角度可以看出SVG响应速度优异。
为了考察SVG补充功率的真实效果，首先投入SVC，然后投入波动无功负荷，并不投入SVG，仿真结果如图 12所示。
在图 12中，QL为负荷无功波动源无功功率，Q1_35 kV为35 kV母线处的无功功率，在4.5 s负荷产生较大的无功波动，而此时只投入SVC，母线处仍有将近18 Mvar的感性无功，这种现象在5.5 s处负荷波动中也有类似的情况，证明在进行电弧炉改造后，SVC系统并不能完全补偿负荷的无功功率; 而在5.3s处，母线有将近20Mvar的容性无功过冲，对比图 11可知，此处的容性无功冲击在SVG并没有呈现，且不是负荷造成，因此是由于SVC不能及时响应造成的，也说明SVC较SVG响应较慢[14-15]。
进一步将SVC和SVG联合注入后，如图 13所示。与图 12相比，在4.5 s和5.5 s处母线的无功缺口功率基本得到了抑制，结合图 11，SVG起到了快速抑制剩余无功的作用。
4. 结论
电弧炉现场无功补偿系统进行SVC和SVG联合应用探索，需要明确两种装置各自的功能空间，即以TCR配套FC为无功功率补偿主体，并特别针对不平衡、负序等工况; 而SVG以其优良的更快速的响应性能，承担剩余无功功率抑制的任务。在这种分工前提下，参考补偿点位置至关重要，本文提出的方案SVC装置是无法捕捉到SVG功率补偿信息，即SVG参考位置在SVC上端，进行剩余无功功率补偿。但在实际中，SVC和SVG不可能同时响应，存在一定延时，因此对于容量有限的SVG，必须做前置限幅来防止过容量补偿指令带来的系统故障保护。处理完善参考点之后进行了PACAD/EMTDC系统仿真，从仿真结果来看，SVC和SVG联合应用可以起到预定的功率补偿效果，为其向广西电弧炉现场应用提供了理论和仿真级参考。
-
图 1 基于共轭SN的MC减方差方法流程图
Figure 1. Flowchart of MC variance reduction method based on adjoint SN
图 3 秦山一期堆芯中平面JSNT网格划分示意图
Figure 3. JSNT mesh for core mid-plane of phase Ⅰ of Qianshan Nuclear Power Plant
图 4 秦山一期前100万无偏和偏倚源粒子分布示意图
Figure 4. First million unbiased and biased source particles for phase Ⅰ of Qianshan Nuclear Power Plant
图 5 CAP1400压力容器内表面快中子注量率JMCT和JSNT计算模型示意图
Figure 5. JMCT and JSNT models for CAP1400 RPV calculation
图 6 CAP1400压力容器内表面JSNT共轭源和JMCT前100万源粒子示意图
Figure 6. Adjoint source of JSNT and first million biased source particles of JMCT
图 7 CAP1400压力容器内表面快中子注量率(E≥1.0 MeV)分布
Figure 7. FNFR distribution (E≥1.0 MeV) on inner RPV of CAP1400
图 8 CAP1400堆腔剂量率JMCT和JSNT计算模型示意图
Figure 8. JMCT and JSNT models for CAP1400 cavity dose rate calculation
表 1 秦山一期快中子注量率(E≥1.0 MeV)计数误差小于0.05时计算总时间和加速比
Table 1. Calculation time and speedup for different cases with fast neutron fluence rate (FNFR) tally error less than 0.05
case geometry type number of mesh fast neutron fluence rate/(109 n·cm-2·s-1) JSNT time/min JMCT time/min total time T/min speedup (T_bias/ T_unbias) measure - - 7.41 - - - - 1(unbiased) - - 6.08 - 19 814.12 19 814.12 1 2 xyz 393 120 6.48 19.66 3 575.56 3 595.22 6 3 xyz 393 120 6.35 19.66 197.92 217.58 91 4 xyz 393 120 6.41 19.66 30.00 49.66 399 5 rtz 402 948 5.92 80.83 100.78 181.61 109 6 rtz 823 368 6.44 161.66 429.62 591.28 34 7 rtz 1 857 492 6.53 397.59 185.07 582.66 34 
下载: 导出CSV
表 2 CAP1400压力容器内表面快中子注量率各个计数误差区间内的网格数目
Table 2. Number of mesh in each tally error bin on inner RPV of CAP1400
tally error number of mesh upper lower unbiased biased 0.00 0.01 0 0 0.01 0.02 0 1 0.02 0.03 0 311 0.03 0.04 0 350 0.04 0.05 343 119 0.05 0.06 322 48 0.06 0.07 114 16 0.07 0.08 50 8 0.08 0.09 21 1 0.09 0.10 7 1 0.10 1.00 3 5
下载: 导出CSV
表 3 CAP1400堆腔剂量率中子和光子剂量率及加速比
Table 3. Neutron dose rate, photon dose rate and speedup for CAP1400 cavity dose rate calculation
dose rate/(mSv·h-1) tally error/10-2 sum of time/min M/10-2 JMCT speedup(M_biased/M_ref) case1: particle type neutron 2.82×106 7.14 10 612.57 1.85 1 biased source (reference) photon 6.82×104 5.35 10 612.57 3.29 1 case2: neutron and photon neutron 2.80×106 1.30 10 384.53 57 31 fluence rate response function photon 6.50×104 3.15 10 384.53 9.70 3 case3: neutron and photon neutron 2.82×106 1.15 10 064.27 75.1 41 dose rate response function photon 6.27×104 1.74 10 064.27 32.8 10
下载: 导出CSV
-
[1] Wagner J C. Acceleration of Monte Carlo shielding calculations with an automated variance reduction technique and parallel processing[R]. University Park, 1997. [2] Wagner J C, Haghighat A. Automated variance reduction of Monte Carlo shielding calculations using the discrete ordinates adjoint function[J]. Nuclear Science and Engineering, 1998, 128: 186-208. doi: 10.13182/NSE98-2 [3] Haghighat A, Wagner J C. Monte Carlo variance reduction with deterministic importance functions[J]. Progress in Nuclear Energy, 2003, 42(1): 25-53. doi: 10.1016/S0149-1970(02)00002-1 [4] Blakeman E D, Peplow D E, Wagner J C, et al. PWR facility dose modeling using MCNP5 and the CADIS/ADVANTG variance-reduction methodology[R]. ORNL/TM-2007/133, 2007. [5] Peplow D E, Evans T M, Wagner J C. Simultaneous optimization of tallies in difficult shielding problems[J]. Nuclear Technology, 2008, 168: 785-792. [6] Wagner J C, Peplow P E, Evans T M. Automated variance reduction applied to nuclear well-logging problems[J]. Nuclear Technology, 2009, 168: 799-809. doi: 10.13182/NT09-A9309 [7] Vasiliev A, Ferroukhi H, Zimmermann M A, et al. Development of a CASMO-4/SIMULATE-3/MCNPX calculation scheme for PWR fast neutron fluence analysis and validation against RPV scraping test data[J]. Annals of Nuclear Energy, 2007, 34: 615-627. doi: 10.1016/j.anucene.2007.02.020 [8] Fero A H. SORCERY user manual[R]. LTR- R EA-06-74, 2006. [9] X-5 Monte Carlo Team, MCNP - Ageneral Monte Carlo N-Particle transport code[R]. LA-UR-03-1987, 2003. 期刊类型引用(2)
1. 刘利，左应红，牛胜利，朱金辉，李夏至. 中子在大气中产生氮俘获γ的蒙特卡罗模拟研究. 强激光与粒子束. 2022(08): 162-168 . 
本站查看2. 王梦琪，郑征，梅其良，黎辉，程汤培. 全局减方差方法在乏燃料干式贮存容器屏蔽计算中的应用. 原子能科学技术. 2019(05): 884-892 . 百度学术其他类型引用(1)
-
下载:
点击查看大图
计量
- 文章访问数: 940
- HTML全文浏览量: 206
- PDF下载量: 91
- 被引次数: 3
