CN105975649A - 一种方形太阳帆表面褶皱形态的建模方法 - Google Patents

Abstract

本发明公开了一种方形太阳帆表面褶皱形态的建模方法，首先，针对方形薄膜太阳帆结构，确定其应力场分布云图，得到褶皱区域；其次，提出了一种新的褶皱形态描述方程，通过引入一个指数衰减函数，建立更接近实际情况的褶皱形态模型。该方法利用Airy应力场分布理论，通过应力叠加准则，并消除边界影响，得到薄膜太阳帆结构的应力场分布云图，进一步确定褶皱区域。在此基础上提出一种更精确的褶皱形态描述方程，基于vonKarman大挠度方程，建立薄膜褶皱参数的一般性方程，采用应力极值法，得到褶皱形态与褶皱方向拉应力的关系模型，求解简单，同时得到的褶皱形态更接近实际情况。

Description

一种方形太阳帆表面褶皱形态的建模方法

技术领域

本发明属于机械工程领域，涉及针对太阳帆结构表面褶皱形态的建模方法，具体是一种方形太阳帆表面褶皱形态的建模方法，用于更准确地描述表面褶皱形态。

背景技术

太阳帆需要通过微弱的太阳光压力获得推进力，所以要求太阳帆具有很大的帆面以及很轻的结构质量，而薄膜材料由于质量轻、收藏体积小、成本低、易折叠等优点，成为理想的帆面材料。然而薄膜的弯曲刚度较小，受到外载荷时容易产生局部不稳定而引起褶皱，由薄膜褶皱引起的光子散射，使得太阳帆产生的有效推力明显降低。因此，准确描述褶皱的变形机理，成为一个热点问题。国内外的学者对此做出大量研究，目前薄膜褶皱的分析主要有两种方法：

一种是基于张力场理论的分析方法。最早是1929年由Wagner提出的，其原理主要是对薄膜的本构关系进行修正，避免薄膜中压应力的出现，结合修正后的本构关系采用数值分析方法进行褶皱分析，得到褶皱形成以后的应力状态。Hongli Ding,Bingen Yang在《International Journal for Numerical Engineering》，2003年10月第12期发表的论文“The modeling and numerical analysis of wrinkled membranes”中提出一种二变参数法，在可变泊松比方法的基础上，同时对弹性模量进行了修正，引入了两个参数，之后L.Zhang,Q.Gao and H.W.Zhang在《International Journal for Numerical Methods inEngineering》，2014年3月第10期发表的论文“Analysis of 2-D bimodular materialsand wrinkled membranes based on the parametric variational principle and co-rotational approach”中给出了进一步研究。采用张力场理论能够判定褶皱区域，但是无法得到褶皱的面外变形信息，如褶皱的波长、幅度和方向等。

另一种是基于稳定性理论的分析方法。该方法主要有两种途径：一是采用非线性屈曲分析方法，通过引入初始缺陷的方式来触发褶皱的形成，Y.W.Wong，S.Pellegrino在《44th AIAA/ASME/ASCE/AHS/ASC Structures,Structural Dynamics,and MaterialsConference》，2003年4月发表的论文“Prediction of Wrinkle Amplitudes in SquareSolar Sails”中结合薄膜褶皱区域的应力平衡关系，得到了褶皱幅度及波长的表达式，并采用非线性屈曲分析的数值方法进行了验证；二是通过分岔分析方法来确定面外变形的分岔路径，Adama Diaby,Anh Le van,Christian Wielgosz在《Finite Elements inAnalysis and Design》，2006年7月第11期发表的论文“Buckling and wrinkling ofprestressed membranes”采用了分岔的分析方法，分析了预应力薄膜的屈曲和褶皱问题。第一种途径可以通过通用的软件进行褶皱的分析，但是引入的缺陷将对分析结果产生影响，第二种途径不需要引入任何缺陷，但在分析过程中需要引入分岔法则，来跟踪分岔路径，求解较为困难。

发明内容

本发明的目的在于避免上述现有技术分析方法的不足，提供一种方形太阳帆薄膜帆面褶皱形态的建模方法，达到求解简单，同时能够更准确地描述出褶皱形态的目的。

本发明的技术方案是：一种方形太阳帆表面褶皱形态的建模方法，它至少包括如下步骤：

(1)引入Airy应力场分布理论，得到更准确地应力圈模型；

(2)在步骤(1)的基础上，通过应力叠加准则，确定应力圈内任意一点处所受的总应力；

(3)基于步骤(2)，针对受等值角拉力方形太阳帆，消除边界影响，分析薄膜帆面的应力场分布，确定褶皱区域；

(4)在步骤(3)的基础上，分析角拉力不等时，薄膜帆面的应力场分布，确定褶皱区域；

(5)基于步骤(2)，步骤(3)和步骤(4)，推导任意正多边形结构受任意角拉力时的应力场分布，确定褶皱区域；

(6)针对方形太阳帆，通过如下的褶皱形态描述方程

其中L为方形太阳帆的边长，A为褶皱幅值，λ是褶皱半波长，h为褶皱长度，引入一个新的坐标系，ξ轴沿褶皱的方向，η轴垂直于褶皱方向，通过引入一个指数衰减函数g(y)，调节常系数k，确定出褶皱区域在Y轴方向的宽度，更准确的描述褶皱形态；

(7)根据步骤(6)中的褶皱形态描述方程，基于von Karman大挠度理论，得到褶皱形态与褶皱方向拉应力的关系模型。

上述步骤(1)中所述的引入Airy应力场分布理论，得到更准确地应力圈模型，其具体方法为：

假设在下半空间无限大的平面结构上的某点O₁处施加沿Y轴向上拉力T，传统的应力场模型中给出的应力等值圈定义为以施力点O₁为圆心，下半平面内任意选取的一点到O₁的距离为半径画圆，得到一组同心圆模型， 该模型只能近似的得到平均应力值，无法准确的描述各点的应力，因此引入Airy应力场分布理论，即针对平面问题，Airy在1862年提出的一个应力函数模型；

根据Airy应力场分布理论对应力等值圈的定义，圆心不再是施力点O₁，而是位于Y轴上的一组点，保证施力点O₁和下半平面内任意选取的一点均落在应力等值圈上，此时得到的应力圈模型也不再是一组同心圆而是都经过点O₁的一组圆；假设其中一个应力圈的圆心位于坐标轴原点，半径为R，应力圈内任意一点A(x,y)与施力点O₁间距离为r，与Y轴间夹角为β，当薄膜厚度为t时，点A处的径向应力σ_r为：

认为不存在环向应力和切应力，则径向应力可以分解到X方向和Y方向，如下：

其中σ_x为X方向应力，σ_y为Y方向应力，τ_xy为切应力，并且存在如下关系式：

将式(1)、(3)代入式(2)中，整理可得：

上述步骤(2)中所述的通过应力叠加准则，确定应力圈内任意一点处所受的总应力，其具体方法为：

同样选取步骤(1)中所述圆心位于坐标轴原点的一个应力圈，当该应力圈上存在两个拉力点时，假设拉力T₁沿Y轴向上，则T₁在应力圈内任意一点A(x,y)点处产生的应力可以由公式(4)得到；假设拉力T₂与拉力T₁间夹角为θ，计算T₂在A点处产生的应力时，可以将该应力圈旋转角θ，使拉力T₂的方向转到沿Y轴向上，在得到的新的坐标系OX′Y′中，点A(x,y)变为点A′(x′,y′)，转换关系式为：

将式(5)代入式(4)，可以得到拉力T₂在A点处产生的应力，与拉力T₁在A点处产生的应力进行叠加，即可得到应力圈内任意一点A处所受的总应力：

上述步骤(3)中所述的针对受等值角拉力方形太阳帆，消除边界影响，分析薄膜帆面的应力场分布，确定褶皱区域，具体方法为：

方形太阳帆展开后，四个角受到支撑梁提供的角拉力,假设分别沿X方向和Y方向且为等值角拉力T；首先定义正方形的外接圆为一个应力圈，该应力圈上任意一点B与正方形四个顶点的连线间的夹角分别为根据圆周角定理可知假设点B与各顶点的连线的方向分别为则点B处所受的沿这四个方向上的总应力是相等的，即：

选取沿方向的总应力其表达式为：

由式(1)可知是一个确定量；

标准的Airy应力场是作用于无限大空间的，然而方形帆面外接圆以外的区域是不存在应力的，因此需要引入一个压应力用来消除边界影响，最终通过简化整理后求得应力圈内各点的总应力：

其中σ_xi(i＝1,2,3,4)，σ_yi(i＝1,2,3,4)，τ_xyi(i＝1,2,3,4)分别为四个角拉力所产生的X方向应力、Y方向应力和切应力；

第一主应力σ₁和第二主应力σ₂分别为：

从式(10)中取出位于正方形内的各点，即可得到方形太阳帆薄膜帆面的应力分布云图；假设薄膜无法承受压应力，则第二主应力小于零的区域即为褶皱区域。

上述步骤(4)中所述的分析角拉力不等时，薄膜帆面的应力场分布，确定褶皱区域的具体方法为：

当方形帆面四个角施加的拉力不相等时，假设Y轴方向受到大小相等方向相反的拉力T₁,X轴方向受到大小相等方向相反的拉力T₂,同样定义正方形的外接圆为一个应力圈，首先计算拉力T₁在应力圈内部各点处的应力，此时应力圈上任意一点与拉力T₁的两个施力点的连线间的夹角为根据式(8)，引入的压应力为：

同理计算T₂在应力圈内部各点处的应力，根据式(8)，消除边界影响， 引入的压应力为：

最后叠加后可以得到应力圈内各点的总应力：

其中σ_x1i(i＝1,2)，σ_y1i(i＝1,2)，τ_xy1i(i＝1,2)为角拉力T₁作用下产生的X方向应力、Y方向应力和切应力，σ_x2i(i＝1,2)，σ_y2i(i＝1,2)，τ_xy2i(i＝1,2)为角拉力T₂作用下产生的X方向应力、Y方向应力和切应力；

根据式(10)计算出第一主应力σ₁和第二主应力σ₂，取出位于正方形内的各点，即可得到方形太阳帆表面的应力分布云图；假设薄膜无法承受压应力，则第二主应力小于零的区域即为褶皱区域。

上述步骤(5)中所述的推导任意正多边形结构受任意角拉力时的应力场分布，确定褶皱区域，具体步骤为：

(a)应力叠加准则

假设多边形边数为n(n≥3)，则相邻两个拉力间的夹角令其中一个拉力方向沿Y轴向上，根据步骤(2)中的应力叠加准则，结构内任意一点处所受的总应力为：

其中T_i(i＝1,2,…,n)为各个顶点施加的角拉力大小；

(b)消除边界影响，得到任意正多边形结构受任意角拉力时的应力场分布，确定褶皱区域，

然后根据步骤(3)和步骤(4)消除边界影响，定义正多边形的外接圆为一个应力圈，假设该应力圈上任意一点B与正多边形顶点的连线方向为根据式(7)可得：

由式(8)，可得任意正多边形外接圆上任意点处所受的沿方向的总应力可以表示为：

其中R为外接圆半径；

引入一个压应力用来消除边界影响，最终通过简化整理后求得应力圈内各点的总应力：

其中σ_xi(i＝1,2,…,n)，σ_yi(i＝1,2,…,n)，τ_xyi(i＝1,2,…,n)为各角拉力T_i(i＝1,2,…,n)作用下产生的X方向应力、Y方向应力和切应力；

最终根据式(10)计算出第一主应力σ₁和第二主应力σ₂，取出位于正多边形内的各点，即可得到正多边形结构内的应力分布云图，第二主应力小于零的区域即为褶皱区域。

上述步骤(6)的具体方法为；

当X轴方向拉力T₂是Y轴方向拉力T₁的4倍的时候，相应的角位移分别为δ₁和δ₂，此时可以根据步骤(4)确定出褶皱区域；假设L为方形太阳帆的边长，A为褶皱幅值，λ是褶皱半波长，λ′是中心处的褶皱半波长，h为褶皱长度，y表示沿坐标轴Y方向，引入一个新的坐标系，ξ轴沿褶皱的方向，η轴垂直于褶皱方向，想要预测对角线区域的褶皱波长及褶皱幅值参数，给出如下形式的褶皱形态描述方程：

w为面外变形，函数用以描述沿褶皱的方向上的褶皱形态，函数用以描述垂直于褶皱的方向上的褶皱形态，同时引入一个指数衰减函数g(y)，通过调节常系数k，可以确定出褶皱区域在Y轴方向的宽度，把结构划分为褶皱区域和非褶皱区域，即将褶皱形态限制在褶皱区域内，更准确的描述褶皱形态。

上述步骤(7)中所述的基于von Karman大挠度理论，得到褶皱形态与褶皱方向拉应力的关系模型的具体方法为：

根据非线性大扰度理论，褶皱区域的平衡方程为：

其中弯曲刚度E为薄膜的杨氏模量，v是泊松比；将式(18)代入式(19)中，整理可得：

在垂直于褶皱方向压应力的作用下，薄膜将从平面状态转变为稳定的褶皱状态，在这个转变过程中，褶皱的半波长将从无限大转变为某个确定值；相应于这个确定的半波长值，垂直于褶皱方向的压应力应当取得以褶皱半波长为自变量函数的极小值；根据函数取得极小值时的条件

将式(20)代入式(21)中，可以得到褶皱的半波长为：

根据式(1)可知由褶皱面外变形引起的垂直于褶皱方向上的收缩量为根据几何关系有：

则褶皱幅值为

求解位移δ₁和δ₂时，为了避免出现应力不收敛的情况，假设正方形薄膜每个角被切去了一个直角三角形，角拉力是均匀分布在一条很小的直边上，令此直边长度为r₁＝R/10，则：

将式(22)、(24)和(25)代入式(18)中，最终得到褶皱形态与褶皱方向拉应力的关系模型。

本发明的有益效果：本发明方法包括如下优点：

1)引入Airy应力场分布理论，通过应力叠加准则，并消除边界影响后，准确地得到方形太阳帆受角拉力时的应力场分布云图，且该方法适用于任意多边形结构，具有普遍性。

2)提出了一种新的褶皱形态描述方程，基于von Karman大挠度方程，采用应力极值法，得到褶皱形态与褶皱方向拉应力的关系模型，求解简单，同时得到的褶皱形态更接近实际情况。

以下将结合附图对本发明做进一步详细说明。

附图说明

图1是方形太阳帆结构图；

图2a是传统应力圈定义；

图2b是Airy应力圈的定义；

图3是一个Airy应力圈示意图；

图4是Airy应力叠加原理图；图4a为旋转前应力圈示意图；图4a为旋转后应力圈示意图；

图5是方形太阳帆结构边界应力图；

图6是正多边形结构受等值角拉力图；

图7是T₂＝4T₁时的褶皱区域分布图；

图8是ξη坐标系示意图；

图9是Ansys仿真得到的正方形结构的一个角区域的第一主应力分布云图；图9a是根据传统应力圈定义标出的应力圈等值线；图9b是根据Airy应力圈定义标出的应力圈等值线；

图10a是正方形受等值角拉力T＝4N时采用本发明方法得到的应力分布云图；图10a1是第一主应力云图；图10a2是第二主应力云图；图10a3是压应力区域；

图10b是正方形受等值角拉力T＝4N时利用Ansys仿真软件得到的应力分布云图；图10b1是第一主应力云图；图10b2是第二主应力云图；

图11是方形太阳帆Y方向拉力T₁＝4N，X方向拉力T₂＝16N时的应力分布云图；图11a是第一主应力云图；图11b是第二主应力云图；图11c是压应力区域；

图12是正六边形受等值角拉力时的应力分布云图，图12a是第一主应力云图；图12b是第二主应力云图；图12c是压应力区域；

图13是方形太阳帆表面褶皱形态图，图13a是采用已有方法得到的褶皱形态图；图13b是采用本发明方法得到的褶皱形态图；图13c是利用Ansys仿真已得到的褶皱形态图；

图14是采用不同方法得到的方形太阳帆表面沿Y方向中心线上的面外变形曲线。

附图标记说明：1、薄膜帆面；2、支撑梁；3、控制块。

具体实施方式

参照图1所示的两级驱动硬盘驱动器的结构图，其主要是由薄膜帆面1、支撑梁2、控制块3三部分组成。

本发明提供了一种方形太阳帆表面褶皱形态的建模方法，其方法步骤及其技术原理如下：

(1)Airy应力场分布理论

假设在下半空间无限大的平面结构上的某点O₁处施加沿Y轴向上拉力T，传统的应力场模型中给出的应力等值圈定义为，以施力点O₁为圆心，下半平面内任意选取的一点到O₁的距离为半径画圆，得到一组同心圆模型，如图2a所示，该模型只能近似的得到平均应力值，无法准确的描述各点的应力，因此引入Airy应力场分布理论；

图2b给出了Airy应力场分布理论对应力等值圈的定义，圆心不再是 施力点O₁，而是位于Y轴上的一组点，保证施力点O₁和下半平面内任意选取的一点均落在应力等值圈上，此时得到的应力圈模型也不再是一组同心圆而是都经过点O₁的一组圆；假设其中一个应力圈的圆心位于坐标轴原点，半径为R，应力圈内任意一点A(x,y)与施力点O₁间距离为r，与Y轴间夹角为β，如图3所示，当薄膜厚度为t时，点A处的径向应力σ_r为：

认为不存在环向应力和切应力，则径向应力可以分解到X方向和Y方向，如下：

其中σ_x为X方向应力，σ_y为Y方向应力，τ_xy为切应力，由图2可知存在如下关系式：

将式(1)、(3)代入式(2)中，整理可得：

(2)应力叠加准则

同样选取圆心位于坐标轴原点的一个应力圈，当该应力圈上存在两个 拉力点时，如图4所示，拉力T₁沿Y轴向上，则T₁在应力圈内任意一点A(x,y)点处产生的应力可以由公式(4)得到；假设拉力T₂与拉力T₁间夹角为θ，计算T₂在A点处产生的应力时，可以将该应力圈旋转角θ，使拉力T₂的方向转到沿Y轴向上，在得到的新的坐标系OX′Y′中，点A(x,y)变为点A′(x′,y′)，转换关系式为：

将式(5)代入式(4)，可以得到拉力T₂在A点处产生的应力，与拉力T₁在A点处产生的应力进行叠加，即可得到应力圈内任意一点A处所受的总应力：

(3)针对受等值角拉力方形太阳帆，消除边界影响，分析薄膜帆面的应力场分布，确定褶皱区域

对于图1所示方形太阳帆，四个角受到支撑梁提供的角拉力,分别沿X方向和Y方向，假设为等值角拉力T；如图5所示，首先定义正方形的外接圆为一个应力圈，该应力圈上任意一点B与正方形四个顶点的连线间的夹角分别为根据圆周角定理可知假设点B与 各顶点的连线的方向分别为则点B处所受的沿这四个方向上的总应力是相等的，即：

选取沿方向的总应力其表达式为：

由式(1)可知是一个确定量；

标准的Airy应力场是作用于无限大空间的，然而方形帆面外接圆以外的区域是不存在应力的，因此需要引入一个压应力用来消除边界影响，最终通过简化整理后求得应力圈内各点的总应力：

其中σ_xi(i＝1,2,3,4)，σ_yi(i＝1,2,3,4)，τ_xyi(i＝1,2,3,4)分别为四个角拉力所产生的X方向应力、Y方向应力和切应力；

第一主应力σ₁和第二主应力σ₂分别为：

从式(10)中取出位于正方形内的各点，即可得到方形太阳帆薄膜帆面的应力分布云图；假设薄膜无法承受压应力，则第二主应力小于零的区 域即为褶皱区域。

(4)角拉力不等时，薄膜帆面的应力场分布，确定褶皱区域

当方形帆面四个角施加的拉力不相等时，假设Y轴方向受到大小相等方向相反的拉力T₁,X轴方向受到大小相等方向相反的拉力T₂,同样定义正方形的外接圆为一个应力圈，首先计算拉力T₁在应力圈内部各点处的应力，此时应力圈上任意一点与拉力T₁的两个施力点的连线间的夹角为根据式(8)，引入的压应力为：

同理计算T₂在应力圈内部各点处的应力，根据式(8)，消除边界影响，引入的压应力为：

最后叠加后可以得到应力圈内各点的总应力：

其中σ_x1i(i＝1,2)，σ_y1i(i＝1,2)，τ_xy1i(i＝1,2)为角拉力T₁作用下产生的X方向应力、Y方向应力和切应力，σ_x2i(i＝1,2)，σ_y2i(i＝1,2)，τ_xy2i(i＝1,2)为角拉力T₂作用下产生的X方向应力、Y方向应力和切应力；

根据式(10)计算出第一主应力σ₁和第二主应力σ₂，取出位于正方形内的各点，即可得到方形太阳帆表面的应力分布云图；假设薄膜无法承受 压应力，则第二主应力小于零的区域即为褶皱区域。

(5)推广到任意正多边形受任意角拉力时的应力场分布，确定褶皱区域

(a)应力叠加准则

假设多边形边数为n(n≥3)，则相邻两个拉力间的夹角令其中一个拉力方向沿Y轴向上，根据应力叠加准则，结构内任意一点处所受的总应力为：

其中T_i(i＝1,2,…,n)为各个顶点施加的角拉力大小；

(b)消除边界影响，得到任意正多边形结构受任意角拉力时的应力场分布，确定褶皱区域，

如图6所示定义正多边形的外接圆为一个应力圈，假设该应力圈上任意一点B与正多边形顶点的连线方向为根据式(7)可得：

由式(8)，可得任意正多边形外接圆上任意点处所受的沿方向的总应力可以表示为：

其中R为外接圆半径；

引入一个压应力用来消除边界影响，最终通过简化整理后求得应力圈内各点的总应力：

其中σ_xi(i＝1,2,…,n)，σ_yi(i＝1,2,…,n)，τ_xyi(i＝1,2,…,n)为各角拉力T_i(i＝1,2,…,n)作用下产生的X方向应力、Y方向应力和切应力；

最终根据式(10)计算出第一主应力σ₁和第二主应力σ₂，取出位于正多边形内的各点，即可得到正多边形结构内的应力分布云图，第二主应力小于零的区域即为褶皱区域。

(6)一种新的褶皱形态描述方程的提出

当X轴方向拉力T₂是Y轴方向拉力T₁的4倍的时候，相应的角位移分别为δ₁和δ₂，此时可以通过式(9)确定出褶皱区域，如图7所示；假设L为方形太阳帆的边长，A为褶皱幅值，λ是褶皱半波长，λ′是中心处的褶皱半波长，h为褶皱长度，y表示沿坐标轴Y方向，引入一个新的坐标系，ξ轴沿褶皱的方向，η轴垂直于褶皱方向，想要预测对角线区域的褶皱波长及褶皱幅值参数，给出如下形式的褶皱形态描述方程：

w为面外变形，函数用以描述沿褶皱的方向上的褶皱形态，函数用以描述垂直于褶皱的方向上的褶皱形态，同时引入一个指数衰减函数g(y)，通过调节常系数k，可以确定出褶皱区域在Y轴方向的宽度，把结构划分为褶皱区域和非褶皱区域，即将褶皱形态限制在褶皱区域内，更准确的描述褶皱形态。

(7)褶皱形态与褶皱方向拉应力的关系模型的建立

根据非线性大扰度理论，褶皱区域的平衡方程为：

其中弯曲刚度E为薄膜的杨氏模量，v是泊松比；将式(18)代入式(19)中，整理可得：

在垂直于褶皱方向压应力的作用下，薄膜将从平面状态转变为稳定的褶皱状态，在这个转变过程中，褶皱的半波长将从无限大转变为某个确定值；相应于这个确定的半波长值，垂直于褶皱方向的压应力应当取得以褶皱半波长为自变量函数的极小值；根据函数取得极小值时的条件

将式(20)代入式(21)中，可以得到褶皱的半波长为：

根据式(1)可知由褶皱面外变形引起的垂直于褶皱方向上的收缩量为根据几何关系有：

则褶皱幅值为

根据Y.W.Wong and S.Pellegrino在《44th AIAA/ASME/ASCE/AHS/ASCStructures,Structural Dynamics,and Materials Conference》2003年发表论文“Prediction of Wrinkle Amplitudes in Square Solar Sails”中对位移δ₁和δ₂的定义，求解位移δ₁和δ₂时，为了避免出现应力不收敛的情况，假设正方形薄膜每个角被切去了一个直角三角形，角拉力是均匀分布在一条很小的直边上，令此直边长度为r₁＝R/10，则：

将式(22)、(24)和(25)代入式(18)中，最终得到褶皱形态与褶皱方向拉应力的关系模型。

下面结合一个典型的方形太阳帆结构说明本发明的具体实施过程，结构参数如表1所示。

表1

模型参数
		边长	L＝500mm
杨氏模量	E＝2.5GPa
		泊松比	v＝0.34
薄膜厚度	t＝25μm

式(18)中所涉及的常系数为k＝0.03，在Ansys软件中选取SHELL181单元模拟薄膜结构进行仿真分析，并通过相应实验验证本发明方法的有效性。

该方法的具体步骤如下：

(1)基于Airy应力场叠加原理，通过消除边界影响，准确地得到在角拉力作用下薄膜结构内部各点处的应力值，可以确定出褶皱区域，同时说明了该方法的普遍适用性，可用来求解任意多边形薄膜结构受任意角拉力时应力值。

(2)给出了一种新的褶皱形态描述方程，基于von Karman大挠度方 程，采用应力极值法，得到褶皱形态与褶皱方向拉应力的关系模型。最后通过Ansys软件进行仿真分析，选取SHELL181单元模拟薄膜结构，同时进行了相关实验验证。

(3)角区域应力圈等值线：图9给出了Ansys仿真得到的正方形结构的一个角区域的第一主应力分布云图，图9a是根据传统应力圈定义标出的应力圈等值线；图9b是根据Airy应力圈定义标出的应力圈等值线；通过比较可以看出，在顶点附近Airy应力圈模型与Ansys仿真结果完全吻合，说明了Airy应力圈模型的准确性；在图10b中的虚线应力圈处与Ansys仿真结果有一定差距是因为其他角拉力产生的应力在此处叠加造成的。

(4)等值角拉力下应力场分布云图：由图10a和图10b所示的方形太阳帆受等值角拉力T＝4N时的应力分布云图，图10a采用本发明方法得到的应力分布云图；图10a1是第一主应力云图；图10a2是第二主应力云图；图10a3是压应力区域。图10b是利用Ansys仿真软件得到的应力分布云图；图10b1是第一主应力云图；图10b2是第二主应力云图。可以看出本发明的方法和仿真结果基本一致，表明了本发明方法的准确性。

(5)本方明方法的推广：图11是方形太阳帆Y方向拉力T₁＝4N，X方向拉力T₂＝16N时的应力分布云图，图11a是第一主应力云图；图11b是第二主应力云图；图11c是压应力区域。图12给出了正六边形受等值角拉力时的应力分布云图，图12a是第一主应力云图；图12b是第二主应力云图；图12c是压应力区域。说明了本发明方法的普遍适用性。

(6)褶皱形态对比：当Y方向拉力T₁＝4N，X方向拉力T₂＝16N时，得到图13所示的方形太阳帆表面褶皱形态图，图13a是采用已有方法得到的褶皱形态图；图13b是采用本发明方法得到的褶皱形态图；图13c是利用Ansys仿真已得到的褶皱形态图。可以看出，本发明中对褶皱形态建模方法得到的褶皱形态更接近实际情况。

(7)面外变形曲线对比：图14给出了方形太阳帆表面沿Y方向中心线上的面外变形曲线，分别为实验得到的曲线、Ansys仿真曲线、采用本发明方法得到曲线和采用已有方法得到的曲线；进一步表明，本发明的建模方法在准确描述褶皱形态方面的有效性。

结果分析：本发明针对太阳帆结构提出了相对有效的褶皱形态建模方法，可以得到薄膜结构受角拉力时的应力场分布云图和褶皱形态，求解简单，且更加准确。

以上例举仅仅是对本发明的举例说明，并不构成对本发明的保护范围的限制，凡是与本发明相同或相似的设计均属于本发明的保护范围之内。

Claims (8)

1.一种方形太阳帆表面褶皱形态的建模方法，其特征是：它至少包括如下步骤：

(1)引入Airy应力场分布理论，得到更准确地应力圈模型；

(2)在步骤(1)的基础上，通过应力叠加准则，确定应力圈内任意一点处所受的总应力；

(3)基于步骤(2)，针对受等值角拉力方形太阳帆，消除边界影响，分析薄膜帆面的应力场分布，确定褶皱区域；

(4)在步骤(3)的基础上，分析角拉力不等时，薄膜帆面的应力场分布，确定褶皱区域；

(5)基于步骤(2)，步骤(3)和步骤(4)，推导任意正多边形结构受任意角拉力时的应力场分布，确定褶皱区域；

(6)针对方形太阳帆，通过如下的褶皱形态描述方程

其中w为面外变形，L为方形太阳帆的边长，A为褶皱幅值，λ是褶皱半波长，λ′是中心处的褶皱半波长，h为褶皱长度，y表示沿坐标轴Y方向，引入一个新的坐标系，ξ轴沿褶皱的方向，η轴垂直于褶皱方向，通过引入一个指数衰减函数g(y)，调节常系数k，确定出褶皱区域在Y轴方向的宽度，更准确的描述褶皱形态；

(7)根据步骤(6)中的褶皱形态描述方程，基于von Karman大挠度理论，得到褶皱形态与褶皱方向拉应力的关系模型。

2.根据权利要求1所述的一种方形太阳帆表面褶皱形态的建模方法，其特征是：步骤(1)中所述的引入Airy应力场分布理论，得到更准确地应力圈模型，其具体方法为：

假设在下半空间无限大的平面结构上的某点O₁处施加沿Y轴向上拉力T，传统的应力场模型中给出的应力等值圈定义为以施力点O₁为圆心，下半平面内任意选取的一点到O₁的距离为半径画圆，得到一组同心圆模型，该模型只能近似的得到平均应力值，无法准确的描述各点的应力，因此引入Airy应力场分布理论；

根据Airy应力场分布理论对应力等值圈的定义，圆心不再是施力点O₁，而是位于Y轴上的一组点，保证施力点O₁和下半平面内任意选取的一点均落在应力等值圈上，此时得到的应力圈模型也不再是一组同心圆而是都经过点O₁的一组圆；假设其中一个应力圈的圆心位于坐标轴原点，半径为R，应力圈内任意一点A(x,y)与施力点O₁间距离为r，与Y轴间夹角为β，当薄膜厚度为t时，点A处的径向应力σ_r为：

${\mathrm{\σ}}_r=\frac{2T}{\mathrm{\π}rt}cos\mathrm{\β}=\frac{T}{\mathrm{\π}Rt}---\left(1\right)$

认为不存在环向应力和切应力，则径向应力可以分解到X方向和Y方向，如下：

$\left\{\begin{array}{c}{\mathrm{\σ}}_x={\mathrm{\σ}}_r{\sin }^2\mathrm{\β}\\ {\mathrm{\σ}}_y={\mathrm{\σ}}_r{\cos }^2\mathrm{\β}\\ {\mathrm{\τ}}_{xy}={\mathrm{\σ}}_{r}sin\mathrm{\β}cos\mathrm{\β}\end{array}\right.---\left(2\right)$

其中σ_x为X方向应力，σ_y为Y方向应力，τ_xy为切应力，并且存在如下关系式：

$sin\mathrm{\β}=\frac{x}{r},cos\mathrm{\β}=\frac{R-y}{r}---\left(3\right)$

将式(1)、(3)代入式(2)中，整理可得：

$\left\{\begin{array}{c}{\mathrm{\σ}}_x=\frac{2T}{\mathrm{\π}t}\left(\frac{(R-y)x^2}{r^4}\right)\\ {\mathrm{\σ}}_y=\frac{2T}{\mathrm{\π}t}\left(\frac{{(R-y)}^3}{r^4}\right)\\ {\mathrm{\τ}}_{xy}=\frac{2T}{\mathrm{\π}t}\left(\frac{{(R-y)}^{2}x}{r^4}\right)\end{array}\right.---\left(4\right).$

3.根据权利要求1所述的一种方形太阳帆表面褶皱形态的建模方法，其特征是：步骤(2)中所述的通过应力叠加准则，确定应力圈内任意一点处所受的总应力，其具体方法为：

同样选取步骤(1)中所述圆心位于坐标轴原点的一个应力圈，当该应力圈上存在两个拉力点时，假设拉力T₁沿Y轴向上，则T₁在应力圈内任意一点A(x,y)点处产生的应力可以由公式(4)得到；假设拉力T₂与拉力T₁间夹角为θ，计算T₂在A点处产生的应力时，可以将该应力圈旋转角θ，使拉力T₂的方向转到沿Y轴向上，在得到的新的坐标系OX′Y′中，点A(x,y)变为点A′(x′,y′)，转换关系式为：

$\left\{\begin{array}{c}{x}^{\′}=xcos\mathrm{\θ}-y\sin \mathrm{\θ}\\ y^{\′}=xsin\mathrm{\θ}+y\cos \mathrm{\θ}\end{array}\right.---\left(5\right)$

将式(5)代入式(4)，可以得到拉力T₂在A点处产生的应力，与拉力T₁在A点处产生的应力进行叠加，即可得到应力圈内任意一点A处所受的总应力：

$\left\{\begin{array}{c}{\mathrm{\σ}}_x=\frac{2T}{\mathrm{\π}t}\left(\frac{\left(R-y\right)x^2}{r^4}\right)+\frac{2T}{\mathrm{\π}t}\left(\frac{\left(R-\left(x\sin \mathrm{\θ}+y\cos \mathrm{\θ}\right)\right){\left(x\cos \mathrm{\θ}-y\sin \mathrm{\θ}\right)}^2}{{\left({\left(x\cos \mathrm{\θ}-y\sin \mathrm{\θ}\right)}^2+{\left(x\sin \mathrm{\θ}+y\cos \mathrm{\θ}\right)}^2\right)}^2}\right)\\ {\mathrm{\σ}}_y=\frac{2T}{\mathrm{\π}t}\left(\frac{{\left(R-y\right)}^3}{r^4}\right)+\frac{2T}{\mathrm{\π}t}\left(\frac{{\left(R-\left(x\sin \mathrm{\θ}+y\cos \mathrm{\θ}\right)\right)}^3}{{\left({\left(x\cos \mathrm{\θ}-y\sin \mathrm{\θ}\right)}^2+{\left(x\sin \mathrm{\θ}+y\cos \mathrm{\θ}\right)}^2\right)}^2}\right)\\ {\mathrm{\τ}}_{xy}=\frac{2T}{\mathrm{\π}t}\left(\frac{{\left(R-y\right)}^{2}x}{r^4}\right)+\frac{2T}{\mathrm{\π}t}\left(\frac{{\left(R-\left(x\sin \mathrm{\θ}+y\cos \mathrm{\θ}\right)\right)}^2\left(x\cos \mathrm{\θ}-y\sin \mathrm{\θ}\right)}{{\left({\left(x\cos \mathrm{\θ}-y\sin \mathrm{\θ}\right)}^2+{\left(x\sin \mathrm{\θ}+y\cos \mathrm{\θ}\right)}^2\right)}^2}\right)\end{array}\right.---\left(6\right).$

4.根据权利要求1所述的一种方形太阳帆表面褶皱形态的建模方法，其特征是：步骤(3)中所述的针对受等值角拉力方形太阳帆，消除边界影响，分析薄膜帆面的应力场分布，确定褶皱区域，具体方法为：

方形太阳帆展开后，四个角受到支撑梁提供的角拉力,假设分别沿X方向和Y方向且为等值角拉力T；首先定义正方形的外接圆为一个应力圈，该应力圈上任意一点B与正方形四个顶点的连线间的夹角分别为根据圆周角定理可知假设点B与各顶点的连线的方向分别为则点B处所受的沿这四个方向上的总应力是相等的，即：

${\mathrm{\σ}}_{r_{\stackrel{\→}{r_1}}}={\mathrm{\σ}}_{r_{\stackrel{\→}{r_2}}}={\mathrm{\σ}}_{r_{\stackrel{\→}{r_3}}}={\mathrm{\σ}}_{r_{\stackrel{\→}{r_4}}}---\left(7\right)$

选取沿方向的总应力其表达式为：

由式(1)可知是一个确定量；

标准的Airy应力场是作用于无限大空间的，然而方形帆面外接圆以外的区域是不存在应力的，因此需要引入一个压应力用来消除边界影响，最终通过简化整理后求得应力圈内各点的总应力：

$\left\{\begin{array}{c}{\mathrm{\σ}}_x=\underset{i=1}{\overset{4}{\Σ}}{\mathrm{\σ}}_{xi}-\frac{2T}{\mathrm{\π}Rt}\\ {\mathrm{\σ}}_y=\underset{i=1}{\overset{4}{\Σ}}{\mathrm{\σ}}_{yi}-\frac{2T}{\mathrm{\π}Rt}\\ {\mathrm{\τ}}_{xy}=\underset{i=1}{\overset{4}{\Σ}}{\mathrm{\τ}}_{xyi}\end{array}\right.---\left(9\right)$

其中σ_xi(i＝1,2,3,4)，σ_yi(i＝1,2,3,4)，τ_xyi(i＝1,2,3,4)分别为四个角拉力所产生的X方向应力、Y方向应力和切应力；

第一主应力σ₁和第二主应力σ₂分别为：

$\left\{\begin{array}{c}{\mathrm{\σ}}_1=\frac{{\mathrm{\σ}}_x+{\mathrm{\σ}}_y}{2}+\sqrt{{\left(\frac{{\mathrm{\σ}}_x+{\mathrm{\σ}}_y}{2}\right)}^2+{\mathrm{\τ}}_{xy}^2}\\ {\mathrm{\σ}}_2=\frac{{\mathrm{\σ}}_x+{\mathrm{\σ}}_y}{2}-\sqrt{{\left(\frac{{\mathrm{\σ}}_x+{\mathrm{\σ}}_y}{2}\right)}^2+{\mathrm{\τ}}_{xy}^2}\end{array}\right.---\left(10\right)$

从式(10)中取出位于正方形内的各点，即可得到方形太阳帆薄膜帆面的应力分布云图；假设薄膜无法承受压应力，则第二主应力小于零的区域即为褶皱区域。

5.根据权利要求1所述的一种方形太阳帆表面褶皱形态的建模方法，其特征是：步骤(4)中所述的分析角拉力不等时，薄膜帆面的应力场分布，确定褶皱区域的具体方法为：

当方形帆面四个角施加的拉力不相等时，假设Y轴方向受到大小相等方向相反的拉力T₁,X轴方向受到大小相等方向相反的拉力T₂,同样定义正方形的外接圆为一个应力圈，首先计算拉力T₁在应力圈内部各点处的应力，此时应力圈上任意一点与拉力T₁的两个施力点的连线间的夹角为根据式(8)，引入的压应力为：

$T_{c1}=-(\mathrm{\σ}+{\mathrm{\σcos}}^2(\mathrm{\π}\left){\cos }^2\right(\frac{\mathrm{\π}}{2})+{\mathrm{\σcos}}^2(\frac{\mathrm{\π}}{2}\left){\cos }^2\right(\frac{\mathrm{\π}}{4}\left)\right)=-\mathrm{\σ}=-\frac{T_1}{\mathrm{\π}Rt}---\left(11\right)$

同理计算T₂在应力圈内部各点处的应力，根据式(8)，消除边界影响，引入的压应力为：

$T_{c2}=-(\mathrm{\σ}+{\mathrm{\σcos}}^2(\mathrm{\π}\left){\cos }^2\right(\frac{\mathrm{\π}}{2})+{\mathrm{\σcos}}^2(\frac{\mathrm{\π}}{2}\left){\cos }^2\right(\frac{\mathrm{\π}}{4}\left)\right)=-\mathrm{\σ}=-\frac{T_2}{\mathrm{\π}Rt}---\left(12\right)$

最后叠加后可以得到应力圈内各点的总应力：

$\left\{\begin{array}{c}{\mathrm{\σ}}_x={\mathrm{\σ}}_{x1}+{\mathrm{\σ}}_{x2}=\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\σ}}_{x1i}+\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\σ}}_{x2i}-\frac{T_1}{\mathrm{\π}Rt}-\frac{T_2}{\mathrm{\π}Rt}\\ {\mathrm{\σ}}_y={\mathrm{\σ}}_{y1}+{\mathrm{\σ}}_{y2}=\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\σ}}_{y1i}+\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\σ}}_{y2i}-\frac{T_1}{\mathrm{\π}Rt}-\frac{T_2}{\mathrm{\π}Rt}\\ {\mathrm{\τ}}_{xy}={\mathrm{\τ}}_{xy1}+{\mathrm{\τ}}_{xy2}=\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\τ}}_{xy1i}+\underset{i=1}{\overset{2}{\Σ}}{\mathrm{\τ}}_{xy2i}\end{array}\right.---\left(13\right)$

其中σ_x1i(i＝1,2)，σ_y1i(i＝1,2)，τ_xy1i(i＝1,2)为角拉力T₁作用下产生的X方向应力、Y方向应力和切应力，σ_x2i(i＝1,2)，σ_y2i(i＝1,2)，τ_xy2i(i＝1,2)为角拉力T₂作用下产生的X方向应力、Y方向应力和切应力；

根据式(10)计算出第一主应力σ₁和第二主应力σ₂，取出位于正方形内的各点，即可得到方形太阳帆表面的应力分布云图；假设薄膜无法承受压应力，则第二主应力小于零的区域即为褶皱区域。

6.根据权利要求1所述的一种方形太阳帆表面褶皱形态的建模方法，其特征是：步骤(5)中所述的推导任意正多边形结构受任意角拉力时的应力场分布，确定褶皱区域，具体步骤为：

(a)应力叠加准则

假设多边形边数为n(n≥3)，则相邻两个拉力间的夹角令其中一个拉力方向沿Y轴向上，根据步骤(2)中的应力叠加准则，结构内任意一点处所受的总应力为：

$\left\{\begin{array}{c}{\mathrm{\σ}}_x=\underset{i=1}{\overset{n}{\Σ}}\frac{2T_i}{\mathrm{\π}t}\left(\frac{\left(R-\left(x\sin \left(\left(i-1\right)\mathrm{\θ}\right)+y\cos \left(\left(i-1\right)\mathrm{\θ}\right)\right)\right){\left(x\cos \left(\left(i-1\right)\mathrm{\θ}\right)-y\sin \left(\left(i-1\right)\mathrm{\θ}\right)\right)}^2}{{\left({\left(x\cos \left(\left(i-1\right)\mathrm{\θ}\right)-y\sin \left(\left(i-1\right)\mathrm{\θ}\right)\right)}^2+{\left(x\sin \left(\left(i-1\right)\mathrm{\θ}\right)+y\cos \left(\left(i-1\right)\mathrm{\θ}\right)\right)}^2\right)}^2}\right)\\ {\mathrm{\σ}}_y=\underset{i=1}{\overset{n}{\Σ}}\frac{2T_i}{\mathrm{\π}t}\left(\frac{{\left(R-\left(x\sin \left(\left(i-1\right)\mathrm{\θ}\right)+y\cos \left(\left(i-1\right)\mathrm{\θ}\right)\right)\right)}^3}{{\left({\left(x\cos \left(\left(i-1\right)\mathrm{\θ}\right)-y\sin \left(\left(i-1\right)\mathrm{\θ}\right)\right)}^2+{\left(x\sin \left(\left(i-1\right)\mathrm{\θ}\right)+y\cos \left(\left(i-1\right)\mathrm{\θ}\right)\right)}^2\right)}^2}\right)\\ {\mathrm{\τ}}_{xy}=\underset{i=1}{\overset{n}{\Σ}}\frac{2T_i}{\mathrm{\π}t}\left(\frac{{\left(R-\left(x\sin \left(\left(i-1\right)\mathrm{\θ}\right)+y\cos \left(\left(i-1\right)\mathrm{\θ}\right)\right)\right)}^2\left(x\cos \left(\left(i-1\right)\mathrm{\θ}\right)-y\sin \left(\left(i-1\right)\mathrm{\θ}\right)\right)}{{\left({\left(x\cos \left(\left(i-1\right)\mathrm{\θ}\right)-y\sin \left(\left(i-1\right)\mathrm{\θ}\right)\right)}^2+{\left(x\sin \left(\left(i-1\right)\mathrm{\θ}\right)+y\cos \left(\left(i-1\right)\mathrm{\θ}\right)\right)}^2\right)}^2}\right)\end{array}\right.---\left(14\right)$

其中T_i(i＝1,2,…,n)为各个顶点施加的角拉力大小；

(b)消除边界影响，得到任意正多边形结构受任意角拉力时的应力场分布，确定褶皱区域，

然后根据步骤(3)和步骤(4)消除边界影响，定义正多边形的外接圆为一个应力圈，假设该应力圈上任意一点B与正多边形顶点的连线方向为根据式(7)可得：

${\mathrm{\σ}}_{r_{\stackrel{\→}{r_1}}}={\mathrm{\σ}}_{r_{\stackrel{\→}{r_2}}}=...={\mathrm{\σ}}_{r_{\stackrel{\→}{r_i}}}---\left(15\right)$

由式(8)，可得任意正多边形外接圆上任意点处所受的沿方向的总应力可以表示为：

其中 R为外接圆半径；

引入一个压应力用来消除边界影响，最终通过简化整理后求得应力圈内各点的总应力：

$\left\{\begin{array}{c}{\mathrm{\σ}}_x=\underset{i=1}{\overset{n}{\Σ}}({\mathrm{\σ}}_{xi}-\frac{{\mathrm{\σ}}_i}{2})\\ {\mathrm{\σ}}_y=\underset{i=1}{\overset{n}{\Σ}}({\mathrm{\σ}}_{yi}-\frac{{\mathrm{\σ}}_i}{2})\\ {\mathrm{\τ}}_{xy}=\underset{i=1}{\overset{n}{\Σ}}{\mathrm{\τ}}_{xyi}\end{array}\right.---\left(17\right)$

其中σ_xi(i＝1,2,…,n)，σ_yi(i＝1,2,…,n)，τ_xyi(i＝1,2,…,n)为各角拉力T_i(i＝1,2,…,n)作用下产生的X方向应力、Y方向应力和切应力；

最终根据式(10)计算出第一主应力σ₁和第二主应力σ₂，取出位于正多边形内的各点，即可得到正多边形结构内的应力分布云图，第二主应力小于零的区域即为褶皱区域。

7.根据权利要求1所述的一种方形太阳帆表面褶皱形态的建模方法，其特征是：步骤(6)的具体方法为；

当X轴方向拉力T₂是Y轴方向拉力T₁的4倍的时候，相应的角位移分别为δ₁和δ₂，此时可以根据步骤(4)确定出褶皱区域；假设L为方形太阳帆的边长，A为褶皱幅值，λ是褶皱半波长，λ′是中心处的褶皱半波长，h为褶皱长度，y表示沿坐标轴Y方向，引入一个新的坐标系，ξ轴沿褶皱的方向，η轴垂直于褶皱方向，想要预测对角线区域的褶皱波长及褶皱幅值参数，给出如下形式的褶皱形态描述方程：

$\left\{\begin{array}{c}w=\left(Ag\right(y\left)\right)\sin \left(\frac{\mathrm{\ξ}}{h}\mathrm{\π}\right)\sin \left(\frac{\mathrm{\η}}{\mathrm{\λ}}\mathrm{\π}\right)\\ g\left(y\right)=\frac{\exp (-k|y-\frac{\sqrt{2}}{2}L\left|\right)}{\exp (-k\frac{{\mathrm{\λ}}^{\′}}{2})}\end{array}\right.---\left(18\right)$

其中w为面外变形，函数用以描述沿褶皱的方向上的褶皱形态，函数用以描述垂直于褶皱的方向上的褶皱形态，同时引入一个指数衰减函数g(y)，通过调节常系数k，可以确定出褶皱区域在Y轴方向的宽度，把结构划分为褶皱区域和非褶皱区域，即将褶皱形态限制在褶皱区域内，更准确的描述褶皱形态。

8.根据权利要求1所述的一种方形太阳帆表面褶皱形态的建模方法，其特征是：步骤(7)中所述的基于von Karman大挠度理论，得到褶皱形态与褶皱方向拉应力的关系模型的具体方法为：

根据非线性大扰度理论，褶皱区域的平衡方程为：

$D{\▿}^{4}w-t\mathrm{\σ}\left(\mathrm{\ξ}\right)\frac{{\∂}^{2}w}{\∂{\mathrm{\ξ}}^2}+t\mathrm{\σ}\left(\mathrm{\η}\right)\frac{{\∂}^{2}w}{\∂{\mathrm{\η}}^2}=0---\left(19\right)$

其中弯曲刚度E为薄膜的杨氏模量，v是泊松比；将式(18)代入式(19)中，整理可得：

$\mathrm{\σ}\left(\mathrm{\η}\right)=\frac{{\mathrm{D\π}}^2}{{\mathrm{t\λ}}^2}+\[\frac{{\mathrm{D\π}}^2+t{\left(h\right)}^2\mathrm{\σ}\left(\mathrm{\ξ}\right)}{t{\left(h\right)}^4}\]{\mathrm{\λ}}^2+\frac{2{\mathrm{D\π}}^2}{t{\left(h\right)}^2}---\left(20\right)$

在垂直于褶皱方向压应力的作用下，薄膜将从平面状态转变为稳定的褶皱状态，在这个转变过程中，褶皱的半波长将从无限大转变为某个确定值；相应于这个确定的半波长值，垂直于褶皱方向的压应力应当取得以褶皱半波长为自变量函数的极小值；根据函数取得极小值时的条件

$\left\{\begin{array}{c}\frac{{\mathrm{d\σ}}_{\mathrm{\η}}\left(\mathrm{\λ}\right)}{d\mathrm{\λ}}=0\\ \frac{d^2{\mathrm{\σ}}_{\mathrm{\η}}\left(\mathrm{\λ}\right)}{{\mathrm{d\λ}}^2}>0\end{array}\right.---\left(21\right)$

将式(20)代入式(21)中，可以得到褶皱的半波长为：

$\mathrm{\λ}={\left(\frac{{\mathrm{D\π}}^2{\left(h\right)}^4}{{\mathrm{D\π}}^2+t{\left(h\right)}^2\mathrm{\σ}\left(\mathrm{\ξ}\right)}\right)}^{1/4}---\left(22\right)$

根据式(1)可知由褶皱面外变形引起的垂直于褶皱方向上的收缩量为根据几何关系有：

$-{\mathrm{\δ}}_2={\mathrm{\δ}}_1-\frac{{\mathrm{\π}}^2{A}^2}{4\mathrm{\λ}}---\left(23\right)$

则褶皱幅值为

$A=\frac{2\sqrt{\mathrm{\λ}({\mathrm{\δ}}_1+{\mathrm{\δ}}_2)}}{\mathrm{\π}}---\left(24\right)$

求解位移δ₁和δ₂时，为了避免出现应力不收敛的情况，假设正方形薄膜每个角被切去了一个直角三角形，角拉力是均匀分布在一条很小的直边上，令此直边长度为r₁＝R/10，则：

$\left\{\begin{array}{c}{\mathrm{\δ}}_1=\frac{T_1}{4Et}\[\mathrm{\π}\ln \left|\frac{\sqrt{2}L}{2r_1}\right|+\frac{2(1-v)}{L^2}\left(2\left(\sqrt{2}-1\right)L^2-\frac{{\left(2-\sqrt{2}\right)}^2{L}^2}{2}\right)\]\\ {\mathrm{\δ}}_2=\frac{T_2}{4Et}\[\mathrm{\π}\ln \left|\left(1-\frac{\sqrt{2}L}{2}\right)\frac{L}{r_1}\right|+\frac{4(1-v)}{\left(2-\sqrt{2}\right)L^2}\left(2\left(\sqrt{2}-1\right)L^2-\frac{{\left(2-\sqrt{2}\right)}^2{L}^2}{2}\right)\]\end{array}\right.---\left(25\right)$

将式(22)、(24)和(25)代入式(18)中，最终得到褶皱形态与褶皱方向拉应力的关系模型。