CN106021644B - 确定混合维模型界面约束方程系数的方法 - Google Patents

Abstract

一种确定混合维模型界面约束方程系数的方法，截取悬臂梁交接面两侧的梁建立静力分析模型，模型边界约束条件同原模型，根据多点约束方程推导所需要的应力分布，进行相关静力分析，并通过有限元软件导出交接面截面处有限元模型节点的应力结果。本发明同时考虑结构总体的振动模态特性和结构局部的精细动态响应计算需要，能够建立多点约束方程的统一求解流程，快速完成多点约束方程求解，方便混合维模型在工程上的应用。本发明形成了应力协调多点约束方程的统一求解步骤，便于统一程序编写，利于在工程应用中的推广，并且采用本发明方法建立的有限元模型与真实模型的前10阶模态频率误差在1％以内，得到的频率响应曲线与真实模型的频率响应曲线变化一致。

Description

确定混合维模型界面约束方程系数的方法

技术领域

本发明属于结构动力学建模领域，具体是一种确定混合维模型界面约束方程系数的方法。

背景技术

对于大型结构，直接对整体结构建立精细的动力学有限元模型进行仿真计算，常常会遇到模型计算规模过大，严重影响有限元分析的效率，甚至使得分析无法进行的情况。采用不同维度单元有针对性的建模，是行之有效的方法，但不同维度单元的连接界面常常需要特殊处理。目前发展的多种方法中，过渡单元方法发展最早，但其推导过程不具有统一性，而且使用中容易锁死，并会产生多余的非物理模态；而近几年发展的使用刚性梁单元完成不同维度单元界面连接的处理的方法，操作简单，但物理意义并不明确，且缺乏理论基础。相较于前两种方法，McCune在International Journal for Numerical Method inEngineering发表的Mixed dimensional coupling in finite element models和Shim在Engineering with Computers发表的Mixed dimensional coupling in finite elementmodels提出的使用交接面功相等原理得到的考虑应力协调的多点约束方程，物理意义明确，各种类型单元连接时方程的推导具有统一性，便于与商用有限元软件对接。但对于非规则结构，或任意截面形式结构，难以获得结构在不同作用力形式下的应力分布解析解，故而难以确定所需的多点约束方程系数。

工程实践表明，结构的振动疲劳破坏都是发生在整体结构的某一个局部区域，即局部结构上，因此要想准确预测结构的振动疲劳寿命，需要提高局部结构动态应力的计算精度。当激励力同时作用在局部结构和周围结构上或只作用在周围结构上时，局部结构的振动响应分析必须计及周围结构振动所传入的振动能量，不可以把局部结构单独分离出来进行建模和计算，必须建立整体结构模型进行分析。此时，采用低维、大尺度单元网格进行周围结构的离散，以最小的计算成本实现振动能量整体传递效应的模拟，而对结构振动疲劳问题所关注的局部结构，则采用局部精细网格或高精度单元进行离散建模，并将周围结构有限元模型和局部结构模型直接采用多点约束方程耦合起来进行同步求解，可以大大减小计算的规模，提高计算的效率。本发明将建立多点约束方程的统一求解流程，快速完成多点约束方程求解，方便混合维模型在工程上的应用。Blanco等在Computer methods inapplied mechanics and engineering发表的A Aariational approach for couplingkinematically incompatible structural models进一步发展了一种变分框架，用于完成运动学不协调单元的耦合。东南大学李兆霞教授在《固体力学学报》发表的《桥梁结构劣化与损伤过程的多尺度分析方法及其应用》中与《东南大学学报》发表的《结构损伤一致多尺度模拟和分析方法》中对多尺度模拟与分析方法进行了研究，通过试验与仿真分析对比说明了子结构方法和多点约束方法建立的多尺度有限元模型可以用来进行局部热点应力分析，但其工作主要集中在子结构中，对于混合维结构研究较少，而混合维结构为本发明的主要研究工作，本发明优势在于用静力分析得到应力分布，不再进行公式计算。检索国内外相关专利，西安电子科技大学的脊柱的三维几何与有限元混合模型的构建方法(专利号：201110114628.0)主要是对脊柱进行三维成像，并为进行动力学建模。广州汽车集团股份有限公司的一种锁止机构的刚柔体混合建模方法(专利号：201410076328.1)，其混合维主要是刚性体和柔性体的混合维建模，但并未涉及有限元中不同自由度单元间的建模问题。

发明内容

为克服现有技术中存在的多尺度有限元模型中对混合维结构研究的不足，本发明提出了一种确定混合维模型界面约束方程系数的方法。

本发明的具体步骤如下：

步骤一：划分管梁的局部结构与周围结构。

所述的管梁为悬臂梁；以该管梁的固支端部分为局部结构；该管梁的其余部分为周围结构。

对所述管梁进行混合维建模；建模时，所述局部结构采用高阶精确单元建立有限元模型，所述周围结构采用低阶减缩单元模拟刚度效应。

将所述有限元模型中发生不同维单元连接的截面定义为交接面A₁。根据有限元模型中的交接面确定管梁的交接面，并使管梁的交接面A₂与有限元模型中的交接面A₁位置相同。

步骤二：以交接面为中心，向两侧分别截取管梁的几何模型。

以交接面A₂为中心截取该交接面两侧的管梁作为几何模型，建立该几何模型的静力分析有限元模型。所述静力分析有限元模型用壳单元建立。

对所述静力分析有限元模型建立边界约束条件。所述静力分析有限元模型的边界约束条件与管梁的边界约束条件相同，将该静力分析有限元模型中的一端视为固支端，静力分析有限元模型中的交接面A₃的位置与管梁的交接面A₂的位置相同。以所述交接面A₃为中心，靠近固支端一侧为局部结构，靠近自由端一侧为周围结构。

所述截取的几何模型的长度需使边界条件不会影响交接面A₂上的应力分布，截取的几何模型的长度>该管梁横截面的最大外形尺寸的3倍。

步骤三：对静力分析有限元模型进行静力分析。

在所述静力分析有限元模型自由端处施加不同方向的任意值的载荷，使该静力分析有限元模型的自由端处分别产生轴力F、弯矩M、扭矩T以及剪力Fs。对所述静力分析有限元模型赋予管梁的材料属性；所述管梁的材料属性包括弹性模量、泊松比和密度。通过有限元分析软件对所述有限元模型进行静力分析。

所述静力分析有限元模型自由端处施加的弯矩M有两个元素M₁与M₂，并且弯矩M₁的方向与坐标系的y向同向；M₂的方向与坐标系的z向同向；所述的剪力Fs有两个元素Fs₁与Fs₂，并且剪力Fs₁的方向与坐标系的y向同向，Fs₂的方向与坐标系的z向同向。

步骤四：分别提取交接面A₃上的正应力与切应力。

采用数值离散方法，获得交接面A₃上在轴力F、弯矩M、扭矩T以及剪力Fs作用下的离散应力分布。

从步骤三得到的静力分析结果中，分别提取位于交接面A₃上静力分析有限元模型的各节点所对应的沿坐标系的x、y、z三轴的正应力与沿该坐标系x、y、z三轴的切应力。

步骤五：确定交接面A₁上梁单元的节点自由度与交接面A₁上壳单元的节点自由度之间的约束方程。

建立梁结构的混合维有限元模型。在建立的混合维有限元模型上划分壳单元，并以该壳单元作为管梁的局部结构；在建立的混合维有限元模型上划分梁单元，并以该梁单元作为管梁的周围结构。

分别建立梁结构混合维有限元模型的梁单元节点的六个自由度的约束条件。

分别将交接面两侧的单元做功显示表示出来，利用交接面A₁两侧做功大小相同原理，通过公式(6)确定所述交接面A₁梁单元的节点自由度与所述交接面A₁壳单元的节点自由度之间的约束方程：

式中，F表示作用在结构上的广义的力，U即为交接面靠近自由端侧梁单元节点位移。

式中，F表示作用在结构上的广义的力，U为交接面靠近自由端侧梁单元节点位移；[N_d]_j为积分子域i上高斯积分点j对应的单元形状函数矩阵取值,j表示积分子域i上选定的高斯积分点的个数，根据实际积分需要的精确度进行合理的选取,下标d无具体含义.σ'表示应力分布函数，x’，y’，z’分别表示局部坐标系的x，y，z轴，[T_d]为界面局部坐标系下的位移向量到结构总体坐标系的位移向量之间的转换关系矩阵，为交接面靠近固支端侧模型壳单元节点位移。

所述确定交接面A₁上梁单元的节点自由度与交接面A₁上壳单元的节点自由度之间的约束方程的具体过程是：

建立整体结构的混合维模型，完成交接面两侧应力做功大小的推导。

对于离散的有限元模型，在局部坐标系下，交接面上应力做功大小表示为：

其中，W_faceA表示为交接面上应力做功大小，faceA表示交接面；N表示交接面上有限元模型单元个数，A_i为不同维单元i交接面上的面积分变量。u'_i(x'_A,y'_A,z'_A)、σ_i'(x'_A,y'_A,z'_A)d中A_i：u_i'和σ_i'分别表示在单元局部坐标系下交接面上对应的位移及相应方向的应力分布函数，x’，y’，z’分别表示局部坐标系的x，y，z轴，下标A表示交接面。

以单元形状函数[N_d]取代公式(1)中的u'_i，即能够将单元局部坐标系下的位移函数用单元节点位移表示：

其中，表示在局部坐标系下界面积分子域i所包含离散结点的位移所构成的位移列向量，下标e无具体含义；[N_d]表示界面积分子域i对应的位移插值形状函数矩阵，下标d无具体含义，有限元中单元的形状函数常在单元局部坐标系下定义。

为了约束方程的添加，需要进一步利用界面局部坐标系下的位移向量到结构总体坐标系的位移向量之间的转换关系矩阵[T_d]，下标d无具体含义。将式(2)中局部坐标系下的位移向量向总体坐标系进行转换：

为了完成式(3)的积分，将已知条件结构交接面上相关节点的应力值代入，利用高斯积分将(3)式的面积分转换为高斯点的应力值叠加：

式中为积分子域i对应的应力做功；j表示积分子域i上选定的高斯积分点的个数，根据实际积分需要的精确度进行合理的选取；[N_d]_j为积分子域i上高斯积分点j对应的单元形状函数矩阵取值；σ'_i(x'_j,y'_j,z'_j)为积分子域i上高斯积分点j对应的具体的应力大小。积分子域内的高斯点的应力值可通过已知节点的应力值插值得到：

σ′_i(x′_j，y′_j，z′_j)＝([N_g]×{σ′_n}_j) (5)

式中插值函数[N_g]的构造与划分结构的单元形状函数无关，根据实际选取的插值节点的个数确定。{σ'_n}_j表示积分子域i上高斯积分点j所包含离散结点的应力所构成的应力列向量，坐标系为界面局部坐标系。所述相关应力值对应为单元局部坐标系方向上的应力值。

分别将交接面两侧的单元做功显示表示出来，利用交接面A₁两侧做功大小相同原理，即得到梁单元离散侧交接面节点自由度与壳单元离散侧交接面节点自由度之间的约束方程。

步骤六：确定交接面的约束方程系数。

通过界面约束方程的通用数值构造方法确定交接面A₁的约束方程系数。

对于沿梁结构混合维有限元模型的梁单元轴向平动自由度的约束条件，将轴力以及对应的节点应力带入公式(5)和(6)中，公式(6)中有限元模型单元的形状函数[N_d]_j以及单元局部坐标系与结构总体坐标系的坐标转换矩阵都能够根据梁结构混合维有限元模型通过有限元分析方法获得，公式(6)中未知量只有U与其中U即为交接面靠近自由端侧梁单元节点位移，即为交接面靠近固支端侧模型壳单元节点位移，即建立了交接面靠近自由端侧梁单元轴向自由度与交接面靠近固支端侧壳单元节点自由度的约束关系。同理，分别计算得到梁单元各方向自由度的约束关系，并确定交接面处各点的混合维模型界面约束方程系数。

本发明提出一种同时考虑结构总体的振动模态特性和结构局部的精细动态响应计算的计算方法，能够建立多点约束方程的统一求解流程，快速完成多点约束方程求解，方便混合维模型在工程上的应用。

本发明中，在划分管梁的局部结构与周围结构时，为了使交接面两侧模型的应力传递关系正确，需要建立应力协调的多点约束方程。方程建立过程需要特定力作用下，交接面所处位置截面上的应力分布情况。对于复杂截面，特定力作用下的应力分布公式往往难以获得。但在线性分析范围内，应力与外力之间总存在有线性关系。因此采用数值离散思路，利用与交接面截面形状相同的结构，在特定力作用下进行静力分析，获得交接面截面上在特定力作用下的离散应力分布。为了获得高的计算效率，并不需要采用原结构进行静力分析，只需截取包括交接面的部分结构即可获得外力与交接面应力的分布关系。对于图1所示梁结构，截取矩形框内结构建立静力分析模型，模型边界约束条件同原模型，根据多点约束方程推导所需要的应力分布，进行相关静力分析，并通过有限元软件导出交接面截面处有限元模型节点的应力结果。

本发明中，界面积分子域的插值形状函数矩阵[N_d]是为了建立已知离散节点参量与任意位置参量取值之间的插值函数关系。[N_d]的具体构造形式可随界面积分子域的形状而变化。如果界面积分子域是由4个离散结点所构成的四边形区域，则可借用四边形平面应力单元构造中类似的形状函数矩阵来类似构造插值函数矩阵。若界面积分子域选择用3个离散结点所围成的三角形区域，则可借鉴平面应力三角元的形状函数矩阵来类似构造插值函数矩阵。如果交接界面退化成了壳元所表征的形式，则界面积分也就相应退化为线积分。因此，原理上，交接界面插值函数矩阵的构造，根据界面积分子域的选择不同进行构造。不过从工程使用的一般性要求出发，基于三角形区域构造面积分插值应该是最通用的。

与现有技术比较，本发明的取得的有益效果如下：

1、针对应力协调多点约束方程推导过程中，复杂截面应力分布公式难以确定的问题，根据线性分析特点，提出采用局部简单结构进行静力分析，通过分析得到的节点应力值近似模拟截面上的应力分布。

2、在数值模拟截面应力分布基础上，进一步提出了交接面内力功的通用数值计算方法，针对任意截面形状，形成了应力协调多点约束方程的统一求解步骤，便于统一程序编写，利于应力协调多点约束方程在工程应用的推广。

3、采用本发明方法建立的有限元模型与真实模型的前10阶模态频率误差在1％以内，如表1所示，并且采用本发明方法建立的有限元模型计算所得的频率响应曲线与真实模型的频率响应曲线变化一致，如图7所示，由此可见，本发明提出的混合维模型界面约束方程的通用数值构造方法可以应用到实际工程应用中。

附图说明

图1是结构静力分析模型建立示意图；

图2是梁式结构的任意混合维单元交接面示意图；

图3是管梁示意图；

图4是梁结构静力分析模型；

图5是梁结构混合维有限元模型；

图6是梁结构整体细化模型。

图7是两种模型的频率响应曲线。

图8是本发明的流程图。图中：

1.三维实体单元；2.梁单元；3.管梁；4.梁结构混合维有限元模型的壳单元；5.混合维模型频率响应曲线；6.全壳模型频率响应曲线；7.梁结构混合维有限元模型的梁单元；8.梁结构静力分析模型中壳单元；V₁.梁结构混合维有限元模型的交接面；V₂.管梁示意图的交接面；V₃.梁结构静力分析模型；

具体实施方式

本发明提供一种兼顾总体动特性时局部结构的精细动态响应计算方法，能够建立多点约束方程的统一求解流程，快速完成多点约束方程求解，方便混合维模型在工程上的应用，主要用于考虑局部结构无法分离时，建立表征周围结构的低维简化网格模型与表征局部结构的高维精细网格模型相一体的混合维动力学模型。

本实施例是采用混合维模型界面约束方程对一个单端固支的管梁3实施通用数值构造的过程，该管梁的截面为环形。管梁中部为交接面A₂，图3中的方框内即为所述交接面A₂。所述管梁的外径为105mm，内径95mm，圆管长度为4000mm。材料属性为：弹性模量210GPa，泊松比0.3，密度7800kg/m³。

具体实施步骤如下：

步骤一：划分管梁的局部结构与周围结构。

对于图3所示管梁，一端固支，另一端自由，中间为交接面。假设动力学分析关心的热点区域为该管梁的后半段区域，定义该区域为局部结构，该管梁的其余部分则定义为周围结构。设管梁3靠近固支端部分为局部结构，其余部分为周围结构。进行混合维建模时，所述局部结构采用高阶精确单元建立有限元模型，所述周围结构采用低阶减缩单元模拟刚度效应，如图5所示。将该有限元模型中发生不同维单元连接的截面定义为交接面A₁。根据有限元模型中的交接面确定管梁的交接面，并使管梁的交接面A₂与有限元模型中的交接面A₁位置相同。

步骤二：以交接面为中心，向两侧分别截取管梁的几何模型。

对于图3所示的管梁3，以交接面A₂为中心截取该交接面两侧的管梁作为几何模型，建立该几何模型的静力分析有限元模型。所述静力分析有限元模型用壳单元8建立，如图4所示。所截取的几何模型的长度需使边界条件不会影响交接面A₂上的应力分布，截取的几何模型的长度>该管梁横截面的最大外形尺寸的3倍。本实施例中，管梁为环形，所述的最大外形尺寸即为该管梁的外径。

为了提高计算效率，在截取几何模型时，该几何模型的长度尽可能短。

对所述静力分析有限元模型建立边界约束条件。所述静力分析有限元模型的边界约束条件与管梁3的边界约束条件相同，如图4所示，将该静力分析有限元模型中的一端视为固支端。图4中，中间方框为静力分析有限元模型中的交接面A₃，所述交接面A₃的位置与管梁的交接面A₂的位置相同。以所述交接面A₃为中心，靠近固支端一侧为局部结构，靠近自由端一侧为周围结构。

本实施例中，截取的几何模型的长度为该管梁长度的1/2，即管梁外径的20倍。

步骤三：对静力分析有限元模型进行静力分析。

如图4所示，在所述静力分析有限元模型自由端处，施加不同方向的任意值的载荷，使该静力分析有限元模型的自由端处分别产生轴力F、弯矩M、扭矩T以及剪力Fs。所述的弯矩M有两个元素M₁与M₂，并且弯矩M₁的方向与坐标系的y向同向，M₂的方向与坐标系的z向同向。所述的剪力Fs有两个元素Fs₁与Fs₂，并且剪力Fs₁的方向与坐标系的y向同向，Fs₂的方向与坐标系的z向同向。对图4所示的有限元模型赋予管梁的材料属性，本实施例中，管梁的材料属性：弹性模量为210GPa，泊松比为0.3，密度为7800kg/m³。通过有限元分析软件对所述有静力分析限元模型进行静力分析。

本实施例中，静力分析时施加的各载荷均为单位载荷，采用的有限元分析软件为nastran，静力分析模块为101求解器。

步骤四：分别提取交接面A₃上的正应力与切应力。

采用数值离散方法，获得交接面A₃上在轴力F、弯矩M、扭矩T以及剪力Fs作用下的离散应力分布。

从步骤3得到的静力分析结果中，分别提取位于交接面A₃上静力分析有限元模型的各节点所对应的沿坐标系的x、y、z三轴的正应力与沿该坐标系x、y、z三轴的切应力。

本实施例中，读取nastran结果文件.f06文件中各节点所对应的沿坐标系的x、y、z三轴正应力与沿该坐标系x、y、z三轴切应力。

步骤五：确定交接面A₁上梁单元的节点自由度与交接面A₁上壳单元的节点自由度之间的约束方程。

建立梁结构的混合维有限元模型。在建立的混合维有限元模型上划分壳单元4，并以该壳单元作为管梁的局部结构；在建立的混合维有限元模型上划分梁单元7，并以该梁单元作为管梁的周围结构。如图5所示。

由于有限元分析方法中假设梁结构混合维有限元模型的梁单元7的节点具有六个方向自由度，即：沿坐标系三个轴的平动以及绕坐标系三个轴的转动。为了能够将采用梁结构混合维有限元模型的梁单元7简化建立的周围结构完整的约束到精确建模的局部结构上，分别建立梁结构混合维有限元模型的梁单元7的节点六个自由度的约束条件。

确定所述交接面A₁梁单元的节点自由度与所述交接面A₁壳单元的节点自由度之间的约束方程。

所述约束方程的确定步骤如下：

建立整体结构的混合维模型，完成交接面两侧应力做功大小的推导。由于有限元分析方法中，单元形状函数是在单元局部坐标系下表示，为了能够方便引用有限元相关知识，交接面上的应力做功大小的表示在单元局部坐标系下完成；所述的单元形状函数是一种数学函数，规定了从节点自由度到单元内所有点处的自由度值的的计算方法。

对于离散的有限元模型，在局部坐标系下，交接面上应力做功大小表示为：

其中，W_faceA表示为交接面上应力做功大小，faceA表示交接面；N表示交接面上有限元模型单元个数，具体见图2，A_i为不同维单元i交接面上的面积分变量。u'_i(x'_A,y'_A,z'_A)、σ_i'(x'_Ay'_A,z'_A)中：u_i'和σ_i'分别表示在单元局部坐标系下交接面上对应的位移及相应方向的应力分布函数，x’，y’，z’分别表示局部坐标系的x，y，z轴，下标A表示交接面。

以单元形状函数[N_d]取代公式(1)中的u'_i，即能够将单元局部坐标系下的位移函数用单元节点位移表示：

其中，表示在局部坐标系下界面积分子域i所包含离散结点的位移所构成的位移列向量，下标e无具体含义；[N_d]表示界面积分子域i对应的位移插值形状函数矩阵，下标d无具体含义，有限元中单元的形状函数常在单元局部坐标系下定义。

为了约束方程的添加，需要进一步利用界面局部坐标系下的位移向量到结构总体坐标系的位移向量之间的转换关系矩阵[T_d]，下标d无具体含义。将式(2)中局部坐标系下的位移向量向总体坐标系进行转换：

所述界面积分子域的插值形状函数矩阵[N_d]仅是为了建立已知离散节点参量与任意位置参量取值之间的插值函数关系。[N_d]的具体构造形式可随界面积分子域的形状而变化。如果界面积分子域是由4个离散结点所构成的四边形区域，则可借用四边形平面应力单元构造中类似的形状函数矩阵来类似构造插值函数矩阵。若界面积分子域选择用3个离散结点所围成的三角形区域，则可借鉴平面应力三角元的形状函数矩阵来类似构造插值函数矩阵。如果交接界面退化成了壳元所表征的形式，则界面积分也就相应退化为线积分。因此，原理上，交接界面插值函数矩阵的构造，可根据界面积分子域的选择不同进行构造。不过从工程使用的一般性要求出发，基于三角形区域构造面积分插值应该是最通用的。

对于任意截面，上述积分中，随界面积分子域的划定，[N_d]，[T_d]都能够显示表达，但任意截面的在特定形式力作用下的应力分布往往难以显示表达。为了完成式(3)的积分，将已知条件：结构交接面上相关节点的应力值代入，利用高斯积分将(3)式的面积分转换为高斯点的应力值叠加：

式中为积分子域i对应的应力做功；j表示积分子域i上选定的高斯积分点的个数，根据实际积分需要的精确度进行合理的选取；本实施例中j选为4；[N_d]_j,为(积x_j分子域i上高斯积分点j对应的单元形状函数矩阵取值；σ'_i(x'_j,y'_j,z'_j)为积分子域i上高斯积分点j对应的具体的应力大小。积分子域内的高斯点的应力值可通过已知节点的应力值插值得到：

σ′_i(x′_j，y′_j，z′_j)＝([N_g]×{σ′_n}_j) (5)

式中插值函数[N_g]的构造与划分结构的单元形状函数无关，根据实际选取的插值节点的个数确定。{σ'_n}_j表示积分子域i上高斯积分点j所包含离散结点的应力所构成的应力列向量，坐标系为界面局部坐标系。所述相关应力值对应为单元局部坐标系方向上的应力值。

分别将交接面两侧的单元做功显示表示出来，利用交接面A₁两侧做功大小相同原理，即得到梁单元离散侧交接面节点自由度与壳单元离散侧交接面节点自由度之间的约束方程。

式中，F表示作用在结构上的广义的力，U表示广义力F对应方向上的广义位移。

步骤六：确定交接面的约束方程系数。

通过界面约束方程的通用数值构造方法确定交接面A₁的约束方程系数。

对于沿梁结构混合维有限元模型的梁单元7轴向平动自由度的约束条件，将轴力以及对应的节点应力带入公式(5)和(6)中，公式(6)中有限元模型单元的形状函数[N_d]_j以及单元局部坐标系与结构总体坐标系的坐标转换矩阵都能够根据梁结构混合维有限元模型通过有限元分析方法获得，公式(6)中未知量只有U与其中U即为交接面靠近自由端侧梁单元节点位移，即为交接面靠近固支端侧模型壳单元节点位移，即建立了交接面靠近自由端侧梁单元轴向自由度与交接面靠近固支端侧壳单元节点自由度的约束关系。同理，分别计算得到梁单元各方向自由度的约束关系，并确定交接面处各点的混合维模型界面约束方程系数。

为了验证本发明所提方法的精确性，将得到的管梁的多点约束方程系数按照节点编号添加进混合维有限元模型中，得到能够兼顾整体传力特性及局部结构应力计算的精确混合维有限元计算模型，具体见图5。这里进行了两种模型计算所得结果的对比分析。对整体悬臂梁结构的细化有限元模型，如图6所示，进行动力学特性分析，并将计算得到的结果作为对比的标准结果。表1为模态频率的结果对比。

表1两种模型模态分析频率结果对比

从表1的频率对比看出，两种方法计算得到模态频率非常接近，误差极小。为了能够方便的从结构响应中读取结构更全面的动力学特性，在模态分析的基础上，对两组模型分别进行频率响应分析。外载为y方向作用力，力的幅值大小为10KN，计算频率范围取为5-500Hz，两组模型的模态阻尼系数同取为0.03；两组模型在相同节点位置频率响应曲线图6。节点具体位置为见图4与图5箭头所指位置，即节点3301。

图7所示两种模型的频率响应曲线说明，在外激励及模态阻尼相同的情况下，两种模型频响分析得到的频响曲线变化一致，混合维模型能够等效实体模型的频率响应分析结果。

故本实施例能够为振动疲劳分析中局部结构的精细应力分析提供一个高精度低阶动力学模型，以减小计算规模，提高计算效率。

Claims (4)

1.一种确定混合维模型界面约束方程系数的方法，其特征在于，具体步骤如下：

步骤一：划分管梁的局部结构与周围结构；

所述的管梁为悬臂梁；以该管梁的固支端部分为局部结构；该管梁的其余部分为周围结构；

对所述管梁进行混合维建模；建模时，所述局部结构采用高阶精确单元建立有限元模型，所述周围结构采用低阶减缩单元模拟刚度效应；

将所述有限元模型中发生不同维单元连接的截面定义为交接面A₁；根据有限元模型中的交接面确定管梁的交接面，并使管梁的交接面A₂与有限元模型中的交接面A₁位置相同；

步骤二：以交接面为中心，向两侧分别截取管梁的几何模型；

以交接面A₂为中心截取该交接面两侧的管梁作为几何模型，建立该几何模型的静力分析有限元模型；所述静力分析有限元模型用壳单元建立；

对所述静力分析有限元模型建立边界约束条件；所述静力分析有限元模型的边界约束条件与管梁的边界约束条件相同，将该静力分析有限元模型中的一端视为固支端，静力分析有限元模型中的交接面A₃的位置与管梁的交接面A₂的位置相同；以所述交接面A₃为中心，靠近固支端一侧为局部结构，靠近自由端一侧为周围结构；

步骤三：对静力分析有限元模型进行静力分析；

在所述静力分析有限元模型自由端处施加不同方向的任意值的载荷，使该静力分析有限元模型的自由端处分别产生轴力F、弯矩M、扭矩T以及剪力Fs；对所述有限元模型赋予管梁的材料属性；所述管梁的材料属性包括弹性模量、泊松比和密度；通过有限元分析软件对所述有限元模型进行静力分析；

步骤四：分别提取交接面A₃上的正应力与切应力；

采用数值离散方法，获得交接面A₃上在轴力F、弯矩M、扭矩T以及剪力Fs作用下的离散应力分布；

从步骤三得到的静力分析结果中，分别提取位于交接面A₃上静力分析有限元模型的各节点所对应的沿坐标系的x、y、z三轴的正应力与沿该坐标系x、y、z三轴的切应力；

步骤五：确定交接面A₁上梁单元的节点自由度与交接面A₁上壳单元的节点自由度之间的约束方程；

建立梁结构的混合维有限元模型；在建立的混合维有限元模型上划分壳单元，并以该壳单元作为管梁的局部结构；在建立的混合维有限元模型上划分梁单元，并以该梁单元作为管梁的周围结构；

分别建立梁结构混合维有限元模型的梁单元节点的六个自由度的约束条件；

分别将交接面两侧的单元做功显示表示出来，利用交接面A₁两侧做功大小相同原理，通过公式(6)确定所述交接面A₁梁单元的节点自由度与所述交接面A₁壳单元的节点自由度之间的约束方程：

式中，F表示作用在结构上的广义的力，U为交接面靠近自由端侧梁单元节点位移；[N_d]_j为积分子域i上高斯积分点j对应的单元形状函数矩阵，j表示积分子域i上选定的高斯积分点的个数，根据实际积分需要的精确度进行合理的选取，下标d无具体含义.σ'表示应力分布函数，x’，y’，z’分别表示局部坐标系的x，y，z轴，[T_d]为界面局部坐标系下的位移向量到结构总体坐标系的位移向量之间的转换关系矩阵，为交接面靠近固支端侧模型壳单元节点位移；

步骤六：确定交接面的约束方程系数；

通过界面约束方程的通用数值构造方法确定交接面A₁的约束方程系数；

积分子域内的高斯点的应力值可通过已知节点的应力值插值得到：

σ'_i(x'_j,y'_j,z'_j)＝([N_g]×{σ'_n}_j) (5)

式中插值函数[N_g]的构造与划分结构的单元形状函数无关，根据实际选取的插值节点的个数确定；{σ'_n}_j表示积分子域i上高斯积分点j所包含离散结点的应力所构成的应力列向量，坐标系为界面局部坐标系；相关应力值对应为单元局部坐标系方向上的应力值；对于沿梁结构混合维有限元模型的梁单元轴向平动自由度的约束条件，将轴力以及对应的节点应力带入公式(5)和(6)中，公式(6)中积分子域i上高斯积分点j对应的单元形状函数矩阵[N_d]_j以及单元局部坐标系与结构总体坐标系的坐标转换矩阵都能够根据梁结构混合维有限元模型通过有限元分析方法获得，公式(6)中未知量只有U与其中U即为交接面靠近自由端侧梁单元节点位移，即为交接面靠近固支端侧模型壳单元节点位移，即建立了交接面靠近自由端侧梁单元轴向自由度与交接面靠近固支端侧壳单元节点自由度的约束关系；同理，分别计算得到梁单元各方向自由度的约束关系，并确定交接面处各点的混合维模型界面约束方程系数。

2.如权利要求1所述确定混合维模型界面约束方程系数的方法，其特征在于，所述截取的几何模型的长度需使边界条件不会影响交接面A₂上的应力分布，截取的几何模型的长度大于该管梁横截面的最大外形尺寸的3倍。

3.如权利要求1所述确定混合维模型界面约束方程系数的方法，其特征在于，所述静力分析有限元模型自由端处施加的弯矩M有两个元素M₁与M₂，并且弯矩M₁的方向与坐标系的y向同向；M₂的方向与坐标系的z向同向；所述的剪力Fs有两个元素Fs₁与Fs₂，并且剪力Fs₁的方向与坐标系的y向同向，Fs₂的方向与坐标系的z向同向。

4.如权利要求1所述确定混合维模型界面约束方程系数的方法，其特征在于，步骤五中确定交接面A₁上梁单元的节点自由度与交接面A₁上壳单元的节点自由度之间的约束方程的具体过程是：

建立整体结构的混合维模型，完成交接面两侧应力做功大小的推导；

对于离散的有限元模型，在局部坐标系下，交接面上应力做功大小表示为：

其中，W_faceA表示为交接面上应力做功大小，faceA表示交接面；N表示交接面上有限元模型单元个数，A_i为不同维单元i交接面上的面积分变量；u'_i(x'_A,y'_A,z'_A)、σ_i'(x'_A,y'_A,z'_A)中：u′_i和σ′_i分别表示在单元局部坐标系下交接面上对应的位移及相应方向的应力分布函数，x’，y’，z’分别表示局部坐标系的x，y，z轴，下标A表示交接面；

以界面积分子域i对应的位移插值形状函数矩阵[N_d]取代公式(1)中的u'_i，即能够将单元局部坐标系下的位移函数用单元节点位移表示：

其中，表示在局部坐标系下界面积分子域i所包含离散结点的位移所构成的位移列向量，下标e无具体含义；[N_d]表示界面积分子域i对应的位移插值形状函数矩阵，下标d无具体含义，有限元中单元的形状函数常在单元局部坐标系下定义；为了约束方程的添加，需要进一步利用界面局部坐标系下的位移向量到结构总体坐标系的位移向量之间的转换关系矩阵[T_d]，下标d无具体含义；将式(2)中局部坐标系下的位移向量向总体坐标系进行转换：

为了完成式(3)的积分，将已知条件结构交接面上相关节点的应力值代入，利用高斯积分将(3)式的面积分转换为高斯点的应力值叠加：

式中为积分子域i对应的应力做功；j表示积分子域i上选定的高斯积分点的个数，根据实际积分需要的精确度选取4；σ'_i(x'_j,y'_j,z'_j)为积分子域i上高斯积分点j对应的具体的应力大小；积分子域内的高斯点的应力值可通过已知节点的应力值插值得到：

σ'_i(x'_j,y'_j,z'_j)＝([N_g]×{σ'_n}_j) (5)

式中插值函数[N_g]的构造与划分结构的单元形状函数无关，根据实际选取的插值节点的个数确定；{σ'_n}_j表示积分子域i上高斯积分点j所包含离散结点的应力所构成的应力列向量，坐标系为界面局部坐标系；相关应力值对应为单元局部坐标系方向上的应力值；

分别将交接面两侧的单元做功显示表示出来，利用交接面A₁两侧做功大小相同原理，即得到梁单元离散侧交接面节点自由度与壳单元离散侧交接面节点自由度之间的约束方程。