核裂变——未尽的探索-无尽站群

核裂变——未尽的探索

|作者:裴俊琛 强雨 乔春源
(北京大学物理学院 核物理与核技术国家重点实验室)
本文选自《物理》2022年第11期
摘要   核裂变的发现深刻地影响了人类社会。核裂变的研究还在不断深入,一方面核裂变有新的应用需求,另一方面核裂变是一个复杂的量子多体动力学过程。近年来,核裂变理论和实验研究有很大进展,人们对核裂变几率、裂变产物和裂变机制都获得了新的认识,这有助于澄清一些唯象模型的经验假设。此外,机器学习的应用为发掘利用不精确不完整的核数据提供了可能。期待未来更精确更自洽的核裂变理论可以更好地支撑应用创新。
关键词  核裂变,先进核能,裂变机制
1. 引 言
1939年2月,Meitner与Frisch首次揭示了铀原子核像液滴一样发生了分裂,并用fission这个词来描述核裂变。更重要的是,他们基于玻尔的液滴模型估算出一次核裂变会释放约200 MeV的能量。实际上fission一词最早是指生物学中的细胞分裂。核裂变释放的能量是如此之巨大,很快就引起了科学家们的极大兴趣。1942年12月,费米在芝加哥大学实现了可控的链式核裂变反应,开启了和平利用原子能的时代。1945年7月,美国成功爆炸了第一颗原子弹,深刻地改变了人类历史。核裂变的发现是一个曲折的传奇故事,一些大科学家曾与之失之交臂,它生动地展现了科学认识积累到一定程度后灵光一现的思想突破。
核裂变的发现至今已经80多年,它深刻地影响了人类社会。人们猜测地球内核就是一个巨大的核裂变反应堆,一直保持着人类生存的温暖。核裂变一方面会释放巨大的能量造福人类,另一方面如果控制不好会带来灾难性的影响。这需要我们进一步研究核裂变,更精准地认识核裂变,更好地利用核裂变。北京大学胡济民先生所著的《核裂变物理》一书对核裂变研究进行了全面系统的阐述。近年来,核裂变的理论和实验研究取得了显著进展,产生了一些新的认识,这为核裂变的应用带来了新的可能。
2. 为什么我们还要研究核裂变?
由于地球上U储量有限,发展先进的可持续、更安全、更清洁的核裂变能将越来越重要。目前主流的核能是压水堆,其中裂变产生的中子经过慢化后变成能量很低的热中子。第四代先进核能(图1(a))的一个主流方向是快中子堆,快中子堆无需中子慢化剂,可以更紧凑。快中子堆可以通过增殖反应将U变成易裂变的Pu,将铀资源的利用率从1%提升到60%。同时快中子堆大幅度地减少了核废料的放射性寿命。相比于压水堆,发展快中子堆需要更精确的、中子能量连续的核裂变数据。而目前国际上主要核数据库的核裂变产物的产额只有热中子、0.5 MeV与14 MeV三个能量点的评价数据。更精确的核模型与核数据也有助于设计更精密紧凑的专用核动力(船用、月基、空间,图1(b))和更好的支撑国防研究。
图片
图1 核裂变的新应用,包括先进核能(a)、空间核动力(b)、同位素电池(c)、天体环境下 R-过程中的核裂变(d)、同位素药物(e)、超重元素的合成(f)、反应堆中微子的研究(g)等
除了利用核裂变释放的巨大能量以外,裂变产物核的循环利用将是一个巨大的机遇。通过核裂变产生的Mo可以获得Tc,Tc是用于核医学诊断的重要同位素,已经有广泛成熟的应用(图1(e))。英国科学家利用核废料长期放射性的特点,通过钻石包裹制成了能够稳定供电两千多年的核电池,创新性地实现了变废为宝。利用反应堆还可以生产Pu,已经将其制成同位素电池用于中国的火星车和月球车,但是其生产还很昂贵(图1(c))。反应堆内的核反应网络十分复杂,这也为实现先进核能提供了新的可能。通过核裂变可以产生数百种同位素核,大部分裂变产物核是不稳定的。通过加速核裂变碎片可以形成放射性束流,国际上新一代放射性束流装置的主要科学目标是研究极端条件下的奇特核物质,如美国的稀有同位素束流装置(FRIB),中国的强流重离子加速器装置(HIAF)等,这将极大地扩展核物理的研究范围。长寿命放射性核素在反应堆内会大量积累起来,对反应堆设计、核废料处理、裂变产物循环利用都十分关键。
核裂变对一些重要的基础问题,比如超重新元素的合成、宇宙中元素的演化过程、反应堆中微子的研究,也不可或缺。实验上熔合反应合成的超重核处于高激发态,它的存活概率取决于中子蒸发与裂变之间的竞争。实验上合成超重核极其困难,往往一年才观测到1—2个事例,需要可靠的理论指引(图1(f))。在双中子星并和与超新星爆发的喷射物中,会发生R-过程快中子俘获反应,从而产生重元素(图1(d))。地球上的铀、钚都起源于天体环境下的R-过程,但是R-过程到了极端丰中子超重核区由于裂变而终止,他们裂变的产物又循环参与R-过程从而显著地影响最终宇宙元素的丰度。此外,在核反应堆中裂变产物的β衰变会产生大量的反中微子(图1(g)),对其能谱的观测将揭示一个基本物理问题,即是否存在第4种中微子——惰性中微子。这些新的应用和基础研究都依赖更可靠的核裂变的几率与产物产额。
图片
图2 核裂变过程示意图。其中横轴是随着裂变核拉长的形变增加,红线表示裂变位垒,曲率表示位势的弯曲程度
核裂变虽然是一个老问题,但是从微观角度看,核裂变是一个极其复杂的非平衡非绝热的量子多体动力学过程,如图2所示。传统的唯象裂变模型通过引入一些参数,对实验数据较多的核区能较好地描述,但是无助于深刻理解核裂变以及外推到实验很难达到的核区。原则上微观核裂变理论可以自洽地描述多种裂变观测量,但是微观模型离应用需求的精度还有一定的距离。近年来,随着超级计算机的发展,科学家对核裂变机制获得了一些新认识。发展能描述多种裂变观测量,包括碎片产额、释放动能、释放γ光子数、释放中子数、裂变几率与裂变截面等观测量的综合可靠的微观裂变理论是一个重要科学目标,可以更深刻地理解核裂变过程,并对很难精确测量的核数据和空白核数据提供关键的补充。此外人工智能与机器学习的应用可以帮助我们更好地模拟核裂变和挖掘核数据。近年来实验上提供了前所未有的精确的裂变观测量,为进一步验证、约束和发展新的裂变理论提供了机遇。
3. 核裂变几率
核裂变的寿命或者裂变几率是一个关键的观测量。重核的自发裂变是裂变位垒的量子隧穿过程,这是一个十分缓慢的过程。裂变位垒是指原子核结合能随着核形状拉长而变化的曲线或多维曲面。原子核的多维集体形状空间由于量子壳效应而呈现复杂的裂变位垒。当原子核处于高激发态,量子效应(对关联、壳效应)逐渐消失,裂变几率可以由统计模型描述,裂变寿命为10—10s。随着激发能增加,裂变寿命先是下降很陡,到了高激发时变成缓慢下降,所以裂变的机制是能量相关的。
传统的玻尔—惠勒统计模型,也叫过渡态理论,在核裂变寿命的计算中有广泛的应用,但依赖较多唯象参数。为了描述裂变位垒的能量相关性,唯象模型通常引入一个参数,来描述壳修正能随激发能增加而指数衰减的因子。而高激发态的裂变位垒可以通过微观的有限温度的能量密度泛函计算给出。微观计算可以自洽地考虑量子效应随温度增加而逐渐消失的过程。德国重离子研究中心(GSI)与日本的理化学研究所(RIKEN)通过冷熔合合成了107至113号元素,俄罗斯的杜布纳联合核子研究所(Dubna)通过热熔合合成了114至118号元素。我们的微观计算结果表明,超重复合核的裂变位垒随激发能增加而减小的因子在不同核区十分不同。热熔合区比冷熔合区的裂变位垒下降要慢,在高激发时仍有一定的位垒,这是唯象模型所没有预料到的。相对于冷熔合所呈现的趋势来说,通过热熔合生成超重元素有很大的截面,这在当时是一个让人困惑和质疑的问题。我们在2009年的理论工作对澄清和推动后续热熔合实验做出了贡献。此外微观计算可以描述位垒形状随能量的变化,热熔合区复合核在高激发时的位谷曲率比一般重核的曲率要小4—5倍,而裂变几率正比于位谷曲率。总之基于微观计算,热熔合超重复合核有显著的存活概率,揭示了114—118号超重元素合成的关键因素。我们理论预言合成119、120号新元素的存活概率与118号相似。此外,我们通过微观计算的复合核裂变位垒可以解释兰州近代物理所的丰质子重核区的新核素实验。
基于微观计算的裂变位垒的能量相关性也可以解释实验观测的裂变产物分布的能量相关性。随着激发能增加,裂变模式从不对称裂变逐渐演化到对称裂变。如铀、钚的第二个裂变位垒的高度在引入反射不对称形变后可以降低2.5 MeV,说明低激发时由不对称裂变主导,但是它们的差别随着激发能增加而逐渐减少。
在热浴环境下,核裂变几率可以用虚自由能法(ImF)来计算。该方法在化学反应中有广泛的应用。我们推导了玻尔—惠勒模型与ImF方法之间的联系,发现它们的主要区别是位垒与位谷的能级密度参数的区别。当处于极高的温度时,量子壳效应消失,玻尔—惠勒模型的裂变几率与虚自由能法只相差一个因子,即
图片
。统计模型十分依赖位谷与位垒的能级密度。原子核的能级密度与核的形状和激发能有关,微观计算能级密度还是一个很大的挑战。随着复合核激发能的增加,可能出现发射中子后的裂变。这一定程度上反映了高温核物质有较大的耗散系数和粘滞性,会延缓裂变动力学过程。通过测量裂变前中子发射多重数可以推测出高激发核的耗散系数在增加。
4. 核裂变产物
核裂变会形成很多不同碎片产物的组合,这些裂变产物核的产额分布是非常重要的裂变观测量。此外,裂变碎片会携带很大的动能,这是核能释放的主要形式。裂变碎片处于激发态,会迅速释放中子而冷却。不同裂变碎片释放的中子数、携带的动能和角动量也不同。往往轻质量碎片释放两个中子,而重碎片释放一个中子,这与我们的预想不一样。裂变碎片会通过β衰变形成最终的累积产额。裂变产物的产额与其他裂变后的观测量是关联在一起的。
理论上描述裂变产物的产额分布主要基于多维裂变位能面的形状演化,比如基于多维朗之万方程求解可以合理地描述裂变碎片的质量分布。与花粉在水中的布朗运动相似,在朗之万方程中,核裂变是原子核在集体形变空间的缓慢演化,而核的单粒子运动作为随机背景在快速的变化,这是一个经典的动力学方程,考虑了涨落—耗散效应。此外基于微观裂变位垒的时间相关的生成坐标法(TD-GCM)也能大致描述裂变产物的分布。TD-GCM主要是基于裂变位能势的驱动,原则上微观计算多维位能面可以更合理地描述裂变产额,但是计算量极大。这些基于静态的裂变位垒的形状演化本质上是绝热动力学,不能自洽考虑碎片激发。基于时间相关的密度泛函(TD-DFT)可以描述裂变的非绝热非平衡的动力学过程。TD-DFT是基于微观的单粒子波函数随时间的自洽演化,不需要计算裂变位能面。TD-DFT能自洽计算多种裂变观测量的平均性质,但TD-DFT长期存在的一个问题是无法给出足够展宽的分布,这是因为TD-DFT缺乏集体自由度的涨落。为了解决这一问题,我们提出了裂变过程中单粒子能级随机跃迁的图像。在TD-DFT中,单粒子运动与集体运动是交织在一起的,随着有效温度增加,随机跃迁的效应增加,经过长时间的演化累积而得到有展宽的分布。此前国际上提出了随机平均场模型来描述产额分布,是基于很大的初始涨落,但是这与自发裂变矛盾。
图片
图3 (a)基于Brosa模型对裂变产物分布的描述;(b)微观计算的Pu核裂变动力学演化路径(其中单位b表示10 m)
唯象的Brosa模型从裂变产物的质量分布出发,认为存在两种不对称的裂变模式(图3(a)),并认为不同的不对称裂变模式的起源是受到裂变位垒的影响。Brosa模型还通过颈部随机断裂来描述裂变观测量的展宽,脖子越长分布越宽。Brosa模型可以合理地解释裂变产额分布、总动能分布、中子发射多重数之间的关联,是物理直觉的很大成功,但是一直缺乏微观理论的支持。我们的结果揭示了动力学涨落效应正是Brosa模型中的S1、S2两种不对称裂变模式的起源。如图3(b)所示,随着涨落增加,长脖子S2裂变道的成分在增加,这与实验是一致的。这两种模式的裂变路径相似,不大可能是静态位垒的影响。
图片
图4 (a)基于贝叶斯机器学习对U的不完整裂变产额的评价;(b)对U裂变碎片Xe的产额—能量关系的评价
近年来人工智能与机器学习在很多学科中的交叉应用获得了很大关注。实验上测量的裂变产物的产额往往是不完整的,有误差或存在分歧。特别是能量相关的裂变产额,在中子入射能量处于2 MeV与14 MeV之间的数据比较稀少。在这种背景下,我们提出基于贝叶斯机器学习来学习补充缺失的裂变产额数据(图4(a)),展示了机器学习的实际应用价值和优势。此外通过输入碎片的电荷奇偶信息,以及在学习中引入负值惩罚等,将物理信息和物理约束与机器学习进行了尝试结合。最近我们提出通过数据融合来更好地评价不完整、有分歧、有误差的核裂变产额。在反应堆中,裂变产物Xe有很大的中子吸收截面,是反应堆的“毒物”,会显著降低反应堆运行功率,其产额的评价很重要。图4(b)是我们基于贝叶斯机器学习对Xe产额的评价。当一种裂变数据在某些能区很稀少时,它与别的数据在其他能区的关联有助于这种数据的推断。数据融合可以考虑潜在的、高维的、非局域的关联,可以给出综合的误差传播,可以发掘出不精确的裂变数据的最大价值,有望形成新的核数据评价方法。
5. 核裂变机制
核裂变是一个极其复杂的非平衡非绝热的量子多体动力学过程,裂变后碎片之间存在量子纠缠。更深入地认识裂变机制有助于发展精确的裂变理论。TD-DFT理论最适合研究核裂变的微观机制。近年来,超级计算机的应用为微观裂变动力学的发展提供了很好的机遇,使我们有可能澄清或者更新一些唯象的裂变模型的经验图像。
在TD-DFT的基础上,对关联是最重要的剩余相互作用,对裂变机制有重要影响。人们认识到对关联相当于裂变的“润滑剂”,可以加速裂变过程。静态对关联也会导致裂变位垒的稍微降低,这会显著减少自发裂变的寿命。当核的密度缓慢变化时,动力学对关联有快速的涨落。当对关联非常大时,核体系形成一个超流的集体态,涨落效应被压制。在TD-DFT计算中,当对关联弱时,可能出现三分裂核裂变,而对关联强时则是二分裂。当对关联很弱时,裂变路径往往走短脖子的S1裂变道,与实验不符。当复合核处于高温激发时,对关联会很快衰减,涨落也更迅速。
核物质的耗散系数或者粘滞系数到底有多大呢?在著名的沥青滴漏实验中,由于沥青有很大的粘滞性,形成一个脖子拉伸很长的液滴,约10年才滴出一滴,如图5(a)所示。我们认为核物质的粘滞性比沥青小,但是比水大。通过把微观动力学裂变路径映射到经典动力学方程,可以提取出形状相关的耗散系数。这个计算需要提取出动力学的裂变位势,相比于静态位垒,动力学的裂变断点更远,碎片间的库仑能更小。这意味着相比于非绝热裂变,绝热计算的总动能会显著偏大,这也验证了非绝热裂变的合理性。我们的结果表明,裂变过程中耗散系数也在发生变化,一般在2×10—4×10s,这有助于约束唯象模型的耗散参数。耗散系数随着激发能增加而增加,随对关联增加而减小。随着激发能增加,由于很强的耗散,裂变动力学演化时间越来越长。这时涨落效应成为裂变的主要驱动机制,这与涨落—耗散定理是一致的。很强的耗散和粘滞性将导致一个拉长的裂变颈部构型(图5(b)),从而使库仑能降低,导致裂变释放的总动能减少。由于能量守恒,总动能减小,导致碎片的激发能增加,从而发射更多中子。微观TD-DFT计算能合理地解释裂变机制的能量相关性。
图片
图5 (a)沥青滴漏实验展示出有拉长脖子的液滴;(b)微观TD-DFT计算给出的Pu裂变的断点构型
裂变分裂成的两个碎片之间存在很强的动力学纠缠,主导着碎片之间的能量分配、核子数分配、角动量关联等。由于断裂的碎片很快飞开,以至于部分的动力学纠缠还来不及塌缩。图6分别展示了不同裂变碎片发射中子数的平均值随碎片质量变化的分布、碎片中的中子与质子的平均比值,以及不同碎片所携带的平均角动量的分布。碎片发射中子的几率主要由碎片的激发能决定。TD-DFT计算得到的轻碎片激发能大于重碎片激发能,且它们的差别随激发能增加而减少。实际上碎片发射中子的多重数分布是一个型的锯齿结构(图6(a)),其能量分配机制还有待微观理论的解释。裂变是颈部随机断裂还是量子纠缠主导呢?唯象Brosa模型认为是随机断裂主导的。基于碎片之间的纠缠,两个碎片的核子数分配和能量分配具有不确定性。微观计算通过粒子数投影也能获得有展宽的分布。近年来法国实验组取得了很大的进展。实验上可以获得所有碎片的产额分布,其中碎片的中子与质子的比值也具有锯齿结构(图6(b)),但与中子多重数的锯齿结构相反,这为进一步认识裂变断裂机制提供了观测量。此外最近实验上获得了裂变碎片的角动量分布(图6(c)),这是一种新的观测量,引起了理论上很大的关注。美国一些理论组很快对裂变碎片角动量分布提出了多个解释。需要指出的是,美国在裂变理论模型方面有长期积累的优势。裂变碎片角动量的获得是与断裂模式(比如脖子的扭曲或弯曲断裂)有关呢,还是断裂后获得的?碎片的角动量分布也有相似的锯齿结构,可能主要由能量分配机制主导。原则上TD-DFT可以描述多种裂变观测量之间的关联,但是还需额外考虑超越平均场的效应,最终形成一个综合、自洽、可靠的微观裂变理论。
图片
图6 (a)实验上观测的裂变碎片的平均发射中子数的分布;(b)碎片中平均的中子/质子比;(c)碎片所携带的平均角动量
6. 总结与展望
核裂变是一团强关联的量子物质分裂成两块的独特的量子动力学过程。核裂变的发现至今已经有80多年,但是核裂变过程非常复杂,对它的认识还有待进一步深入。裂变过程中既有单粒子自由度,也有集体自由度、集团自由度,还有涨落—耗散效应等交织在一起。裂变断裂前的脖子构型对裂变后观测量有重要影响。裂变断点既有随机性,也存在碎片之间的量子纠缠。断裂前体系的裂变位垒、能级密度、耗散系数也具有能量相关性,此外对关联对裂变机制有重要影响。我们看到微观的TD-DFT可以成功地解释裂变机制,有助于澄清或更新一些唯象的裂变模型的图像。基于BBGKY框架,进一步考虑更高阶的关联动力学,可以更现实地描述裂变。随着计算能力的增加,考虑高阶关联的裂变动力学将是一个重要方向。目前微观核裂变理论基于有效核力,存在一定的不确定性,包括裂变位垒的预言也存在误差。从现实核力出发,发展从头计算(ab initio)核结构是核物理的前沿方向。基于ab initio计算重核裂变过程还很遥远,但是可以为发展更精确的有效核力和有效哈密顿量提供指引。
无疑,核裂变的研究有很强的应用背景。为了应对气候变化,先进核能的发展将受到更大的重视。发展更精确可靠的核裂变理论,对升级核能应用十分重要,对一些重大的基础研究也很关键,比如超重元素合成、天体环境中的R-过程、中微子研究等。核裂变的研究既是迷人的量子多体问题,也有很强的交叉应用需求。随着超级计算、机器学习、量子计算的应用以及实验装置的发展,核裂变的基础研究迎来了新的机遇,将为核裂变应用提供新的线索。近年来,美国、法国等在核裂变基础研究方面取得系列进展,中国在核裂变的基础研究方面还有很大的发展空间和前景。
参考文献
[1] Meitner L,Frisch O R. Nature(London),1939,143:239
[2] 胡济民. 核裂变物理学. 北京:北京大学出版社,2014
[3] Kolos K,Sobes V,Vogt R et al. Phys. Rev. Research,2022,4:021001
[4] 核物理与等离子体物理—学科前沿及发展战略 . 北京:科学出版社,2017
[5] http://hiaf.impcas.ac.cn/
[6] 周善贵. 物理,2014,43(12):817
[7] Eichler M,Arcones A,Kelic A et al. Astrophys. J.,2015,808:30
[8] Mention G,Fechner M,Lasserre Th et al. Phys. Rev. D,2011,83:073006
[9] Bender M,Bernard R,Bertsch G et al. J. Phys. G,2020,47:113002
[10] Bohr N,Wheeler J A. Phys. Rev.,1939,56:426
[11] Pei J C,Nazarewicz W,Sheikh J A et al. Phys. Rev. Lett.,2009,102:192501
[12] Zhu Y,Pei J C. Phys. Rev. C,2016,94:024329
[13] Qiao C Y,Pei J C. Phys. Rev. C,2022,106:014608
[14] Ma L et al. Phys. Rev. C,2022,106:034316
[15] Pei J C,Zhu Y. Nucl. Phys. Rev.,2017,34(1):87
[16] Liu L L,Wu X Z,Chen Y J et al. Phys. Rev. C,2019,99:044614
[17] Randrup J,Moller P,Sierk A J. Phys. Rev. C,2011,84:034613
[18] Regnier D,Dubray N,Schunck N et al. Phys. Rev. C,2016,93:054611
[19] Zhao J,Nikšić T,Vretenar D et al. Phys. Rev. C,2019,99:014618
[20] Lu B N,Zhao J,Zhao E G et al. Phys. Rev. C,2014,89:014323
[21] Negele J W,Koonin S E,Möller P et al. Phys. Rev. C,1978,17:1098
[22] Bulgac A,Jin S,Roche K J et al. Phys. Rev. C,2019,100:034615
[23] Qiang Y,Pei J C,Stevenson P D. Phys. Rev. C,2021,103:L031304
[24] Lacroix D,Ayik S. Eur. Phys. J. A,2014,50:95
[25] Brosa U,Grossmann S,Müller A. Phys. Rep.,1990,197:167
[26] Wang Z A,Pei J C,Liu Y et al. Phys. Rev. Lett.,2019,123:122501
[27] Qiao C Y,Pei J C,Wang Z A et al. Phys. Rev. C,2021,103:034621
[28] Wang Z A,Pei J C. Phys. Rev. C,2021,104:064608
[29] Wang Z A,Pei J C,Chen Y J et al. Phys. Rev. C,2022,106:L021304
[30] Qiang Y,Pei J C. Phys. Rev. C,2021,104:054604
[31] Schmidt K H,Estienne M,Fallot M et al. EPJ Conf.,2021,256:00015
[32] Martin J F et al. Phys. Rev. C,2021,104:044602
[33] Ramos D et al. Phys. Rev. Lett.,2019,123:092503
[34] Wilson J N et al. Nature,2021,590:566
[35] Verriere M,Schunck N,Regnier D. Phys. Rev. C,2021,103:054602
转载内容仅代表作者观点
不代表中科院物理所立场
如需转载请联系原公众号
来源:中国物理学会期刊网
编辑:草莓熊

相关内容推荐

欧阳波
唇之下
以岭药业
张宝胜
电影本能在线观看
来秀
吴国胜
一亩良田
高中老师电影
云涧
99大香蕉视频
韵婷
铁翅
朋友的丈夫电影
sifang
虎道
香港三级黄色视频
张传波
隋侯珠
瑞萨半导体
上海物贸
采星
月照
皓彩
糊涂酒业
叶文玉
李正太
曼卡
美女互操
林跃平
益丰大药房高毅
铁中快运官网
中拓互联
禁忌第一季
杨超超
吴中杰
华美整形美容院
解放军总后勤部
吴江农商行
黄色a片免费观看
人邮
朱拥军
凶手还未睡完整版
郑振宇
陈雪良
搜我
综合色惰
关庆维简介
浮屠镇
源之安达速运
韩国片年轻的妈妈
沈昌
董君雅
沈深
点微
林敏聪
金流
瑞丽服饰
卡特彼勒公司
祁军
中华英才半月刊
谭跃
王泓钦
金惠玉
刘启雄
赵国珍
香元
午夜福利导航视频
大乐透走势图综合版
西红柿影院
家居城
西安隆基
黄片免费电影
林德雄
宁波物流公司
深能集团
金兰香
王会林
李道亮
刘洪宇
新特药
安永会计事务所
友聚
wanyu
王梓均
金飞达
冰野
免费jiZZ
林丹青
陈红娟
迷人的保姆免费
苏易登
鲜果粒
成人免费手机电影
方微
山东核电有限公司
重庆考试教育院
胡婷
吴世宏
鑫达
芯瑞达
我屋
一丁集团
柠檬之恋
大漠长河
王屋山黑加仑
5号特工组2
色在线看
情人岛度假村
洪韬
斯丹德
米米七月
桃花影
久草国产在线
办公室恋情韩国
在线午夜福利电影
知风
徐艳秋
众润多
鹿梦
李星雨
松野
谭冲
文宇
金陵体育
影视先锋色
贵州省委组织部
马晓雪
西牛
魔装学园动漫
妈妈的朋友19
好棋牌
铸金
崔金龙
远东资信
李婵娟
申露
饿了么公司
小姿
徐志雷
张立平
紫金龙
束越新
闽侯四中
北京12315
陈一兵
高宗祥
龙岩三杰
金保华
杨蕤
四零影视
胡茂林
刘学兵
全新好
诺心
硕德
江永忠
电影处女
林伟鹏
国林网
李社教
金融资产投资公司
赛斯
新易盛官网
黄启
陆婷
和德华
我要看电视剧全集
蒋妍
刘伟奇
尚美装饰
李宏远
易鑫公司
货捕头
青藤网
周韬
董建刚
壹梦
仙人桥镇
雪山红
欧美做爱a片
林美玲
李金泽
仙人桥镇
晋城城建
平顶山人才市场
鸿达
安伏
宗昌
美丽之旅
爱欲网
史金龙
陈燕辉
高婷
朱丁
裕华集团
江能
神马伦理网
南通二建集团
亚洲人成人

合作伙伴

无尽站群

sdsrjt.com
www.dgsenshuo.com
bxdLk.com
www.suqinbeiye.com
www.sdsrjt.com
sdsrjt.com
xys-piano.com
www.fxnihao.com
dgsenshuo.com
www.ylph.com
www.yxcandle.com
www.fxnihao.com
bxdLk.com
www.hfgsdb.com
tzystec.com
www.zbyuzhiyuan.com
www.zhujiaohui.com
www.xys-piano.com
tzystec.com
www.yzsdyxh.com