[发明专利]一种异型蜂窝结构的有限元建模方法

权利要求书

1.一种异型蜂窝结构的有限元建模方法，其特征在于：该方法具体步骤如下：

步骤一：建立异型蜂窝结构的截面几何模型；

首先，根据给定的异型蜂窝结构的边界特性，选择一个垂直于蜂窝厚度方向的截面，建立该蜂窝结构的截面几何模型，蜂窝结构的截面几何模型建立在笛卡尔坐标系的x-y面上，蜂窝结构的厚度方向位于z轴方向；

然后，将截面几何模型划分成多个形状规则的四边形分区，这里需进行多次点是否位于平面几何模型区域内的判断，以及蜂窝结构厚度方向的插值计算，这两种计算都基于四边形平面几何区域，故根据模型实际情况对截面几何模型进行划分；

最后，根据六边形蜂窝的几何特性进行有限元建模，假设蜂窝结构粘接方向位于y轴方向，若需要蜂窝粘接方向位于其它方向，则要将蜂窝结构的截面几何模型在坐标系中旋转相应角度进行有限元建模，待建模完成后，将有限元模型旋转回原角度，其中，以一点P₀(x₀，y₀)为旋转中心点，将点P(x，y)逆时针旋转θ角得到点P_t(x_t，y_t)的坐标值的运算用矩阵表示为：

xtyt=x0y0+xy-x0y0cosθsinθ-sinθcosθ---(1)]]>

步骤二：生成第K(0≤K≤Kp)环蜂窝结构的蜂窝中心点；

选择截面几何模型上的一个蜂窝结构中心点C(x_C，y_C)作为基准点，向外生成第K(0≤K≤Kp)环蜂窝结构的蜂窝中心点，这里以C点为中心点的蜂窝视为第0环，向外依次为第1，2，...，Kp环，Kp为本发明中预设的蜂窝环数上限；

当K＝0时，第k环上的蜂窝中心点为C(x_C，y_C)，当K≠0时，按下述方法生成第k环蜂窝中心点：

首先，生成第K(1≤K≤Kp)环上蜂窝中心点和C点的连线与x轴夹角为(i＝1，2...，6)的蜂窝结构中心点，这样的蜂窝中心点每环有六个，其坐标为：

xk,i=3Khcos[π6+(i-1)π3]+xCyk,i=3Khsin[π6+(i-1)π3]+yCi=1,2,..,6---(2)]]>

其中，h为蜂窝的边长；

然后，考虑到蜂窝在y轴方向上存在粘接厚度，将生成的上述中心点与基准点C间的y轴方向距离扩大Ap倍，即y′_k，i＝y_C+Ap(y_k，i-y_C)，其中，Ap为考虑粘接厚度的蜂窝宽度比，t为两蜂窝粘接处的厚度，即两蜂窝粘接处的实际厚度与蜂窝胞壁厚度之差；

最后，生成第K(2≤K≤Kp)环上其余的蜂窝中心点，这些蜂窝中心点位于每对相邻的已生成中心点连线的等分点上，选取已生成的两相邻中心点C_k，i(x_k，i，y′_k，i)和C_k，i+1(x_k，i+1，y′_k，i+1)，通过这两点生成二者之间的蜂窝中心点C_k，j的坐标值为：

xk,j=xk,i(1-q)+xk,i+1qyk,j=y′k,i(1-q)+y′k,i+1qj=1,...,K-1---(3)]]>

其中，q＝j/K，这样的操作总共进行六次，将第k环上其余的蜂窝中心点生成；

步骤三：在第K(0≤K≤Kp)环蜂窝结构中筛选出有角点位于截面几何模型区域内的蜂窝结构即域内蜂窝，记第k环域内蜂窝的数量为nf；

首先，对第K(0≤K≤Kp)环的每个蜂窝中心点C′(x_i，y_i)，写出其对应的蜂窝结构的第j(j＝1，2，...，6)个角点的坐标：

xi,j=xi+hcos(j-1)π3yi,j=yi+hsin(j-1)π3j=1,2,...,6---(4)]]>

其中，h为蜂窝的边长；

然后，对该蜂窝的第j(j＝1，2，...，6)个角点(x_i，j，y_i，j)，判断其是否位于截面几何模型区域内，若该蜂窝有角点位于截面几何模型区域内，则将该蜂窝的中心点坐标记录到一个矩阵cmf中，最终得到第k环的域内蜂窝个数，即矩阵cmf的行数nf；

判断一点P(x_P，y_P)是否位于几何模型区域内使用等参元逆变换插值法或者三角形面积法，两种方法所适用的几何模型都是平面四边形区域，故步骤一中将异型蜂窝结构的截面几何模型划分为多个四边形区域，针对这多个四边形区域分别进行点是否位于其域内的判断，若截面模型中存在这样一个区域使点在其之内，则说明点位于截面几何模型区域内；

其中，等参元逆变换插值法是根据四边形区域的解析性质对四边形区域和点P(x_P，y_P)进行坐标变换，坐标变换之后四边形区域变为D{(ξ，η)|-1≤ξ≤1，-1≤η≤1}，而点的坐标变为P′(ξ_P，η_P)，若点P(ξ_P，η_P)位于在区域D{(ξ，η)|-1≤ξ≤1，-1≤η≤1}之内，则说明点P(x_P，y_P)位于原四边形区域内；三角形面积法是以点P(x_P，y_P)为顶点和待判断四边形区域Ω的四条边构造四个三角形，若四个三角形的面积之和小于或等于四边形区域Ω的面积则说明点P(x_P，y_P)位于区域Ω内；

步骤四：若nf＜1，即第K(0≤K≤Kp)环不存在域内蜂窝，说明已生成的蜂窝中心点完全覆盖了异型蜂窝结构的截面几何模型区域，进而执行步骤十一，否则继续执行步骤五，本步骤的作用是避免计算资源的浪费；

步骤五：生成第K(0≤K≤Kp)环蜂窝结构的第i(i＝1，...，nf)个域内蜂窝在截面几何模型上的节点，并记此蜂窝位于几何模型区域内的节点个数为nin1；

本步骤生成第K(0≤K≤Kp)环蜂窝结构的第i(i＝1，...，nf)个域内蜂窝在截面几何模型上的节点，并将这个域内蜂窝的各节点信息写入一个矩阵NInf，NInf是一个6nd×5的矩阵，nd是蜂窝边上的单元密度，NInf的第一列到第五列所要记录的信息分别为：节点序号、节点是否位于蜂窝结构的角点、节点是否位于几何模型区域内、节点的横坐标、节点的纵坐标；

按以下方法生成该蜂窝结构的节点：

首先，对于中心点坐标为(x_i，y_i)的蜂窝，其对应的各节点坐标值为：

xi,j=sinπ3hsin(2π3-α)cosθ+xiyi,j=sinπ3hsin(2π3-α)sinθ+yij=1,2,...,6nd---(5)]]>

其中，h为蜂窝结构的边长，(j＝1，2，...，6nd)，

并且，由于蜂窝在y轴方向上存在粘接厚度，故需将各节点与中心点(x_i，y_i)在y轴方向上的距离扩大为原来的Ap倍，即y′_i，j＝y_i+(y_i，j-y_i)Ap，其中，Ap为考虑粘接厚度的蜂窝宽度比，h为蜂窝的边长，t为两蜂窝粘接处的粘接厚度；

然后，向矩阵NInf中写入各个节点的信息：

将该蜂窝结构的节点的坐标值写入矩阵NInf的第四列和第五列，将该蜂窝结构的节点序号即上述j值写入矩阵NInf的第一列，将该蜂窝结构的节点在蜂窝结构中的位置属性写入矩阵NInf的第二列，若节点位于蜂窝角点则元素值为1，若节点位于蜂窝边上则元素值为0，采用式m＝mod(j-1，nd)来进行节点在蜂窝结构中的位置属性的判断，将该蜂窝结构的节点在几何模型中的位置属性写入矩阵NInf的第三列，若节点位于几何模型区域内则元素值为1，若节点位于几何模型区域外则元素值为0；

第i(i＝1，...，nf)个域内蜂窝的节点全部生成之后，找出该蜂窝结构中位于几何模型区域内的节点个数，即矩阵NInf中第三列元素值为1的行数，记作nin1；

步骤六：对第K(0≤K≤Kp)环蜂窝结构的第i(i＝1，...，nf)个域内蜂窝进行处理，并记经过处理后此蜂窝位于几何模型区域内的节点个数为nin2；

本步骤旨在处理位于截面几何模型边界上的不完整的蜂窝结构，对其采取特殊的节点排布方法，以达到有限元模型与几何模型的匹配，对第K(0≤K≤Kp)环蜂窝结构的第i(i＝1，...，nf)个域内蜂窝进行处理之后，该蜂窝结构的信息矩阵NInf会有相应的更新，对蜂窝结构的节点进行处理的方法将另外详述；

经过处理之后，第i(i＝1，...，nf)个域内蜂窝的节点属性将发生一定的变化，找出此时该蜂窝结构位于几何模型区域内的节点个数，即矩阵NInf中第三列元素值为1的行数，记作nin2；

步骤七：判断nin2是否大于2，若大于2，则继续执行步骤八，否则直接执行步骤九；

nin2≤2说明第i(i＝1，...，nf)个域内蜂窝可近似地看做只有一条蜂窝边位于截面几何模型边界上，在实际情况中该蜂窝边并不存在，故摒弃此蜂窝，进而处理第K(0≤K≤Kp)环的第i+1个域内蜂窝；

步骤八：划分第K(0≤K≤Kp)环蜂窝结构的第i(i＝1，...，nf)个域内蜂窝的线单元；

将蜂窝结构的域内节点按顺序依次写入矩阵nodes，从矩阵nodes中选取适当的节点来构造蜂窝结构的线单元，将构成线单元的节点号写入一个矩阵elems，对于完整的域内蜂窝结构，其蜂窝结构中每两个相邻节点生成一个线单元，对于位于截面几何模型边界上的不完整的蜂窝结构，每两个相邻的域内节点生成一个线单元，但第一个域内节点和最后一个域内节点之间不生成线单元，即蜂窝结构被模型边界切割处没有胞壁结构；

由于蜂窝粘接处存在粘接厚度，进行有限元分析时粘接处壳单元需要赋给另外的厚度值，故需将粘接处的壳单元归类为一个集合，使用矩阵sets来记录这些单元，假设y轴方向为粘接方向，则将每个蜂窝结构中垂直于y轴的蜂窝边上的线单元号写入矩阵sets，对完整蜂窝来说，按其节点生成的顺序，蜂窝中第n个单元为粘接处单元的条件是：将(n-1)/nd向零取整值为1或4，其中，nd为蜂窝边上的单元密度；

当截面几何模型上的线单元全部生成之后，矩阵nodes储存截面几何模型上的全部节点信息，矩阵elems储存截面几何模型上的全部线单元信息，矩阵sets储存截面几何模型上蜂窝结构粘接处的单元号；

步骤九：i＝i+1，并判断i值是否大于nf，若小于或等于nf，则返回执行步骤五，若大于nf，则执行步骤十；

i＞nf说明第K(0≤K≤Kp)环的域内蜂窝结构的有限元模型已建立完成，进而处理第K+1环蜂窝结构；

步骤十：k＝k+1，并判断k值是否大于Kp，若小于或等于Kp，则返回执行步骤二，若大于Kp，规定的最大环数蜂窝结构的有限元模型已全部生成，继续执行步骤十一；

步骤十一：输出壳单元信息；

得到截面几何模型上线单元总数ElemNum和节点总数NodeNum，假设在z轴方向上有NLay层壳单元；

对于第iLay(iLay＝1，...，NLay)层上的第j(j＝1，...，ElemNum)个壳单元，其节点号为n1，n2，n1+NodeNum和n2+NodeNum，显然，iLay＝1时，n1和n2即矩阵elems中的对应线单元节点号，iLay＞1时，n1和n2即第iLay-1层壳单元的两个较大值的节点号；

按照输出要求，将单元所包含的节点号按行写入单元信息输出文件即可；

步骤十二：输出节点信息；

对于位于截面几何模型上的节点，其竖坐标为0，其它节点的竖坐标根据异型蜂窝结构几何模型特性利用厚度方向的插值方法得到，对于竖坐标值不为0的节点，为其构造一个节点坐标向量coord3，coord3的前两个元素是该节点对应在矩阵nodes中所记录的横坐标和纵坐标值，coord3的第三个元素，即该节点的竖坐标值，采用基于几何模型四边形分区的四个角点坐标进行平面插值的方法得到，最终形成的有限元模型在蜂窝结构的顶面将会是分片的平面，或者更精确的是基于节点在x-y面上的投影所在的四边形分区选用相应的插值函数进行曲面插值，将在蜂窝结构的顶面得到与实际模型更加精确匹配的分片曲面；

将每个节点的坐标经过角度旋转的坐标变换，得到原角度位置下的节点坐标，按输出要求输出，这里需要注意节点信息的输出顺序，位于截面几何模型上的节点应按照矩阵nodes的存储的顺序，其它每层节点都应对应于矩阵nodes的存储顺序，这样才能使节点信息与单元信息输相对应；

步骤十三：输出蜂窝结构粘接处的单元集合信息；

根据输出要求，将蜂窝结构粘接处的单元号写入集合信息输出文件即可。

2.根据权利要求1所述的一种异型蜂窝结构的有限元建模方法，其特征在于：步骤六中所述的对蜂窝结构的节点进行处理的方法，其具体实施步骤如下：

步骤一：找出该蜂窝结构位于截面几何模型区域内的节点个数，即该蜂窝结构的节点信息矩阵NInf中第三列元素值大于零的行数，记作nin；

步骤二：判断nin是否等于6rd，若nin＝6nd，说明该蜂窝结构是一个完整的域内蜂窝，直接执行步骤十三，若nin≠6nd，则说明蜂窝结构被截面几何模型边界截断，失去了完整性，执行步骤三；

步骤三：对此蜂窝结构的节点进行排序；

本步骤是使此蜂窝的节点信息矩阵NInf中从上到下的各行所代表的节点在蜂窝结构中成逆时针顺序排列，并且，矩阵NInf中的第一行代表按逆时针顺序排列的截面几何模型区域内的第一点；将矩阵NInf的第一行挪至最后一行，其余各行分别向上移一行，这种操作进行数次直至矩阵NInf符合本步骤的要求；

步骤四：将矩阵NInf的最后一行的信息储存在一个行向量pend中，显然，pend表示按逆时针顺序将要进入截面几何模型区域的节点信息；

步骤五：规定一个误差范围herr；

步骤六：处理该蜂窝结构位于截面几何模型区域内按逆时针顺序排列的最后一个节点处；

首先，将截面几何模型区域内最后一个节点记作P1即对应矩阵NInf的第nin行，P₁的下一个节点记作P₂即对应矩阵NInf的第nin+1行，求得P₁和P₂连线与截面几何模型边界的交点pin1；

然后，考察P₁与pin1之间的距离，若距离大于herr，则将区域外的节点移到截面几何模型边界上，即将点pin1的坐标值赋给矩阵NInf的第nin+1行的第四列和第五列，同时，矩阵NInf的第nin+1行所代表的节点的属性改为“不位于蜂窝角点”和“位于截面几何模型区域内”，即将矩阵NInf的第nin+1行的第二列元素值改为0，第三列元素值改为1；

若P₁与pin1之间的距离小于或等于herr，考虑nin为1的情况，这说明这个蜂窝只有一个节点在区域内，且界内点P₁与交点pin1的距离在误差范围之内，则将这个蜂窝舍弃，返回到权利要求1中的步骤十二，若nin不为1，则将区域内最后一个节点移到截面几何模型边界上，即将点pin1的坐标值赋给矩阵NInf的第nin行的第四列和第五列，同时，矩阵NInf的第nin行所代表的节点的属性改为“不位于蜂窝角点”和“位于截面几何模型区域内”，即将矩阵NInf的第nin行的第二列元素值改为0，第三列元素值改为1；

步骤七：找出矩阵NInf的第三列元素值大于0的行数，并重新赋给nin，这表示这个蜂窝经过以上步骤的操作，该蜂窝结构位于截面几何模型区域内的节点个数nin进行了更新；

步骤八：处理该蜂窝结构位于截面几何模型区域内按逆时针顺序排列的第一个节点处；

首先，将行向量pend所记录的节点记作P₁，将截面几何模型区域内第一点记作P₂即对应矩阵NInf的第1行，求得蜂窝结构与截面几何模型边界的交点为pin2；

然后，考察P₂与pin2之间的距离，若距离大于herr，考虑nin等于矩阵NInf的行数的情况，此时该蜂窝经过步骤六的操作所有节点都变为域内点，需要在矩阵NInf的最后添加一行，即增加第6nd+1个节点，这一行的内容可暂时复制矩阵NInf的最后一行的内容，若nin不等于矩阵NInf的行数，则不需要做增加行的操作；

将矩阵NInf的最后一行挪到第一行，其余各行分别向后移一行，将此时矩阵NInf的第一行代表域内第一个节点的信息，将pin2的坐标值赋给矩阵NInf第一行的第四列和第五列，同时，矩阵NInf第一行所代表的节点属性改为“不位于蜂窝角点”和“位于截面几何模型区域内”，即将矩阵NInf的第一行的第二列元素值改为0，第三列元素值改为1；

若P₂与pin2之间的距离小于或等于herr，则将区域内的第一个节点移到截面几何模型边界上，即将pin2的坐标值赋给矩阵NInf第一行的第四列和第五列，同时，矩阵NInf第一行所代表的节点属性改为“不位于蜂窝角点”和“位于截面几何模型区域内”，即将矩阵NInf的第一行的第二列元素值改为0，第三列元素值改为1；

步骤九：找出矩阵NInf的第三列元素值大于0的行数，并重新赋给nin，这表示这个蜂窝经过上述步骤的操作，该蜂窝结构位于截面几何模型区域内的节点个数nin进行了更新；

步骤十：若蜂窝结构位于截面几何模型区域内的第一个节点和最后一个节点距离小于herr，则进行节点合并；

考察矩阵NInf的第1行所表示的节点与第nin行所表示的节点之间的距离值，若距离小于herr，则将矩阵NInf的第1行所表示的节点移到两节点的中点处，即将矩阵NInf的第一行的第四列和第五列元素值改为两节点的中点坐标值，同时，将第nin个节点的属性改为“位于截面几何模型区域外”，即将第nin行的第三列的元素值改为0；

步骤十一：找出该蜂窝结构中位于截面几何模型区域内且位于蜂窝角点上的节点个数，即矩阵NInf中第二列和第三列元素值都为1的行数，记作ncd，若ncd＜1，则说明该蜂窝没有角点位于截面几何模型区域内，返回权利要求1中的步骤十二，若ncd≥1，则找出矩阵NInf的第三列元素值大于0的行数，并重新赋给nin；

步骤十二：判断这个蜂窝的样式并进行蜂窝节点排序；

首先，记该蜂窝在截面几何模型区域内的第一个位于蜂窝角点上的节点在矩阵NInf中位于第pos₁行，并记d₁＝pos₁-1，记该蜂窝在截面几何模型区域内的最后一个位于蜂窝角点上的节点在矩阵NInf中位于第pos₂行，并记d₂＝nin-pos₂；

然后，若d₁＞d₂，则需将该蜂窝的节点变为如下顺序：原域内的最后一点作为域内第一点即对应矩阵NInf的第一行，并且各节点按顺时针顺序排列；

步骤十三：找出矩阵NInf的第三列元素值大于0的行数，并重新赋给nin，这表示这个蜂窝经过上述步骤的操作，该蜂窝结构位于截面几何模型区域内的节点个数nin进行了更新；

对蜂窝进行处理时，需要求四点构成的两线段的交点坐标，具体方法如下：

首先，假设两条线段和的端点坐标分别为P₁(x₁，y₁)、P₂(x₂，y₂)、P₃(x₃，y₃)、P₄(x₄，y₄)；

然后，构造下列参数：

d1=-(x3-x4)(y1-y2)+(x1-x2)(y3-y4)d2=-(x3-x4)(y1-y2)+(x1-x2)(y3-y4)---(6)]]>

若d₁d₂＝0，则两线段平行，无交点，否则，两线段的交点或延长线的交点P₀(x，y)的坐标值为：

x=x1[x4(-y2+y3)+x3(y2-y4)]+x2[x4(y1-y3)+x3(-y1+y4)]d1y=x4(y1-y2)y3+x1y2(y3-y4)+x3(-y1+y2)y4+x2y1(-y3+y4)d2---(7).]]>