蜂窝结构非线性本构关系的建立方法、介质和设备与流程

本发明涉及蜂窝结构技术领域，特别涉及一种蜂窝结构非线性本构关系的建立方法、介质和设备。

背景技术：

蜂窝结构由于具有比强度大、比刚度高、隔热隔音性能好以及良好的能量吸收特性，是一种具有良好力学特性的轻量化结构。伴随着现代工业技术的发展，蜂窝结构由于其诸多优点而广泛应用于航空、航天、机械及交通等各个工程领域。典型的复合材料夹芯结构通常由上下两层薄面板和中间的一层芯体材料胶结而成。上下面板通常有较高的强度和刚度，分为纤维板、铝合金板、胶合板、玻璃钢板等，相比于蜂窝芯体结构，其通常较薄。蜂窝芯的网格形式有六边形，正方形，菱形，圆形，椭圆形等；根据材料的不同又可以分为纸、铝合金、nomex等蜂窝芯体结构。随着工业制造技术的不断发展，蜂窝结构的生存成本也在大幅下降，也使得蜂窝结构得以广泛应用于民用建筑和包装中。

蜂窝芯体的集合构型决定其具有各项异性的力学性能，在工程设计和分析中，复杂的几何形式和各项异性的力学性能给蜂窝芯体结构的设计与分析带来很大的苦难，通常的试验方法，耗时费力又不经济。随着计算科学的发展，采用通用有限元程序对大型复杂结构进行结构建模及力学仿真已成为现实，然而，这种有限元方法在结构系统建模过程中，详细的建模将大大增加模型规模和难度，使得计算效率显著降低。同时复杂多样的结构及边界条件，也使得模型具有一定的局限性。因此，工程技术人员往往希望将非连续的蜂窝芯体等效替代为连续的正交各向异性体，从而用材料等效模型和板壳单元来模拟蜂窝结构，达到降低建模难度并提高计算效率，但在精度方面有所降低。

技术实现要素：

本发明的第一目的在于克服现有技术的缺点与不足，提供一种蜂窝结构非线性本构关系的建立方法，该方法能够方便地建立非线性本构关系，所建立的非线性本构关系能够良好地体现蜂窝结构的非线性，有利于快速准确地分析蜂窝结构的力学特性。

本发明的第二目的在于提供一种存储介质。

本发明的第三目的在于提供一种计算设备。

本发明的第一目的通过下述技术方案实现：一种蜂窝结构非线性本构关系的建立方法，包括如下步骤：

针对于在不同方向上承受不同荷载和不同变形情况下的蜂窝结构代表胞元，根据大变形理论计算出对应的等效应力和等效应变；

根据等效应力和等效应变建立等效应力-等效应变关系，得到非线性等效材料参数；

计算出代表胞元的线性等效材料参数，根据非线性等效材料参数和线性等效材料参数，计算两者的商并作为非线性修正因子；

根据非线性修正因子和线性等效材料参数，建立同一类蜂窝结构的非线性本构关系。

优选的，针对于在不同方向上承受不同荷载和不同变形情况下的六边形代表胞元，根据大变形理论计算出对应的等效应力和等效应变，具体如下：

针对于在x左右方向、y上下方向上承受不同荷载和不同变形情况下的代表胞元，分别计算蜂窝结构的无量纲荷载ζ、在x方向上对应的位移量和以及在y方向上对应的位移量和

根据无量纲荷载ζ、在x方向上的位移量和以及在y方向上的位移量和计算出在x方向上承受荷载时，x方向上和y方向上的等效应力σ′x、σ′y和等效应变ε′x、ε′y，以及在y方向上承受荷载时，x方向上和y方向上的等效应力σ″x、σ″y和等效应变ε″x、ε″y，计算公式具体如下：

等效应力计算公式：

其中，p是施加在代表胞元上的等效荷载；b是蜂窝结构在z前后方向上的深度；l是代表胞元斜壁板的长度；t是蜂窝结构的壁板厚度；h是代表胞元竖壁板的长度；es是蜂窝结构的材料杨氏模量；

等效应变计算公式：

其中，φ是代表胞元斜壁板与y方向的夹角；是代表胞元在y方向上的投影长度，as是蜂窝结构壁板的横截面积；

根据等效应力和等效应变建立等效应力-等效应变关系，具体如下；

非线性等效材料参数为等效应力和等效应变的商，非线性等效材料参数包括：代表胞元在x方向上受荷载时的横向非线性等效杨氏模量ex和横向非线性等效泊松比νxy、在y方向上受荷载时的纵向非线性等效杨氏模量ey和纵向非线性等效泊松比νyx。

更进一步的，针对于在x、y方向上承受不同荷载和不同变形情况下的六边形代表胞元，分别计算蜂窝结构的无量纲荷载ζ、在x方向上对应的位移量和以及在y方向上对应的位移量和具体过程如下：

步骤1、针对于在x、y方向上承受不同荷载和变形情况下的代表胞元，根据梁的弯曲理论建立对应的半斜壁板的控制方程；

步骤2、将控制方程转化为无量纲方程，并对无量纲方程进行简化；

步骤3、针对于简化后的无量纲方程，通过第一类椭圆积分f(β)对无量纲荷载ζ进行表示；

通过第一类椭圆积分f(β)对蜂窝结构的斜壁板沿x方向的位移量进行表示；

通过第一类椭圆积分f(β)和第二类椭圆积分e(β)对蜂窝结构的斜壁板沿y方向的位移量进行表示。

更进一步的，在计算蜂窝结构的无量纲荷载ζ、在x方向上对应的位移量和以及在y方向上对应的位移量和的过程，还包括：判断斜壁板处于弹性变形阶段或者是塑性变形阶段，具体过程如下：

计算出蜂窝结构的斜壁板两个端点所承受的当前弯矩；

计算代表胞元斜壁板两个端点承受的极限弯矩mu：

其中，σs是蜂窝材料屈服强度；

判断当前弯矩与极限弯矩的大小，若当前弯矩小于极限弯矩，则未出现塑性转角β0，半斜壁板处于弹性变形阶段；

若当前弯矩大于极限弯矩，则出现塑性转角β0，半斜壁板处于塑性变形阶段；并在构建半斜壁板的控制方程和计算无量纲荷载和位移量时，将当前弯矩等于极限弯矩。

更进一步的，当蜂窝结构在y方向上承受荷载时：

(1)当荷载为拉荷载时，针对于半斜壁板的弹性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；i是蜂窝结构壁板的截面惯性矩；α是弧坐标s处的转角；

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；β是斜壁板中点的转角；η0是当弧坐标s＝0，即弧坐标s处的转角α＝0时的η值；

(2)当荷载为拉荷载时，针对于半斜壁板的塑性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；i是蜂窝结构壁板的截面惯性矩；α是弧坐标s处的转角；β0是塑性变形阶段所存在的塑性转角，β0由斜壁板端点处的弯矩mh确定：

其中，β是斜壁板中点的转角；

控制方程的边界条件为：

其中，mo是代表胞元斜壁板中点的弯矩；

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；η′0是当弧坐标s＝0，即弧坐标s处的转角α＝β0时的η值；

(3)当荷载为压荷载时，针对于半斜壁板的弹性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；i是蜂窝结构壁板的截面惯性矩；α是弧坐标s处的转角；

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；β是斜壁板中点的转角；η0是当弧坐标s＝0，即弧坐标s处的转角α＝0时的η值；

(4)当荷载为压荷载时，针对于半斜壁板的塑性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；i是蜂窝结构壁板的截面惯性矩；α是弧坐标s处的转角；β0是塑性变形阶段所存在的塑性转角，β0由斜壁板端点处的弯矩mh确定：

其中，β是斜壁板中点的转角；

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；η′0是当弧坐标s＝0，即弧坐标s处的转角α＝β0时的η值。

更进一步的，当蜂窝结构在x方向上承受荷载时：

(1)当荷载为拉荷载时，针对于半斜壁板的弹性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；i是蜂窝结构壁板的截面惯性矩；θ是代表胞元斜壁板与x方向的夹角，θ＝π/2-φ；α是弧坐标s处的转角；

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；β是斜壁板中点的转角；η0是当弧坐标s＝0，即弧坐标s处的转角α＝0时的η值；

(2)当荷载为拉荷载时，针对于半斜壁板的塑性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；i是蜂窝结构壁板的截面惯性矩；θ是代表胞元斜壁板与x方向的夹角，θ＝π/2-φ；α是弧坐标s处的转角，β0是塑性变形阶段所存在的塑性转角，β0由斜壁板端点处的弯矩mh确定：

其中，β是斜壁板中点的转角；

控制方程的边界条件为：

其中，mo是代表胞元斜壁板中点的弯矩；

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；η′0是当弧坐标s＝0，即弧坐标s处的转角α＝β0时的η值；

(3)当荷载为压荷载时，针对于半斜壁板的弹性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；i是蜂窝结构壁板的截面惯性矩；θ是代表胞元斜壁板与x方向的夹角，θ＝π/2-φ；α是弧坐标s处的转角；

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；β是斜壁板中点的转角；η0是当弧坐标s＝0，即弧坐标s处的转角α＝0时的η值；

(4)当荷载为压荷载时，针对于半斜壁板的塑性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；i是蜂窝结构壁板的截面惯性矩；θ是代表胞元斜壁板与x方向的夹角，θ＝π/2-φ；α是弧坐标s处的转角；β0是塑性变形阶段所存在的塑性转角，β0由斜壁板端点处的弯矩mh确定：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；β是斜壁板中点的转角；

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；η′0是当弧坐标s＝0，即弧坐标s处的转角α＝β0时的η值。

优选的，根据ginson公式计算出代表胞元的线性等效材料参数，公式具体如下：

其中，是代表胞元在x左右方向上承受荷载时的横向等效杨氏模量，是代表胞元在x方向上承受荷载时的横向等效泊松比，是代表胞元在y上下方向上承受荷载时的纵向等效杨氏模量，是代表胞元在在y方向上承受荷载时的纵向等效泊松比；

es是蜂窝结构的材料杨氏模量，t是蜂窝结构的壁板厚度，l是代表胞元斜壁板的长度，h是代表胞元竖壁板的长度，θ是代表胞元斜壁板与x方向的夹角；

根据非线性等效材料参数和线性等效材料参数计算两者的商，并将所计算出的商作为非线性修正因子，非线性修正因子包括：代表胞元在x方向上承受荷载时的等效杨氏模量修正因子λx和泊松比修正因子ωx、在y方向上承受荷载时的等效杨氏模量修正因子λy和泊松比修正因子ωy，具体如下：

其中，ζ为无量纲荷载，φ是代表胞元斜壁板与y方向的夹角，和分别为在x方向上承受荷载时沿x和y方向的位移量，和分别为在y方向上承受荷载时沿x和y方向的位移量，是代表胞元在y方向上的投影长度，当竖壁板在y方向上承受拉荷载时变形后，长度当竖壁板在y方向上承受压荷载时变形后，长度t是蜂窝结构的壁板厚度；

根据非线性修正因子和线性等效材料参数，建立蜂窝结构的非线性本构关系，所建立的非线性本构关系具体如下：

优选的，蜂窝结构的代表胞元的竖壁板垂直于代表胞元在x左右方向上所承受的荷载，竖壁板平行于代表胞元在y上下方向上所承受的荷载。

本发明的第二目的通过下述技术方案实现：一种存储介质，存储有程序，所述程序被处理器执行时，实现本发明第一目的所述的蜂窝结构非线性本构关系的建立方法。

本发明的第三目的通过下述技术方案实现：一种计算设备，包括处理器以及用于存储处理器可执行程序的存储器，所述处理器执行存储器存储的程序时，实现本发明第一目的所述的蜂窝结构非线性本构关系的建立方法。

本发明相对于现有技术具有如下的优点及效果：

(1)本发明蜂窝结构非线性本构关系的建立方法，包括步骤：首先对于在不同方向上承受不同荷载和不同变形情况下的蜂窝结构代表胞元，根据大变形理论计算出对应的等效应力和等效应变；根据等效应力和等效应变建立等效应力-等效应变关系，得到非线性等效材料参数；计算出代表胞元的线性等效材料参数，并计算非线性等效材料参数和线性等效材料参数的商，得到非线性修正因子；最终根据非线性修正因子和线性等效材料参数，建立起同一类蜂窝结构的非线性本构关系。本发明方法通过非线性等效材料参数和线性等效材料参数得到非线性修正因子，能够方便地建立起非线性本构关系，所建立的非线性本构关系考虑到了复杂多样的结构参数，能够良好地体现蜂窝结构的非线性，适用于模拟蜂窝结构进行分析，有利于快速准确地分析蜂窝结构的力学特性。

(2)本发明蜂窝结构非线性本构关系的建立方法中，分别针对于六边形代表胞元在x、y方向上承受压或者拉荷载时的弹性和塑性变形阶段进行非线性本构关系的构建，使得非线性本构关系能够充分考虑材料的弹塑性力学行为，即在不同方向上承受不同荷载和不同变形下的情况，使用非常方便，同时具有计算精度高和适用范围广的优点。

附图说明

图1是本发明蜂窝结构的代表胞元的示意图。

图2(a)是本发明蜂窝结构在y方向上承受拉荷载时的整体受力图。

图2(b)是图2(a)中代表胞元的受力示意图。

图2(c)是图2(b)中斜壁板的变形示意图。

图2(d)是图2(c)中的斜壁板不考虑塑性时的变形分析图。

图2(e)是图2(c)中的斜壁板考虑塑性时的变形分析图。

图3(a)是本发明蜂窝结构在y方向上承受压荷载时的整体受力图。

图3(b)是图3(a)中代表胞元的受力示意图。

图3(c)是图3(b)中斜壁板的变形示意图。

图3(d)是图3(c)中的斜壁板不考虑塑性时的变形分析图。

图3(e)是图3(c)中的斜壁板考虑塑性时的变形分析图。

图4是本发明等效应力和等效应变的关系示意图。

图5是本发明等效应变与非线性修正因子的关系示意图。

图6是本发明等效应力与非线性修正因子的关系示意图。

图7(a)至图7(d)是在不同大小的夹角φ下的非线性等效材料参数的结果示意图。

图8(a)至图8(d)是不同细长比t/l下的非线性等效材料参数的结果示意图。

图9(a)至图9(d)是不同蜂窝材料屈服强度σs下的非线性等效材料参数的结果示意图。

图10是本发明方法与有限元方法的等效应变对比图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

实施例1

本实施例公开了一种蜂窝结构非线性本构关系的建立方法，其中，蜂窝结构作为一种具有周期特性的各向异性材料，由一系列的周期网格组成，周期网格形状具体可以是六边形、正方形、菱形、圆形或者椭圆形等，周期网格的各个壁板厚度相同。蜂窝结构的代表胞元等效为蜂窝结构中的任意一个网格，非线性本构关系的建立方法包括步骤如下：

针对于在不同方向上承受不同荷载和不同变形情况下的蜂窝结构代表胞元，根据大变形理论计算出对应的等效应力和等效应变。

根据等效应力和等效应变建立等效应力-等效应变关系，得到非线性等效材料参数。

计算出代表胞元的线性等效材料参数，根据非线性等效材料参数和线性等效材料参数，计算两者的商并作为非线性修正因子。

根据非线性修正因子和线性等效材料参数，建立同一类蜂窝结构的非线性本构关系。

其中，如图1所示，本实施例的周期网格形状为六边形。如图2(a)和图3(a)所示，建立全局直角坐标系xyz，x轴方向为左右方向，y轴方向为上下方向，z轴方向为前后方向，本实施例的不同受力方向具体分为x方向和y方向。代表胞元承受的不同荷载分为拉荷载和压荷载，受力变形分为弹性变形阶段和塑性变形阶段。当代表胞元在x方向上承受荷载时，代表胞元的竖壁板hg垂直于代表胞元在x方向上承受的荷载；当代表胞元在y方向上承受荷载时，代表胞元的竖壁板hg平行于代表胞元在y方向上承受的荷载。

由于六边形代表胞元的形状和承受的荷载具有对称性，因此将代表胞元的竖壁板厚度设置为斜壁板厚度的一半，且针对于代表胞元的受力变形分析可转化为对半斜壁板ho或者bo的受力变形分析。

针对于在不同方向上承受不同荷载和不同变形情况下的六边形代表胞元，根据大变形理论计算出对应的等效应力和等效应变，具体如下：

步骤1、针对于在x、y方向上承受不同荷载和不同变形情况下的代表胞元，分别计算蜂窝结构的无量纲荷载ζ、在x方向上对应的位移量和以及在y方向上对应的位移量和具体过程如下：

步骤11、针对于在x、y方向上承受不同荷载和变形情况下的代表胞元，根据梁的弯曲理论建立对应的半斜壁板的控制方程；

步骤12、将控制方程转化为无量纲方程，并对无量纲方程进行简化；

步骤13、针对于简化后的无量纲方程，通过第一类椭圆积分f(β)对无量纲荷载ζ进行表示；

通过第一类椭圆积分f(β)对蜂窝结构的斜壁板沿x方向的位移量进行表示；

通过第一类椭圆积分f(β)和第二类椭圆积分e(β)对蜂窝结构的斜壁板沿y方向的位移量进行表示。

在上述计算无量纲荷载和位移量的过程中，还包括：判断斜壁板处于弹性变形阶段或者是塑性变形阶段，具体过程如下：

计算出蜂窝结构的斜壁板两个端点h点、b点承受的当前弯矩；

计算代表胞元斜壁板两个端点承受的极限弯矩mu：

其中，b是蜂窝结构在z轴方向上的深度，σs是蜂窝材料屈服强度；t是蜂窝结构的壁板厚度；

判断当前弯矩与极限弯矩的大小，若当前弯矩小于极限弯矩，则未出现塑性转角β0，斜壁板处于弹性变形阶段；

若当前弯矩大于极限弯矩，则说明斜壁板所承受的荷载已经超过极限，导致出现塑性转角β0，斜壁板处于塑性变形阶段。

在塑性变形阶段构建半斜壁板的控制方程和计算无量纲荷载和位移量时，再将当前弯矩调整等于极限弯矩。

所构建的半斜壁板的控制方程和计算得到的无量纲荷载和位移量分别如下：

当蜂窝结构在y方向上承受荷载时：

(1)当荷载为拉荷载时，如图1、图2(a)至图2(d)所示，针对于半斜壁板的弹性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；es是蜂窝结构的材料杨氏模量；i是蜂窝结构壁板的截面惯性矩；p是施加在代表胞元上的等效荷载；φ是代表胞元斜壁板与y方向的夹角；α是弧坐标s处的转角。

接着，将控制方程转化为无量纲方程：

其中，ζ是无量纲荷载，pcr是斜壁板弯曲时的临界荷载，s是无量纲弧坐标，s＝s/l(0≤s≤0.5)；

控制方程的边界条件为：mo＝0，mo是蜂窝胞元斜板中点的弯矩。

将无量纲方程进行简化，得到：

其中，β是斜壁板中点o点的转角。

再根据简化后的无量纲方程，得到无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；β是斜壁板中点o点的转角；η0是当弧坐标s＝0，即弧坐标s处的转角α＝0时的η值。

(2)当荷载为拉荷载时，如图1、图2(a)至图2(c)和图2(e)所示，针对于半斜壁板的塑性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；es是蜂窝结构的材料杨氏模量；i是蜂窝结构壁板的截面惯性矩；p是施加在代表胞元上的等效荷载；φ是代表胞元斜壁板与y方向的夹角；α是弧坐标s处的转角；β0是塑性变形阶段所存在的塑性转角，β0由斜壁板端点处的弯矩mh确定：

其中，β是斜壁板中点o点的转角；

控制方程的边界条件为：

其中，mo是代表胞元斜壁板中点o点的弯矩。

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；η′0是当弧坐标s＝0，即弧坐标s处的转角α＝β0时的η值。

(3)当荷载为压荷载时，如图1、图3(a)至图3(d)所示，针对于半斜壁板的弹性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；es是蜂窝结构的材料杨氏模量；i是蜂窝结构壁板的截面惯性矩；p是施加在代表胞元上的等效荷载；φ是代表胞元斜壁板与y方向的夹角；α是弧坐标s处的转角。

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；β是斜壁板中点o点的转角；η0是当弧坐标s＝0，即弧坐标s处的转角α＝0时的η值。

(4)当荷载为压荷载时，如图1、图3(a)至图3(c)和图3(e)所示，针对于半斜壁板的塑性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；es是蜂窝结构的材料杨氏模量；i是蜂窝结构壁板的截面惯性矩；p是施加在代表胞元上的等效荷载；φ是代表胞元斜壁板与y方向的夹角；α是弧坐标s处的转角；β0是塑性变形阶段所存在的塑性转角，β0由斜壁板端点处的弯矩mh确定：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；β是斜壁板中点o点的转角。

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；η′0是当弧坐标s＝0，即弧坐标s处的转角α＝β0时的η值。

当蜂窝结构在x方向上承受荷载时：

(1)当荷载为拉荷载时，如图1所示，针对于半斜壁板的弹性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；es是蜂窝结构的材料杨氏模量；i是蜂窝结构壁板的截面惯性矩；p是施加在代表胞元上的等效荷载；θ是代表胞元斜壁板与x方向的夹角，θ＝π/2-φ；α是弧坐标s处的转角。

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；β是斜壁板中点o点的转角；η0是当弧坐标s＝0，即弧坐标s处的转角α＝0时的η值。

(2)当荷载为拉荷载时，如图1所示，针对于半斜壁板的塑性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；es是蜂窝结构的材料杨氏模量；i是蜂窝结构壁板的截面惯性矩；p是施加在代表胞元上的等效荷载；θ是代表胞元斜壁板与x方向的夹角，θ＝π/2-φ；α是弧坐标s处的转角，β0是塑性变形阶段所存在的塑性转角，β0由斜壁板端点处的弯矩mh确定：

其中，β是斜壁板中点o点的转角；

控制方程的边界条件为：

其中，mo是代表胞元斜壁板中点o点的弯矩。

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；η′0是当弧坐标s＝0，即弧坐标s处的转角α＝β0时的η值。

(3)当荷载为压荷载时，如图1所示，针对于半斜壁板的弹性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；es是蜂窝结构的材料杨氏模量；i是蜂窝结构壁板的截面惯性矩；p是施加在代表胞元上的等效荷载；θ是代表胞元斜壁板与x方向的夹角，θ＝π/2-φ；α是弧坐标s处的转角。

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；β是斜壁板中点o点的转角；η0是当弧坐标s＝0，即弧坐标s处的转角α＝0时的η值。

(4)当荷载为压荷载时，如图1所示，针对于半斜壁板的塑性变形阶段，所建立的半斜壁板的控制方程为：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；es是蜂窝结构的材料杨氏模量；i是蜂窝结构壁板的截面惯性矩；p是施加在代表胞元上的等效荷载；θ是代表胞元斜壁板与x方向的夹角，θ＝π/2-φ；α是弧坐标s处的转角；β0是塑性变形阶段所存在的塑性转角，β0由斜壁板端点处的弯矩mh确定：

其中，esi是蜂窝结构斜壁板截面的抗弯刚度；β是斜壁板中点o点的转角。

无量纲荷载ζ、斜壁板沿x方向的位移量和y方向的位移量分别如下：

其中，

f(β)是第一类椭圆积分；e(β)是第二类椭圆积分；η是根据坐标变换公式产生的新变量；η′0是当弧坐标s＝0，即弧坐标s处的转角α＝β0时的η值。

步骤2、根据无量纲荷载ζ、在x方向上的位移量和以及在y方向上的位移量和计算出在x方向上承受荷载时，x方向上和y方向上的等效应力σ′x、σ′y和等效应变ε′x、ε′y，以及在y方向上承受荷载时，x方向上和y方向上的等效应力σ″x、σ″y和等效应变ε″x、ε″y，计算公式具体如下：

等效应力计算公式：

其中，l是代表胞元斜壁板的长度；h是代表胞元竖壁板的长度；

等效应变计算公式：

其中，是代表胞元在y方向上的投影长度，as是蜂窝壁板的横截面积。

如图4所示，本实施例六边形代表胞元的等效应力-等效应变关系具体如下；

由上述等效应力-等效应变关系可得，非线性等效材料参数为等效应力和等效应变的商，其中，非线性等效材料参数包括：代表胞元在x方向上承受荷载时的横向非线性等效杨氏模量ex和横向非线性等效泊松比νxy、在y方向上承受荷载时的纵向非线性等效杨氏模量ey和纵向非线性等效泊松比νyx。

在本实施例中，代表胞元的线性等效材料参数具体是根据ginson公式计算得到，ginson公式具体如下：

其中，是代表胞元在x方向上承受荷载时的横向等效杨氏模量，是代表胞元在x方向上承受荷载时的横向等效泊松比，是代表胞元在y方向上承受荷载时的纵向等效杨氏模量，是代表胞元在在y方向上承受荷载时的纵向等效泊松比。

根据非线性等效材料参数和线性等效材料参数计算得到的非线性修正因子λx、ωx、λy、ωy具体如下：

其中，λx是代表胞元在x方向上承受荷载时的等效杨氏模量修正因子；ωx是泊松比修正因子；λy是在y方向上承受荷载时的等效杨氏模量修正因子；ωy是泊松比修正因子。

当竖壁板在y方向上承受拉荷载时变形后，投影长度当竖壁板在y方向上承受压荷载时变形后，投影长度

根据非线性修正因子和线性等效材料参数，建立同一类蜂窝结构的非线性本构关系，所建立的非线性本构关系具体如下：

在本实施例的蜂窝结构中，壁板厚度t＝1mm；材料杨氏模量es＝12mpa，材料屈服强度σs＝10mpa；蜂窝结构在z轴方向上的深度b＝1mm。如图1所示，在本实施例的六边形代表胞元中，斜壁板hb与y方向的夹角φ＝60°；斜壁板的长度l＝10mm，斜壁板的厚度为1mm；竖壁板hg的长度h＝10mm，竖壁板的厚度为0.5mm。

如图4所示，正应变表示的是蜂窝结构承受的拉荷载，负应变表示蜂窝结构承受的压荷载。当正应变逐渐增大到b’点或者d’点时，蜂窝结构从弹性变形阶段开始进入塑性变形阶段；当负应变逐渐增大到a’点或者c’点时，蜂窝结构从弹性变形阶段开始进入塑性变形阶段。因此，点a’、b’、c’和d’均表示蜂窝结构开始进入塑性变形，也即是出现了塑性转角。从曲线的斜率可以看出，随着塑性转角的出现，斜率变小，说明相同应变增量下所增加的应力变小，因此在塑性变形阶段，蜂窝结构的承载能力有所下降，故图4至图6中的曲线符合实际情况，非线性修正因子能够体现等效应力和等效应变之间的关系。

图7(a)至图7(d)表示的是在不同大小的夹角φ下的非线性等效材料参数随着转角β增大的变化情况，其中，夹角φ分别为：30°、45°、60°、75°；转角β的增大也即是代表胞元变形程度的变大。从图7(a)至图7(d)可以看出：

1、夹角φ越大，x方向杨氏模量ex越大，y方向杨氏模量ey越小，且x方向杨氏模量ex和y方向杨氏模量ey均随着代表胞元变形程度的变大而变大，x方向杨氏模量ex和y方向杨氏模量ey均表现出极强的非线性。

2、夹角φ越大，x方向泊松比νxy越大，y方向泊松比νyx越小，且x方向泊松比νxy和y方向泊松比νyx均随着代表胞元变形程度的变大而变大，x方向泊松比νxy和y方向泊松比νyx均表现出极强的非线性。

图8(a)至图8(d)表示的是不同细长比t/l下的非线性等效材料参数随着转角β增大的变化情况，其中，t/l分别为：0.05、0.1、0.2、0.4、0.6，图8(c)和图8(d)中的小曲线图为虚线框的放大图。从图8(a)至图8(d)可以看出：

1、细长比t/l越大，x方向杨氏模量ex越大，y方向杨氏模量ey越大，且两者均随着代表胞元变形程度的变大而变大。

2、细长比t/l越大，x方向泊松比νxy和y方向泊松比νyx均先变小随后变大，且两者均随着代表胞元变形程度的变大而变大。

图9(a)至图9(d)表示的是不同蜂窝材料屈服强度σs下的非线性等效材料参数随着转角β增大的变化情况，其中，σs分别为：3mpa，5mpa，8mpa，10mpa和15mpa，图9(d)中的小曲线图为虚线框的放大图。从结果可以看出：

1、在弹性阶段，蜂窝材料屈服强度σs对结构参数没有影响；

2、蜂窝材料屈服强度σs越大，代表胞元越晚出现塑性变形，导致x方向杨氏模量ex和y方向杨氏模量ey均变大。

3、蜂窝材料屈服强度σs越大，x方向泊松比νxy和y方向泊松比νyx均变小，且两者均随着代表胞元变形程度的变大而变大，且y方向的泊松比变大效果明显强于x方向的泊松比变大效果。

图10为在本实施例方法与有限元方法下的等效应变随着转角β增大的变化情况，从图10可以看出，本实施例方法的等效应变曲线和有限元方法的等效应变曲线大体一致，说明本实施例中的蜂窝结构与相对精细的有限元蜂窝结构吻合良好。

综合来说，从图4至图10可得，本实施例方法所构建的非线性本构关系考虑到了复杂多样的结构参数，能够良好地体现蜂窝结构的非线性，所等效的代表胞元具有良好的代表性，能够用于模拟蜂窝结构进行分析。

实施例2

本实施例公开了一种存储介质，存储有程序，所述程序被处理器执行时，实现实施例1中所述的蜂窝结构非线性本构关系的建立方法，具体如下：

针对于在不同方向上承受不同荷载和不同变形情况下的蜂窝结构代表胞元，根据大变形理论计算出对应的等效应力和等效应变；

根据等效应力和等效应变建立等效应力-等效应变关系，得到非线性等效材料参数；

计算出代表胞元的线性等效材料参数，根据非线性等效材料参数和线性等效材料参数，计算两者的商并作为非线性修正因子；

根据非线性修正因子和线性等效材料参数，建立同一类蜂窝结构的非线性本构关系。

其中，代表胞元的形状具体可以是六边形、正方形、菱形、圆形或者椭圆形等，代表胞元可以等效为蜂窝结构中的任意一个网格。

本实施例中的存储介质可以是磁盘、光盘、计算机存储器、只读存储器(rom，read-onlymemory)、随机存取存储器(ram，randomaccessmemory)、u盘、移动硬盘等介质。

实施例3

本实施例公开了一种计算设备，包括处理器以及用于存储处理器可执行程序的存储器，所述处理器执行存储器存储的程序时，实现实施例1中所述的蜂窝结构非线性本构关系的建立方法，具体如下：

针对于在不同方向上承受不同荷载和不同变形情况下的蜂窝结构代表胞元，根据大变形理论计算出对应的等效应力和等效应变；

根据等效应力和等效应变建立等效应力-等效应变关系，得到非线性等效材料参数；

计算出代表胞元的线性等效材料参数，根据非线性等效材料参数和线性等效材料参数，计算两者的商并作为非线性修正因子；

根据非线性修正因子和线性等效材料参数，建立同一类蜂窝结构的非线性本构关系。

其中，代表胞元的形状具体可以是六边形、正方形、菱形、圆形或者椭圆形等，代表胞元可以等效为蜂窝结构中的任意一个网格。

本实施例中所述的计算设备可以是台式电脑、笔记本电脑、智能手机、pda手持终端、平板电脑或其他具有处理器功能的终端设备。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。