李胜强 , 谭铭 , 张展博
清华大学 核能与新能源技术研究院, 北京 100084

收稿日期：2017-12-08
作者简介：李胜强(1974—), 男, 副研究员。E-mail:sqli@tsinghua.edu.cn


摘要：为优化加速器驱动次临界系统（ADS）靶件束窗形状设计，该文研究了一类含黏性力的最速降线问题。采用变分法进行分析，最终需数值求解一组非线性积分方程。针对Newton迭代法在高黏性系数情况下难以获得收敛解且无法获得近似解的缺陷，提出一种最优化求解思路。通过数值试验比较了4种启发式最优化算法的计算效率，并选择粒子群算法与Newton迭代法作进一步对比。结果表明：在高黏性系数情况下粒子群算法计算效率优于Newton迭代法，并且粒子群算法在该问题上近似线性收敛。用最速降线取代半椭圆的束窗形状后，流动分离更晚且流道滞止区更小，有利于提升换热效率。
关键词：加速器驱动次临界系统(ADS)最速降线黏性力粒子群优化算法
An optimization method of brachistochrone problem with viscous friction and its application in ADS design
LI Shengqiang, TAN Ming, ZHANG Zhanbo
Institute of Nuclear and New Energy Technology, Tsinghua University, Beijing 100084, China


Abstract: The brachistochrone problem with viscous friction was studied to optimize the beam window (BW) shape of an accelerator driven sub-critical system (ADS) target. A set of nonlinear integral equations was derived using the variational method. The Newton iteration method could not get a convergent solution or an approximate solution for the highly viscous conditions, so an optimization method was developed. The particle swarm optimization (PSO) algorithm was found to be more efficient for the highly viscous conditions than other three heuristic algorithms with approximately linear convergence. The flow separation is later and the stagnation region is smaller for the brachistochrone BW instead of the semielliptical BW, which enhances the heat transfer.
Key words: accelerator driven sub-critical system (ADS)brachistochroneviscous frictionparticle swarm optimization (PSO) algorithm
在核电领域，加速器驱动次临界系统(accelerator driven sub-critical system，ADS)^[1]由于其固有安全性和强大的核废料嬗变能力成为目前的研究热点。ADS由加速器、散裂靶、次临界堆芯3部分组成，其中散裂靶是ADS的关键组件，其结构如图 1所示。

图 1 ADS散裂靶结构示意图

图选项

加速器产生的高能质子束流轰击散裂靶靶件(通常是铅或铅铋合金^[2])，释放中子维持次临界堆芯中的链式反应。其中的质子束窗用于隔离真空束流管和散裂靶靶件，其表面存在很大的热负荷，需及时冷却。如何设计质子束窗的形状以提高冷却效率是ADS设计的一个技术难题。相较于经典的圆形束窗，Sugawara等^[3]提出采用两段圆弧组合的束窗形状以降低局部最高温度；徐长江等^[4]提出采用半椭圆形的束窗(长轴位于径向位置)能够减少冷却剂(通常是液态的铅或铅铋合金^[2])流动分离，提高换热效率；秦雪^[5]进一步比较了长轴置于径向位置的半椭圆、半圆及长轴置于轴向位置的半椭圆这3种束窗形状，发现流动滞止区依次增大，换热效率依次降低，据此认为冷却剂通过束窗的平均速度是影响换热效率的关键。对流换热的Newton冷却公式为

$q = h\Delta T.$

(1)

其中：q为热流密度，h为对流换热系数，ΔT为温差。在温差一定时，可以通过提高对流换热系数来提高热流密度。对流换热系数与流速、流体物性、流动起因、流动形态、相变等诸多因素有关^[6]，在控制其他因素不变的情况下，对流换热系数与流速正相关。因此，本文提出一种改进的束窗设计：将长轴置于径向位置的椭圆更换为最速降线，通过最小化相同流量冷却剂通过束窗的时间使平均速度最大化，从而提高换热效率。
最速降线问题可以表述为：给定起点和终点，寻找一条使得小球下落时间最短的曲线。1744年Euler建立了变分法，获得了最速降线问题的一般解法^[7]。随着研究的不断深入，最优控制理论和离散点近似法也被用于最速降线求解，如Golubev等^[8-11]基于最优控制理论的极大值原理分别获得了无阻力、含黏性力、含摩擦力的最速降线的解析式，Szarkowicz^[12]则采用多阶段的Monte Carlo法获得一系列离散点来近似经典最速降线。但是，极大值原理只是最优解的必要不充分条件，其实际求解效果往往不如变分法，而离散点近似法则计算量过大，且无法估计误差。因此，变分法仍是求解最速降线问题的首选方法。
近年来，变分法被****们广泛应用于各种存在阻力或约束的最速降线问题。史友进等^[13]研究了受摩擦阻力的最速降线问题；周嘉炜等^[14]研究了受Stokes黏性阻力的最速降线问题；Smith等^[15]研究了在地球内部受地心引力作用的最速降线问题；Singh等^[16]研究了在地球外部受地心引力作用的最速降线问题；Vratanar等^[17]研究了受非保守力作用的最速降线问题；Stork等^[18]研究了在重力作用下终点速度方向被约束的最速降线问题；Cherkasov等^[19]研究了固定运动时间和终点速度方向且存在高阶黏性阻力时的最速降线问题。
将最速降线运用于ADS靶件束窗形状设计时，由于冷却剂为液态金属，需要考虑黏性力的影响。采用变分法分析含黏性力的最速降线问题，最终必须求解一组非线性积分方程^[17]，若采用Newton迭代法求解，在高黏性系数情况下难以获取收敛解，且无法获取近似解。因此，本文将方程组求解问题转化为最优化问题，比较并选择合适的算法对目标函数进行最优化计算。数值试验结果表明，在高黏性系数情况下，最优化方法较Newton迭代法计算效率高，且能够获得一个较好的近似解，为最速降线应用于ADS靶件束窗形状设计提供了可能。
1 变分法分析变分法是一种用于求泛函极值的方法，在稳定性分析、结构分析等方面有广泛应用，下面简述利用变分法分析含黏性力最速降线问题的过程。
如图 2所示，小球同时受重力mg、支持力F_N、黏性力F_r作用。假设黏性力大小正比于速度大小，即F_r=kmv (k为等效阻力系数，单位为s^-1，对于Re＜0.5的低速运动，黏性力可用Stokes公式表示^[14]，此时k=6πηr/m，η为流体黏性系数，r为小球半径，m为小球质量)，由切向力平衡可得

$\frac{{{\rm{d}}v}}{{{\rm{d}}t}} - g\sin \varphi - kv = 0.$

(2)

图 2 最速降线问题的小球受力示意图

图选项

假设小球经过的弧长为s(φ)，可将式(1)改写为

$vv' - gs'\sin \varphi - kvs' = 0.$

(3)

其中(·)′=d(·)/dφ。式(3)为力学约束，另外存在以下几何约束：

$\left\{ \begin{array}{l}{\rm{d}}x = {\rm{d}}s\cos \varphi ,\\{\rm{d}}y = {\rm{d}}s\sin \varphi .\end{array} \right.$

(4)

式(4)两边均除以dφ后可化为

$\left\{ \begin{array}{l}x' - s'\cos \varphi = 0,\\y' - s'\sin \varphi = 0.\end{array} \right.$

(5)

以(·)₀、(·)₁分别代表起点和终点处的变量值，小球下滑的总时间为

${t_{\rm{s}}} = \int_0^{{s_1}} {\frac{{{\rm{d}}s}}{v}} = \int_{{\varphi _0}}^{{\varphi _1}} {\frac{{s'}}{v}{\rm{d}}\varphi } .$

(6)

引入Lagrange乘子κ(φ)、λ(φ)、μ(φ)，问题转化为求式(7)泛函的最小值，

$F = \int_{{\varphi _0}}^{{\varphi _1}} {G\left( {\varphi ,\mathit{\boldsymbol{q}},\mathit{\boldsymbol{q'}}} \right){\rm{d}}\varphi } .$

(7)

其中：变量q =(v，s，x，y，κ，λ，μ)，核函数G(φ，q，q′)为

$\begin{array}{*{20}{c}}{G\left( {s,\mathit{\boldsymbol{q}},\mathit{\boldsymbol{q'}}} \right) = \frac{{s'}}{v} + \kappa \left( {x' - s'\cos \varphi } \right) + }\\{\lambda \left( {y' - s'\sin \varphi } \right) + \mu \left( {vv' - gs'\sin \varphi - kvs'} \right).}\end{array}$

(8)

式(7)对应的Euler方程^[20]为

$\frac{{\partial G}}{{\partial \mathit{\boldsymbol{q}}}} - \frac{{\rm{d}}}{{{\rm{d}}\varphi }}\left( {\frac{{\partial G}}{{\partial \mathit{\boldsymbol{q'}}}}} \right) = 0.$

(9)

式(9)对应的自然边界条件^[20]为

$\begin{array}{l} - \frac{{\partial G}}{{\partial \mathit{\boldsymbol{q'}}}}\left| {_{\varphi = {\varphi _0}}} \right.{\rm{d}}{\mathit{\boldsymbol{q}}_0} + \frac{{\partial G}}{{\partial \mathit{\boldsymbol{q'}}}}\left| {_{\varphi = {\varphi _1}}} \right.{\rm{d}}{\mathit{\boldsymbol{q}}_1} + \\\;\;\;\left( {G - \mathit{\boldsymbol{q'}}\frac{{\partial G}}{{\partial \mathit{\boldsymbol{q'}}}}} \right)\left| {_{\varphi = {\varphi _1}}} \right.{\rm{d}}{\varphi _1} = 0.\end{array}$

(10)

利用初始条件x₀、y₀、x₁、y₁、v₀和式(10)可求解式(9)，得x(φ)、y(φ)、s(φ)的表达式，其中所含的未知量φ₁、v₁则可通过求解式(11)积分方程组获得，

$\left\{ \begin{array}{l}M\left( {{\varphi _1},{v_1}} \right) = \int_{{\varphi _0}}^{{\varphi _1}} {s'\cos \varphi \ {\rm{d}}\varphi } - {x_1} + {x_0} = 0,\\N\left( {{\varphi _1},{v_1}} \right) = \int_{{\varphi _0}}^{{\varphi _1}} {s'\sin \varphi \ {\rm{d}}\varphi } - {y_1} + {y_0} = 0.\end{array} \right.$

(11)

2 最优化数值求解方法文[17]采用Newton迭代法对方程组(11)进行数值求解，需要选取迭代初值(φ₁⁽⁰⁾，v₁⁽⁰⁾)。若初值选取合适，迭代可获得任意精度的数值解；若初值选取不合适，迭代会直接发散。实际工程中发散是需要严格避免的，因为此时的计算没有任何价值，尤其是在高黏性系数情况下，由于很难选取到可收敛的迭代初值，会产生大量的无效计算，导致时间成本不可控，所以Newton迭代法不适合在实际工程中使用。本文将方程组(11)的求解转化为一个最优化问题，使用启发式搜索算法来获取近似解，从而避免了无效计算问题。
任取初值点(φ₁^*，v₁^*)，均可获得一条最速降线，其终点(x₁^*，y₁^*)由式(12)确定，

$\left\{ \begin{array}{l}x_1^ * = \int_{{\varphi _0}}^{\varphi _1^ * } {s'\cos \varphi \ {\rm{d}}\varphi } + {x_0},\\y_1^ * = \int_{{\varphi _0}}^{\varphi _1^ * } {s'\sin \varphi \ {\rm{d}}\varphi } + {y_0}.\end{array} \right.$

(12)

计算当前终点与目标终点间的Euclid距离，

$d = \sqrt {{{\left( {{x_1} - x_1^ * } \right)}^2} + {{\left( {{y_1} - y_1^ * } \right)}^2}} .$

(13)

现只需最小化d，采用启发式最优化算法不断更新(φ₁^*，v₁^*)，经过有限步后可获得满足精度要求的近似解。
本文考察了4种启发式最优化算法的计算效果。其中：a代表(φ₁，v₁)，ε代表收敛误差。
2.1 模拟退火算法模拟退火算法^[21-22]将最优化过程类比为固体物质的退火过程，通过Metropolis准则概率接受新解，避免陷入局部最优解，通过降温过程保证收敛性。计算步骤如下：
步骤1??随机产生初始解a₀，选取初温T₀。
步骤2??按照式(12)和(13)计算d，若d＜ε，则停止迭代输出结果，否则转步骤3。
步骤3??在当前解a_i邻域内随机产生新解a_j并计算d_j，若d_j＜d_i，则将a_j作为新的当前解，否则继续判断exp[-(d_j-d_i)/T_m]＞random(0，1)成立与否，若成立，则接受a_j作为新的当前解。
步骤4??判断内循环收敛准则是否满足，若满足则转步骤5，否则转步骤3。
步骤5??温度下降至T_m+1，转步骤2。
2.2 禁忌搜索算法禁忌搜索算法^[22]是局部邻域搜索的一种扩展，通过引入禁忌准则来避免陷入局部最优解，并通过藐视准则来赦免一些被禁忌的优良解，保证收敛性。计算步骤如下：
步骤1??随机产生初始解a₀，置空禁忌表。
步骤2??按照式(12)和(13)计算d，若d＜ε，则停止迭代输出结果，否则转步骤3。
步骤3??在当前解a_i邻域内选取若干候选解a_j1，…，a_jn，计算d_j1，…，d_jn。
步骤4??若存在d_jm小于历史最小值d_best，则令a_i=a_jm，转步骤2，否则转步骤5。
步骤5??选择非禁忌域内的最优候选解作为当前解。
步骤6??若禁忌解数量大于禁忌长度l，则用当前解a_i替换最早加入禁忌表的解，否则直接将当前解a_i加入禁忌表，转步骤2。
2.3 遗传算法遗传算法^[21-22]将最优化过程类比为基因的适者生存过程，通过选择、交叉过程保证收敛性，通过变异过程避免陷入局部最优解。计算步骤如下：
步骤1??随机产生一个初始种群a₁，…，a_n，按照式(12)和(13)计算其适应值d₁，…，d_n。
步骤2??若存在d_j＜ε，则停止迭代输出结果，否则转步骤3。
步骤3??依概率${{p}_{j}}={\left( 1/{{d}_{j}} \right)}/{\sum\limits_{j=1}^{n}{\left( 1/{{d}_{j}} \right)}}\;$，j=1, 2, …, n, 选取n对个体。
步骤4??按交叉概率p_c执行交叉操作。
步骤5??按变异概率p_v执行变异操作。
步骤6??计算新种群的个体适应值d_j，j=1, 2, …, n, 转步骤2。
2.4 粒子群算法粒子群算法^[23]是一种群智能算法，通过群体间的信息共享实现最优化过程。每一次迭代粒子均向自身最优解和群体最优解的加权中心运动，因此在保证收敛性的同时也具有一定的跳出局部最优解的能力。粒子的速度和位置按式(14)更新，

$\left\{ \begin{array}{l}{\mathit{\boldsymbol{v}}^{\left( {m + 1} \right)}} = {c_0}{\mathit{\boldsymbol{v}}^{\left( m \right)}} + {c_1}{r_1}\left( {\mathit{\boldsymbol{a}}_{\rm{p}}^{\left( m \right)} - {\mathit{\boldsymbol{a}}^{\left( m \right)}}} \right) + {c_2}{r_2}\left( {\mathit{\boldsymbol{a}}_{\rm{g}}^{\left( m \right)} - {\mathit{\boldsymbol{a}}^{\left( m \right)}}} \right),\\{\mathit{\boldsymbol{a}}^{\left( {m + 1} \right)}} = {\mathit{\boldsymbol{a}}^{\left( m \right)}} + {\mathit{\boldsymbol{v}}^{\left( {m + 1} \right)}}.\end{array} \right.$

(14)

其中：c₀、c₁、c₂为学习因子，r₁、r₂为(0，1)间的随机数，a_p、a_g分别为粒子自身最优解和群体最优解。计算步骤如下：
步骤1??随机产生一个初始种群a₁，…，a_n和对应的初始速度v₁，…，v_n，按照式(12)和(13)计算d₁，…，d_n。
步骤2??若存在d_j＜ε，则停止迭代输出结果，否则转步骤3。
步骤3??更新粒子自身最优解和群体最优解。
步骤4??按照式(14)更新粒子的位置与速度，并计算d_j，j=1, 2, …, n, 转步骤2。
3 数值试验3.1 4种最优化算法计算效率对比本文以一个具体的高黏性力最速降线问题为例比较4种最优化算法的计算效率，问题具体参数如下：起点x₀=0 m、y₀=3.5 m，终点x₁=6 m、y₁=0.5 m，小球初速度v₀=0 m·s^-1，重力加速度g=9.8 m·s^-2，对于高黏性阻力情况k=1.5 s^-1，计算最速降线。
设置搜索区域^[17]A={a|φ₁∈(-π/2，π/2)，v₁∈(0，min(g(1+sinφ₁)/2k，$\sqrt{2g\left( {{y}_{0}}-{{y}_{1}} \right)}$))}，收敛误差ε=10^-6 m。试验在一台CPU主频为2.3 GHz的计算机上进行，考虑到最优化算法的随机性和计算能力的限制，重复进行100次试验，每次试验最多迭代计算10 000次，若无法收敛则输出误差d最小的近似解。最优化方法的参数设置如下：
模拟退火算法：初温T₀=1，当前解a_i的邻域S={a_j||φ_1j-φ_1i|＜0.1π，|v_1j-v₁|＜0.1·min(g(1+sinφ₁)/2k，$\sqrt{2g\left( {{y}_{0}}-{{y}_{1}} \right)}$)}，内循环收敛准则为当前解连续10次保持不变，退温函数T_m=T₀/ln(m+2)。
禁忌搜索算法：候选解个数n=8，邻域同上，禁忌域为P={a_j||φ_1j-φ_1i|＜0.1π/3，|v_1j-v₁|＜0.1·min(g(1+sinφ₁)/2k，$\sqrt{2g\left( {{y}_{0}}-{{y}_{1}} \right)}$)/3}，禁忌长度l=3。
遗传算法：种群个体数n=10，基因为十进制个位至十万分位的数位值，共6位，交叉概率p_c=1，交叉操作为基因取平均数，变异概率p_v=0.1，变异操作为基因取随机数。
粒子群算法：粒子数n=10，学习因子c₀=0.68，c₁=c₂=2。
试验结果列于表 1。由表 1可知，粒子群算法远优于其他3种算法，只有粒子群算法能获得收敛解，这是因为随着搜索的进行，Δd=d^(m+1)-d^(m)逐渐减小，相应地，Δa=a^(m+1)-a^(m)也逐渐减小，由式(14)可知，v也逐渐减小，相当于粒子群算法能自适应地调节搜索邻域的大小，而其他3种算法均没有此特性，因此粒子群算法是本问题最适合的算法。
表 1 4种最优化算法计算效率对比

最优化算法	收敛概率	近似解平均误差/m
模拟退火算法	0	0.448
禁忌搜索算法	0	0.602
遗传算法	0	0.835
粒子群算法	0.7	0.797

表选项






3.2 粒子群算法与Newton迭代法计算效率对比将粒子群算法与Newton迭代法的计算效率进行对比。粒子群算法参数设置同3.1节；对于Newton迭代法，初值不收敛时记录其对应的误差d，并从中选出最优近似解，此时等效于使用Monte Carlo法进行最优化计算，从而能与粒子群算法进行近似解误差对比。共对比了k取0.25、0.50、0.75、1.00、1.25、1.50、1.75、2.00 s^-1这8种情况，计算结果如图 3所示。

图 3 粒子群算法与Newton迭代法计算效率对比

图选项

由图 3可知，在k=1.5 s^-1时，粒子群算法计算效率高于Newton迭代法，由于计算次数的限制，在k=1.75 s^-1和k=2 s^-1时，两者均无法获得收敛解，但粒子群算法获得的近似解误差较小，因此可以认为在高黏性系数情况下，粒子群算法优于Newton迭代法。但是，在低黏性系数情况下，粒子群算法的计算效率不如Newton迭代法，这是因为粒子群算法的收敛速度慢于Newton迭代法。
3.3 粒子群算法收敛速度分析任取一次收敛试验，作误差d关于计算次数N的图并进行最小二乘法线性拟合，如图 4所示。

图 4 粒子群算法收敛速度分析

图选项

收敛速度判定准则(C为常数，K)：
次线性收敛：N≥C/d^K，即lgN~O(lg(1/d))；
线性收敛：N≥C lg(1/d)，即N~O(lg(1/d))；
二阶收敛：N≥C lg lg(1/d)，即N~O(lg lg(1/d))。
由图 4可知，拟合相关系数R²最大的是线性收敛拟合，因此粒子群算法在此问题上是近似线性收敛的，弱于二阶收敛的Newton迭代法。
3.4 最速降线计算结果k取0.25、0.50、0.75、1.00、1.25、1.50 s^-1时的6条最速降线计算结果如图 5所示。

图 5 不同等效阻力系数下含黏性力最速降线计算结果

图选项

由图 5可知，随着k的增大最速降线逐渐趋近于一条直线。这是由于平均时间由运动距离和平均速度的比值决定，在k较小时，驱动力(重力)占主导地位，粒子容易加速，从而使最速降线在前期表现出更大的斜率；在k较大时，黏性力占主导地位，粒子难以加速，从而使最速降线趋于运动距离较小的直线。
采用图 5中k=0.25 s^-1情况下的最速降线设计靶件束窗(径向60 mm，轴向30 mm)，并与等长的半椭圆束窗进行数值模拟对比，参数设置同文[1]，即工质为20 ℃的水，向下流动，入口流速为1 m·s^-1，入口宽10 mm，出口宽20 mm，出入口压强均为0.1 MPa，边界绝热。二者的速度场和局部流场如图 6所示，图中L为流动分离点距出口中心水平距离，S为流道滞止区(流速小于0.145 m·s^-1)面积。

图 6 (网络版彩图)最速降线束窗与椭圆束窗的速度场及局部流场对比

图选项

由图 6可知，最速降线束窗相较于半椭圆束窗，推迟了流动分离，减小了滞止区面积，有利于提升换热效率。本文作者推测这是由于最速降线存在一定的凹陷，减弱了回流对主流的削弱作用所致。
4 结论采用变分法分析含黏性力的最速降线问题，最终可获得一组非线性积分方程。若使用Newton迭代法数值求解，在高黏性系数情况下很难获得收敛解，因而存在大量无效计算。本文将方程组的求解问题转化为一个最优化问题，在无法获得收敛解时可确保获得一个较好的近似解，从而避免了无效计算。通过比较模拟退火算法、禁忌搜索算法、遗传算法、粒子群算法这4种启发式最优化算法，最终选择近似线性收敛的粒子群算法进行计算，在高黏性系数情况下其计算效率优于Newton迭代法。数值模拟结果表明：获得的最速降线用于ADS散裂靶靶件束窗设计时，相较于半椭圆束窗，推迟了流动分离，减小了滞止区面积，有利于提升换热效率。

参考文献

[1]	陈海燕, 徐长江. ADS靶区流场的数值模拟[J]. 核科学与工程, 2002, 22(4): 361–364. CHEN H Y, XU C J. Numerical simulations of flow field in the target region of accelerator-driven subcritical reactor system[J]. Chinese Journal of Nuclear Science and Engineering, 2002, 22(4): 361–364. (in Chinese)
[2]	徐敬尧. 先进核反应堆用铅铋合金性能及纯净化技术研究[D]. 合肥: 中国科学技术大学, 2013. XU J Y. Study on performance and purification technology of lead-bismuth alloy for advanced nuclear reactor[D]. Hefei: University of Science and Technology of China, 2013. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10358-1014121213.htm
[3]	SUGAWARA T, NISHIHARA K, OBAYASHI H, et al. Conceptual design study of beam window for accelerator-driven system[J]. Journal of Nuclear Science and Technology, 2010, 47(10): 953–962. DOI:10.1080/18811248.2010.9720974
[4]	徐长江, 陈海燕. 基于流场数值模拟的ADS靶件结构的设计优化[J]. 原子能科学技术, 2003, 37(2): 111–115. XU C J, CHEN H Y. Design optimization of ADS target based on numerical simulations of flow field[J]. Atomic Energy Science and Technology, 2003, 37(2): 111–115. (in Chinese)
[5]	秦雪. 加速器驱动次临界系统中液态有窗和无窗靶的中子物理及热工水力特性[D]. 重庆: 重庆大学, 2015. QIN X. Neutron physics and thermal-hydraulic characteristics of liquid window and windowless target in accelerator driven sub-critical systems[D]. Chongqing: Chongqing University, 2015. (in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10611-1015969053.htm
[6]	章熙民, 任泽霈, 梅飞鸣. 传热学[M]. 5版, 北京: 中国建筑工业出版社, 2007. ZHANG X M, REN Z P, MEI F M. Principles of heat transfer[M]. 5th ed, Beijing: China Architecture & Building Press, 2007. (in Chinese)
[7]	吴佩萱. "最速降线"问题——数学史上最激动人心的一次公开挑战[J]. 数学通报, 2006, 45(5): 42–44. WU P X. Brachistochrone:The most stirring challenge in mathematics history[J]. Bulletin Des Sciences Mathematics, 2006, 45(5): 42–44. (in Chinese)
[8]	GOLUBEV Y F. Brachistochrone with friction[J]. Journal of Computer and Systems Sciences International, 2010, 49(5): 719–730. DOI:10.1134/S1064230710050060
[9]	GOLUBEV Y F. Brachistochrone with dry and arbitrary viscous friction[J]. Journal of Computer and Systems Sciences International, 2012, 51(1): 22–37. DOI:10.1134/S1064230712010078
[10]	VONDRUKHOV A S, GOLUBEV Y F. Brachistochrone with an accelerating force[J]. Journal of Computer and Systems Sciences International, 2014, 53(6): 824–838. DOI:10.1134/S1064230714060124
[11]	VONDRUKHOV A S, GOLUBEV Y F. Optimal trajectories in the brachistochrone problem with an accelerating force[J]. Journal of Computer and Systems Sciences International, 2015, 54(4): 514–524. DOI:10.1134/S1064230715040139
[12]	SZARKOWICZ D S. Investigating the brachistochrone with a multistage Monte Carlo method[J]. International Journal of Systems Science, 1995, 26(2): 233–243. DOI:10.1080/00207729508929033
[13]	史友进, 俞晓明. 库仑摩擦最速降曲线问题的讨论[J]. 盐城工学院学报(自然科学版), 2012, 25(2): 1–4. SHI Y J, YU X M. On Brachistochrone with Coulomb friction[J]. Journal of Yancheng Institute of Technology (Natural Science Edition), 2012, 25(2): 1–4. (in Chinese)
[14]	周嘉炜, 曹炳阳, 过增元. 具有黏滞阻力时的最速降线[J]. 力学与实践, 2011, 33(3): 11–14. ZHOU J W, CAO B Y, GUO Z Y. Brachistochrone with viscous resistance[J]. Mechanics in Engineering, 2011, 33(3): 11–14. (in Chinese)
[15]	SMITH D R. Variational methods in optimization[M]. New York: Dover Publications, 1974.
[16]	SINGH B, KUMAR R. Brachistochrone problem in nonuniform gravity[J]. Indian Journal of Pure and Applied Mathematics, 1988, 19(6): 575–585.
[17]	VRATANAR B, SAJE M. On the analytical solution of the brachistochrone problem in a non-conservative field[J]. International Journal of Non-Linear Mechanics, 1998, 33(3): 489–505. DOI:10.1016/S0020-7462(97)00026-7
[18]	STORK D G, YANG J X. The general unrestrained brachistochrone[J]. American Journal of Physics, 1988, 56(1): 22–26. DOI:10.1119/1.15423
[19]	CHERKASOV O Y, ZARODNYUK A V. Brachistochrone problem with linear and quadratic drag[J]. AIP Conference Proceedings, 2014, 1637: 195–200. DOI:10.1063/1.4904579
[20]	老大中. 变分法基础[M]. 2版, 北京: 国防工业出版社, 2007. LAO D Z. Fundamentals of the calculus of variations[M]. 2nd ed, Beijing: National Defense Industry Press, 2007. (in Chinese)
[21]	蒋金山, 何春雄, 潘少华. 最优化计算方法[M]. 广州: 华南理工大学出版社, 2007. JIANG J S, HE C X, PAN S H. Optimization computational method[M]. Guangzhou: South China University of Technology Press, 2007. (in Chinese)
[22]	王凌. 智能优化算法及其应用[M]. 北京: 清华大学出版社, 2001. WANG L. Intelligent optimization algorithms with applications[M]. Beijing: Tsinghua University Press, 2001. (in Chinese)
[23]	高尚, 杨静宇. 群智能算法及其应用[M]. 北京: 中国水利水电出版社, 2006. GAO S, YANG J Y. Swarm intelligence algorithms and applications[M]. Beijing: China Water & Power Press, 2006. (in Chinese)