离散微分几何在膜结构设计中的运用
发布时间:2020年4月11日 点击数:3513
0 引言
离散微分几何是一门新兴的数学分支, 多年来众多学者致力于该领域的研究。在工程设计中, 一般采用有限元思想, 将连续的结构离散成为有限个单元。特别地, 在膜结构设计中, 将膜曲面离散为三角形或四边形单元, 将索离散为直线或抛物面单元。在本文的讨论中, 分别为三角形单元和直线单元。
膜结构设计步骤包括形态分析、荷载分析、裁剪分析和成形分析, 它们构成了结构分析的全过程。在膜结构设计中, 多个分析过程需要求解曲面的曲面曲率, 法向量, 主方向
1 非线性有限元法
在膜结构的形态分析和荷载分析中, 一般采用非线性有限元法, 形态分析求解一定的膜索初应力和边界条件下达到的平衡形态, 而荷载分析求解膜结构在预应力和外荷载作用下达到的平衡形态。从初始状态到最终的几何形状膜材在预应力和外荷载的作用下发生位移, 结构势必要发生很大的变形, 体现出很强的几何非线性。通常都采用非线性有限元的增量法来进行非线性分析, 本文采用UL格式, 所有变量以时间T位形作为参考位形。
膜结构中, 膜采用常应变平面三角形单元, 索采用两节点直线单元。索只考虑其轴力, 位移可以很大, 并且是空间的, 索和膜之间无相对位移。为了描述单元运动, 需建立位移模式, 即单元内任一点位移与节点位移间的关系。
其中, h
考虑几何非线性, 在应变和位移的关系中计及二次项:
材料均为线弹性材料, 应力应变关系为:
式中, σ0为预应力;D为材料本构矩阵。
由虚功原理和变分原理得, 在t+Δt时刻局部坐标系下膜单元的平衡方程为:
其中:ttKL结构线性应变增量刚度矩阵:
ttKNL结构非线性应变增量刚度矩阵:
t+ΔtR荷载列向量, 初始形态分析只为内压, 荷载分析包括内压和其他荷载;
ttF:T时刻单元应力的等效节点力向量。
式中, ttB
在局部坐标下求得每个单元的平衡方程后, 进行坐标转化和整体叠加, 得到整个结构整体坐标下的平衡方程:
采用Newton-Raphson法求解上述得到的非线性方程组, 计算收敛准则采用不平衡力准则。计算收敛后得到的形态即为外荷载作用下的平衡形态。
2 离散曲面曲率求解
由于膜结构设计的曲面一般为复杂曲面, 在形态和荷载分析后得到的曲面由一系列离散点表示, 没有显式的曲面方程, 用曲面拟合的方法计算量大且引入误差。基于分析得到的离散曲面求得膜曲面的各种曲率 (高斯曲率, 平均曲率, 主曲率) 、主方向矢量和法矢量, 能够方便地跟踪膜曲面在预应力和外荷载作用下的曲率变化情况。
首先定义三角形离散曲面中任意点的Voronoi面积。
如果三角形i (j-1) j为非钝角三角形, 点O为三角形的外接圆心, 如图1所示。
在三角形中的Voronoi面积为:
如果三角形xixj-1xj为钝角三角形, 点O为钝角对应边的中点, 若∠xi为钝角, 点xi在三角形中的Voronoi面积为:
若环绕点xi的三角形都是非钝角三角形, 那么点xi的Voronoi面积为:
其中, αj, βj为边xixj所对的两个角;‖xi-xj‖为边ij的长度, Ne (i) 为所有和点xi相连点的集合;
可得点xi的平均曲率:
高斯曲率:
其中, n为点xi的曲面法向量, 得点xi的切平面方程为:
由平均曲率, 高斯曲率可得点xi的主曲率:
其中, Δ=κH2-κG
3 离散曲面测地线求解
将测地线作为膜裁剪分析中的裁剪线。由微分几何知, 测地线上每点的主法线向量与曲面在这点的法线向量平行出发。基于这一性质提出一种新的测地线生成算法, 用测地线在某点的密切平面和所在曲面的交线近似测地线。该算法能够通过简单迭代生成任意两端点间的测地线, 收敛稳定, 易于实现。
算法首先确定测地线上各点的密切平面, 测地线上某点的主法线向量与曲面在这点的法线向量平行, 而某点处的曲面法线向量N为:
式中:m为共有该点的网格平面总数;Nitria为第i个网格平面的法线向量。
测地线在各点的切线向量:
当网格划分足够密时这种近似是合理的, 引入的误差很小且处理方便。
于是过点A的密切平面的法线向量Np:
由式 (16) 求出测地线和边线 (图2中边线MN) 的交点B, 连接AB即为密切平面和初始平衡曲面的交线。
按上述方法确定点B处的切线向量和主法线向量即可确定密切平面, 和初始平衡曲面相交得到测地线的下一个延伸段BC。以此类推不断地确定密切平面及其和初始平衡曲面的交线, 直至曲面边界。
已知初始点及该点切线向量运用上述方法就可求得测地线轨迹。测地线初始点A, 初始切线向量Ti, 对应的测地线终点Ai′有以下关系:
式中, 映射关系F为测地线延伸算法, 并规定函数值:两端点确定的竖直面P把三维空间分为两个子空间, 在左子空间F (A, Ti) -Ai′<0, 在右子空间F (A, Ti) -Ai′>0。
已知两端点A和A′求测地线的问题, 转化为搜索点Ai′逼近点A′的初始切线向量T′, 再由点A和切向量T′得到最终的测地线。本文采用二分迭代法搜索切线向量Ti, 当满足收敛准则, 迭代过程结束, T′即为所求的初始切线向量。
二分迭代法的实施步骤如下:
(1) 在点A处的测地线从切平面内, 设切向量T的搜索区间为[D0, E0];
(2) 取D0和E0的角平分向量T0, 运用密切平面相交法得到F (A, T0) 。若F (A, T0) -C=0, 则T0确定的切向量就是待求的切向量;
(3) 否则检查不等式 (F (A, E0) -C) · (F (A, T0) -C) <0是否成立。如果成立, 则T0在右子空间, 取D1=D0, E1=T0;否则取D1=T0, E1=E0;
(4) 对收缩的求解区间[D1, E1], 取D1和E1的角平分向量T1, 返回②进行搜索。
由上述算法即可求出任意给定两端点间的测地线。
4 算例
算例1:某膜结构看台, 取其1/6为例说明膜曲面曲率设计。分别取:膜面应力2.5kN/m, 边索30kN, 面内索2kN进行形态分析, 得到的初始平衡形态见图3, 曲面离散点的法矢、主方向见图4。
4-10-1, 10-11-12, 3-12-6是索段, 四周均为固定边界, 6个控制点的坐标分别为:1 (207.72, -3.88, 9.15) , 2 (217.64, -4.06, 7.98) , 3 (225.29, -4.21, 0.14) , 4 (207.72, 3.88, 9.15) , 5 (217.64, 4.06, 7.98) , 6 (225.29, 4.21, 0.14) 。
张拉膜面曲率的分布情况:膜面高斯曲率为负值, 除局部曲面外, 一般都比较小, 曲拱上膜面点的高斯曲率为正值。膜面的平均曲率为极小正值, 可近似为零, 和薄壳的无矩理论吻合。索上膜面点和拱上点的平均曲率比膜面大几个数量级。
荷载取值:恒载 (DL) 膜重1.05kg/m2, 索重由程序自动计算得到;索膜预应力 (PL) 膜两向预应力2.5kN/m, 边索预拉力取30kN, 膜面压索取2kN;风吸 (WX) , 1245区域取1170N/m2, 2356区域取900N/m2;风压 (WY) 和风吸方向相反;雪载 (SN) , 1245区域取630N/m2, 2356区域取504N/m2。荷载组合工况:①DL+PL;②DL+PL+WX;③DL+PL+WY;④DL+PL+SN。
利用有限元法编制的荷载分析程序进行以上工况分析, 求出各个工况下曲面的高斯曲率, 平均曲率, 下表给出了其中六个点的高斯曲率, 平均曲率。
初始平衡曲面及各工况下曲面的曲率变化表 导出到EXCEL
荷载组合 |
初始平衡曲面 |
工况1 | 工况2 | 工况3 | 工况4 | |||||
Gauss (E-3) |
Mean (E-3) |
Gauss (E-3) |
Mean (E-3) |
Gauss (E-3) |
Mean (E-3) |
Gauss (E-3) |
Mean (E-3) |
Gauss (E-3) |
Mean (E-3) |
|
点7 |
-4.81 | 51.25 | -4.84 | 0.284 | -3.91 | 0.109 | -5.71 | 7.51 | -5.55 | 8.495 |
点8 |
-2.18 | 75.26 | -2.205 | 0.366 | -1.01 | 14.83 | -2.823 | 6.636 | -2.796 | 6.686 |
点9 |
-17.68 | 96.27 | -17.71 | 0.022 | -17.08 | 0.4863 | -17.99 | 3.140 | -18.46 | 2.065 |
点10 |
-3.144 | 257.19 | -3.11 | 7.176 | -7.061 | 37.14 | -2.387 | 35.10 | -2.289 | 36.74 |
点11 |
-4.647 | 379.03 | -4.49 | 19.34 | -12.83 | 79.20 | -3.293 | 26.57 | -3.342 | 25.49 |
点12 |
-11.50 | 189.31 | -11.41 | 2.337 | -12.89 | 4.229 | -11.10 | 15.98 | -10.74 | 4.668 |
算例2:如图5锥形膜结构, 顶部半径是a=0.7 m, 底部固定点的半径是b=3.8 m, 高h=2.4m。形面为中心对称结构, 顶部的边界固定, 底部为五边形的顶点固定, 顶点之间有索连接, 图5所示 (b) 为膜曲面形状的正视图。
考虑结构对称性和美观, 膜曲面上测地线布置如图6, 以测地线为裁剪线, 将整个曲面分为15片。将每一曲面片依次展开, 得到如图7的平面展开图。展开前整个结构曲面的面积为30.535m2;展开后15片曲面片的面积之和为:30.539m2。面积误差很小, 仅为0.02%。
5 结语
数学作为一门基础学科, 其发展必然推动工程设计研究。本文将离散微分几何和膜结构设计相结合, 首先利用非线性有限元法进行膜结构的形态分析和荷载分析, 再基于有限元分析得到的离散模型, 求解曲面曲率, 法向量及测地线, 进行结构分析设计。算例表明此方法快速简便, 避免了曲面拟合带来的计算成本和误差。