CN100587696C - 一种参数化的静态超单元构造方法 - Google Patents

Abstract

本发明公开了一种参数化的静态超单元构造方法，该单元利用子结构的方法，降低模型的自由度，然后对结构不同尺寸进行参数化；对刚度矩阵K、应变矩阵B进行拟合，构造出超单元，该超单元结点数量少，可大量节省计算时间，易于优化，且计算结果能保证精度要求。基于这种超单元的构造方法，可构造出计算量大大降低的、各种各样的超单元。

Description

一种参数化的静态超单元构造方法

技术领域

本发明涉及数值模拟和CAE领域，尤其是一种参数化的静态超单元构造方法。

背景技术

随着现代计算机技术的发展，数值建模技术已经成为描述结构特性的主要手段，目前在工程领域内常用有限单元法的数值模拟方法来进行数值建模。有限单元法的基本思想是将问题的求解域划分为一系列的单元，单元之间仅靠结点相连。单元内部的待求量可由单元结点通过选定的函数关系插值得到。由于单元形状简单，易于平衡关系和能量关系建立结点量的方程式，然后将各单元方程组集成总体代数方程组，计入边界条件后可对方程求解。

有限元理论自出现以来，许多学者在新型有限单元、有效的数值求解方法、结构模型网格划分技术以及有限元计算软件等方面做了大量的研究工作。但是，随着现代工程结构逐渐大型化、复杂化的趋势，为了确保得到比较精确的计算模型，基于有限元方法建立起来的模型往往规模过于庞大。由于有限元模型过于复杂，会导致计算容量过大，计算所需时间过长等等。针对这些问题，很多的模型简化方法被提出，对于存在重复结构的静态问题求解，可以利用子结构的方法，实现对复杂有限元模型的简化处理。

对于一些大型复杂结构中具有大量连续平板中心开孔的结构，如图1所示，平板的两侧受到约束、载荷作用，中心孔处存在应力集中，对其进行有限元分析时，为了得到开孔处较为精确的结果，需要划分成较密的网格进行计算，工作量大，计算时间长，而且不容易对其进行优化。

发明内容

针对现有技术存在的问题，本发明的目的在于提供一种建立模型简单，结点数量少，易于优化，且能保证计算结果精度要求的参数化的静态超单元构造方法。

为实现上述目的，本发明一种参数化的静态超单元构造方法，具体为：

1)针对各种有重复结构的构件，将其重复结构的特征尺寸设为参数；

2)利用有限元软件对该重复结构建立不同特征尺寸的有限元模型，划分网格，分别读出每个结构的刚度矩阵；

3)利用有限元中子结构的方法，仅把关心部位的结点作为外部结点保留，做成超单元，得到超单元的刚度矩阵K；

4)对3)中所得到各种模型的刚度矩阵K进行拟合，得到刚度矩阵K和其特征参数相关的函数；

5)对应变矩阵B进行拟合，得到应变矩阵B和其特征参数相关的函数；利用已求得的结点位移，进而计算出单元高斯点的应力值，再通过外插法得到结点的应力值。

本发明一种关于结构参数化的静态超单元的构造方法，该单元利用子结构的方法对重复结构进行简化，建立模型简单，结点数量少，可大量节省计算时间，易于优化，并以平板开孔结构为例，将平板的中心孔的长短半径设定为参数，对不同情况进行拟合后，计算出位移、应力的结果能保证精度要求。

附图说明

图1为连续平板中心开孔的结构；

图2为两侧边界上受外载荷和约束的超单元；

图3为图2中超单元的1/4有限元模型；

图4为图2中超单元的有限元模型；

图5为物理空间中任意的四边形单元；

图6为在参考坐标系ζη下的图5中四边形单元；

图7为实例中超单元的有限元模型；

图8为图7中超单元的总刚点号；

图9为3个连续平板中心开孔结构的有限元模型；

图10为图9中超单元的总刚点号。

图11为超单元中心孔的单刚点号；

图12为实施例1中超单元的边界结点位移的ANSYS计算结果和超单元计算结果；

图13为实施例1中超单元的中心孔位移的ANSYS计算结果和超单元计算结果；

图14为实施例1中超单元的中心孔应力的ANSYS计算结果和超单元计算结果。

具体实施方式

本发明中所构造的超单元为平板中心开孔结构，只在板的两侧边界上受外载荷和约束，建立超单元如图2所示，每个超单元只有边界结点。每侧17个，共34个结点。

其中，超单元的长和宽尺寸固定0.3×0.2m；椭圆孔的长半径r1＝0.03～0.11m；短半径r2＝0.02～0.06m。

在ANSYS软件中建立有限元模型，如图3所示，先建立1/4有限元模型，单元划分时画成放射性的网格，比率是0.5，中心处网格较密，边界网格较疏松，这样可以保证内外网格尺寸变化不大，然后再镜像出整个模型(如图4所示)。

利用ANSYS内部程序读出整体单元刚度矩阵K，凝聚内部自由度，将所读出的刚度矩阵按照一定的顺序重新排列(逆时针)：外部结点(外部两侧边界结点，中间孔结点，中间孔第二圈结点)，内部结点。

$\left[\begin{array}{cc}{K}_{11}& K_{12}\\ K_{21}& K_{22}\end{array}\right]\left[\begin{array}{c}{u}_1\\ u_2\end{array}\right]=\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]---\left(1\right)$

其中：u₁为外部两侧边界结点，中间孔结点，中间孔第二圈结点的位移；

u₂为其它内部结点位移；

f₁为外部载荷；

f₂为内部载荷0。

K₁₁×u₁+K₁₂×u₂＝f₁                        (2)

K₂₁u₁+K₂₂u2＝0                             (3)

由方程(3)可得到内部位移由外部位移表示的关系式

$u_2=-K_{22}^{-1}{K}_{21}u1---\left(4\right)$

将(4)式带入(2)式，可得由外部结点(外部两侧边界结点，中间孔结点，中间孔第二圈结点)刚度矩阵、位移和载荷之间的关系。

$(K_{11}-K_{12}{K}_{22}^{-1}{K}_{21})u_1=f_1---\left(5\right)$

令 $K0=(K_{11}-K_{12}{K}_{22}^{-1}{K}_{21})$ 为外部结点(外部两侧边界结点，中间孔结点，中间孔第二圈结点)的刚度矩阵。

将K0再进行划分

$\left[\begin{array}{cc}{\mathrm{KK}}_{11}& {\mathrm{KK}}_{12}\\ {\mathrm{KK}}_{21}& {\mathrm{KK}}_{22}\end{array}\right]\left[\begin{array}{c}{\mathrm{uu}}_1\\ {\mathrm{uu}}_2\end{array}\right]=\left[\begin{array}{c}{\mathrm{ff}}_1\\ {\mathrm{ff}}_2\end{array}\right]---\left(6\right)$

其中：uu₁为外部两侧边界结点；

uu₂为中间孔结点中间孔第二圈结点的位移(计算应力时使用)；

ff₁为外部载荷；

ff₂为内部载荷0。

KK₁₁uu₁+KK₁₂uu₂＝ff₁                      (7)

KK₂₁uu₁+KK₂₂uu₂＝0                        (8)

由方程(7)可得到内部位移由外部位移表示

${\mathrm{uu}}_2=-{\mathrm{KK}}_{22}^{-1}{\mathrm{KK}}_{21}\mathrm{uu}1---\left(9\right)$

即由外部结点位移到内部结点位移的转换矩阵为

$T\_u=-{\mathrm{KK}}_{22}^{-1}{\mathrm{KK}}_{21}---\left(10\right)$

将(9)式带入(7)式，可得由外部结点(外部两侧边界结点)刚度矩阵、位移和载荷之间的关系。

$({\mathrm{KK}}_{11}-{\mathrm{KK}}_{12}{\mathrm{KK}}_{22}^{-1}{\mathrm{KK}}_{21}){\mathrm{uu}}_1={\mathrm{ff}}_1---\left(11\right)$

令 $\mathrm{KK}0=({\mathrm{KK}}_{11}-{\mathrm{KK}}_{12}{\mathrm{KK}}_{22}^{-1}{\mathrm{KK}}_{21})$ 为外部结点(外部两侧边界结点)的刚度矩阵。经过自由度凝聚后得到外部结点刚度矩阵、位移、载荷的关系为

KK0uu₁＝ff₁                (12)

将上述模型左边约束，另一边结点分别均施加x、y方向结点载荷。

由公式(12)计算出外部结点的位移值，再代入公式(9)中可得到内部中心孔处单元结点的位移值。

在ANSYS中分别计算中心椭圆孔长半径为0.03～0.11m；短半径0.02～0.06m的情况，根据情况选择较好的拟合公式，将计算结果中的T_u、KK0分别进行拟合，得到相应系数，计算过程如下：

$\left[\begin{array}{ccccc}a1& a2& a3& a4& a5\end{array}\right]\left[\begin{array}{c}{r1}^2\·{r2}^2\\ r1\·r2\\ r1\\ r2\\ 1\end{array}\right]=\mathrm{KK}0---\left(13\right)$

$\left[\begin{array}{ccccccc}b1& b2& b3& b4& b5& b6& b7\end{array}\right]\left[\begin{array}{c}{r1}^2\·{r2}^2\\ r1\·r2\\ r1\\ r2\\ 1\\ r1\·r1\\ r2\·r2\end{array}\right]=T\_u---\left(14\right)$

可拟合出矩阵a1、a2、a3、a4、a5；b1、b2、b3、b4、b5、b6、b7。

因此，在r1＝0.03～0.11m；r2＝0.02～0.06m范围内，只需计算将r1、r2代入公式(13)、(14)中，即可计算出外部结点的刚度矩阵KK0和T_u，通过施加外部约束和载荷，计算出外部结点位移及内部结点位移。

在ANSYS中读取中心孔处一圈的单元号、每个单元的结点号及结点坐标，并按逆时针排列。编写四边形平面应力单元的应变矩阵(B矩阵)，选择积分阶数为2，代入高斯点±0.577350269189626，权系数为1，计算出应变矩阵B。

任意四边形单元应变矩阵具体计算过程如下：

如图5所示为一个物理空间中任意的四边形单元。在它上面建立一个参考坐标系ζη，他们不需要互相正交，也不需要与整体坐标系平行。单元的边被ξ和η平分，它们的方程分别是ξ＝±1和η＝±1，在参考坐标系ζη下，该四边形变换成一个边长为2的正方形(如图6所示)。

根据等参元的含义，它的位移坐标都采用相同的形函数表示，即

$u=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_i(\mathrm{\ξ},\mathrm{\η})u_i,$ $v=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_i(\mathrm{\ξ},\mathrm{\η})v_i---\left(15\right)$

$x=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_i(\mathrm{\ξ},\mathrm{\η})x_i$ $y=\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_i(\mathrm{\ξ},\mathrm{\η})y_i---\left(16\right)$

式中u_i和v_i(i＝1，2，3，4)是整体坐标系xy下结点的位移分量，x_i和y_i是结点的整体坐标系，N_i(ξ，η)是由参考坐标系ξ和η表示的形函数。

N_i的表达式可写成：

$N_i=\frac{1}{4}(1+{\mathrm{\ξ}}_i\mathrm{\ξ})(1+{\mathrm{\η}}_i\mathrm{\η})(i=\mathrm{1,2,3,4})---\left(17\right)$

式中ξ_i和η_i是4个结点的参考坐标值，他们的值为±1。参考坐标系ξ和η在整体坐标系xy下的方向是由单元的结点编号决定的。如图5所示结点编号是1-2-3-4。

将位移表达式带入几何方程，变得到应变分量的计算公式

$\mathrm{\ϵ}=\left[\begin{array}{c}\frac{\∂u}{\∂x}\\ \frac{\∂v}{\∂y}\\ \frac{\∂u}{\∂y}+\frac{\∂v}{\∂x}\end{array}\right]={\mathrm{B\δ}}^e=\left[\begin{array}{cccc}{B}_1& B_2& B_3& B_4\end{array}\right]{\mathrm{\δ}}^e---\left(18\right)$

式中 ${\mathrm{\δ}}^e={\left[\begin{array}{cccc}{\mathrm{\δ}}_1^T& {\mathrm{\δ}}_2^T& {\mathrm{\δ}}_3^T& {\mathrm{\δ}}_4^T\end{array}\right]}^T$ 是单元结点的位移列阵，以及

$B_i=\left[\begin{array}{cc}{N}_{i^{\′}x}& 0\\ 0& N_{i^{\′}y}\\ N_{i^{\′}y}& N_{i^{\′}x}\end{array}\right],$ ${\mathrm{\δ}}_i=\left[\begin{array}{c}{u}_i\\ v_i\end{array}\right](i=\mathrm{1,2,3,4})---\left(19\right)$

这里下标中的逗号表示求导数，即 $N_{i,x}={\∂N}_i/\∂x,$ $N_{i,y}={\∂N}_i/\∂y.$ 由于形函数是由参考坐标ξ和η给出的，这两个导数一般不能显示给出。

根据复合函数求导规则，它们可以按下列步骤计算：

$\left[\begin{array}{c}{N}_{i^{\′}\mathrm{\ξ}}\\ N_{i^{\′}\mathrm{\η}}\end{array}\right]=\left[\begin{array}{cc}x_{\′\mathrm{\ξ}}& y_{\′\mathrm{\ξ}}\\ x_{\′\mathrm{\η}}& y_{\′\mathrm{\η}}\end{array}\right]\left[\begin{array}{c}{N}_{i^{\′}x}\\ N_{i^{\′}y}\end{array}\right]=J\left[\begin{array}{c}{N}_{i^{\′}x}\\ N_{i^{\′}y}\end{array}\right]---\left(20\right)$

式中的矩阵J称为雅克比矩阵，即

$J=\left[\begin{array}{cc}x_{\′\mathrm{\ξ}}& y_{\′\mathrm{\ξ}}\\ x_{\′\mathrm{\η}}& y_{\′\mathrm{\η}}\end{array}\right]=\left[\begin{array}{cc}\underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_{i^{\′}\mathrm{\ξ}}{x}_i& \underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_{i^{\′}\mathrm{\ξ}}{y}_i\\ \underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_{i^{\′\mathrm{\η}}}{x}_i& \underset{i=1}{\overset{4}{\mathrm{\Σ}}}{N}_{i^{\′\mathrm{\η}}}{y}_i\end{array}\right]---\left(21\right)$

则形函数对x和y的导数为

$\left[\begin{array}{c}{N}_{i^{\′}x}\\ N_{i^{\′}y}\end{array}\right]=J^{-1}\left[\begin{array}{c}{N}_{i^{\′}\mathrm{\ξ}}\\ N_{i^{\′}\mathrm{\η}}\end{array}\right]---\left(22\right)$

由此可得应变子矩阵Bi

$B_i=J^{-1}\left[\begin{array}{cc}{N}_{i^{\′}\mathrm{\ξ}}& 0\\ 0& N_{i^{\′}\mathrm{\η}}\\ N_{i^{\′}\mathrm{\η}}& N_{i^{\′}\mathrm{\ξ}}\end{array}\right]---\left(23\right)$

在ANSYS中分别计算椭圆孔长半径为0.03～0.11m；短半径0.02～0.06m的情况，对于半径不同的r1、r2的有限元模型分别计算出中心孔一圈单元的B矩阵，选择较好的拟合公式，对B矩阵进行拟合，计算过程如下：

$\left[\begin{array}{ccccccccc}c1& c2& c3& c4& c5& c6& c7& c8& c9\end{array}\right]\left[\begin{array}{c}{r1}^2\·{r2}^2\\ r1\·r2\\ r1\\ r2\\ 1\\ {r1}^2\\ {r2}^2\\ {r2}^2\·r1\\ {r1}^2\·r2\end{array}\right]=B---\left(24\right)$

可拟合出矩阵c1、c2、c3、c4、c5、c6、c7、c8、c9。

因此，在r1＝0.03～0.11m；r2＝0.02～0.06m范围内，只需计算将r1、r2代入公式(12)中，即可计算出中心孔处一圈的应变矩阵B，

再计算出中心孔的应力值：

σ＝D×B×u                (25)

其中：u为位移(计算中得到的中心孔单元的各结点位移值)，σ为单元高斯点处的应力值。

平面应力问题中 $D=\frac{E}{1-\mathrm{\μ}}\left[\begin{array}{ccc}1& \mathrm{\μ}& 0\\ \mathrm{\μ}& 1& 0\\ 0& 0& \frac{1-\mathrm{\μ}}{2}\end{array}\right]---\left(26\right)$

其中，E为弹性模量，μ为泊松比。

分别计算出中心孔单元每个单元的4个高斯点的应力值，再用外插法计算出4个结点的应力值。

实施例1：

1个超单元结构

确定超单元的基本尺寸为0.3×0.2m；中心开孔长半径为r1＝0.07m，短半径为r2＝0.03m，采用solid182号单元；材料属性：弹性模量E＝2e11Pa，泊松比Pr＝0.3。建立有限元模型如图7所示，超单元的总刚点号如图8所示，超单元中心孔的单刚点号如图11所示。其中结点1～17全约束，结点18～23施加x方向载荷800N，结点24～34施加x方向载荷1000N；结点18～34施加y方向载荷1000N。

实施例2：

3个连续超单元结构

采用solid182号单元。材料属性：弹性模量E＝2e11Pa，泊松比Pr＝0.3。

中心椭圆孔半径从左向右：

第一个单元r1＝0.06m，r2＝0.04m；

第二个单元r1＝0.07m，r2＝0.05m；

第三个单元r1＝0.10m，r2＝0.06m。

有限元模型如图9所示，整体总刚点号如图10所示，超单元中心孔的单刚点号如图11所示，其中结点1～17全约束，结点52～60施加x方向载荷800N，结点61～68施加x方向载荷1000N，结点52～60施加y方向载荷800N，结点61～68施加y方向载荷1000N。

ANSYS计算结果和超单元计算结果如图12、图13、图14所示。

由结果分析可知，结果分析：使用常规算法，先要计算512个四边形单元的刚度矩阵，每个四边形单元是8×8的矩阵，得到每个刚度矩阵需要计算100次左右，然后再将512个单元刚度矩阵组合成总体刚度矩阵，共576个结点，需要计算G^TKG，每个单刚扩展成总刚的大小1152×1152的矩阵，再相加得到总体刚度矩阵，总体计算量大约10⁵的数量级。而使用超单元计算量大大减少，只需计算544次乘法，4次加法可以得到刚度矩阵；计算816次乘法，6次加法可以得到外部位移和内部位移的关系矩阵；计算1224次乘法，8次加法可以得到应变矩阵B，总计算量在10³数量级。通过以上计算量的分析，使用超单元的方法可以使速度提高100倍左右。

分别用ANSYS直接建立有限元模型计算和用超单元计算，并比较了计算误差(两者的相对误差和相对于最大值的误差)。计算结果可知，超单元的边界结点位移值、中心孔结点位移值的相对误差和相对最大值的误差都小于1％；中心孔处在应力值较大的地方相对误差较小，应力值较小的地方相对误差较大，相对于最大值的误差小于4％。这些误差值在误差要求的范围之内，满足精度的要求。

Claims (1)

1、一种参数化的静态超单元构造方法，具体为：

1)针对各种有重复结构的构件，将其重复结构的特征尺寸设为参数；

2)利用有限元软件对该重复结构建立不同特征尺寸的有限元模型，划分网格，分别读出每个结构的刚度矩阵；

3)利用有限元中子结构的方法，仅把关心部位的结点作为外部结点保留，做成超单元，得到超单元的刚度矩阵K；

4)对步骤3)中所得到各种模型的刚度矩阵K进行拟合，得到刚度矩阵K和其特征参数相关的函数；

5)对应变矩阵B进行拟合，得到应变矩阵B和其特征参数相关的函数；利用已求得的结点位移，进而计算出单元高斯点的应力值，再通过外插法得到结点的应力值。

Images