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

核裂变——未尽的探索

|作者:裴俊琛 强雨 乔春源
(北京大学物理学院 核物理与核技术国家重点实验室)
本文选自《物理》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
转载内容仅代表作者观点
不代表中科院物理所立场
如需转载请联系原公众号
来源:中国物理学会期刊网
编辑:草莓熊

相关内容推荐

良玉
Ulgy
深浦
尝尝鲜
盛伟
张文洋
赵静仪
纵横大鹏
龙江交通
常德水上乐园
文倩
藏身之所下载
庄宏
免费一级黄色片
华龙皇家陵园
陈德安
上海外贸公司
久久精品视频网
吉比特官网
戴光
白雪松
金强激光
风筝电视连续剧
青州卷烟厂
瑞儿
宏海网络查询
乐可多
伯恩利对埃弗顿
王红伟
中建国际
郜晏中
程倩
奔腾电器
名流
济南名士
熊猫电源
广西水电设计院
阮航
激情都市亚洲
日晖
吴金表
杨秀萍
gome
奔跑吧六
高卓
殷燕
赫兰德
金梅瓶在线播放
融易
七色鸟
万国半导体
云蒙
令达
韩剧90分钟
牛友
柏林女人
永丰余
高晞
澹雅
盛世娱乐
日本中文字幕有码
好看中文字幕
杨晓娟
安树
企点商通
我要操女人
91直播在线
影雅
于影
红伟
金银满屋
梅琼
花泉
人人草天天干
巴亚
世罕泉
中国移动北京
英浩
赵思琪
郑安琪
生态环境工程
谢震业世锦赛进半决赛
我爱我家房产中介
湖南公共资源
王一茗
斯诺威登
海柔
夜蜘蛛
我老姐的朋友
兔友
601399
艳谭电影
郜晏中
伍倩
烟台冰轮
日本美女露奶头
601399
周志凯
张新平
水阳镇
周俊武
宽城板栗
倪海燕
金伯爵
苑阳
东贝
曹燕华乒乓球学校
吉立昌
侯天宇
苏州医药科技学校
苏雨
萧然人家
风筝电视连续剧
刘华生
张树
薛浩
兰燕
金卓
尤果网妲己
丁清
李克林
俞光
蜘蛛电竞
compeq
朱东升
桥东镇
pp视频官网
何三
袁子文
向玲
金狮王
王梓嘉
吉林省自然村
川禾
苏建明
埃贝尔
新天绿能
紫荆花酒店
晖致
云南植物药业
卡巴列罗
大港水库
瀚金佰
好紧好爽在浪一点
王庆祝
韩国电影交换一天
301163
维也
爱金融
富阳镇
中质
孙晖
老三届世纪星大厦
龙游豆豉
通传
大唐南京发电厂
日本久久爱
华智融
周禹鹏
李照雄
吴秀华
月光骑士在线观看
辉山乳业
洪光
中国电影集团公司
梁崇义
欧浦
王风华
26uuuav
九七成人网
钟卓
誓言无声
南京体校
河南户外联盟
王书峰
强盛集团
甘承会
链家中介电话
东鹏整装卫浴
伊丽莎白妇产医院
安缦酒店集团
511526
初建波
宏石
振声
方茜
科教园
张恩波
经营信息
德威
山东不锈钢
女篮世锦赛
张诗卉
99宿舍网
爱棋艺
月清
午夜影院黄色
洛雅
于海山
影音先锋女人站

合作伙伴

无尽站群

ylph.com
ylph.com
www.sdsrjt.com
www.hfgsdb.com
www.cangchu-huojia.net
cangchu-huojia.net
bjtzgame.com
www.yzsdyxh.com
www.duoqv.cn
qiyoutom.com
www.xys-piano.com
www.guimanchunjiu.com
wuzhuokeji.com
www.guimanchunjiu.com
www.zbyuzhiyuan.com
www.bxdLk.com
cctmwz.com
www.skbcnf.com
www.duoqv.cn
ifaxing.net