张拉膜结构是一种柔性结构, 需要确定结构的初始平衡形状[1], 即对结构进行找形分析。当膜结构符合膜面等应力条件时, 其膜面找形分析等价于数学上的极小曲面求解[2], 即寻求满足给定边界曲线条件下的面积最小的曲面。

膜结构的极小曲面找形分析是一个高度非线性求解过程, 用数值方法求解时会由于解曲面的表述不唯一, 导致问题求解困难。常用解法如非线性有限元法[3]、动力松弛法[4]、力密度法[5]和URS (Updated Reference Strategy) 法[6]等, 在一定程度上减小了问题求解的难度, 但各自还存在一定的缺点。例如非线性有限元法的求解迭代步长要求比较小, 收敛比较慢;动力松弛法由于初始解曲面的假设和最终曲面差别较大, 导致了收敛速度比较慢;力密度法只是达到了结构静力平衡状态, 曲面上的应力分布不均匀[7];URS法会造成找形后曲面的网格出现畸变[8]。此外, 对于某些特殊问题 (如平面找形问题) , 多数现有数值方法都因解答不唯一存在求解困难。因此, 寻求更加高效、简单、可靠的膜结构极小曲面找形方法是很有必要的。

有限元线法 (FEMOL) 是一种基于常微分方程 (ODE) 求解的半解析方法[9], 用有限元线法进行膜结构的找形分析具有很多优点:膜曲面在结线方向不仅光滑而且具有高阶连续性, 便于进行后续的剪裁分析;非线性常微分方程组的求解由ODE求解器COLSYS[10]自动完成, 无需人为控制求解过程;由于对问题只进行半离散, 求解精度高。

文献[11]通过对曲面面积公式的近似简化, 将原非线性问题近似为一线性问题, 用FEMOL求解该线性问题, 得到极小曲面的高质量的近似解。本文将这一线性解作为极小曲面找形分析的初始解, 从膜曲面面积的表达式出发, 在其变分方程中, 用初始解对面积微元和坐标基矢量夹角进行限定, 有效地排除了解答不唯一的干扰。

本文首先介绍张拉膜结构极小曲面找形分析的变分方程的推导, 然后提出克服解答不唯一的限定方法, 并提出一套基于非线性FEMOL的迭代求解算法。大量数值算例表明, 本法效力很强, 可以求解各类刁难问题, 且收敛迅速, 精度和效率均相当出色。

图1所示为张拉膜曲面单元示意图, 在参数坐标系 (ξ, η) 下, 膜曲面Se的面积[12]为:

其中g1、g2分别为曲面上任意一点P (x, y, z) 的协变基矢量:

式中:(⋅)ξ=∂(⋅)∂ξ(⋅)ξ=∂(⋅)∂ξ, (⋅)η=∂(⋅)∂η(⋅)η=∂(⋅)∂η;e1、e2、e3分别为x、y、z方向的单位矢量。被积函数R为面积微元的Jacobi式, 可表示为:

R=g11⋅g22−g212−−−−−−−−−−−√ (3)R=g11⋅g22-g122(3)

其中

集成所有单元面积后得整体曲面面积A=∑eAe, 而等应力膜曲面的找形分析等价于数学上极小曲面问题:在给定曲面边界坐标 (xˉxˉ, yˉyˉ, zˉzˉ) 的条件下, 寻求一个满足边界条件并使曲面面积A最小的曲面, 即要求:

上式即为极小曲面找形分析的变分方程, 由此可以导出相应的非线性微分方程。

用常规有限元线法[13]对式 (5) 进行离散求解时, 由于解可能不唯一, 会导致问题求解困难[13]。以图2所示的极小曲面求解为例, 在参数坐标系 (ξ, η) 下, 用两个线性FEMOL单元求解, 两个单元之间的公共结线位置待定, 曲面的面积为两个单元面积A1、A2之和。对于图2 (a) 的情况, 若结线按照图示虚线改变位置, 则两个单元面积不变, 整个曲面的面积亦不变, 即两种结线布置方式都是极小曲面的解答。对于图2 (b) 的情况, 若结线如图示改变位置, 两个单元的面积A1和A2发生变化, 但其和保持不变, 整个曲面的面积仍然不变, 即两种结线布置方式也都是极小曲面的解答。在这两种情况中, 结线位置的确定都不是唯一的, 从而导致问题求解困难。

为了克服极小曲面数值解答的不唯一性, 本文提出一种限定求解方案, 其思路如下。分析图2的两个示例可以看出:图2 (a) 中单元的面积保持不变, 但是两个单元的基矢量夹角ϕ发生了变化;而图2 (b) 中两个单元的基矢量夹角ϕ0保持不变, 两个单元的各自面积发生了变化。因此, 同时限定单元的面积和基矢量夹角, 单元之间的结线可以得到唯一的解答。

如前所述, 本文利用文献[11]方法得到一个高质量的近似极小曲面, 将其作为本文非线性迭代的初始解, 用以克服后续迭代中出现的求解困难。设初始解曲面上任意一点为P0 (X, Y, Z) , 则相应的协变基矢量表达式为:

初始解曲面上面积微元的Jacobi式和基矢量夹角的余切分别为:

Rˉˉˉ=|G1×G2|=G11⋅G22−G212−−−−−−−−−−−−√Rˉ=|G1×G2|=G11⋅G22-G122, cotφ=G12/Rˉˉˉ (7)cotφ=G12/Rˉ(7)

其中G11、G22、G12的计算与式 (4) 类似, 只是要用初始解曲面上的坐标 (X, Y, Z) 。本文用初始解曲面面积微元的RˉˉˉRˉ来限定控制方程 (5) 中极小曲面的R, 用初始解曲面的基矢量夹角余切G12/RˉˉˉG12/Rˉ来限定极小曲面的g12/R。这样, 式 (5) 的变分方程修正为:

式 (8) 是本文提出的极小曲面迭代求解的变分方程, 由此导出的控制方程仍为非线性微分方程, 但其非线性程度较弱, 更容易求解。

为了检验式 (8) 代替式 (5) 所带来的误差, 记ε1=R−Rˉˉˉε1=R-Rˉ和ε2=g12-G12。由式 (8) 可以看出Ae00e=∫AeR0dξdη, 其中被积函数R0可取为

R0=g22⋅g11−2G12⋅g12+G22⋅G112Rˉˉˉ (9)R0=g22⋅g11-2G12⋅g12+G22⋅G112Rˉ(9)

对式 (9) 稍加变形, 则不难得到:

R0=R+ε21+ε222Rˉˉˉ (10)R0=R+ε12+ε222Rˉ(10)

可见, R0和R之间的误差为ε1和ε2的二阶微量, 即式 (8) 的解以二阶速度收敛于式 (5) 的解。

根据以上讨论, 可得到一套完整的迭代求解算法如下:

(1) 令k=0, 给定曲面误差限tol, 用文献[11]的方法得到一个高精度的初始解曲面S (0) (X, Y, Z) , 并计算其面积A (0) ;

(2) 令k≡k+1, 用S (k-1) (X, Y, Z) 构建式 (8) 导出的非线性ODEs, 求解得新的曲面S (k) (x, y, z) 并计算其面积A (k) ;

(3) 检验A (k+1) -A (k) ≤tol· (1+A (k) ) 是否满足, 若不满足, 用 (x, y, z) 更新 (X, Y, Z) , 即用当前曲面作为初始曲面, 返回到步2;

(4) 若满足, 得到极小曲面S (k) (x, y, z) , 停止计算。

图3所示为FEMOL的膜单元映射, 其中参数坐标ξ和η的定义域均为[-1, 1]。单元结线上的坐标向量为:

其中p为单元的次数。膜单元内任意一点的坐标 (x, y, z) 采用Lagrange插值得到:

其中

由式 (12) 不难得到坐标x、y和z对ξ、η的各阶偏导数。

对于FEMOL膜单元, 式 (8) 的变分式可写为:

由式 (4) 可以得到:

将式 (15) 代入式 (14) 中, 形式上完成ξ方向的积分, 并在η方向上作必要的分部积分, 注意到{δx}e、{δy}e和{δz}e在η=±1的边界上均为零, 有

其中

式中各系数矩阵分别为

将所有单元集成得到δA0, 又由δA0=0以及{δx}、{δy}和{δz}的任意性, 可以导出一组二阶非线性常微分方程组 (ODEs) :

Ax′′+Bx′+Cx={0}Ay′′+By′+Cy={0}‚ −1<η<1Az′′+Bz′+Cz={0} (20)Ax″+Bx′+Cx={0}Ay″+By′+Cy={0}‚-1<η<1Az″+Bz′+Cz={0}(20)

其中各系数矩阵分别由单元系数矩阵集成。边界条件 (BCs) 为:

其中xˉˉxˉ、yˉyˉ和zˉzˉ为给定的边界值。对于以上导出的非线性常微分方程组, 本文采用ODE求解器COLSYS[10]求解, 求解时简单地将x、y和z的初始解取为X、Y和Z即可顺利求解。

本节给出若干典型的数值算例用以展示本文方法的可靠性和有效性。算例中, COLSYS的误差限统一取为10-5, 膜曲面面积误差限取为10-3, 初始解曲面用文献[11]方法线性求解得到。

Scherk曲面是为数不多的具有解析表达式的极小曲面之一, 其表达式为[14]

z=1k[ln(cos(ky)−ln(cos(kx))]+c (22)z=1k[ln(cos(ky)-ln(cos(kx))]+c(22)

本例中取k=1.0, c=5.0, -5π/12≤x (或y) ≤5π/12, 4条边界曲线给定。本文用一个8次FEMOL元求解, 结线沿y轴方向布置。用文献[11]方法得到的初始解曲面如图4 (a) 所示, 迭代求解1次后, 得到收敛的极小曲面如图4 (b) 所示。

表1给出了坐标z的计算结果与精确解的比较, 可见本文结果精度很高。表2给出不同次数单元的计算结果比较, 可见单元次数越高结果越好。

旋转悬链面亦是有解析表达式的极小曲面。本例求解下述悬链面[13]:

z=a[ln(x2+y2−−−−−−√+x2+y2−a2−−−−−−−−−−√)−lna] (23)z=a[ln(x2+y2+x2+y2-a2)-lna](23)

其中a≤x (或y) ≤5a, 0≤z≤2.29243167a, 并取参数a=10, 四条边界给定。本文采用一个6次FEMOL单元计算, 结线沿z轴方向布置。用文献[11]方法得到的初始解曲面如图5 (a) 所示, 迭代求解1次后, 得到收敛的的悬链面如图5 (b) 所示。

表3给出了坐标z的计算结果与精确解的比较, 表4给出不同单元网格划分下得到的面积比较, 可见, 单元次数越高, 结果越好, 通常一个高次元的精度比两个低次元的精度高。

Scherk-like曲面也是极小曲面找形的经典算例。取边长为1的立方体的边界, 本算例用3个5次FEMOL单元, 结线沿穿洞的方向布置。用文献[11]方法得到的初始解曲面如图6 (a) 所示, 迭代求解4次后, 得到收敛的极小曲面如图6 (b) 所示。本例得到的曲面面积为2.4676, 文献[8]中得到的曲面面积分别为2.4710 (四边形单元) 和2.4704 (三角形单元) 。从极小曲面的属性来看, 本例得到的面积更小, 精度更高, 收敛也快。

如图7所示, 3条边界组成扇形平面区域, 半径r=10, 对应圆心角为π/3的圆弧。用一个4次FEMOL单元求解, 结线沿径向布置, 成为退化边单元。图7 (a) 为初始解曲面, 经过1次迭代求解, 得到收敛的曲面如图7 (b) 所示。本文找形后的曲面面积为52.3602, 精确解为52.3600。本例表明, 本文方法完全适用于平面找形, 即便有退化边单元也无困难。

本例旨在检验本文方法对三角形区域的适用性。考虑锥形曲面:顶点坐标是 (0, 0, 10) , 底边为半径为10的半圆。本文用2个6次FEMOL元求解, 单元对称分布, 结线沿z方向布置。用文献[11]方法得到的初始解曲面如图8 (a) 所示。用该解作为初始解, 迭代求解2次, 得到收敛的极小曲面如图8 (b) 所示。本文找形后的曲面面积为205.54。值得一提的是, 对于本例, ANSYS未能给出收敛的结果。

文献[11]中将极小曲面问题近似为一线性问题, 用有限元线法求解该线性问题, 无需迭代求解, 一步即可得到逼近度很高的初始解。本文对有限元线法分析极小曲面找形问题提出了一种限定求解方案, 有效地克服了解答不唯一导致的求解困难, 进而提出一套基于非线性有限元线法的迭代求解算法。算例表明, 本法效力强、适用广、收敛快、精度高, 是一种可靠、实用、高效的找形方法。这一限定方案有望在有限元法中也能得以百富策略白菜网。