CN101520814A - 一种采用耦合的变截面梁等效铰链波纹管的设计方法 - Google Patents

Abstract

一种采用耦合的变截面梁等效铰链波纹管的设计方法，用于简化含有铰链波纹管的复杂管路系统设计的有限元分析工作。本发明首先基于悬臂梁小挠度假设，推导了任意界面为矩形的变截面梁的自由端施加横向力作用下的力和位移关系，考虑变截面梁的弹性和弹塑性变形；其次提出了简化铰链波纹管模型的方法和利用有限元计算获得简化后的铰链波纹管固定部分两个方向刚度数据的方法；再次提出了利用混合遗传算法进行变截面梁参数的优化的方法，并提出针对混合遗传算法的空间渐缩方法，进而得到两根变截面梁的几何和材料参数；最后提出通过约束方程实现两根变截面梁在自由端的耦合，并通过建立辅助圆盘实现变截面梁根部与各自法兰的连接。

Description

一种采用耦合的变截面梁等效铰链波纹管的设计方法

技术领域

本发明涉及在管路系统有限元分析中采用耦合的变截面梁来代替铰链波纹管的设计方法。利用等效出的变截面梁可以简化含有铰链波纹管的复杂管路系统的有限元分析工作。

背景技术

在大跨度复杂管路系统中，铰链波纹管通常被用来吸收管路接头处的角位移，但也不可避免地要承受一定的横向位移。铰链波纹管结构较为复杂，其包括两端的法兰、波纹、圆环、单耳片、双耳片、铆钉、垫片、接头等零件，在铰链波纹管结构中存在10处接触，如圆环与耳片之间，铆钉与耳片之间等，如果利用有限元方法分析整个管路系统的在载荷作用下的响应时不将铰链波纹管进行简化或等效，那么整个模型将包含管壁的壳单元，即Shell elements，和铰链波纹管上的复杂实体单元，即Solid elements，又由于众多接触对的存在，大变形的考虑和结构弹塑性变形的存在，这很可能导致包含实际铰链波纹管的管路系统的有限元分析计算无法进行下去，即使可以进行，有限元调试的难度也很大，调试周期会很长，而这些对于急于把握系统的变形特性是很不利的，更使得管路有限元分析工作变得异常艰难。

发明内容

本发明提出一种将真实的铰链波纹管等效为两根耦合的变截面梁的设计方法，就很好地解决了含铰链波纹管管路系统有限元计算调试难度大的问题。在实际使用铰链波纹管的管路系统中，铰链波纹管所承担的主要是垂直于其轴线的横向位移，而承受的轴向的位移较小，所以本发明针对的主要是铰链波纹管横向刚度的等效。通常铰链波纹管在管路系统中安装后其横向位移方向和铰链波纹管中铆钉轴向方向平行或垂直，本发明只考虑这种安装情形。

等效的主要思路是：将铰链波纹管简化为固定部分和移动部分通过圆环连接的结构，而固定部分和移动部分在结构上除了法兰不同外其他零件都相同，见装配图10和11，从装配角度看，固定部分绕其轴线选择90°即为移动部分，设固定部分上与铆钉轴线平行的方向为Y向，与铆钉轴线垂直并与法兰端面平行的方向为Z向，所以如果得到在得到的在Y向和Z向与固定部分的结构刚度相同或相近的两根梁，再将两根梁在自由端耦合，这样耦合的变截面梁在横向与原始铰链波纹管是等效的，同时耦合的变截面梁也可以吸收角位移，达到在横向位移和角位移上等效实际铰链波纹管的目的，这里的梁采用变截面梁。

耦合的变截面梁具有以下四方面的特征：

a)变截面梁自由端横向力和位移的关系基于小变形假设，变截面梁的任意截面都是矩形，变截面梁的变形包括弹性变形和弹塑性变形，变截面梁采用双线性材料模型；

b)等效需要利用有限元方法获得简化的铰链波纹管固定端的两个相互垂直方向上的刚度数据，首先简化铰链波纹管得到固定部分模型：去除外部波纹、圆环和垫片，将铆钉与单而片和双耳片固连，双耳片端面与法兰连接，简化接头、铆钉、单而片、双耳片和法兰盘组成简化铰链波纹管的固定部分，再利用常用CAD软件建立简化铰链波纹管固定部分的几何模型；其次对简化的铰链波纹管固定部分进行数值试验：建立参考点RP，该参考点位于两个铆钉的中间空间部位，外力施加在参考点RP上，参考点与两个铆钉固连在一起，计算时法兰端面固定，载荷F施加在RP上，F与法兰端面平行，为了获得两个方向上的刚度数据，两次施加的F方向分别与铆钉轴线垂直和平行，与铆钉轴线垂直且与法兰端面平行的方向为Z向，与铆钉轴线平行的方向为Y向，通过计算获得Y方向和Z方向的力和位移数据；

c)铰链波纹管的等效采用混合遗传算法，采用基于前次搜索得到的最优点的空间渐缩方法不断调整遗传算法的搜索空间，最终得到两根等刚度的变截面梁。

d)通过约束方程将两根变截面梁在自由端进行耦合，并通过建立辅助圆盘实现变截面梁根部与各自法兰的连接，取代铰链波纹管，最终将等效的变截面梁应用在管路系统中。

变截面梁中间部位的矩形截面的尺寸为根部和自由端的截面尺寸的线性过度，其中根部截面的宽度为B0，高度为H0，自由端截面的宽度为B1，高度为H1；双线性材料模型需要参数为E、K和σ_s，其中E为材料弹性模量，K为屈服后应力应变增量之比，σ_s为材料屈服应力，最终遗传算法优化的7个参数为E，K，σ_s，B0，H0，B1，H1，得到求解变截面梁自由端横向刚度计算流程。变截面梁自由端横向刚度计算的流程如下：

8)输入变截面梁的几何参数：B0、B1、H0、H1和L，材料参数：E、K和σ_s，输入外力F；

9)计算结构的弹性极限载荷F_e，弹性极限挠度U_e，判断如果F≤F_e，输出U＝U_e*F/F_e，停止；否则进行下一步；

10)y为中性层挠度，x为水平坐标系，与中性层平行，计算弹塑性区的边界，L1，L2；

11)对于[0，L1]弹性段，带入边界条件 $\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=0}=0,$ y|_x＝0＝0，计算y(L1)，y′(L1)；

12)对于[L1，L2]弹塑性段，带入边界条件 $\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=L1}=y\′\left(L1\right),$ 计算y(L2)，y′(L2)；

13)对于[L2，L]弹性段，带入边界条件 $\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=L2}=y\′\left(L2\right),$ y|_x＝L2＝y(L2)，计算y(L)，y′(L)；

14)输出自由端挠度U＝y(L)。

上述流程中提及的F_e与U_e的确定方法如下：

首先求解一元三次方程

At³+BBt²+Ct+D＝0                        (1)

其中 $\left\{\begin{array}{c}A=2B{H}^2\\ \mathrm{BB}=-2\mathrm{BH}0H-H^{2}B0-3B{H}^2\\ C=4\mathrm{BH}0H+2H^{2}B0\\ D=-\mathrm{BH}{0}^2-2\mathrm{HB}0H0+B0H{0}^2\\ B=B0-B1,H=H0-H1\end{array}\right.$

该方程确定变截面梁最大应变的位置。由一元三次方程至少有一个实根，对于方程(1)求解出的结果，只取[0，1]范围内的解，如果求解出的实根在[0，1]之外，则真实的解为0。假设求解处最大值应变位置为t_max，令x_max＝t_maxL，最终可得极限弹性载荷为

$F_e=\frac{{\mathrm{\ϵ}}_s\mathrm{EB}\left(x_{\max }\right)H^2\left(x_{\max }\right)}{6(L-x_{\max })}---\left(2\right)$

对于F<F_e和F≥F_e，挠度计算方法如下：

i)横向力F<Fe，变截面梁处于弹性变形阶段

令 $d0=\frac{B0}{B0-B1},$ $d1=\frac{H0}{H0-H1},$ $G=\frac{-12L}{E(B1-B0){(H1-H0)}^3},$ 归一化令 $t=\frac{x}{L}.$

如果d0≠d1，则得截面转角和中性层挠度计算公式

$\frac{\mathrm{dy}}{\mathrm{dx}}\frac{1}{\mathrm{FGL}}=\frac{(1-d0)}{{(d0-d1)}^3}\ln \left(\frac{t-d1}{t-d0}\right)+\frac{d0-1}{{(d0-d1)}^2(t-d1)}+\frac{d1-1}{2(d0-d1){(t-d1)}^2}+C0---\left(3\right)$

     $\frac{y}{{\mathrm{FGL}}^2}=\frac{(1-d0)}{{(d0-d1)}^3}\left\{(t-d1)\ln \right|t-d1|-(t-d0)\ln |t-d0\left|\right\}$

                                      (4)

     $+\frac{d0-1}{{(d0-d1)}^2}\ln |t-d1|-\frac{d1-1}{2(d0-d1)(t-d1)}+C0t+C1$

如果d0＝d1，则有

$\frac{\mathrm{dy}}{\mathrm{dx}}=-\mathrm{GFL}[\frac{1}{2{(d0-t)}^2}+\frac{1-d0}{3{(d0-t)}^3}]+C0---\left(5\right)$

$y=-{\mathrm{GFL}}^2[\frac{1}{2(d0-t)}+\frac{1-d0}{6{(d0-t)}^2}]+C0\mathrm{Lt}+C1---\left(6\right)$

其中C0，C1为积分常数，通过边界条件确定，将(2)式带入(4)或(6)，并令t＝1求得弹性极限挠度U_e；

ii)横向力F≥Fe，变截面梁处于弹塑性变形阶段

需要求解如下微分方程

$\frac{d^{2}y}{d{x}^2}=\frac{3F(L-x)-3B\left(x\right)({\mathrm{\&sigma;}}_s-K{\mathrm{\&epsiv;}}_s)(H^2\left(x\right)/4-a^2)}{2B\left(x\right)[(E-K)a^3+K{H}^3\left(x\right)/8]}---\left(7\right)$

其中 ${\hat{z}}_s=\frac{\mathrm{EI}\left(x\right)}{F(L-x)}{\mathrm{\&epsiv;}}_s,$ a通过下式求解

$a^3+a(\frac{3F}{{\mathrm{\&epsiv;}}_s(E-K)}\frac{(L-x)}{B\left(x\right)}-3h^2)+2h^3\frac{K}{K-E}=0---\left(8\right)$

方程(7)的求解采用如下数值方法

1)设有边界条件 $\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=L1}=C2,$ y|_x＝L1＝C1，将区间[L1，L2]等分为n份，即

[x_i，x_i+1]，i＝1，...，n，x_i＝L1+(i-1)(L2-L1)/n，其中n自定，分别得到(7)式在各个区间上的数值积分 ${\mathrm{INT}2}_i=C2+{\&Integral;}_{x_i}^{x_{i+1}}\frac{d^{2}y}{d{x}^2}\mathrm{dx},$ i＝1，...，n，于是得到 ${\&Integral;}_{L1}^{L2}\frac{d^{2}y}{d{x}^2}\mathrm{dx}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{INT}2}_i=\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=L2};$

2)令INT2₀＝C2，可得在n+1个点处的值，即INT2_i，i＝0，...，n，这样便可插值得到的表达式，这里采用三次样条插值，插值得到n段三次多项式函数直接积分得到n段函数

在各自区间上的积分 ${\mathrm{INT}1}_i={\&Integral;}_{x_i}^{x_{i+1}}{f}_i^{\&prime;}\mathrm{dx},$ 则 ${\&Integral;}_{L1}^{L2}\frac{\mathrm{dy}}{\mathrm{dx}}\mathrm{dx}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{INT}1}_i=y{|}_{x=L2};$

c)对于L1，L2的确定方法如下

令 $Q=\frac{6\mathrm{FL}}{{\mathrm{\&epsiv;}}_s{\mathrm{EBH}}^2},$ 求解如下方程

t³+(-2d1-d0)t²+(d1²+2d0d1-Q)t+Q-d0d1²＝0，0≤t≤1        (9)

如果(9)式只求解出一个[0，1]间的实根t2，则必有t1＝0；如果求解出2个或3个实根，则只取在[0，1]间的根t1，t2，t1≤t2；则L1＝t1*L，L2＝t2*L。

混合遗传算法采用MATLAB遗传算法工具箱中的ga函数，设获得固定部分的刚度数据FY-UY和FZ-UZ后，从中挑选出在Y或Z方向的m个特征力Fi，i＝1，...，m作用下的位移分别为Uri，令列向量x＝(E，K，σ_s，B0，H0，B1，H1)^T，则变截面梁等效问题对应的优化问题表述为

$\min \mathrm{Res}\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(\mathrm{URi}-\mathrm{Ui})}^2=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(\mathrm{URi}-U(F_i,x^T))}^2---\left(10\right)$

s.t.xLB≤x≤xUB

其中[xLB，xUB]为x的取值范围，第一次搜索范围确定后，以后各次搜索范围通过空间渐缩方法确定。空间渐缩方法是指：首先指定一个较大的搜索空间S_i＝{x|xLB_i≤x≤xUB_i}，i＝0，该空间应该尽量大以包含理论上的全局最小值点，随后利用混合遗传算法在S_i内搜索，算法返回最优解x_i*及误差Res_i，如果Res_i大于设计的误差容限Err，则新空间S_i+i是以x_i*为中心，各维长度是原来的Scal倍，其中Scal<1，即

S_i+1＝{x|xLB_i+1≤x≤xUB_i+1}

其中 $\left\{\begin{array}{c}{\mathrm{xLB}}_{i+1}=x_i^{*}-\mathrm{Scalg}(x{\mathrm{UB}}_i-{\mathrm{xLB}}_i)/2\\ {\mathrm{xUB}}_{i+1}=x_i^{*}+\mathrm{Scalg}({\mathrm{xUB}}_i-{\mathrm{xLB}}_i)/2\end{array}\right.---\left(11\right)$

关于变截面梁的耦合及其与法兰的连接方法如下：

a)变截面梁的耦合：采用有限元方法需要将变截面梁划分成实体单元，而两个变截面梁在自由端需要耦合以模拟铰链波纹管上法兰盘之间的类似万向节的转动结构；设两根等刚度的变截面梁中一个变截面梁自由端的4个角节点编号分别为Ni，i＝1，2，3，4，另一个变截面梁自由端的4个角节点编号分别为Ni，i＝5，6，7，8，则建立耦合的需要建立如下的约束方程：

$\left\{\begin{array}{c}0=\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{\mathrm{ux}}_{\mathrm{Ni}}-\underset{i=5}{\overset{8}{\mathrm{\&Sigma;}}}{\mathrm{ux}}_{\mathrm{Ni}}\\ 0=\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{\mathrm{uy}}_{\mathrm{Ni}}-\underset{i=5}{\overset{8}{\mathrm{\&Sigma;}}}{\mathrm{uy}}_{\mathrm{Ni}}\\ 0=\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{\mathrm{uz}}_{\mathrm{Ni}}-\underset{i=5}{\overset{8}{\mathrm{\&Sigma;}}}{\mathrm{uz}}_{\mathrm{Ni}}\end{array}\right.---\left(12\right)$

其中ux_Ni，uy_Ni，uz_Ni，分别表示Ni节点的x，y，z方向的位移。

b)变截面梁与法兰的连接：在有限元软件中，对于整个管路的几何模型，在两个变截面梁的根部分别建立两个辅助圆盘与相应的法兰连接，辅助圆盘的直径与法兰的内径相同，并将辅助圆盘与法兰的内缘固连，辅助圆盘与变截面梁的根部贴合并固连，辅助圆盘与变截面梁同轴。辅助圆盘采用弹性材料模型，材料的弹性模量为管路系统管壁材料弹性模量的10倍，辅助圆盘厚度与管壁厚度相同，泊松比和线膨胀系数与管壁材料相同。

最终，通过约束方程建立变截面梁之间以及变截面梁与法兰的连接之后，便可进行整个管路系统的有限元分析。采用等效的变截面梁方法，省略了复杂的实体铰链波纹管的建模和网格划分工作，使得在整体分析阶段不用考虑波纹管中众多的接触问题，这样有限元计算工作可以顺利进行下去，有限元调试的难度也大大降低了，调试周期也大大缩短了，因此，利用简化后的变截面弹塑性可以较快地把握系统的变形特性；同时，如果其他管路系统中使用到同尺寸的铰链波纹管，那么前期得到的变截面梁的参数同样可以使用，所以简化等效工作带来的有限元分析工作的简化是非常可观的，本发明具有较好工业应用价值。

附图说明

图1为采用耦合的变截面梁等效铰链波纹管的实施流程图。

图2为变截面梁的5个几何参数B0，H0，B1，H1，L示意图。

图3为变截面梁的3个材料参数E，K，σ_s示意图，采用双线性材料模型。

图4(a)为变截面梁中性层挠度示意图，含坐标系；(b)为在弹性变形范围内变截面梁任意截面的应力分布示意图。

图5为弹塑性区边界L1，L2的示意图以及屈服区域，此时F>F_e。其中下凹的虚线为假设Q截面仍满足弹性应力应变关系前提下在

方向应变达到εs的屈服范围的下界高度

，下凹的实线为Q截面屈服后真实的应力达到σ_s的点的集合，屈服范围的真实下界高度a。

图6为Q假设满足弹性应力应变关系和满足真实应力应变关系应力分布，及

与a的关系图。

图7为屈服范围的真实下界高度a的关系图。

图8为变截面梁挠度计算流程图。

图9为真实铰链波纹管的CAD模型，其中忽略了法兰盘上的螺纹孔。

图10为铰链波纹管不含外部波纹的装配图。

图11为铰链波纹管半剖图。

图12为采用有限元方法计算的铰链波纹管的固定部分刚度的边界示意图，图示为位移边界和施加在参考点RP上的Z向外力FZ，其中参考点RP和铆钉固连在一起。

图13为采用有限元方法计算的铰链波纹管的固定部分刚度的边界示意图，图示为位移边界和施加在参考点RP上的Y向外力FY，其中参考点RP和铆钉固连在一起。

图14为混合遗传算法所采用的空间渐缩方法的流程图。

图15与图12对应，为在参考点上的Z方向施加外力后结构的有限元计算变形图。

图16与图13对应，为在参考点上的Y方向施加外力后结构的有限元计算变形图。

图17为采用耦合的变截面梁等效铰链波纹管的实施后的效果图，(a)为简化后的铰链波纹管，(b)为耦合的等效变截面梁，图示的F方向为横向，即等效所在的方向。

其中：1-移动端法兰，2-双而片，3-接头，4-铆钉，5-波纹(见图12)，6-移动端法兰，7-圆环，8-单耳片。

具体实施方式

如附图1所示的发明实施步骤，本发明相应包括如下4部分内容：1)推导悬臂的变截面梁自由端横向力和位移关系的表达式；2)对铰链波纹管进行简化，通过有限元计算获得真实铰链波纹管两部分刚度数据；3)利用混合遗传算法进行各个变截面梁参数的优化得到两根在横向等刚度的变截面梁；4)给出变截面梁的耦合及其与法兰的连接方法。下面分四部分分别介绍。

1 小挠度变截面梁的横向力和位移关系

1.1 变截面梁几何参数

说明书附图2显示了变截面梁的5个几何参数：固定端截面的宽度B0，高度H0，自由端截面的宽度B1，高度H1，长度L，由于截面的宽度和高度线性变化，则变截面梁任意截面处的宽度B(x)、高度H(x)为

$\left\{\begin{array}{c}B\left(x\right)=B1-\frac{x}{L}(B1-B0)\\ H\left(x\right)=H1-\frac{x}{L}(H1-H0)\end{array}\right.0\&le;x\&le;L---\left(14\right)$

任意截面惯性矩 $I\left(x\right)=\frac{1}{12}B\left(x\right)H^3\left(x\right),0\&le;x\&le;L.$

1.2 材料模型及参数

如图3所示，这里采用双线性材料模型，涉及到3个参数：材料弹性模量E，屈服应力σs，屈服后材料应力应变增量之比K，相应的应力应变关系如下

$\mathrm{\&sigma;}=\left\{\begin{array}{cc}\mathrm{E\&epsiv;}& \mathrm{\&epsiv;}\&le;{\mathrm{\&epsiv;}}_s\\ {\mathrm{\&sigma;}}_s+K(\mathrm{\&epsiv;}-{\mathrm{\&epsiv;}}_s)& \mathrm{\&epsiv;}>{\mathrm{\&epsiv;}}_s\end{array}\right.---\left(15\right)$

其中ε_s＝σ_s/E为屈服应变。

图4显示了在F作用下，任意截面上的应力分布，与中性层相垂直的轴为

轴。在小挠度假设下，中性层曲率 $k=\frac{1}{\mathrm{\&rho;}}\&ap;\frac{d^{2}y}{d{x}^2},$ 则截面应变 ${\mathrm{\&epsiv;}}_x\left(\hat{z}\right)=-\frac{\hat{z}}{\mathrm{\&rho;}}=-\hat{z}\frac{d^{2}y}{d{x}^2}.$

1.3 弹性变形下自由端的挠度

假设使得整个变截面梁进入塑性变形的最小载荷为F_e，弹性变形及施加载荷F≤F_e。在小挠度下，任意截面的弯矩M(x)＝-F(L-x)。

已知变截面梁的小挠度微分方程

$\frac{d^{2}y}{d{x}^2}=-\frac{M\left(x\right)}{\mathrm{EI}\left(x\right)}---\left(16\right)$

方程(16)可以显示地积分出来：

令 $d0=\frac{B0}{B0-B1},$ $d1=\frac{H0}{H0-H1},$ $G=\frac{-12L}{E(B1-B0){(H1-H0)}^3},$ 归一化令 $t=\frac{x}{L}.$

如果d0≠d1，则得截面转角和中性层挠度计算公式

$\frac{\mathrm{dy}}{\mathrm{dx}}\frac{1}{\mathrm{FGL}}=\frac{(1-d0)}{{(d0-d1)}^3}\ln \left(\frac{t-d1}{t-d0}\right)+\frac{d0-1}{{(d0-d1)}^2(t-d1)}+\frac{d1-1}{2(d0-d1){(t-d1)}^2}+C0---\left(17\right)$

     $\frac{y}{{\mathrm{FGL}}^2}=\frac{(1-d0)}{{(d0-d1)}^3}\left\{(t-d1)\ln \right|t-d1|-(t-d0)\ln |t-d0\left|\right\}$

                                      (18)

     $+\frac{d0-1}{{(d0-d1)}^2}\ln |t-d1|-\frac{d1-1}{2(d0-d1)(t-d1)}+C0t+C1$

如果d0＝d1，则有

$\frac{\mathrm{dy}}{\mathrm{dx}}=-\mathrm{GFL}[\frac{1}{2{(d0-t)}^2}+\frac{1-d0}{3{(d0-t)}^3}]+C0---\left(19\right)$

$y=-{\mathrm{GFL}}^2[\frac{1}{2(d0-t)}+\frac{1-d0}{6{(d0-t)}^2}]+C0\mathrm{Lt}+C1---\left(20\right)$

其中C0，C1为积分常数，通过边界条件确定。

1.4 弹塑性变形下自由端的挠度

对于给定的载荷F，首先需要判断F与F_e的关系，下一小节主要解决这个问题。

1.4.1 极限弹性载荷的F_e确定

由同一截面上应变分布的对称性可知，变截面梁的应力应变分布关于中性层对称，同时易知同一截面上最大应力应变一定出现在最外侧，所以判断极限载荷就考虑变截面梁的最外侧表面是否屈服。考虑图2所示变截面梁上侧最外表面M的应变可表达为：

${\mathrm{\&epsiv;}}_o={\mathrm{\&epsiv;}}_x\left(H\left(x\right)\right)=\frac{-H\left(x\right)}{2}\frac{d^{2}y}{d{x}^2}=\frac{-M\left(x\right)H\left(x\right)}{2\mathrm{EI}\left(x\right)}---\left(21\right)$

上式是关于x的函数，首先要找出最大应变的位置，得到该处的最大应变ε_omax，然后令ε_omax＝ε_s便可得到F_e。细化ε_o得

${\mathrm{\&epsiv;}}_o=\frac{6F(L-x)}{\mathrm{EB}\left(x\right)H^2\left(x\right)}=\frac{6\mathrm{FL}}{E}\frac{1-\frac{x}{L}}{(B0-\frac{x}{L}(B0-B1)){(H0-\frac{x}{L}(H0-H1))}^2}0\&le;x\&le;L---\left(22\right)$

令B＝B0-B1，H＝H0-H1，问题转化为求下式的最小值

$y=\frac{(B0-\mathrm{tB}){(H0-\mathrm{tH})}^2}{1-t}0\&le;t\&le;1---\left(23\right)$

将(23)两边对x求导得y′，并令y′＝0，问题转化为求前面的一元三次方程

At³+BBt²+Ct+D＝0                (24)

其中 $\left\{\begin{array}{c}A=2B{H}^2\\ \mathrm{BB}=-2\mathrm{BH}0H-H^{2}B0-3B{H}^2\\ C=4\mathrm{BH}0H+2H^{2}B0\\ D=-\mathrm{BH}{0}^2-2\mathrm{HB}0H0+B0H{0}^2\\ B=B0-B1,H=H0-H1\end{array}\right.$

由一元三次方程至少有一个实根，对于方程(1)求解出的结果，只取[0，1]范围内的解，如果求解出的实根在[0，1]之外，则真实的解为0，原因是实际上(23)式是一个有边界的函数最值问题，如果(24)式求出的驻点位于[0，1]之外则实际最值一定在边界上取得，而1肯定不是符合实际的真实解，因为该处的弯矩永远为0，不可能屈服，所以真实解只能为0。

假设取得通过(22)-(24)式综合求出的最大值应变位置为t_max，令x_max＝L×t_max，则有最大应变

${\mathrm{\&epsiv;}}_{o\max }=\frac{6F(L-x_{\max })}{\mathrm{EB}\left(x_{\max }\right)H^2\left(x_{\max }\right)}---\left(25\right)$

最终可得极限弹性载荷F_e

$F_e=\frac{{\mathrm{\&epsiv;}}_s\mathrm{EB}\left(x_{\max }\right)H^2\left(x_{\max }\right)}{6(L-x_{\max })}---\left(26\right)$

令t＝1，再将F_e分别带入公式(18)和(20)便可得到自由端F_e对应下的弹性极限挠度U_e。

1.4.2 弹塑性区边界L1，L2的计算

当F>F_e时，首要任务是确定变截面梁上屈服区域的范围，即图5中的[L1，L2]段。与确定极限载荷的思路类似，只需确定出变截面梁的外表面哪些点的应变达到了屈服应变，只需求解下面的方程：

${\mathrm{\&epsiv;}}_s=\frac{6F(L-x)}{\mathrm{EB}\left(x\right)H^2\left(x\right)},0\&le;x\&le;L---\left(27\right)$

细化上述方程：令 $Q=\frac{6\mathrm{FL}}{{\mathrm{\&epsiv;}}_s{\mathrm{EBH}}^2},$ 方程(27)转化为

t³+(-2d1-d0)t²+(d1²+2d0d1-Q)t+Q-d0d1²＝0，0≤t≤1        (28)

(28)同样是一个一元三次方程，也存在着根的取舍问题。同样的道理，变截面梁在自由端不可能屈服，如果(28)只求解出一个[0，1]间的实根t2，则必有t1＝0；如果求解出2个或3个实根，则只取在[0，1]间的根t1，t2，t1≤t2。则L1＝t1*L，L2＝t2*L。

1.4.3 屈服区域的计算

通过(28)已经确定了弹塑性区边界，即屈服区域在x方向的范围，而屈服区域在

方向的分布尚未确定，这将在本节讨论。

如图6所示，当F>F_e，假设Q截面仍满足弹性应力应变关系，(图5中的虚线为变截面梁上应力达到σ_s的点)，这时Q截面上的弯矩(即图6中Q截面上ODS三角形内应力产生的弯矩的两倍)为

$M_Q=2\&times;\frac{1}{2}h{\mathrm{\&sigma;}}_{A}B\left(x\right)\frac{2}{3}h=2\&times;\frac{1}{3}B\left(x\right)h^2{\mathrm{\&sigma;}}_{s}h/{\hat{z}}_s---\left(29\right)$

其中 ${\hat{z}}_s=\frac{\mathrm{EI}\left(x\right)}{F(L-x)}{\mathrm{\&epsiv;}}_s$ 表示假设Q截面仍满足弹性应力应变关系前提下在

方向应变达到ε_s的屈服范围的下界高度，在

范围内应变大于ε_s。

真实的双线性材料的应力在

方向分布为阴影部分，其弯矩为

$M_Q=2[\frac{1}{2}a{\mathrm{\&sigma;}}_{s}B\left(x\right)\frac{2}{3}a+((h-a)/2+a)(h-a){\mathrm{\&sigma;}}_{s}B\left(x\right)$

     $+B\left(x\right)({\mathrm{\&sigma;}}_p-{\mathrm{\&sigma;}}_s)(a+2(h-a)/3)(h-a)/2]---\left(30\right)$

   $=2B\left(x\right)[a^2{\mathrm{\&sigma;}}_s/3+(h^2-a^2){\mathrm{\&sigma;}}_s/2+({\mathrm{\&sigma;}}_p-{\mathrm{\&sigma;}}_s)(2h+a)(h-a)/6]$

联立(29)与(30)得到

${\mathrm{\&sigma;}}_s{h}^3/\left(3{\hat{z}}_s\right)=a^2{\mathrm{\&sigma;}}_s/3+(h^2-a^2){\mathrm{\&sigma;}}_s/2+({\mathrm{\&sigma;}}_p-{\mathrm{\&sigma;}}_s)(2h+a)(h-a)/6----\left(31\right)$

其中h＝H(x)/2，a为变截面梁在方向应力达到σ_s的点，由(31)还不能确定a，因为σ_p未知，还需要添加下列关系。

因为C、D在同一个截面上，其挠度y的y＂相同，而C点的应力应变关系为

-ay＂E＝σ_s                         (32)

C、D点的应力应变关系为

(-hy＂+ay＂)K＝σ_p-σ_s               (33)

(25)与(26)联立得

(h-a)K/aE＝(σ_p-σ_s)/σ_s              (34)

(24)与(27)联立可以求得a与

的关系为

$a^3+a(\frac{3F}{{\mathrm{\&epsiv;}}_s(E-K)}\frac{(L-x)}{B\left(x\right)}-3h^2)+2h^3\frac{K}{K-E}=0---\left(35\right)$

1.4.4 弹塑性区微分方程的建立

如图7所示，积分区域 $\&PartialD;\mathrm{\&Omega;}=A_e+A_p,$ 其中A_e为弹性变形区，A_P为塑性变形区。则弹塑性区截面弯矩为

$M\left(x\right)=-\underset{\&PartialD;\mathrm{\&Omega;}}{\&Integral;}\mathrm{\&sigma;}\left(\hat{z}\right)\hat{z}\mathrm{dA}=-\underset{A_e}{\&Integral;}\mathrm{\&sigma;}\left(\hat{z}\right)\hat{z}\mathrm{dA}-\underset{A_p}{\&Integral;}\mathrm{\&sigma;}\left(\hat{z}\right)\hat{z}\mathrm{dA}$

    $=-2({\&Integral;}_0^{a}E{\mathrm{\&epsiv;}}_x\left(\hat{z}\right)B\left(x\right)\hat{z}d\hat{z}+{\&Integral;}_a^{H\left(x\right)/2}({\mathrm{\&sigma;}}_s+K({\mathrm{\&epsiv;}}_x\left(\hat{z}\right)-{\mathrm{\&epsiv;}}_s))B\left(x\right)\hat{z}d\hat{z})---\left(36\right)$

    $=-\frac{2}{3}B\left(x\right)((E-K)a^3+{\mathrm{KH}}^3\left(x\right)/8)\frac{d^{2}y}{d{x}^2}-B\left(x\right)({\mathrm{\&sigma;}}_s-{\mathrm{K\&epsiv;}}_s)(H^2\left(x\right)/4-a^2)$

由M(x)＝-F(L-x)得最终微分方程为

$\frac{d^{2}y}{d{x}^2}=\frac{3F(L-x)-3B\left(x\right)({\mathrm{\&sigma;}}_s-K{\mathrm{\&epsiv;}}_s)(H^2\left(x\right)/4-a^2)}{2B\left(x\right)[(E-K)a^3+K{H}^3\left(x\right)/8]}---\left(37\right)$

(37)无法显式求解，这里采用了下列数值求解的办法。

1)设有边界条件 $\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=L1}=C2,$ y|_z＝L1＝C1，将区间[L1，L2]等分为n份，即

[x_i，x_i+1]，i＝1，..，n，x_i＝L1+(i-1)(L2-L1)/n，其中n自定，分别得到(37)式在各个区间上的数值积分 ${\mathrm{INT}2}_i=C2+{\&Integral;}_{x_i}^{x_{i+1}}\frac{d^{2}y}{d{x}^2}\mathrm{dx},$ i＝1，...，n，积分方法任选，比如辛普森求积Simpson，高斯积分Gauss，龙贝格积分Romberg等，于是得到 ${\&Integral;}_{L1}^{L2}\frac{d^{2}y}{d{x}^2}\mathrm{dx}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{INT}2}_i=\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=L2};$

2)令INT2₀＝C2，可得

在n+1个点处的值，即INT2_i，i＝0，...，n，这样便可插值得到

的表达式，这里采用三次样条插值，插值得到n段三次多项式函数

，三次多项式函数可以直接积分，得到n段函数

在各自区间上的积分 ${\mathrm{INT}1}_i={\&Integral;}_{x_i}^{x_{i+1}}{f}_i^{\&prime;}\mathrm{dx},$ 于是得到

${\&Integral;}_{L1}^{L2}\frac{dy}{dx}\mathrm{dx}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{INT}1}_i={y|}_{x=L2};$

1.4.5 变截面梁挠度计算流程图

图8为整个变截面梁自由端挠度计算的流程，基本内容如下：

1)输入变截面梁的几何参数B0，B1，H0，H1，L和材料参数E，K，σ_s，输入外力F；

2)计算结构的弹性极限载荷F_e，弹性极限挠度U_e，判断如果F≤F_e，输出U＝Ue*F/F_e，停止；否则进行下一步；

3)计算弹塑性区的边界，L1，L2；

4)对于[0，L1]弹性段，带入边界条件 $\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=0}=0,$ y|_x＝0＝0，通过计算(4)或(6)计算y(L1)，通过(3)或(5)计算y′(L1)；

5)对于[L1，L2]弹塑性段，带入边界条件 $\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=L1}=y\&prime;\left(L1\right),$ y|_x＝L1＝y(L1)，通过(7)计算y(L2)，y′(L2)；

6)对于[L2，L]弹性段，带入边界条件 $\frac{\mathrm{dy}}{\mathrm{dx}}{|}_{x=L2}=y\&prime;\left(L2\right),$ y|_x＝L2＝y(L2)，通过计算(4)或(6)计算y(L)，通过(3)或(5)计算y′(L)；

7)输出自由端挠度U＝y(L)；

综上，可以得到一个非完全解析的函数U＝U(F，E，K，σ_s，L，B0，H0，B1，H1)，表示在9个参数给定条件下便可惟一确定自由端的挠度U。

2 获得真实铰链波纹管两部分刚度数据

本发明提出的利用两个耦合的变截面梁代替实际铰链波纹管的方法主要考虑实际铰链波纹管横向刚度的等效，所以等效前必须获得实际铰链波纹管的刚度数据，本发明提出了数值试验方法：即采用有限元方法，建立好铰链波纹管的固定部分的几何模型，固定法兰端面，并在参考点RP上施加横向力进行计算，最终获得在参考点RP横向载荷和位移关系数据。

简化铰链波纹管：由于铰链波纹管结构复杂，本发明提出在进行物理试验或有限元数值模拟时将需要试验的结构进行简化。图9为实际铰链波纹管的CAD模型，其中两个法兰盘的端面的螺钉孔被去除了，实际铰链波纹管外部有多层波纹，但经过计算单层波纹的轴向刚度是相同材料同厚度圆筒轴向刚度1/1000左右，所以本发明在进行刚度计算时忽略了外部波纹，将外部的波纹去掉，同时将垫片与圆环固连，将铆钉与单而片和双耳片固连，忽略了接头外凸部分只保留接头与单耳片和双耳片焊接的部分，即一个较短的圆柱筒壁，称之为简化接头。最终图10为简化铰链波纹管，图11是简化铰链波纹管的半剖图。

简化铰链波纹管固定部分数值试验：将简化铰链波纹管分成两部分，即固定部分和移动部分，固定部分和移动部分分别由简化接头、铆钉、单而片、双耳片和法兰盘连接而成，这两部分通过圆环连接在一起，固定部分和移动部分在结构上非常相似，只有各自的法兰盘外径不同，固定端稍大，但是这对两部分的刚度试验没有影响，所以本发明数值试验方法采用同一种结构，只是施加载荷的方向不同。图12和图13是数值试验的简化铰链波纹管的固定部分，但是这里的固定部分的接头结构也被简化了，与图10相比忽略了外部不受力的接头外凸部分，采用有限元分析软件进行固定部分刚度的计算工作，外力施加在参考点RP上，该点位于两个铆钉的中间部位，参考点与两个铆钉固连在一起。计算时法兰端面固定，为了得到固定端在Y和Z方向的刚度，载荷FY和FZ分别施加在RP上，FY与FZ和法兰端面平行，为了获得两个方向上的刚度数据，FY方向与铆钉轴线平行，即Y向，FZ方向与铆钉轴线垂直并与法兰端面平行，即Z向。通过计算分别获得两个方向上的FY-UR和FZ-UR数据，UR为参考点在对应方向上位移。

需要注意的是，推导变截面梁考虑的是小变形假设，所以FY与FZ不需要施加太大，以产生的位移不要太大，一般要求位移不超过变截面梁长度的10％。

3变截面梁参数的优化

综合上面两节可以看出，可以得到一个非完全解析的函数U＝U(F，E，K，σ_s，L，B0，H0，B1，H1)，该函数表示在对于变截面梁，在材料参数、结构参数确定的条件下，给定外力F便可确定自由端的挠度U。对于实际的某个铰链波纹管而言，其长度L是确定的，所以就不作为优化参数，而F只是外部作用力，不是变截面梁本身的属性，所以也不作为优化参数，这样优化问题的设计参数减少为7个：E，K，σ_s，B0，H0，B1，H1，则非完全解析的函数变为U＝U(F，E，K，σ_s，B0，H0，B1，H1)。

优化的目标是得到与实际铰链波纹管的两部分刚度分别相同的两个变截面梁，评价变截面梁的刚度与实际波纹管刚度差距的方法是计算相同力作用的自由端位移差，以Y方向为例，为了减少计算量首先从FY-UR中提取m个特征力Fi，i＝1，...，m作用下的位移分别URi，而变截面梁在Fi作用的位移为Ui，则优化目标函数为

$\mathrm{Res}=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(\mathrm{URi}-\mathrm{Ui})}^2=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(\mathrm{URi}-U(F_i,E,K,{\mathrm{\&sigma;}}_s,B0,H0,B1,H1))}^2---\left(38\right)$

可以把(38)式看作Res是(E，K，σ_s，B0，H0，B1，HI)的函数，令列向量x＝(E，K，σ_s，B0，H0，B1，H1)^T，同时设x的取值范围为[xLB，xUB]，则优化问题表述为

$\min \mathrm{Res}\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(\mathrm{URi}-\mathrm{Ui})}^2=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(\mathrm{URi}-U(F_i,x^T))}^2---\left(39\right)$

s.t.xLB≤x≤xUB

下一个问题就是xLB和xUB应该如何取，这应该和混合遗传算法结合起来考虑。

本发明采用混合遗传算法来解决铰链波纹管刚度等效的问题，这里所谓的混合即在遗传算法运行结束后，采用传统的局部优化算法以遗传算法返回的最优个体为初值进行局部解的改进，以提高解的精度，弥补遗传算法局部搜索精度的不足。目前科学计算语言MATALB具有非常强大的遗传算法工具箱，该工具箱可以实现经典遗传算法与多种局部的混合，比如单纯形搜索，高斯-牛顿法等方法，这里的混合遗传算法采用MATLAB遗传算法工具箱中的ga函数。由于ga函数需要提供变量取值范围，而对于等效变截面梁问题本身而言，一开始并不能准确给定各个变量的取值范围，本发明提出空间渐缩的方法。

空间渐缩，是一个取值范围不断缩减的方法。由于搜索空间的大小对遗传算法搜索成功率影响较大，所以需要尽量控制搜索空间。图14为空间渐缩方法的流程图，空间渐缩，首先指定一个较大的搜索空间S_i＝{x|xLB_i≤x≤xUB_i}，i＝0，该空间应该尽量大以包含理论上的全局最小值点，随后利用混合遗传算法在S_i内搜索，算法返回最优解x_i*及误差Res_i，如果Res_i大于设计的误差容限且x_i*优于上次遗传算法返回的解，则新空间S_i+i是以x_i*为中心，各维长度是原来的Scal倍，Scal<1，即

S_i+1＝{x|xLB_i+1≤x≤xUB_i+1}

其中 $\left\{\begin{array}{c}{\mathrm{xLB}}_{i+1}=x_i^{*}-\mathrm{Scalg}(x{\mathrm{UB}}_i-{\mathrm{xLB}}_i)/2\\ {\mathrm{xUB}}_{i+1}=x_i^{*}+\mathrm{Scalg}({\mathrm{xUB}}_i-{\mathrm{xLB}}_i)/2\end{array}\right.---\left(40\right)$

由于通过一次混合遗传算法优化通常不能得到理想的结果，所以需要不断地在新的小空间内搜索才能保证较好的优化结果。

4 变截面梁的耦合及其与法兰的连接

在有限元分析中，两根变截面梁采用具有弹塑性分析能力的实体单元。在利用混合遗传算法优化获得两根不同尺寸的变截面梁后并划分完单元后，需要将这两根变截面梁进行耦合，并将各根梁与各自的法兰进行连接才能与原始铰链波纹管在横向刚度等效。下面分两部分说明变截面梁的耦合及其与法兰的连接。

1)变截面梁的耦合：采用有限元方法需要将变截面梁划分成实体单元Solid Element，而两个变截面梁在自由端需要耦合以模拟铰链波纹管上法兰盘之间的类似万向节的转动结构。设两个变截面梁中一个梁自由端的4个角节点编号分别为Ni，i＝1，2，3，4，另一个梁自由端的4个角节点编号分别为Ni，i＝5，6，7，8，则建立耦合的需要建立式(12)的约束方程。其中ux_Ni，uy_Ni，uz_Ni，分别表示Ni节点的x，y，z方向的位移。该约束方程的作用是耦合两个自由端的矩形截面的形心的在空间三个方向的自由度，即相当于一个球铰链的作用，在图17(b)中的两个变截面梁接头处的球即表示耦合。

2)变截面梁与法兰的连接：实际的铰链波纹管的法兰是与双耳片焊接在一起的，而采用本发明把铰链波纹管的中间多个零件等效成了两根变截面梁，所以变截面梁与法兰之间的连接关系同样是焊接关系。在有限元软件中，对于整个管路的几何模型，在两个变截面梁的根部分别建立两个辅助圆盘与相应的法兰连接，辅助圆盘的直径与法兰的内径相同，并将辅助圆盘与法兰的内缘固连，辅助圆盘与变截面梁的根部贴合并固连，辅助圆盘与变截面梁同轴。辅助圆盘采用弹性材料模型，材料的弹性模量为管路系统管壁材料弹性模量的10倍，辅助圆盘厚度与管壁厚度相同，泊松比和线膨胀系数与管壁材料相同。

实现变截面梁的耦合及其与法兰的连接，就可以将本发明提出的变截面梁应用在整个管路系统的有限元分析中，从而大大简化由铰链波纹管带来的强非线性分析工作。

下面结合图1说明本发明的实施流程。

本发明的实施分四步：步骤1、编写挠度计算程序；步骤2、刚度数据准备；步骤3、利用遗传算法优化；步骤4、建立约束方程和辅助圆盘。下面详细说明这四步。

步骤1、编写挠度计算程序

该步骤的目的是编程实现推导出的变截面梁自由端横向力和位移关系表达式。根据图8的变截面梁自由端挠度计算流程图，编写挠度计算程序，程序的最终的基本格式可以为

U＝U(F，E，K，σ_s，L，B0，H0，B1，H1)，表示输入变截面梁的几何参数(L，B0，H0，B1，H1)、材料参数(E，K，σ_s)、外力F后返回自由端的挠度U。

步骤2、刚度数据准备

利用有限元数值模拟的方法获得固定部分的FY-UR和FZ-UR数据。图12和图13显示的是有限元计算时采用的几何模型和模型的边界，其中法兰盘端面固定，外力FY和FZ分别施加在参考点RP上，最终得到FY，FZ与RP点在Y，Z方向上的位移实验数据UR，需要注意的是，推到变截面梁考虑的是小变形假设，所以无论物理试验还是数值试验都不需要施加太大的外力，以产生的位移不要太大，一般要求位移不超过变截面梁长度的10％。

从获得的FY-UR，FZ-UR数据中各提出m个特征数据点作为结构的真实刚度数据，m个特征力Fi，i＝1，...，m作用下的位移分别为URi，这些数据将应用在遗传算法优化中。特征数据即可以反应整个试验数据特征的数据。

步骤3、利用遗传算法优化

优化需要分两部分进行，最终得到两种尺寸的变截面梁。下面以施加FY方向力的变截面梁的优化为例，说明如何利用混合遗传算法。

通过步骤1和2已经得到了在Y和Z方向的挠度计算流程，可以总体表示如下：

U＝U(F，E，K，σ_s，L，B0，H0，B1，H1)                (41)

其中对于实际的某种型号铰链波纹管而言，L为常数，表示固定端和移动端法兰接头间的距离的一半。

已经获得FY-UR和FZ-UR试验数据中的m个特征数据，Fi，i＝1，...，m，然后利用公式(41)结合流程图14的空间渐缩方法，则可完成在Y和Z方向上整个变截面梁7个参数E，K，σ_s，B0，H0，B1，H1的优化。

需要注意的是，空间渐缩方法需要提供初始搜索空间S₀，即各个参数的取值范围，这里S₀要足够大以包含所需的最优解。终止误差Err，缩减因子Scal需要根据实际问题而定。

步骤4、建立约束方程和辅助圆盘

为了将优化得到的变截面梁应用在实际管路系统的有限元分析中，需要将这两根变截面梁进行耦合，并将各根变截面梁与各自的法兰进行连接。

变截面梁的耦合：获取两个变截面梁中一个梁自由端的4个角节点编号：Ni，i＝1，2，3，4，另一个梁自由端的4个角节点编号：Ni，i＝5，6，7，8，然后在有限元分析中建立(12)式的约束方程。

变截面梁与法兰的连接：在两个变截面梁的根部分别建立两个辅助圆盘与相应的法兰连接，辅助圆盘的直径与法兰的内径相同，并将辅助圆盘与法兰的内缘固连，辅助圆盘与变截面梁的根部贴合并固连，辅助圆盘与变截面梁同轴。辅助圆盘采用弹性材料模型，材料的弹性模量为管路系统管壁材料弹性模量的10倍，辅助圆盘厚度与管壁厚度相同，泊松比和线膨胀系数与管壁材料相同。

工业上的可利用性

本发明主要为简化含铰链波纹管的管路系统有限元分析工作。获得铰链波纹管两部分的刚度数据后，编程实现变截面梁自由端挠度计算，同时利用现有的混合遗传算法，便可得到等刚度的变截面梁，通过约束方程建立变截面梁之间，通过建立辅助圆盘实现变截面梁与法兰的连接，进而简化整个管路的有限元分析工作，本发明具有较好工业应用价值。

Claims (7)

1.一种采用耦合的变截面梁等效铰链波纹管的设计方法，其特征在于：

a)变截面梁自由端横向力和位移的关系基于小变形假设，变截面梁的任意截面都是矩形，变截面梁的变形包括弹性变形和弹塑性变形，变截面梁采用双线性材料模型；

b)等效需要利用有限元方法获得简化的铰链波纹管固定端的两个相互垂直方向上的刚度数据，首先简化铰链波纹管得到固定部分模型：去除外部波纹、圆环和垫片，将铆钉与单而片和双耳片固连，双耳片端面与法兰连接，简化接头、铆钉、单而片、双耳片和法兰盘组成简化铰链波纹管的固定部分，再利用常用CAD软件建立简化铰链波纹管固定部分的几何模型；其次对简化的铰链波纹管固定部分进行数值试验：建立参考点RP，该参考点位于两个铆钉的中间空间部位，外力施加在参考点RP上，参考点与两个铆钉固连在一起，计算时法兰端面固定，载荷F施加在RP上，F与法兰端面平行，为了获得两个方向上的刚度数据，两次施加的F方向分别与铆钉轴线垂直和平行，与铆钉轴线垂直且与法兰端面平行的方向为Z向，与铆钉轴线平行的方向为Y向，通过计算获得两个方向上刚度数据：Y方向和Z方向的力和位移数据；

c)铰链波纹管的等效采用混合遗传算法，采用基于前次搜索得到的最优点的空间渐缩方法不断调整遗传算法的搜索空间，最终得到两根等刚度的变截面梁；

d)通过约束方程将两根变截面梁在自由端进行耦合，并通过建立辅助圆盘实现变截面梁

根部与各自法兰的连接，取代铰链波纹管，最终将等效的变截面梁应用在管路系统中。

2.根据权利要求所述1的设计方法，其特征在于：变截面梁中间部位的矩形截面的尺寸为根部和自由端的截面尺寸的线性过度，其中根部截面的宽度为B0，高度为H0，自由端截面的宽度为B1，高度为H1；双线性材料模型需要参数为E，K，σ_s，其中E为材料弹性模量，K为屈服后应力应变增量之比，σ_s为材料屈服应力，最终遗传算法优化的7个参数为E、K、σ_s、B0、H0、B1和H1，得到求解变截面梁自由端横向刚度计算流程。

3.根据权利要求2所述的设计方法，其特征在于：变截面梁自由端横向刚度计算的流程如下

1)输入变截面梁的几何参数：B0、B1、H0、H1和L，材料参数：E、K和σ_s，输入外力F；

2)计算结构的弹性极限载荷Fe，弹性极限挠度Ue，判断如果F≤Fe，输出U＝Ue*F/Fe，停止；否则进行下一步；

3)y为中性层挠度，x为水平坐标系，与中性层平行，计算弹塑性区的边界，L1，L2；

4)对于[0，L1]弹性段，带入边界条件 ${\frac{\mathrm{dy}}{\mathrm{dx}}|}_{x=0}=0,$ y|x＝0＝0，计算y(L1)，y′(L1)；

5)对于[L1，L2]弹塑性段，带入边界条件 ${\frac{\mathrm{dy}}{\mathrm{dx}}|}_{x=L1}=y\&prime;\left(L1\right),$ 计算y(L2)，y′(L2)；

6)对于[L2，L]弹性段，带入边界条件 ${\frac{\mathrm{dy}}{\mathrm{dx}}|}_{x=L2}=y\&prime;\left(L2\right),$ y|_x＝L2＝y(L2)，计算y(L)，y′(L)；

7)输出自由端挠度U＝y(L)。

4.根据权利要求3所述的设计方法，其特征在于：

a)对于Fe的确定方法如下

首先求解一元三次方程

At³+BBt²+Ct+D＝0           (1)

其中 $\left\{\begin{array}{c}A={2\mathrm{BH}}^2\\ \mathrm{BB}=-2\mathrm{BH}0H-H^{2}B0-{3\mathrm{BH}}^2\\ C=4\mathrm{BH}0H+{2H}^{2}B0\\ D=-{\mathrm{BH}0}^2-2\mathrm{HB}0H0+B0{H0}^2\\ B=B0-B1,H=H0-H1\end{array}\right.$

该方程确定变截面梁最大应变的位置；由一元三次方程至少有一个实根，对于方程(1)求解出的结果，只取[0，1]范围内的解，如果求解出的实根在[0，1]之外，则真实的解为0；假设求解处最大值应变位置为t_max，令x_max＝t_maxL，最终得极限弹性载荷为

$F_e=\frac{{\mathrm{\&epsiv;}}_s\mathrm{EB}\left(x_{\max }\right)H^2\left(x_{\max }\right)}{6(L-x_{\max })}---\left(2\right)$

b)对于F<Fe和F≥Fe，挠度计算方法如下：

i)横向力F<Fe，变截面梁处于弹性变形阶段

令 $d0=\frac{B0}{B0-B1},$ $d1=\frac{H0}{H0-H1},$ $G=\frac{-12L}{E(B1-B0){(H1-H0)}^3},$ 归一化令 $t=\frac{x}{L};$

如果d0≠d1，则得截面转角和中性层挠度计算公式

$\frac{\mathrm{dy}}{\mathrm{dx}}\frac{1}{\mathrm{FGL}}=\frac{(1-d0)}{{(d0-d1)}^3}\ln \left(\frac{t-d1}{t-d0}\right)+\frac{d0-1}{{(d0-d1)}^2(t-d1)}+\frac{d1-1}{2(d0-d1){(t-d1)}^2}+C0---\left(3\right)$

$\frac{y}{{\mathrm{FGL}}^2}=\frac{(1-d0)}{{(d0-d1)}^3}\{(t-d1)\ln |t-d1|-(t-d0)\ln |t-d0|\}$

                                           (4)

$+\frac{d0-1}{{(d0-d1)}^2}\ln |t-d1|-\frac{d1-1}{2(d0-d1)(t-d1)}+C0t+C1$

如果d0＝d1，则有

$\frac{\mathrm{dy}}{\mathrm{dx}}=-\mathrm{GFL}[\frac{1}{2{(d0-t)}^2}+\frac{1-d0}{3{(d0-t)}^3}]+C0---\left(5\right)$

$y=-{\mathrm{GFL}}^2[\frac{1}{2(d0-t)}+\frac{1-d0}{6{(d0-t)}^2}]+C0\mathrm{Lt}+C1---\left(6\right)$

其中C0，C1为积分常数，通过边界条件确定，将(2)式带入(4)或(6)，并令t＝1求得弹性极限挠度Ue；

ii)横向力F≥Fe，变截面梁处于弹塑性变形阶段

需要求解如下微分方程

$\frac{d^{2}y}{{\mathrm{dx}}^2}=\frac{3F(L-x)-3B\left(x\right)({\mathrm{\&sigma;}}_s-{\mathrm{K\&epsiv;}}_s)(H^2\left(x\right)/4-a^2)}{2B\left(x\right)[(E-K)a^3+{\mathrm{KH}}^3\left(x\right)/8]}---\left(7\right)$

其中 ${\hat{z}}_s=\frac{\mathrm{EI}\left(x\right)}{F(L-x)}{\mathrm{\&epsiv;}}_s,$ a通过下式求解

$a^3+a(\frac{3F}{{\mathrm{\&epsiv;}}_s(E-K)}\frac{(L-x)}{B\left(x\right)}-{3h}^2)+{2h}^3\frac{K}{K-E}=0---\left(8\right)$

方程(7)的求解采用如下数值方法

1)设有边界条件 ${\frac{\mathrm{dy}}{\mathrm{dx}}|}_{x=L1}=C2,$ y|_x＝L1＝C1，将区间[L1，L2]等分为n份，即[x_i，x_i+1]，i＝1，...，n，x_i＝L1+(i-1)(L2-L1)/n，其中n自定，分别得到(7)式在各个区间上的数值积分 ${\mathrm{INT}2}_i=C2+{\&Integral;}_{x_i}^{x_{i+1}}\frac{d^{2}y}{{\mathrm{dx}}^2}\mathrm{dx},$ i＝1，...，n，积分方法任选，比如辛普森求积Simpson，高斯积分Gauss，龙贝格积分Romberg等，于是得到 ${\&Integral;}_{L1}^{L2}\frac{d^{2}y}{{\mathrm{dx}}^2}\mathrm{dx}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{INT}2}_i={\frac{\mathrm{dy}}{\mathrm{dx}}|}_{x=L2};$

2)令INT2₀＝C2，可得

在n+1个点处的值，即INT2_i，i＝0，...，n，这样便可插值得到

的表达式，这里采用三次样条插值，插值得到n段三次多项式函数

三次多项式函数可以直接积分，得到n段函数

在各自区间上的积分 ${\mathrm{INT}1}_i={\&Integral;}_{x_i}^{x_{i+1}}{f}_i^{\&prime;}\mathrm{dx},$ 于是得到 ${\&Integral;}_{L1}^{L2}\frac{dy}{\mathrm{dx}}\mathrm{dx}=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{\mathrm{INT}1}_i={y|}_{x=L2};$

c)对于L1，L2的确定方法如下：

令 $Q=\frac{6\mathrm{FL}}{{\mathrm{\&epsiv;}}_s{\mathrm{EBH}}^2},$ 求解如下方程

t³+(-2d1-d0)t²+(d1²+2d0d1-Q)t+Q-d0d1²＝0，0≤t≤1          (9)

如果(9)式只求解出一个[0，1]间的实根t2，则必有t1＝0；如果求解出2个或3个实根，则只取在[0，1]间的根t1，t2，t1≤t2；则L1＝t1*L，L2＝t2*L。

5.根据权利要求1所述的设计方法，其特征在于：步骤c)中混合遗传算法采用MATLAB遗传算法工具箱中的ga函数，设获得固定部分的刚度数据FY-UY和FZ-UZ后，从中挑选出在Y或Z方向的m个特征力Fi，i＝1，...，m作用下的位移分别为Uri，令列向量x＝(E，K，σ_s，B0，H0，B1，H1)^T，则变截面梁等效问题对应的优化问题表述为

$\min \mathrm{Res}\left(x\right)=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(\mathrm{URi}-\mathrm{Ui})}^2=\underset{i=1}{\overset{m}{\mathrm{\&Sigma;}}}{(\mathrm{URi}-U(F_i,x^T))}^2---\left(10\right)$

s.t.xLB≤x≤xUB

其中[xLB，xUB]为x的取值范围，第一次搜索范围确定后，以后各次搜索范围通过空间渐缩方法确定。

6.根据权利要求5或1所述的设计方法，其特征在于，空间渐缩方法如下：

首先指定一个较大的搜索空间S_i＝{x|xLB_i≤x≤xUB_i}，i＝0，该空间应该尽量大以包含理论上的全局最小值点，随后利用混合遗传算法在S_i内搜索，算法返回最优解x_i ^*及误差Res_i，如果Res_i大于设计的误差容限且x_i ^*优于上次遗传算法返回的解，则新空间S_i+i是以x_i ^*为中心，各维长度是原来的Scal倍，其中Scal<1，即

S_i+1＝{x|xLB_i+1≤x≤xUB_i+1}

其中 $\left\{\begin{array}{c}{\mathrm{xLB}}_{i+1}=x_i^{*}-\mathrm{Scalg}({\mathrm{xUB}}_i-{\mathrm{xLB}}_i)/2\\ {\mathrm{xUB}}_{i+1}=x_i^{*}+\mathrm{Scalg}({\mathrm{xUB}}_i-{\mathrm{xLB}}_i)/2\end{array}\right.---\left(11\right)$

7.根据权利要求1所述的设计方法，其特征在于，步骤d)中：

a)变截面梁的耦合：采用有限元方法需要将变截面梁划分成实体单元，而两个变截面梁在自由端需要耦合以模拟铰链波纹管上法兰盘之间的类似万向节的转动结构；设两根等刚度的变截面梁中一个变截面梁自由端的4个角节点编号分别为Ni，i＝1，2，3，4，另一个变截面梁自由端的4个角节点编号分别为Ni，i＝5，6，7，8，则建立耦合的需要建立如下的约束方程：

$\left\{\begin{array}{c}0=\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{\mathrm{ux}}_{\mathrm{Ni}}-\underset{i=5}{\overset{8}{\mathrm{\&Sigma;}}}{\mathrm{ux}}_{\mathrm{Ni}}\\ 0=\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{\mathrm{uy}}_{\mathrm{Ni}}-\underset{i=5}{\overset{8}{\mathrm{\&Sigma;}}}{\mathrm{uy}}_{\mathrm{Ni}}\\ 0=\underset{i=1}{\overset{4}{\mathrm{\&Sigma;}}}{\mathrm{uz}}_{\mathrm{Ni}}-\underset{i=5}{\overset{8}{\mathrm{\&Sigma;}}}{\mathrm{uz}}_{\mathrm{Ni}}\end{array}\right.---\left(12\right)$

其中ux_Ni，uy_Ni，uz_Ni，分别表示Ni节点的x，y，z方向的位移；

b)变截面梁与法兰的连接：在有限元软件中，对于整个管路的几何模型，在两个变截面梁的根部分别建立两个辅助圆盘与相应的法兰连接，辅助圆盘的直径与法兰的内径相同，并将辅助圆盘与法兰的内缘固连，辅助圆盘与变截面梁的根部贴合并固连，辅助圆盘与变截面梁同轴。辅助圆盘采用弹性材料模型，材料的弹性模量为管路系统管壁材料弹性模量的10倍，辅助圆盘厚度与管壁厚度相同，泊松比和线膨胀系数与管壁材料相同。

Images