膜结构的找形分析是其设计过程中极为基本也是较难解决的问题，等应力膜曲面是其中一类重要的找形曲面，其膜面的找形分析等价于数学上求解极小曲面问题[1]。随着计算机技术的发展，数值求解方法已成为膜结构找形分析的主要方法，常见的主要有三类：力密度法[2－4]、动力松弛法[5－8]和非线性有限元法[9－11]。与力密度法和动力松弛法相比，非线性有限元法是一种更为精确的数值计算方法，不足之处是非线性有限元法需要形成刚度矩阵并求解方程组，网格较密时占用内存较大，计算速度较慢。
有限元法(finite element method，FEM)自适应分析是数值计算方法研究的热点，自适应有限元法能够根据用户设定的误差限对网格自适应地进行调整，最终得到优化的网格，并给出满足误差限的解答。对于膜结构极小曲面找形分析问题，当用户对解答有一定的精度要求时，相比于传统有限元“盲目”地细分网格，自适应有限元能够针对性地仅在需要调整的区域进行调整，在保证求解质量的同时，大大提高求解的效率，这在一定程度上能够缓解了有限元法相对占用内存较大，计算速度较慢的问题。
本文基于单元能量投影(element energy projection，EEP)法，提出了一种二维非线性有限元自适应方法，可用于求解膜结构极小曲面找形分析问题。EEP法是袁驷教授等从结构力学的矩阵位移法得到启发[12]，通过在单元内创造性地运用有限元基本数学理论中的投影定理[13]，提出的一种后处理超收敛算法，已在一维有限元中得到充分运用和系统的开发[14－19]，并被成功推广到二维[20－22]、三维[23]问题中去。该法简便直观，行之有效，充分利用其优良特性，用EEP超收敛解代替真解对有限元解逐点按最大模进行误差估计，并自适应地进行网格细分，从而形成了一整套基于 EEP法的有限元自适应分析方法，已在一维有限元及二维线性有限元应用中获得了一致的成功。
基于 EEP法的一维非线性有限元自适应分析已经取得成功[16－18]，如何利用二维线性有限元的EEP超收敛解和自适应方法，实现技术的转移，将其应用到二维非线性自适应求解中去、将非线性有限元的自适应分析从一维问题拓展到二维问题是一个值得研究的课题。本文对此进行研究，并取得了良好进展，成功将之应用于求解膜结构极小曲面找形分析问题。下文将对所提出的方法进行介绍，包括方法的基本思想，相应的计算公式以及算法。文末给出的数值算例表明，本法可靠，效果颇佳。
直角坐标系下，一膜曲面u(x,y)面积可表示为关于u(x,y)的泛函π(u)：
其中下标x,y表示对其求偏导。求膜面面积最小值，对式(1)取一阶变分，由其驻值条件可导出相应的Euler方程：
其中，
；N为非线性算子。此外，还应满足给定的边界条件：
这里，
为给定的边界值。
式(2)为一非线性偏微分方程，难以直接对其求解，通常解法是将其线性化，并进行迭代求解。非线性迭代方法众多，二维非线性问题也存在着多种线性化路径。Newton法是一种广为采用的非线性迭代方法，具有局部二阶收敛特性。本文采用基于弱形式的Newton法将原非线性模型问题线性化，与此同时，建立起相应的Newton迭代格式；思路自然直接、简明清晰，推导过程如下。
首先写出与式(2)对应的Galerkin法弱形式为：
其中，v为检验函数。
以下统一规定，字母下标k均表示该量非线性迭代第k步的值。令
，将F在
处作Taylor级数展开，舍去二阶及以上的微量有：
将式(4)代入式(3)，定义双线性型和线性型为：
整理得：
上式中线性化参数具体表达式如下：
由式(7)可以导出第k步 后关于u的线性微分方程：
边界条件为：
这样，给定第k步的解答，由式(7)计算线性化参数，求解式(8)可得第k+1步的解答，反复迭代直至收敛。值得一提的是，式(7)中的线性化荷载项
包含uk的二阶导数，可以通过分部积分避免，这也是本文考虑基于弱形式推导的原因之一。
考虑式(6)的有限元求解。本小节只对单元模式和网格划分做一简要介绍，其他方面遵循通常的有限元方法。本文采用双m次四边形单元，单元上的试函数借助双m次形函数对结点位移插值得到：
其中：Ni和Nj是m次 Lagrange或 Hierarchical形函数；dij为相应的结点位移。单元上的检验函数vh采用同uh相同的形函数插值形式。记
为满足位移边界条件、直至一阶导数均平方可积的函数空间，
为通常的有限元试探函数空间，有限元的求解归结于求uh∈Sh使其满足如下的虚功方程：
由于二维线性有限元的 EEP超收敛计算要求有限元网格是“拟线法网格”，本文的有限元网格也须是“拟线法网格”：先利用有限元线法(FEMOL)的离散方式用一组结线对求解区域进行离散，然后再对其结线方向进一步离散，得到的网格既是有限元线法常微分方程组(ordinary differential equations，ODEs)的广义一维有限元网格，也是原问题的二维有限元网格，如图1所示。
采用本小节的有限元模式在当前网格上按上述的Newton迭代格式做非线性迭代，直至收敛，得到当前网格上的非线性有限元解。
图1 二维有限元网格逐维离散示意图
Fig. 1 D-by-D discretization of 2D problems
上节得到非线性有限元解是当前网格下最接近原非线性问题的有限元解，如果不细分网格，继续迭代下去并不能改进当前的有限元解，亦即得到了一个等价线性问题[18]。此时，可以引入基于EEP法的二维线性有限元超收敛计算和自适应技术，对等价线性问题的有限元解进行误差估计和网格细分，再在新的网格下重新进行非线性迭代求解，如此反复直至满足用户设定的误差限。随着网格加密和重复非线性迭代，在设定的误差限意义下，等价线性问题最终逼近到原非线性问题，二者同解，继续迭代或细分网格无益，求解结束。
这里仅简单介绍二维线性有限元EEP超收敛计算思路，并针对本文模型问题简要给出必要的表达式。基于EEP法的二维线性有限元超收敛计算采用的是“逐维恢复”的方案，借助有限元线法[24]作为桥梁，分两步利用单元能量投影定理导出单元上任一点的超收敛计算公式，先是对线法ODEs问题的广义一维有限元解进行超收敛计算得到“拟线法解”[21]，用之代替线法解，继而将问题看作是有限元线法二维超收敛计算问题[20]，套用线法的 EEP超收敛公式对全域进行精度恢复，得到二维有限元的超收敛解。EEP超收敛计算思路简要总结如下：
1) 二维有限元解。求得有限元解，并认定图1(c)中某方向为“拟线法”方向。
2) 拟线法解。将有限元解看作是“拟结线”方向线法ODEs的广义一维有限元解，采用线法ODEs的EEP超收敛计算公式(11)求得“拟结线”上的超收敛解，作为“拟线法”解：
3) 全域超收敛解。将“拟线法”解看作是线法解，采用式(12)求得全域的二维有限元超收敛解：
自适应求解的两大关键：可靠的误差估计和高效的网格细分技术。EEP超收敛解相比于有限元解具有逐点超收敛性，可将其替代真解对有限元解按最大模进行误差估计。有了可靠的误差估计后，便可以进一步考虑网格细分的问题。基于单元边线解答的均差法[22]被大量算例验证了是一种高效稳健的网格细分技术，该法生成的网格分布合理、误差冗余度小。对于二维有限元，基于单元边线解答的均差法主要包含两个步骤，即选择加密方向和确定加密位置，具体实施方法如下。
1) 选择加密方向。遍历所有有限元单元，若某单元误差超限，找到误差最大的边，确定并标记该单元二分的方向。遍历完所有单元后，统计所有待加密的单元及加密方向，以便下一步网格细分。
2) 确定加密位置。对于二维有限元网格的细分方式采用通条插线的方式，即对误差超限的单元所在的行(或列)单元通条插入网格线进行二分，二分的位置按照误差的平方在二分后的两个新单元上均等的原则来确定。
对于非线性问题的自适应求解，其停机准则为同时满足以下两个条件：1) 非线性迭代收敛；2) 网格无需再细分。基于以上的分析讨论，可以大致将基于 EEP法的二维非线性有限元的自适应求解总结为“三步走”策略，如下所述。
1) 有限元解。在新网格上(初始网格人为给定)进行有限元迭代计算直至收敛，得到非线性有限元解。
2) EEP解。计算上一步得到的有限元解的EEP超收敛解。
3) 网格细分。用 EEP超收敛解替代真解对有限元解进行误差估计，并采用基于单元边线解答的均差法将单元细分得到新的网格，返回第一步；如果所有单元均满足设定的误差限则停机。
本文算法已经编制成 Fortran90程序。为了检验算法有效性，本节的数值算例均有解析解可比较，统一采用最大模绝对误差进行度量，设定非线性有限元迭代收敛的误差限为10-8，自适应误差限tol=10-4，并定义误差比为(u-uh)/tol，全域误差比在[-1, 1]表明自适应求解圆满完成。初始网格均为 1个单元。初始解通过求解满足边界条件式(8b)的Poisson方程获得，这在程序实施上十分简便，只需将式(8a)中
赋值为零。
例1. 悬链面
悬链面的解析表达式为[25]：
本例取求解域为
x,y∈R+,a=0.1}，如图2所示，边界AB、CD采用精确的圆弧曲线并给定位移值，AB、CD为自由边。计算结果如表1所示，其中m为单元次数，自适应最终网格的单元数用两个方向的单元数相乘来表示，下同。以二次元为例，最终得到自适应有限元网格划分如图2所示，其有限元近似解如图3所示，同精确解比其全域误差比如图4所示。靠近圆弧AB区域单元较密集且单元径向尺寸都很小，最小尺寸仅有 3.81× 1 0-7，这揭示了该问题的特性所在，不难从解析解看出沿径向的导数值在圆弧AB处为无穷大，使得在该点及其邻域的求解难度较大，数值实验显示自适应过程不断地在该区域细分网格直至有限元解满足设定的误差限。
例2. Scherk第一曲面
案例I：Scherk第一曲面(first surface)的解析表达式为[25]：
本例取求解域为
边界条件给定如图5所示。计算结果如表2所示，以二次元为例，最终得到自适应有限元网格划分如图5所示，其有限元近似解如图6所示，同精确解比其全域误差比如图7所示。
表1 悬链面计算结果
Table 1 Calculation results of Catenoid
单元次数 最终网格 自适应步数maxhu-utolm=2 1×21 10 0.85m=3 1×11 6 0.84m=4 1×9 5 0.69
图2m=2，悬链面最终自适应网格划分
Fig. 2m=2, Final mesh of Catenoid
图3m=2，悬链面自适应有限元解示意图
Fig. 3m=2, The results of Catenoid
图4m=2，悬链面自适应有限元解误差比图
Fig. 4m=2,Erroru-uhof Catenoid
表2 Scherk第一曲面计算结果
Table 2 Calculation results of Scherk’s first surface
单元次数 最终网格 自适应步数maxhu-utol m=2 17×16 8 0.96m=3 7×8 6 0.98m=4 4×4 4 0.29
图5m=2，Scherk第一曲面最终自适应网格划分
Fig. 5m=2, Final mesh of Scherk’s first surface
图6m=2, Scherk第一曲面自适应有限元解示意图
Fig. 6m=2, The results of Scherk’s first surface
图7m=2, Scherk第一曲面自适应有限元解误差比图
Fig. 7m=2, Erroru-uhof Scherk’s first surface
案例II：为了验证算法对于更广泛几何区域的适用性，就Scherk第一曲面问题构造另一非规则区域算例，图8给出求解域示意图和4个顶点的坐标，四边给定位移边界条件。计算结果如表3所示，以二次元为例，最终得到自适应有限元网格划分如图8所示，其有限元近似解如图9所示，同精确解比其全域误差比如图10所示。
表3 Scherk第一曲面(非规则区域)计算结果
Table 3 Calculation results of Scherk’s first surface(Irregular domain)
单元次数 最终网格 自适应步数maxhu-utolm=2 21×30 7 0.93m=3 8×12 6 0.69m=4 5×6 3 0.64
图8m=2, Scherk第一曲面(不规则区域)最终自适应网格划分
Fig. 8m=2, Final mesh of Scherk’s first surface(Irregular domain)
例3.Scherk第五曲面
Scherk第五曲面的解析表达式为[25]：
本例取求解域为
四边给定位移边界条件。计算结果如表4所示，以二次元为例，最终得到自适应有限元网格划分如图11所示，其有限元近似解如图12所示，同精确解比其全域误差比如图13所示。
图9m=2, Scherk第一曲面(不规则区域)自适应有限元解示意图
Fig. 9m=2, The results of Scherk’s first surface(Irregular domain)
表4 Scherk第五曲面计算结果
Table 4 Calculation results of Scherk’s fifth surface
单元次数 最终网格 自适应步数maxhu-utolm=2 16×20 7 0.72m=3 8×8 5 0.36m=4 6×6 5 0.25
图10m=2, Scherk第一曲面(不规则区域)自适应有限元解误差比图
Fig. 10m=2, Erroru-uhof Scherk’s first surface(Irregular domain)
图11m=2，Scherk第五曲面最终自适应网格划分
Fig. 11m=2, Final mesh of Scherk’s fifth surface
图12m=2，Scherk第五曲面自适应有限元解示意图
Fig. 12m=2, The results of Scherk’s fifth surface
本文提出的基于 EEP法的二维非线性有限元自适应方法思路清晰。算法简明，将非线性迭代和线性问题自适应算法有机结合起来，其核心是充分利用EEP超收敛解的优良特性。文中给出的三个算例均为经典的极小曲面问题，从数值试验结果来看，可以大致总结如下：首先是基于弱形式的Newton法迭代的求解效率较高，就目前的数值试验而言迭代次数均在6以内；其次是本文所提出的自适应算法的可靠性很高，最终的有限元解与精确解相比最大模误差均在设定误差限内，且误差分布较为均匀，冗余度较小。这表明本文算法的适用性、可靠性以及求解效率和解答质量都颇具竞争力，对于膜结构的极小曲面找形分析是一种高性能的自适应数值算法。
图13m=2，Scherk第五曲面自适应有限元解误差比图
Fig.13m=2,Erroru-uhof Scherk’s fifth surface
[1]Hildebrandt S, Tromba A. Mathematics and optimal form[J]. College Mathematics Journal, 1985, 18(1): 84.
[2]Schek H J. The force density method for form finding and computation of general networks [J]. Computer Methods in Applied Mechanics ＆ Engineering, 1974,3(1): 115―134.
[3]Gründig L. Minimal surfaces for finding forms of structural membranes [J]. Computers ＆ Structures, 1988,30(3): 679―683.
[4]Maurin B, Motro R. The surface stress density method as a form-finding tool for tensile membranes [J].Engineering Structures, 1998, 20(8): 712―719.
[5]Day A S. An introduction to dynamic relaxation [J]. The Engineer, 1965, 29(1): 218―221.
[6]Topping B H V. The application of dynamic relaxation to the design of modular space structures [D]. London: City University London, 1978.
[7]Lewis W J, Jones M S, Rushton K R. Dynamic relaxation analysis of the non-linear static response of pretensioned cable roofs [J]. Computers ＆ Structures, 1984, 18(6):989―997.
[8]Barnes M. Form stress engineering of tension structures[J]. Structural Engineering Review, 1994, 6(3): 175―202.
[9]Haug E, Powell G H. Finite element analysis of nonlinear membrane structures [C]// Tension and Space Structures(v.2): Proceedings of the 1971 IASS Pacific Symposium.Tokyo and Kyoto, 1972: 165―175.
[10]Argyris J H, Balmer H, Kleiber M, et al. Natural description of large inelastic deformations for shells of arbitrary shape-application of trump element [J].Computer Methods in Applied Mechanics ＆ Engineering, 1980, 22(3): 361―389.
[11]Nishimura T, Tosaka N, Honma T. Membrane structure analysis using the finite element technique [C]. IASS Symposium, Osaka, 1986: 9―16.
[12]袁驷. 从矩阵位移法看有限元应力精度的损失与恢复[J]. 力学与实践，1998, 20(4): 1―6.Yuan Si. The loss and recovery of stress accuracy in FEM as seen from matrix displacement method [J].Mechanics in Engineering, 1998, 20(4): 1―6. (in Chinese)
[13]Strang G, Fix G J. An analysis of the finite element method [M]. Englewood Cliffs, NJ: Prentice-Hall, 1973.
[14]袁驷, 王枚. 一维有限元后处理超收敛解答计算的EEP法[J]. 工程力学, 2004, 21(2): 1―9.Yuan Si, Wang Mei. An element-energy-projection method for post-computation of super-convergent solutions in one-dimensional FEM [J]. Engineering Mechanics, 2004, 21(2): 1―9. (in Chinese)
[15]袁驷, 王旭, 邢沁妍, 等. 具有最佳超收敛阶的EEP法计算格式: I算法公式[J]. 工程力学, 2007, 24(10): 1―5.Yuan Si, Wang Xu, Xing Qinyan, et al. Ascheme with optimal order of super-convergence based on the element energy projection method – I formulation [J].Engineering Mechanics, 2007, 24(10): 1―5. (in Chinese)
[16]袁驷, 杜炎, 邢沁妍, 等. 一维 EEP自适应技术新进展: 从线性到非线性[J]. 工程力学, 2012, 29(增刊Ⅱ):1―8.Yuan Si, Du Yan, Xing Qinyan, et al. New progress in self-adaptive analysis of 1D problems: from linear to nonlinear [J]. Engineering Mechanics, 2012, 29 (Suppl II): 1―8. (in Chinese)
[17]Yuan Si, Du Yan, Xing Qinyan, et al. Self-adaptive one-dimensional nonlinear finite element method based on element energy projection method [J]. Applied Mathematics and Mechanics (English Edition), 2014,35(10): 1223―1232.
[18]杜炎. 基于EEP法的一维非线性有限元法自适应分析[D]. 北京: 清华大学, 2012.Du Yan. Adaptive analysis of 1D nonlinear FEM based on EEP super-convergent method [D]. Beijing: Tsinghua University, 2012. (in Chinese)
[19]袁驷, 刘泽洲, 邢沁妍. 一维变分不等式问题的自适应有限元分析新探[J]. 工程力学, 2015, 32(7): 11―16.Yuan Si, Liu Zezhou, Xing Qinyan. A new approach to self-adaptive FEM for one-dimensional variational inequality problems [J]. Engineering Mechanics, 2015,32(7): 11―16. (in Chinese)
[20]袁驷, 王枚, 王旭. 二维有限元线法超收敛解答计算的EEP法[J]. 工程力学, 2007, 24(1): 1―10.Yuan Si, Wang Mei, Wang Xu. An element-energy projection method for super-convergence solutions in two-dimensional finite element method of lines [J].Engineering Mechanics, 2007, 24(1): 1―10. (in Chinese)
[21]袁驷, 肖嘉, 叶康生. 线法二阶常微分方程组有限元分析的 EEP超收敛计算[J]. 工程力学, 2009, 26(11):1―9, 22.Yuan Si, Xiao Jia, Ye Kangsheng. EEP super-convergent computation in FEM analysis of FEMOL second order ODEs [J]. Engineering Mechanics, 2009, 26(11): 1―9,22. (in Chinese)
[22]袁驷, 徐俊杰, 叶康生, 等. 二维自适应技术新进展:从有限元线法到有限元法[J]. 工程力学, 2011, 28(增刊II): 1―10.Yuan Si, Xu Junjie, Ye Kangsheng, et al. New progress in self-adaptive analysis of 2D problems: from FEMOL to FEM [J]. Engineering Mechanics, 2011, 28 (Suppl II):1―10. (in Chinese)
[23]袁驷, 吴越, 徐俊杰, 等. 基于EEP法的三维有限元超收敛计算初探[J]. 工程力学, 2016, 33(9): 15―20.Yuan Si, Wu Yue, Xu Junjie, et al. Exploration on super-convergent solutions of 3D FEM based on EEP method [J]. Engineering Mechanics, 2016, 33(9): 15―20. (in Chinese)
[24]Yuan Si. The finite element method of lines [M].Beijing-New York: Science Press, 1993.
[25]Dierkes U, Hildebrandt S, Wohlrab O. Minimal surfaces I[J]. Grundlehren Der Mathematischen Wissenschaften,1992, 295(11): 13―40.
A NEW ADAPTIVE FEM FOR MINIMAL SURFACES FORM-FINDING OF MEMBRANE STRUCTURES