面向货车编队的水泥混凝土路面结构拓扑优化设计及验证

摘要

货车编队的渠化交通荷载对路面结构提出了高承载力及精细化设计需求。结合拓扑优化理论，基于变密度法、拓扑梯度法与移动渐近线法，提出了面向水泥混凝土路面的结构拓扑优化算法；分析货车编队轮迹偏移特征，利用数值计算设计了面向货车编队的水泥混凝土路面拓扑优化原型结构；从力学性能角度通过室内缩尺试验探究了拓扑优化原型结构的合理性。结果表明，对于轮迹零偏移和15 cm偏移范围的两种荷载分布，拓扑优化原型结构最大疲劳应力分别为传统结构的45.0%和37.5%，同时其质量分别为传统结构的34.3%和30.5%；试验测得的拓扑优化原型结构的最大疲劳应力与理论计算值相近，且相较于传统结构，其极限承载力提升了18.3%，证明了设计方法的有效性。

关键词

道路工程; 货车编队; 路面结构; 水泥混凝土路面; 拓扑优化

货车编队是指多辆货车通过车联网和自动驾驶技术以较小车距和稳定轨迹运行的编队行驶模式，该行驶模式可提高运输效率、扩大路网容量、实现节能减排（车辆节省燃油约8%~13%^［1］），已成为各国落实网联汽车发展的普遍战略。相较于传统货车，编队运行货车具有高精度的定位和操控系统，其轮迹横向分布较为集中，交通渠化现象突出。研究表明，自动驾驶等级不同，轮迹横向分布标准差可相差1倍以上^［2］。在轮迹横向“零”偏移时，即极度渠化的交通荷载下，路面结构累计疲劳损伤会提高近2倍^［3］。因此，货车编队荷载的高频和渠化作用对路面承载能力和寿命提出了更高要求。同时，渠化的交通荷载造成的路面结构局部受荷现象突出，使路面结构不同部位力学特性要求的差异性更为显著，因此对路面结构进行非均匀性的精细化设计变得尤为重要。

路面结构的非均匀性精细化设计也为路面智能化与功能化发展提供了更多的可能性。根据智能铺面的概念与技术发展路径^［4-5］，路面在实现智能化与功能化过程中，如传感器埋设、排水路面以及管线埋设等，均对路面结构提出了埋设空间的需求，而非均匀性精细化设计可以在保证路面结构承载能力与使用寿命的前提下拓展更多的结构内部空间。同时，对于以装配式技术^［6］为主的智能铺面建造方式，非均匀性精细化设计也可以使路面结构实现轻量化，有效提高建造效率，降低建造成本。

目前，传统路面结构设计方法均采用尺寸优化与形状优化方法，无法适应货车编队下路面结构的非均匀性精细化设计需求。而拓扑优化方法作为结构优化方法之一，可以在设定的区域内，根据给定的荷载、边界条件与性能指标，对区域内材料分布进行优化，在路面结构精细化设计方面具有广阔的应用前景。

拓扑优化理论最早可追溯至1904年Michell研究的离散结构拓扑优化中典型的Michell桁架理论^［7］。随着拓扑优化方法的发展，根据研究对象不同，拓扑优化可分为离散结构拓扑优化与连续体结构拓扑优化，目前结构拓扑优化主要研究对象多为连续体结构。自从1988年Bendsoe与Kikuchi提出了连续体结构拓扑优化的均匀化方法^［8］，拓扑优化理论进入蓬勃发展时期。在均匀化方法的基础上，变密度法（SIMP）于1989年由Bendsoe^［9］提出，该方法的核心是引入一种假想的相对密度在0~1之间可变的材料，将结构优化问题转变为材料本构模型的插值问题，通过改变微观尺度上的材料本构关系来实现宏观尺度的结构优化，是目前使用最为广泛的拓扑优化方法之一^［10-13］。除此之外，拓扑优化的主流方法还包括由Sethian等提出的水平集方法^［14-16］，以及由谢亿民和Steven提出的渐进结构拓扑优化法^［17-18］与双向渐进拓扑优化方法^［19-20］等。目前，随着拓扑优化理论的成熟，该方法已被广泛应用于结构设计，尤其是航空器^［21-22］、车辆的结构轻量化设计^［23-24］，以及建筑结构轻量化设计^［25-26］。

本文面向货车编队，以水泥混凝土路面结构为研究对象，基于拓扑优化理论和方法，构建了水泥混凝土路面结构拓扑优化模型并提出了优化算法；基于此，在分析货车编队轮迹横向分布特征的基础上，利用数值计算方法设计了面向货车编队的水泥混凝土路面拓扑优化原型结构；最后通过3D打印方法制作了原型结构的缩尺模型，开展力学性能试验以分析拓扑优化设计方法的合理性。

1 水泥混凝土路面结构拓扑优化模型和算法

1.1　拓扑优化问题模型

水泥混凝土路面结构拓扑优化设计首先需保证优化结构具有足够的承载能力，即水泥混凝土路面拓扑优化结构应以设计基准期内行车荷载与温度梯度综合作用产生的面层板疲劳断裂作为设计标准。同时，考虑到水泥混凝土路面结构的功能化以及轻量化需求，在拓扑优化设计算法中，问题模型的建立应当以行车荷载与温度荷载作用下设计域内任一点的疲劳应力不超过水泥混凝土材料的极限拉应力作为约束条件，并以结构质量最小化为优化目标。将优化目标与约束条件建模即可得水泥混凝土路面结构拓扑优化设计的问题模型，如式（1）所示。

$$\begin{array}{l}\mathrm{m}\mathrm{i}\mathrm{n}m\left(\rho \right)={\int }_{\Omega }\rho d\Omega \\ \mathrm{s}.\mathrm{t}\left\{\begin{array}{l}a\left(\mathit{u},\mathit{v}\right)=l\left(\mathit{v}\right)\\ {\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\le f_{\mathrm{r}}\end{array}\right.\end{array}$$

（1）

式中：ρ为材料密度，是表征结构内材料分布状态的变量；Ω为表征结构区域的有界开集；$a\left(\mathit{u},\mathit{v}\right)=l\left(\mathit{v}\right)$为结构内外力平衡方程的弱形式，u、v分别为实际位移场与虚拟位移场；σ_max为外部荷载与温度梯度综合作用下产生的最大疲劳应力；f_r为水泥混凝土材料的抗拉强度。

在上述模型基础上，本文选用变密度法（solid isotropic microstructure with penalization，SIMP）作为结构内材料分布的表述方法，对上述问题模型进行改进。变密度法的原理为通过离散初始设计域，并以离散后每个单元内材料密度作为设计变量来表征结构中材料分布状态，并通过改变不同位置的单元内材料密度表征结构变化，如图1所示。

图 1  SIMP方法

Fig. 1  SIMP method

由于SIMP方法会导致大量的灰度材料出现，即单元内材料密度介于0与1之间，不利于工程建造。为避免优化结果出现此类情况，将拓扑优化问题模型（式（1））中的密度ρ替代为体积插值函数$m_{\mathrm{V}}\left(\mathit{y}\right)$，可用式（2）表示。

$$\begin{array}{l}\underset{z\in \left[\mathrm{0,1}\right]}{\mathrm{m}\mathrm{i}\mathrm{n}}f\left(\mathit{z}\right)=\frac{{\mathit{A}}^{\mathrm{T}}{m}_{\mathrm{V}}\left(\mathit{y}\right)}{{\mathit{A}}^{\mathrm{T}}I}\\ \mathrm{s}.\mathrm{t}\left\{\begin{array}{l}a\left(\mathit{u},\mathit{v}\right)=l\left(\mathit{v}\right)\\ g_j\left(\mathit{z},\mathit{u}\right)=m_{\mathrm{E}}\left(y_j\right){\Lambda }_j\left({\Lambda }_j^{\mathrm{2}}+\mathrm{1}\right)\le \mathrm{0}\\ j=\mathrm{1},\cdots ,N\end{array}\right.\end{array}$$

（2）

其中

$$m_{\mathrm{V}}\left(y_j\right)=\frac{\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\beta \eta \right)+\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\beta \left(y_j-\eta \right)\right)}{\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\beta \eta \right)+\mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\beta \left(\mathrm{1}-\eta \right)\right)}$$

（3）

$$m_{\mathrm{E}}\left(y_j\right)=\epsilon +\left(\mathrm{1}-\epsilon \right){\left[m_{\mathrm{V}}\left(y_j\right)\right]}^p$$

（4）

$${\Lambda }_j=\frac{{\sigma }_j^{\mathrm{m}\mathrm{a}\mathrm{x}}}{f_{\mathrm{r}}}-\mathrm{1}$$

（5）

式（2）~（5）中：z为结构设计变量，其值决定了单元密度的大小；${\mathit{A}}^{\mathrm{T}}$为设计域内各单位面积的向量；$I$为与其维度大小相同的单位向量；单元密度矩阵$\mathit{y}=\mathit{P}\mathit{z}$，P为灵敏度过滤矩阵，作用是减缓结构边界出现的锯齿化现象；体积插值函数$m_{\mathrm{V}}\left(\mathit{y}\right)$中，β参数控制函数的平滑程度，而η参数代表大于该阈值的数将被投射为1；$m_{\mathrm{E}}\left(y_j\right)$为材料插值函数，主要作用为关联单元材料密度$m_{\mathrm{V}}\left(y_j\right)$与模量ε；y_j为单元密度矩阵y的子项，j为单元编号；$g_j\left(\mathit{z},\mathit{u}\right)$为式（1）中应力约束${\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}\le f_{\mathrm{r}}$的变体，可以对超出最大应力限制的局部应力进行有效约束；${\Lambda }_j$为单元j的应力相对于约束值的偏离程度；${\sigma }_j^{\mathrm{m}\mathrm{a}\mathrm{x}}$为单元j的最大疲劳应力。

为寻找优化问题的最优解，利用增广拉格朗日法（augmented Lagrange method，ALM）将有约束优化问题模型（式（2））转变为无约束优化问题模型。

$$\mathrm{m}\mathrm{i}\mathrm{n}J\left(\mathit{z},\mathit{u}\right)=f\left(\mathit{z}\right)+\frac{\mathrm{1}}{N}p\left(\mathit{z},\mathit{u}\right)$$

（6）

其中

$$p\left(\mathit{z},\mathit{u}\right)={\displaystyle \sum _{j=\mathrm{1}}^N}\left[{\lambda }_j{h}_j\left(\mathit{z},\mathit{u}\right)+\frac{\mu }{\mathrm{2}}{h}_j{\left(\mathit{z},\mathit{u}\right)}^{\mathrm{2}}\right]$$

（7）

$$h_j\left(\mathit{z},\mathit{u}\right)=\mathrm{m}\mathrm{a}\mathrm{x}\left[g_j\left(\mathit{z},\mathit{u}\right),-\frac{{\lambda }_j}{\mu }\right]$$

（8）

式（6）~（8）中：N为单元总数；$p\left(\mathit{z},\mathit{u}\right)$为惩罚项；λ_j为单元j的拉格朗日乘子；μ为罚函数。其中，拉格朗日乘子与罚函数随着计算步的迭代而更新，其更新准则为

$${\mu }^{k+\mathrm{1}}=\mathrm{m}\mathrm{i}\mathrm{n}\left[\alpha {\mu }^k,{\mu }_{\mathrm{m}\mathrm{a}\mathrm{x}}\right]$$

（9）

$${\lambda }_j^{k+\mathrm{1}}={\lambda }_j^k+{\mu }^k{h}_j\left(z^k,\mathit{u}\right)$$

（10）

式（9）、（10）中：k为计算步；α为大于1的更新系数；μ_max为罚函数上限值。

通过逐步求解式（6）即可获得设计变量z，从而得到单元密度ρ并形成最后的优化结果。

1.2　灵敏度求解

在求解式（6）的过程中，需要对优化模型的灵敏度进行分析。灵敏度是指设计变量的改变对目标方程产生的影响，是拓扑优化计算过程中指导设计变量改变方向的重要参数。本文选用拓扑梯度法获取灵敏度。将式（6）关于设计变量$\mathit{z}$求导即可得本文优化模型的灵敏度。

$$\frac{\mathrm{d}J}{\mathrm{d}\mathit{z}}=\frac{\partial \mathit{E}}{\partial \mathit{z}}\left(\frac{\partial f}{\partial \mathit{E}}+\frac{\mathrm{1}}{N}\frac{\partial p}{\partial \mathit{E}}\right)+\frac{\partial \mathit{V}}{\partial \mathit{z}}\left(\frac{\partial f}{\partial \mathit{V}}+\frac{\mathrm{1}}{N}\frac{\partial p}{\partial \mathit{V}}\right)$$

（11）

式中：$\mathit{E}=m_{\mathrm{E}}(\mathit{P}\mathit{z})$；$\mathit{V}=m_{\mathrm{V}}(\mathit{P}\mathit{z})$。

材料刚度与密度关于设计变量的灵敏度为

$$\frac{\partial \mathit{E}}{\partial \mathit{z}}={\mathit{P}}^{\mathrm{T}}{\mathit{J}}_{m_{\mathrm{E}}}\left(\mathit{P}\mathit{z}\right)$$

（12）

$$\frac{\partial \mathit{V}}{\partial \mathit{z}}={\mathit{P}}^{\mathrm{T}}{\mathit{J}}_{m_{\mathrm{V}}}\left(\mathit{P}\mathit{z}\right)$$

（13）

式（12）、（13）中：${\mathit{J}}_{m_{\mathrm{E}}}$与${\mathit{J}}_{m_{\mathrm{V}}}$分别为$m_{\mathrm{E}}\left(y_j\right)\mathrm{\text{'}}$与$m_{\mathrm{V}}\left(y_j\right)\mathrm{\text{'}}$的对角矩阵。

对于其他项，有

$$\frac{\partial f}{\partial \mathit{E}}=\mathrm{0}$$

（14）

$$\frac{\partial f}{\partial \mathit{V}}=\frac{\mathit{A}}{{\mathit{A}}^{\mathrm{T}}I}$$

（15）

$$\frac{\partial p}{\partial \mathit{V}}=\mathrm{0}$$

（16）

$$\frac{\partial p}{\partial \mathit{E}}={\displaystyle \sum _{j=\mathrm{1}}^N}\left[{\lambda }_j+\frac{\mu }{\mathrm{2}}{h}_j\left(\mathit{z},\mathit{u}\right)\right]\left[\frac{\partial h_j\left(\mathit{z},\mathit{u}\right)}{\partial \mathit{E}}+\frac{\partial h_j\left(\mathit{z},\mathit{u}\right)}{\partial \mathit{u}}\frac{\partial \mathit{u}}{\partial \mathit{E}}\right]$$

（17）

由于式（17）中$\partial \mathit{u}/\partial \mathit{E}$项在实际计算过程中对计算资源占用过高，故采用伴随方程法避免该项计算。将平衡方程与伴随变量$\mathit{\xi }$引入式（17），有

$$\frac{\partial p}{\partial \mathit{E}}={\displaystyle \sum _{j=\mathrm{1}}^N}\left[{\lambda }_j+\frac{\mu }{\mathrm{2}}{h}_j\left(\mathit{z},\mathit{u}\right)\right]\left[\frac{\partial h_j\left(\mathit{z},\mathit{u}\right)}{\partial \mathit{E}}+\frac{\partial h_j\left(\mathit{z},\mathit{u}\right)}{\partial \mathit{u}}\frac{\partial \mathit{u}}{\partial \mathit{E}}\right]+\frac{\partial a\left(\mathit{u},\mathit{\xi }\right)}{\partial \mathit{u}}\frac{\partial \mathit{u}}{\partial \mathit{E}}-\frac{\partial l\left(\mathit{\xi }\right)}{\partial \mathit{u}}\frac{\partial \mathit{u}}{\partial \mathit{E}}$$

（18）

令包含$\partial \mathit{u}/\partial \mathit{E}$的项之和为零，利用有限元方法即可求解伴随变量$\mathit{\xi }$为

$$\mathit{K}\mathit{\xi }=-{\displaystyle \sum _{j=\mathrm{1}}^N}\left[{\lambda }_j+\frac{\mu }{\mathrm{2}}{h}_j\left(\mathit{z},\mathit{u}\right)\right]\frac{\partial h_j\left(\mathit{z},\mathit{u}\right)}{\partial \mathit{u}}$$

（19）

式中：$\mathit{K}$为刚度矩阵。将伴随变量带入式（11）即可得本文拓扑优化问题的灵敏度为

$$\begin{array}{l}\frac{\mathrm{d}J}{\mathrm{d}\mathit{z}}=\frac{\partial \mathit{E}}{\partial \mathit{z}}\left(\frac{\partial f}{\partial \mathit{E}}+\frac{\mathrm{1}}{N}\frac{\partial p}{\partial \mathit{E}}+\frac{\partial a\left(\mathit{u},\mathit{\xi }\right)}{\partial \mathit{E}}\right)+\\ \frac{\partial \mathit{V}}{\partial \mathit{z}}\left(\frac{\partial f}{\partial \mathit{V}}+\frac{\mathrm{1}}{N}\frac{\partial p}{\partial \mathit{V}}\right)\end{array}$$

（20）

1.3　设计变量更新

在获取优化模型关于设计变量的灵敏度后，需要利用适合的更新算法对设计变量进行逐步迭代更新，实现结构的拓扑改变。本文选用移动渐近线法（method of moving asymptotes，MMA）对设计变量$\mathit{z}$进行更新，其形式可以表现为

$$\begin{array}{l}\underset{\mathit{z}\in {\left[\mathrm{0,1}\right]}^N}{\mathrm{m}\mathrm{i}\mathrm{n}}{\tilde{J}}_k\left(\mathit{z}\right)=J^k+{\displaystyle \sum _{\mathrm{𝓁}=\mathrm{1}}^N}\left(\frac{p_{\mathrm{𝓁}}^k}{U_{\mathrm{𝓁}}^k-z_{\mathrm{𝓁}}}+\frac{q_{\mathrm{𝓁}}^k}{z_{\mathrm{𝓁}}-L_{\mathrm{𝓁}}^k}\right)\\ \mathrm{s}.\mathrm{t}{\alpha }_{\mathrm{𝓁}}^k\le z_{\mathrm{𝓁}}\le {\beta }_{\mathrm{𝓁}}^k,\mathrm{𝓁}=\mathrm{1},\dots ,N\end{array}$$

（21）

其中

$$p_{\mathrm{𝓁}}^k={\left(U_{\mathrm{𝓁}}^k-z_{\mathrm{𝓁}}^k\right)}^{\mathrm{2}}\left[\mathrm{m}\mathrm{a}\mathrm{x}\left(\frac{\partial J}{\partial z_{\mathrm{𝓁}}},\mathrm{0}\right)+\tau \left|\frac{\partial J}{\partial z_{\mathrm{𝓁}}}\right|+\frac{\theta }{U_{\mathrm{𝓁}}^k-L_{\mathrm{𝓁}}^k}\right]$$

（22）

$$\begin{array}{l}{q}_{\mathrm{𝓁}}^k={\left(z_{\mathrm{𝓁}}^k-L_{\mathrm{𝓁}}^k\right)}^{\mathrm{2}}\left[-\mathrm{m}\mathrm{i}\mathrm{n}\left(\frac{\partial J}{\partial z_{\mathrm{𝓁}}},\mathrm{0}\right)+\tau \left|\frac{\partial J}{\partial z_{\mathrm{𝓁}}}\right|+\frac{\theta }{U_{\mathrm{𝓁}}^k-L_{\mathrm{𝓁}}^k}\right]\\ \end{array}$$

（23）

$$L_{\mathrm{𝓁}}^k=z_{\mathrm{𝓁}}^k-s_{\mathrm{𝓁}}^k\left(z_{\mathrm{𝓁}}^{k-\mathrm{1}}-L_{\mathrm{𝓁}}^{k-\mathrm{1}}\right)$$

（24）

$$U_{\mathrm{𝓁}}^k=z_{\mathrm{𝓁}}^k+s_{\mathrm{𝓁}}^k\left(U_{\mathrm{𝓁}}^{k-\mathrm{1}}-z_{\mathrm{𝓁}}^{k-\mathrm{1}}\right)$$

（25）

$${\alpha }_{\mathrm{𝓁}}^k=\mathrm{m}\mathrm{a}\mathrm{x}\left[z_{\mathrm{𝓁}}^{\mathrm{m}\mathrm{i}\mathrm{n}},\mathrm{0.9}{L}_{\mathrm{𝓁}}^k+\mathrm{0.1}{z}_{\mathrm{𝓁}}^k\right]$$

（26）

$${\beta }_{\mathrm{𝓁}}^k=\mathrm{m}\mathrm{i}\mathrm{n}\left[z_{\mathrm{𝓁}}^{\mathrm{m}\mathrm{a}\mathrm{x}},\mathrm{0.9}{U}_{\mathrm{𝓁}}^k+\mathrm{0.1}{z}_{\mathrm{𝓁}}^k\right]$$

（27）

$$z_{\mathrm{𝓁}}^{\mathrm{m}\mathrm{i}\mathrm{n}}=\mathrm{m}\mathrm{a}\mathrm{x}\left[\mathrm{0},z^k-\nu \right]$$

（28）

$$z_{\mathrm{𝓁}}^{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{m}\mathrm{i}\mathrm{n}\left[\mathrm{1},z^k+\nu \right]$$

（29）

式（21）~（29）中：k为迭代步数；l为单元编号；$z_{\mathrm{𝓁}}$为单元的设计变量；$\partial J/\partial z_{\mathrm{𝓁}}$为式（20）计算得到的灵敏度；τ与θ为算法参数，默认分别取10^-3与10^-6；U与L分别对应MMA方法中设计变量的上界与下界；s为上下界移动步的步长；$\nu$为设计变量更新步长。

通过构造移动近似子问题（式（21））逐步逼近原问题（式（6）），并随着迭代步更新设计变量的上下界值，最终可得到式（6）的近似解，即最终的优化结果。

2 面向货车编队的水泥混凝土路面拓扑优化原型结构

货车编队荷载是水泥混凝土路面拓扑优化的重要参数之一，在水泥混凝土路面结构拓扑优化模型和算法的基础上，本文以EcoTwin货车编队测试项目的开源数据^［27］为基础，分析货车编队运行中的轮迹横向偏移特征，以此作为代表性荷载开展拓扑优化原型结构设计。

2.1　货车编队的轮迹偏移特征

EcoTwin项目起始于2018年，测试路段为位于荷兰海尔蒙德市（Helmond）至埃因霍温市（Eindhoven）之间的公路，货车编队数为1队，队列中车辆为2辆，运行过程中以居中式编队为主要队列形态，其测试数据包括每辆货车的位置坐标、速度、传感器信息、车辆交互信息等。图2和图3分别展示了2019年4月9日测试的编队货车速度状态和跟车状态，其中A270路段车辆间距全程保持在21~22 m，车头时距稳定保持在0.7~0.8 s。由此可知，货车编队在A270路段行驶状态较为稳定，因此本文选择该路段的行车数据作为分析对象。

图 2  编队货车速度状态

Fig. 2  Speed status of trucks in platoon

图 3  编队测试过程中的车辆间距与车头时距

Fig. 3  Location and time distance of trucks in test program

提取A270路段的64 000条车辆位置数据，与开源地图数据库OpenStreetMap（OSM）提供的道路数据进行比对，计算车辆中轴线相对于车道的横向偏移分布状态。车道宽度为3.75 m，图4为以1 cm为区间绘制所得的车辆中轴线对于车道的横向偏移分布状态，由图4可知，货车在编队运行过程中车辆中轴线相对于车道中线的偏移量约为15 cm。虽然其运行模式设定为严格沿车道中线行驶，但实际情况中3~10 cm以内的偏移出现概率较高。

图 4  编队测试中车辆中轴线横向偏移分布

Fig.4  Distribution of truck centerline in test program

根据国际自动机工程师学会（Society of Automotive Engineers， SAE）制定的5层级货车编队等级^［28］，EcoTwin中采用的货车编队等级仅为Level 2级，在行驶过程中会产生一定程度的轮迹偏移现象。随着控制器与自动驾驶技术的发展，未来在更高等级的编队中，轮迹偏移量会逐渐减小甚至达到零偏移状态^［3］。因此，本文选用面向未来货车编队的轮迹零偏移和实测数据得到的轮迹15 cm偏移两种场景，作为代表性荷载特征开展后续分析。

2.2　拓扑优化原型结构数值计算

首先，选取传统水泥混凝土路面结构的横截面为拓扑优化初始设计域，尺寸为3.5 m×0.3 m，边界条件为无传力杆构件的四边自由板。为改善刚性混凝土结构整体的抗开裂能力，同时也考虑到在路面板下方拓展空间的需求，本文基于平衡刚度设计理念，将路面结构底部支撑设置为刚性梁基础^［29］。为简化分析，本文选取单轴单轮作为分析荷载模型，轴载大小取为100 kN，轮胎宽度为20 cm，基于荷载零偏移和15 cm偏移两种场景，在轮迹分布范围内均设置均布荷载，分别如图5、图6所示。路面设计场景选取为位于公路自然区划Ⅱ区的三级公路，交通荷载等级为中等。

图 5  初始模型：零偏移的交通荷载

Fig.5  Initial model: drainage traffic loading

图 6  初始模型：15 cm偏移的交通荷载

Fig. 6  Initial model: traffic loading with an offset of 15 cm

在拓扑优化模型中，路面结构材料为水泥混凝土，弹性模量取3.15×10⁴ MPa，泊松比取0.15，线膨胀系数α_c=1×10^-5 ℃^-1，温度梯度取为88 ℃·m^-1；可靠度系数γ=1.04，荷载安全系数k_c=1.0，应力折减系数k_j=1.00，荷载应力疲劳系数由下式确定：

$$k_{\mathrm{f}}=\frac{\mathrm{1}}{a}{\left(N_{\mathrm{e}}\right)}^b$$

（30）

其中标准轴载累计作用次数N_e取值为$\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$，回归系数为a=1.0，b=0.057。温度疲劳应力系数计算公式为

$$k_{\mathrm{t}}=\frac{f_{\mathrm{r}}}{{\sigma }_{\mathrm{c},\mathrm{m}\mathrm{a}\mathrm{x}}}\left[a_{\mathrm{t}}{\left(\frac{{\sigma }_{\mathrm{c},\mathrm{m}\mathrm{a}\mathrm{x}}}{f}\right)}^{b_{\mathrm{t}}}-c_{\mathrm{t}}\right]$$

（31）

其中混凝土弯拉强度f_r=4.5 MPa，回归系数分别为a_t=0.828、b_t=0.041、c_t=1.323。

基于本文提出的拓扑优化算法，利用MATLAB软件进行有限单元法与拓扑优化的编程与计算。考虑到计算成本，本文的拓扑优化设计主要在二维空间展开。将初始设计区域划分为10 500个单元，即延宽度方向350个单元，厚度方向30个单元，单元尺寸为1 cm×1 cm，单元类型为双线性矩形单元。在计算过程中，单元应变与单元应力的计算点位选取为单元中心点处，中心点应力值通过对高斯积分点应力进行差值得到。为方便比较，采用单元应力值与设定的应力约束值之比作为结果。

对于零偏移的渠化交通荷载场景，路面结构拓扑优化数值计算的材料分布与应力云图如图7所示。在图7中，左图黑色部分的材料密度为1，是主要承力结构，白色区域为无材料区域，可自行根据功能和建造需求进行拓展。根据优化结果，应力约束在优化过程中有较为明显的振荡，但总体满足约束，收敛时结构最大疲劳应力为相同条件下初始结构的45.0%，约束值的68.0%，如图8所示。此外，优化后的结构材料用量为原结构的34.3%。而对于轮迹15 cm偏移的交通荷载场景，其材料分布状态与应力云图如图9所示。根据优化结果，该结构最大疲劳应力为相同荷载与边界条件下原结构的37.5%，约束值的75.0%，如图10所示。此外，优化后的结构材料用量为原结构的30.5%。

图 7  优化后材料分布与疲劳应力云图（轮迹零偏移交通荷载）

Fig. 7  Optimized material distribution and fatigue stress (drainage traffic loading)

图 8  最大疲劳应力水平迭代过程（轮迹零偏移交通荷载）

Fig. 8  Iteration courses of fatigue stress (drainage traffic loading)

图 9  优化后材料分布与疲劳应力云图（轮迹15 cm偏移交通荷载）

Fig. 9  Optimized material distribution and fatigue stress (traffic loading with an offset of 15 cm)

图 10  最大疲劳应力水平迭代过程（轮迹15 cm偏移的交通荷载）

Fig. 10  Iteration courses of fatigue stress (traffic loading with an offset of 15 cm)

3 拓扑优化原型结构性能验证

为探究拓扑优化设计方法的合理性，本文采用3D打印技术制作零偏移的交通荷载作用下拓扑优化路面结构的缩尺原型结构，并利用静载加载试验开展力学性能分析。

3.1　拓扑优化原型结构建造

利用有限单元法对初始结构进行了离散化计算，所以最终的优化结果为数值离散状态，无法直接依据计算结果进行建造。本文主要利用AutoCAD中的NURBS曲线对计算结果进行可行性优化。由于应力计算点位为单元中心点处，故以单元中心点为控制点构造NURBS曲线，对原离散结构连续化，使原型结构具有清晰明确的边界，如图11、12所示。

图 11  离散结构连续化

Fig. 11  Continuation of discrete structure

图 12  优化前后结构对比

Fig. 12  Comparison of structure before and after optimization

根据相似第二定理中的量纲分析方法，将原型结构（截面尺寸为3.5 m×0.3 m）等比例缩小为截面尺寸为23.3 cm×2.0 cm的板结构，如图13所示。由于目前水泥混凝土结构制造技术尚无法精确成型大尺寸异形板结构，本文选用尼龙材料替代水泥混凝土材料，使用3D打印技术制作缩尺结构，比较不同荷载等级下优化前后结构的最大疲劳应力值，以探究本文提出的拓扑优化设计方法的合理性。尼龙材料参数如表1所示，最终的缩尺结构如图14所示。

图 13  原型结构尺寸（单位：cm）

Fig. 13  Size of prototype structure(unit:cm)

表 1  尼龙材料主要参数

Tab. 1  Material properties of nylon

弹性模量/MPa	泊松比	屈服强度/MPa	抗拉强度/MPa	抗压强度/MPa
2650	0.2	36	44	20

图 14  缩尺试验试件

Fig. 14  Test specimen of scale model

3.2　试验方法

试验采用电液式压力试验机（最大加载压力为2 000 kN）进行加载，并利用DIC（digital image correlation）设备对加载过程中的形变进行图像捕捉，如图15所示，最后利用MATLAB软件对结构形变进行分析。加载试件由应力扩散板、原型结构以及底部支撑板组成，如图14所示。其中，应力扩散板用于将顶部加载头的圆形均布荷载转变为作用于主要受力结构上方的条形均布荷载，底部支撑板与缩尺结构利用结构胶黏接，作用为模拟边界条件。

图 15  单点加载试验设备布置

Fig. 15  Equipment of single point loading experiment

试验加载过程分为2个阶段，分别为20 kN至90 kN逐级（10 kN为一级）加载阶段，以及加载至破坏阶段。试验过程中，DIC设备会在未加载的初始阶段捕捉图片，作为后期变形识别的参照标准，同时在每一级加载时对试件进行拍摄，通过后期比对加载时图片与未加载时的图片的像素位置，获取结构的应变状态，即可计算并绘制结构的疲劳应力云图以及最大疲劳应力值，如图16所示。

图 16  疲劳应力云图

Fig. 16  Map of fatigue stress

3.3　结果及分析

利用DIC设备，分别捕捉了1张初始未加载状态下的结构图、50张加载过程中结构图以及10张破坏时刻结构图，通过图像处理与后期计算，可得结构在每级加载时内部产生的最大疲劳应力值。同时，本文计算了采用尼龙材料的传统结构与拓扑优化后的原型结构在相同等级荷载下产生的最大疲劳应力，并与试验观测的最大疲劳应力值进行比较，结果如图17所示。

图 17  最大疲劳应力计算值与实测值对比

Fig. 17  Comparison of calculated and measured values of maximum equivalent fatigue stress

试验结果表明，在每级荷载作用下，优化结构最大疲劳应力值均约为传统结构的45.0%，符合理论计算结果。优化结构的实测值普遍大于计算值，是由于在拓扑优化计算过程中，边界条件均为理想状态，故实际测量值与计算值有所不同。但总体而言，实测的最大疲劳应力值与理论计算值差距较小，且均低于传统结构的最大疲劳应力值。

为探究结构的最大承载极限，在逐级加载后，继续提升荷载作用水平直至结构破坏。当荷载加载至104 kN时，结构出现了过大的变形，并伴随拓扑优化结构的开裂，如图18所示。根据计算可得，尼龙材料建造的传统结构在同样边界条件下，当荷载加载至88 kN时会在底部支撑处产生超过材料抗拉强度的应力，因而可知优化后结构的承载能力相较于传统结构提升约18.3%。

图 18  试件破坏

Fig. 18  Failure of scale model specimen

由上述试验可知，经过拓扑优化设计后的结构可以在有效提高结构承载能力的同时降低结构质量，表明本文提出的拓扑优化设计方法有效可行。

4 结论

（1）综合考虑荷载疲劳应力与温度疲劳应力，构建了基于变密度法的水泥混凝土路面结构拓扑优化模型。采用拓扑梯度法，对拓扑优化问题的灵敏度求解通式进行推导，并结合移动渐近线更新算法，提出了适用于水泥混凝土路面结构的拓扑优化算法。

（2）数值分析结果表明，在轮迹零偏移的交通荷载场景下，拓扑优化结构的质量降低至传统结构的34.3%，产生的最大疲劳应力降低至传统结构的45.0%；在轮迹15 cm偏移的交通荷载场景下，其质量降低至传统结构的30.5%，最大疲劳应力降低至传统结构的37.5%。

（3）缩尺试验结果表明，拓扑优化原型结构在不同等级荷载下产生的最大疲劳应力值均小于相同荷载与边界条件下的传统结构，与理论值较为接近；同时，拓扑优化原型结构的极限承载能力相较于传统结构提升了18.3%，表明本文提出的拓扑优化设计方法有效可行。

拓扑优化设计方法在道路领域的实际应用较少，并且非均质的拓扑优化路面结构的建造具有一定技术难度，建议深入研究拓扑优化算法的适用性与可靠性，对拓扑优化设计的原型结构进行验算与路用性能分析，并探讨包括水泥混凝土3D打印技术在内的适用于轻量化装配式路面结构的建造方式。

作者贡献声明

高达辰：论文撰写，提出算法以及试验方案。

赵鸿铎：论文审定，论文整体结构与研究方法指导。

马鲁宽：研究思路和技术指导。

陈昊昱：算法公式推导与参数分析。

成 可：数据采集与分析。