CN106354954A - 一种基于叠层基函数的三维力学模态仿真模拟方法 - Google Patents

Abstract

本发明属于三维力学振动分析数值求解技术领域，涉及一种基于叠层基函数的三维力学模态仿真模拟方法。本发明首先对目标电子器件进行建模，引入位移或者应力边界条件建立对应的几何结构模型，并采用四面体网格剖分求解域，然后通过选择标量叠层基函数，将位移在所有网格内用标量叠层基函数展开，并运用标准变分原理得到目标电子器件结构的有限元方程，最后求解该方程获得特征值和特征向量，并进行后处理获得振动模态频率和振动振型，从而实现了快速构建高阶基函数，且获得高精度的数值计算结果。

Description

一种基于叠层基函数的三维力学模态仿真模拟方法

技术领域

本发明属于三维力学振动分析数值求解技术领域，涉及一种基于叠层基函数的三维力学模态仿真模拟方法。

背景技术

电子器件的使用环境往往十分恶劣，例如在崎岖上路上运输时的振动、飞机起飞、坦克高速行进、卫星和导弹上升阶段的重力加速等对电子器件的机械强度提出了十分严格的要求。机械性能却又是电子器件的可靠性和稳定性的重要组成部分，这直接影响到器件能否正常工作。因此对电子器件的机械性能进行优化设计是有必要的，而模态分析可以获得电子器件的振动特性，是其机械性能设计的重要环节，因此模态分析中高精度的获得器件的振动特性具有极其重要的意义。

目前，利用各种力学计算方法对电子器件结构模态进行仿真分析时，都是采用的有限元本征分析方法。有限元分析一般包括，单元划分、单元分析、系统综合、引入条件、求解方程组、后处理等几个步骤。在单元分析中针对单元划分得到的每个单元，假设待求函数的近似式，选择若干个结点，将近似函数表示成节点上插值函数的形式，然后利用单元上的插值函数，通过物理分析或数学分析，建立单元节点的待求函数与外界条件之间满足的方程组。在这一过程中，插值函数一般是一阶、二阶、三阶甚至更高阶数，目前各种力学计算方法都采用的是[有限元法原理简明教程，110-117页，作者：廖日东]一文中提到的画面法来形成各阶基函数，相对于高阶基函数来说，这种方法实现起来比价困难，而且不利于后期加速算法加载。当几何机构比较规则，结构比较简单采用低阶基函数来进行有限元分析，其精度是可以接受的。随着分析的结构复杂化，低阶基函数的精度已经不能满足设计者的要求，因此需要构造高阶基函数来实现有限元分析。

发明内容

针对上述存在问题或不足，为解决构造有限元模态分析中的高阶插值基函数，从而获得高精度的数值模拟结果；本发明提供了一种基于叠层基函数的三维力学模态仿真模拟方法。

该基于叠层基函数的三维力学模态仿真模拟方法，包括以下步骤：

A.将目标电子器件结构进行建模，引入位移边界条件或者应力边界条件建立对应的几何结构模型；

B.建立电子器件结构的线性弹性力学基本边值问题矩阵形式；

C.采用四面体网格剖分求解域；

D.选择标量叠层基函数，将位移在所有网格内用标量叠层基函数展开，并运用标准变分原理得到目标电子器件结构的有限元方程；

E.引入结构的惯性力，得到结构的自由振动有限元广义本征方程；

F.求E步骤所获得的本征方程，获得特征值λ和对应的特征向量即振幅向量；

G.对F步骤获得的特征值和对应特征向量进行后处理获得振动模态频率和对应振动振型。

进一步优选，所述步骤D中，基函数为标量叠层基函数，基函数根据所求问题的复杂度和计算结果的精度要求按照如下叠层规则进行确定：

$W^p=\left\{{\mathrm{\λ}}_1^i{\mathrm{\λ}}_2^j{\mathrm{\λ}}_3^k{\mathrm{\λ}}_4^l|\begin{array}{c}i=0,1,\mathrm{...},p\\ j=0,1,\mathrm{...},p-i\\ k=0,1,\mathrm{...},p-i-j\\ l=0,1,\mathrm{...},p-i-j-k\end{array}\right\}---\left(1\right)$

$Dim\left(W^p\right)=\frac{(p+1)(p+2)(p+3)}{6}---\left(2\right)$

式(1)(2)中，Wp表示所有基函数的集合，p表示选取基函数的阶数，Dim(W^p)表示基函数的个数，i,j,k,l表示上标。λ₁，λ₂，λ₃，λ₄为四面体网格最基本的标量基函数，其表达式和推导过程是一种公知，这里不再阐述。

一般基函数会附着在点、边、面、体上则基函数集合表示成

$W^p=W_n^p\⊕W_e^p\⊕W_f^p\⊕W_v^p---\left(3\right)$

其中n,e,f,v表示点、边、面、体。

附着在点上的基函数为

$\begin{array}{c}{W}_n^p=\{{\mathrm{\λ}}_1,{\mathrm{\λ}}_2,{\mathrm{\λ}}_3,{\mathrm{\λ}}_4\}\\ Dim\left(W_n^p\right)=4,(p=1)\end{array}---\left(4\right)$

附着在边上的基函数为

$\begin{array}{c}{W}_e^p=\{e_{1,2}^p,e_{1,3}^p,e_{1,4}^p,e_{2,3}^p,e_{2,4}^p,e_{3,4}^p\}\\ Dim\left(W_e^p\right)=6(p-1),(p=2,3)\end{array}---\left(5\right)$

其中

附着在面上的基函数为

$\begin{array}{c}{W}_f^p=\{f_{1,2,3}^p,f_{1,2,4}^p,f_{1,3,4}^p,f_{2,3,4}^p\}\\ Dim\left(W_f^p\right)=4\frac{(p-2)(p-1)}{2}\end{array}---\left(6\right)$

其中

附着在体上的基函数为

$\begin{array}{c}{W}_v^p=\left\{v_{1,2,3,4}^p\right\}\\ Dim\left(W_v^p\right)=\frac{(p-3)(p-2)(p-1)}{6}\end{array}---\left(7\right)$

其中

进一步优选，所述步骤G中后处理具体为：运用处理特征值λ，处理对应特征向量u为位移矢量，f₀为振动模态频率，N为插值基函数的矩阵形式。

与现有技术相比，本发明的有益效果：

本发明提出的基于叠层基函数的三维力学模态仿真模拟方法可以快速的构建高阶基函数，且获得高精度的数值计算结果。

附图说明

图1是本发明基于叠层基函数的三维力学模态仿真模拟方法的流程图；

图2是实施例四节点四面体单元示意图；

图3是基函数附着在点上示意图；

图4是基函数附着在边上示意图；

图5是基函数附着在面上示意图；

图6是基函数附着在体上示意图。

具体实施方式

下面结合附图和实施例来详细说明本发明的技术方案。

参照图1，一种基于叠层基函数的三维力学模态仿真模拟方法，包括以下步骤：

A.将目标电子器件结构进行建模，引入位移边界条件或者应力边界条件建立对应的几何结构模型。

根据电子器件的特性，引入位移边界条件建立对应的几何结构模型仿真整个结构的振动特性。

B.建立电子器件结构的线性弹性力学基本边值问题矩阵形式。

对于空间边值问题，在结构(弹性体)内部我们要考虑静力学、几何学、物理学三方面条件，分别建立三套方程；并给定约束或面力的边界上，建立位移边界条件或应力边界条件。具体如下：

平衡微分方程

$\begin{array}{c}\frac{\∂{\mathrm{\σ}}_x}{\∂x}+\frac{\∂{\mathrm{\τ}}_{yx}}{\∂y}+\frac{\∂{\mathrm{\τ}}_{zx}}{\∂z}+f_x=0\\ \frac{\∂{\mathrm{\σ}}_y}{\∂y}+\frac{\∂{\mathrm{\τ}}_{zy}}{\∂z}+\frac{\∂{\mathrm{\τ}}_{xy}}{\∂x}+f_y=0\\ \frac{\∂{\mathrm{\σ}}_z}{\∂z}+\frac{\∂{\mathrm{\τ}}_{xz}}{\∂x}+\frac{\∂{\mathrm{\τ}}_{yz}}{\∂y}+f_z=0\end{array}---\left(8\right)$

几何方程

$\begin{array}{c}{\mathrm{\ϵ}}_x=\frac{\∂u}{\∂x},{\mathrm{\ϵ}}_y=\frac{\∂v}{\∂y},{\mathrm{\ϵ}}_z=\frac{\∂w}{\∂z}\\ {\mathrm{\γ}}_{xy}=\frac{\∂v}{\∂x}+\frac{\∂u}{\∂y},{\mathrm{\γ}}_{yz}=\frac{\∂w}{\∂y}+\frac{\∂v}{\∂z},{\mathrm{\γ}}_{zx}=\frac{\∂w}{\∂x}+\frac{\∂u}{\∂z}\end{array}---\left(9\right)$

物理方程

$\begin{array}{c}{\mathrm{\σ}}_x=\frac{E}{1+\mathrm{\μ}}{\mathrm{\ϵ}}_x+\frac{\mathrm{\μ}E}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}({\mathrm{\ϵ}}_x+{\mathrm{\ϵ}}_y+{\mathrm{\ϵ}}_z)\\ {\mathrm{\σ}}_y=\frac{E}{1+\mathrm{\μ}}{\mathrm{\ϵ}}_y+\frac{\mathrm{\μ}E}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}({\mathrm{\ϵ}}_x+{\mathrm{\ϵ}}_y+{\mathrm{\ϵ}}_z)\\ {\mathrm{\σ}}_z=\frac{E}{1+\mathrm{\μ}}{\mathrm{\ϵ}}_z+\frac{\mathrm{\μ}E}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}({\mathrm{\ϵ}}_x+{\mathrm{\ϵ}}_y+{\mathrm{\ϵ}}_z)\\ {\mathrm{\τ}}_{xy}=\frac{E}{2(1+\mathrm{\μ})}{\mathrm{\γ}}_{xy}\\ {\mathrm{\τ}}_{yz}=\frac{E}{2(1+\mathrm{\μ})}{\mathrm{\γ}}_{yz}\\ {\mathrm{\τ}}_{zx}=\frac{E}{2(1+\mathrm{\μ})}{\mathrm{\γ}}_{zx}\end{array}---\left(10\right)$

位移边界条件

${\left(u\right)}_{s_u}=\stackrel{\‾}{u},{\left(v\right)}_{s_u}=\stackrel{\‾}{v},{\left(w\right)}_{s_u}=\stackrel{\‾}{w}---\left(11\right)$

应力边界条件

$\begin{array}{c}{({\mathrm{l\σ}}_x+{\mathrm{m\τ}}_{yx}+{\mathrm{n\τ}}_{zx})}_{s_{\mathrm{\σ}}}={\stackrel{\‾}{p}}_x\\ {({\mathrm{m\σ}}_y+{\mathrm{n\τ}}_{zy}+{\mathrm{l\τ}}_{xy})}_{s_{\mathrm{\σ}}}={\stackrel{\‾}{p}}_y\\ {({\mathrm{n\σ}}_z+{\mathrm{l\τ}}_{xz}+{\mathrm{m\τ}}_{yz})}_{s_{\mathrm{\σ}}}={\stackrel{\‾}{p}}_z\end{array}---\left(12\right)$

上述(8)(9)(10)(11)(12)式中，σ_x,σ_y,σ_z,τ_xy＝τ_yx,τ_yz＝τ_zy,τ_zx＝τ_xz表示求解区域中6个应力分量，ε_x,ε_y,ε_z,γ_xy,γ_yz,γ_zx表示求解区域中6个形变应力分量，u,v,w表示求解区域中3个位移分量。E是求解区域中结构的杨氏弹性模量，μ是求解区域中结构的泊松比，S_u表示位移边界面，S_σ表示应力边界面。l＝cos(n′,x),m＝cos(n′,y),n＝cos(n′,z)，表示应力边界面S_σ上的方向余弦，其中n′为应力边界面S_σ的外法线，x,y,z为应力边界面S_σ上三个方向的坐标值。为位移边界面S_u上的位移值，f_x,f_y,f_z为求解区域内结构受到的各个方向的体力，为应力边界面S_σ上受到的各个方向的面力，具体推导过程为一种公知过程，这里不再阐述。

在笛卡尔坐标系oxyz中，我们假设弹性体由空间区域Ω定义，任意点(x,y,z)的位移矢量记为

u＝[u v w]^T (13)

体力矢量记为

f＝[f_x f_y f_z]^T (14)

面力矢量记为

$\stackrel{\‾}{p}={\left[\begin{array}{ccc}{\stackrel{\‾}{p}}_x& {\stackrel{\‾}{p}}_y& {\stackrel{\‾}{p}}_z\end{array}\right]}^T---\left(15\right)$

应变张量记为

ε＝[ε_x ε_x ε_x γ_xy γ_yz γ_zx]^T (16)

应力张量记为

σ＝[σ_x σ_x σ_x τ_xy τ_yz τ_zx]^T (17)

需要强调的是，上述体力矢量，面力矢量，应变张量，应力张量均为空间坐标函数，式中为简洁起见未明确写出。

令

$A=\left[\begin{array}{cccccc}\frac{\∂}{\∂x}& 0& 0& \frac{\∂}{\∂y}& 0& \frac{\∂}{\∂z}\\ 0& \frac{\∂}{\∂y}& 0& \frac{\∂}{\∂x}& \frac{\∂}{\∂z}& 0\\ 0& 0& \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂y}& \frac{\∂}{\∂x}\end{array}\right]---\left(18\right)$

$L=\left[\begin{array}{ccc}\frac{\∂}{\∂x}& 0& 0\\ 0& \frac{\∂}{\∂y}& 0\\ 0& 0& \frac{\∂}{\∂z}\\ \frac{\∂}{\∂y}& \frac{\∂}{\∂x}& 0\\ 0& \frac{\∂}{\∂z}& \frac{\∂}{\∂y}\\ \frac{\∂}{\∂z}& 0& \frac{\∂}{\∂x}\end{array}\right]=A^T---\left(19\right)$

$D=\frac{E(1-\mathrm{\μ})}{(1+\mathrm{\μ})(1-2\mathrm{\μ})}\left[\begin{array}{cccccc}1& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 0& 0& 0\\ \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 1& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 0& 0& 0\\ \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& \frac{\mathrm{\μ}}{1-\mathrm{\μ}}& 1& 0& 0& 0\\ 0& 0& 0& \frac{(1-2\mathrm{\μ})}{2(1-\mathrm{\μ})}& 0& 0\\ 0& 0& 0& 0& \frac{(1-2\mathrm{\μ})}{2(1-\mathrm{\μ})}& 0\\ 0& 0& 0& 0& 0& \frac{(1-2\mathrm{\μ})}{2(1-\mathrm{\μ})}\end{array}\right]---\left(20\right)$

$n=\left[\begin{array}{cccccc}{n}_x& 0& 0& n_y& 0& n_z\\ 0& n_y& 0& n_x& n_z& 0\\ 0& 0& n_z& 0& n_y& n_x\end{array}\right]---\left(21\right)$

其中n_x＝l,n_y＝m,n_z＝n

则平衡微分方程(8)、几何方程(9)、物理方程(10)位移边界(11)、应力边界(12)可分别写成矩阵形式为

Aσ+f＝0 (22)

ε＝Lu (23)

σ＝Dε (24)

C.采用四面体网格剖分求解域。

采用四面体网格剖分求解域是有限元方法中的一种公知过程，因此本步骤不再详细描述。剖分后的求解域被人为分割为多个三维四面体网格，从而将连续的几何结构空间转化为离散的网格空间。

D.选择叠层基函数，将位移在所有网格内用标量叠层基函数展开，并运用标准变分原理得到电子器件结构的有限元方程。

如果B步骤中的位移矢量u是我们的目标函数，那么对于位移矢量u我们用插值基函数展开如下形式

$\begin{array}{c}u={\left[\begin{array}{ccc}\underset{i=1}{\overset{n}{\Σ}}{N}_i{u}_i& \underset{i=1}{\overset{n}{\Σ}}{N}_i{v}_i& \underset{i=1}{\overset{n}{\Σ}}{N}_i{w}_i\end{array}\right]}^T\\ =\left[\begin{array}{cccccccccc}{N}_1& 0& 0& N_2& 0& 0& \mathrm{...}& N_n& 0& 0\\ 0& N_1& 0& 0& N_2& 0& \mathrm{...}& 0& N_n& 0\\ 0& 0& N_1& 0& 0& N_2& \mathrm{...}& 0& 0& N_n\end{array}\right]\\ \·{\left[\begin{array}{cccccccccc}{u}_1& v_1& w_1& u_2& v_2& w_2& \mathrm{...}& u_n& v_n& w_n\end{array}\right]}^T\end{array}---\left(27\right)$

令

$N_i=\left[\begin{array}{ccc}{N}_i& 0& 0\\ 0& N_i& 0\\ 0& 0& N_i\end{array}\right]---\left(28\right)$

N＝[N₁ N₂ … N_n] (29)

α＝[u₁ v₁ w₁ u₂ v₂ w₂ … u_n v_n w_n]^T (30)

则式(27)可以记为

u＝Nα (31)

其中N_i为我们要选择的插值基函数。

如图2所示四面体单元中i,j,k,l代表四个顶点的编号，我们首先得到四个最基本的基函数，具体推导过程为一种公知过程，这里不再阐述：

${\mathrm{\λ}}_1=N_i=\frac{a_i+b_{i}x+c_{i}y+d_{i}z}{6V}---\left(32\right)$

${\mathrm{\λ}}_2=N_j=\frac{a_j+b_{j}x+c_{j}y+d_{j}z}{6V}---\left(33\right)$

${\mathrm{\λ}}_3=N_k=\frac{a_k+b_{k}x+c_{k}y+d_{k}z}{6V}---\left(34\right)$

${\mathrm{\λ}}_4=N_l=\frac{a_l+b_{l}x+c_{l}y+d_{l}z}{6V}---\left(35\right)$

式中

$a_i=\left|\begin{array}{ccc}{x}_j& y_j& z_j\\ x_k& y_k& z_k\\ x_l& y_l& z_l\end{array}\right|---\left(36\right)$

$b_i=\left|\begin{array}{ccc}1& y_j& z_j\\ 1& y_k& z_k\\ 1& y_l& z_l\end{array}\right|---\left(37\right)$

$c_i=\left|\begin{array}{ccc}1& x_j& z_j\\ 1& x_k& z_k\\ 1& x_l& z_l\end{array}\right|---\left(38\right)$

$d_i=\left|\begin{array}{ccc}1& x_j& y_j\\ 1& x_k& y_k\\ 1& x_l& y_l\end{array}\right|---\left(39\right)$

$V=\frac{1}{6}\left|\begin{array}{cccc}1& x_i& y_i& z_i\\ 1& x_j& y_j& z_j\\ 1& x_k& y_k& z_k\\ 1& x_l& y_l& z_l\end{array}\right|---\left(40\right)$

将(36)式、(37)式、(38)式、(39)式中的i,j,k,l轮换，得到a_j,a_k,a_l，b_j,b_k,b_l，c_j,c_k,c_l，d_j,d_k,d_l。V为四面体的体积。为了使四面体的体积不为负值，单元节点编号i,j,k,l必须依照一定的顺序。在右手坐标系中，当i,j,k的方向转动时，右手螺旋应指向l的方向。

对于标量叠层高阶基函数的选择有如下规则叠层规则

$W^p=\left\{{\mathrm{\λ}}_1^i{\mathrm{\λ}}_2^j{\mathrm{\λ}}_3^k{\mathrm{\λ}}_4^l|\begin{array}{c}i=0,1,\mathrm{...},p\\ j=0,1,\mathrm{...},p-i\\ k=0,1,\mathrm{...},p-i-j\\ l=0,1,\mathrm{...},p-i-j-k\end{array}\right\}---\left(41\right)$

$Dim\left(W^p\right)=\frac{(p+1)(p+2)(p+3)}{6}---\left(42\right)$

式(41)(42)中，Wp表示所有基函数的集合，p表示选取基函数的阶数，Dim(W^p)表示基函数的个数，i,j,k,l表示上标。

一般基函数会附着在点、边、面、体上则基函数集合表示成

$W^p=W_n^p\⊕W_e^p\⊕W_f^p\⊕W_v^p---\left(43\right)$

其中n,e,f,v表示点、边、面、体。

如图3所示附着在点上的基函数为

$\begin{array}{c}{W}_n^p=\{{\mathrm{\λ}}_1,{\mathrm{\λ}}_2,{\mathrm{\λ}}_3,{\mathrm{\λ}}_4\}\\ Dim\left(W_n^p\right)=4,(p=1)\end{array}---\left(44\right)$

如图4所示附着在边上的基函数为

$\begin{array}{c}{W}_e^p=\{e_{1,2}^p,e_{1,3}^p,e_{1,4}^p,e_{2,3}^p,e_{2,4}^p,e_{3,4}^p\}\\ Dim\left(W_e^p\right)=6(p-1),(p=2,3)\end{array}---\left(45\right)$

其中

如图5(p＝3)所示附着在面上的基函数为

$\begin{array}{c}{W}_f^p=\{f_{1,2,3}^p,f_{1,2,4}^p,f_{1,3,4}^p,f_{2,3,4}^p\}\\ Dim\left(W_f^p\right)=4\frac{(p-2)(p-1)}{2}\end{array}---\left(46\right)$

其中

如图6(p＝4)所示附着在体上的基函数为

$\begin{array}{c}{W}_v^p=\left\{v_{1,2,3,4}^p\right\}\\ Dim\left(W_v^p\right)=\frac{(p-3)(p-2)(p-1)}{6}\end{array}---\left(47\right)$

其中

定义弹性体的全部形变势能密度U₁为，具体推导为公知过程这里不再阐述

$U_1=\frac{1}{2}({\mathrm{\σ}}_x{\mathrm{\ϵ}}_x+{\mathrm{\σ}}_y{\mathrm{\ϵ}}_y+{\mathrm{\σ}}_z{\mathrm{\ϵ}}_z+{\mathrm{\τ}}_{xy}{\mathrm{\γ}}_{xy}+{\mathrm{\τ}}_{yz}{\mathrm{\γ}}_{yz}+{\mathrm{\τ}}_{zx}{\mathrm{\γ}}_{zx})---\left(48\right)$

用矩阵表示为

$U_1=\frac{1}{2}{\mathrm{\ϵ}}^T\mathrm{\σ}---\left(49\right)$

形变势能可以全部用位移分量来表示。为此，利用物理方程矩阵式(24)和几何方程矩阵式(23)代入上式得

$U_1=\frac{1}{2}{\left(Lu\right)}^T\left(DLu\right)---\left(50\right)$

则得弹性体全部形变势能为

$U=\∫\∫{\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(Lu\right)}^T\left(DLu\right)dV---\left(51\right)$

若弹性体受体力和面力作用，空间区域Ω内的体力分量为f＝[f_x f_y f_z]^T，s_σ边界上的面力分量为则外力在实际位移上所做的功称为外力功

$W=\∫\∫{\∫}_{\mathrm{\Ω}}{u}^{T}fdV+\∫{\∫}_{s_{\mathrm{\σ}}}{u}^T\stackrel{\‾}{p}dS---\left(52\right)$

由于外力做了功，消耗了外力势能，因此在发生实际位移时，弹性体的外力势能V是

$\begin{array}{c}V=-W\\ =-\∫\∫{\∫}_{\mathrm{\Ω}}{u}^{T}fdV-\∫{\∫}_{s_{\mathrm{\σ}}}{u}^T\stackrel{\‾}{p}dS\end{array}---\left(53\right)$

弹性体的形变势能与外力势能之和，即为弹性体的总势能E_P

E_P＝U+V (54)

将式(51)、式(31)和式(53)代入式(54)得

$\begin{array}{c}{E}_P=\∫\∫{\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(Lu\right)}^T\left(DLu\right)dV-\∫\∫{\∫}_{\mathrm{\Ω}}{u}^{T}fdV-\∫{\∫}_{s_{\mathrm{\σ}}}{u}^T\stackrel{\‾}{p}dS\\ =\∫\∫{\∫}_{\mathrm{\Ω}}\frac{1}{2}{\left(LN\mathrm{\α}\right)}^T\left(DLN\mathrm{\α}\right)dV-\∫\∫{\∫}_{\mathrm{\Ω}}{\left(N\mathrm{\α}\right)}^{T}fdV-\∫{\∫}_{s_{\mathrm{\σ}}}{\left(N\mathrm{\α}\right)}^T\stackrel{\‾}{p}dS\\ =\∫\∫{\∫}_{\mathrm{\Ω}}\frac{1}{2}{\mathrm{\α}}^T\left(N^T{L}^{T}DLN\right)\mathrm{\α}dV-\∫\∫{\∫}_{\mathrm{\Ω}}{\mathrm{\α}}^T{N}^{T}fdV-\∫{\∫}_{s_{\mathrm{\σ}}}{\mathrm{\α}}^T{N}^T\stackrel{\‾}{p}dS\\ =\frac{1}{2}{\mathrm{\α}}^T\[\∫\∫{\∫}_{\mathrm{\Ω}}\left(N^T{L}^{T}DLN\right)dV\]\mathrm{\α}-{\mathrm{\α}}^T\∫\∫{\∫}_{\mathrm{\Ω}}{N}^{T}fdV-{\mathrm{\α}}^T\∫{\∫}_{s_{\mathrm{\σ}}}{N}^T\stackrel{\‾}{p}dS\end{array}---\left(55\right)$

根据最小势能原理，即变分原理，弹性力学微分方程定解问题等价于要求总势能E_P的最小值，则有即

$(\∫\∫{\∫}_{\mathrm{\Ω}}{N}^T{L}^{T}DLNdV)\mathrm{\α}-\∫\∫{\∫}_{\mathrm{\Ω}}{N}^{T}fdV-\∫{\∫}_{s_{\mathrm{\σ}}}{N}^T\stackrel{\‾}{p}dS=0---\left(56\right)$

也即

$(\∫\∫{\∫}_{\mathrm{\Ω}}{N}^T{L}^{T}DLNdV)\mathrm{\α}=\∫\∫{\∫}_{\mathrm{\Ω}}{N}^{T}fdV+\∫{\∫}_{s_{\mathrm{\σ}}}{N}^T\stackrel{\‾}{p}dS---\left(57\right)$

令

K＝∫∫∫_ΩN^TL^TDLNdV (58)

$F=\∫\∫{\∫}_{\mathrm{\Ω}}{N}^{T}fdV+\∫{\∫}_{s_{\mathrm{\σ}}}{N}^T\stackrel{\‾}{p}dS---\left(59\right)$

则方程也可简记为

Kα＝F (60)

其中K为弹性体的刚度矩阵，α为位移向量，F为外载荷量。

E.引入结构的惯性力，得到结构的自由振动有限元广义本征方程。

当研究结构振动问题时，上述E步骤的α位移向量为时间的函数，我们重新定义

α(t)＝[u₁(t) v₁(t) w₁(t) u₂(t) v₂(t) w₂(t) … u_n(t) v_n(t) w_n(t)]^T (61)

则根据E步骤讨论得到的有限元方程(60)，引入物体的惯性力得到

$M\stackrel{\·\·}{\mathrm{\α}}\left(t\right)+K\mathrm{\α}\left(t\right)=F\left(t\right)---\left(62\right)$

其中

M＝∫∫∫_ΩρN^TNdΩ (63)

M为质量矩阵，ρ为求解区域物体的密度，为对时间的二阶导数。

当物体自由振动时，此时F(t)＝0方程(62)退化为

$M\stackrel{\·\·}{\mathrm{\α}}\left(t\right)+K\mathrm{\α}\left(t\right)=0---\left(64\right)$

其振动形式叫做自由振动，该方程有解的形式

$\mathrm{\α}\left(t\right)=\widehat{\mathrm{\α}}{e}^{j\mathrm{\ω}t}---\left(65\right)$

这是简谐振动的形式，其中ω为常数，为振幅向量，将其代入式(64)中，有

$(-{\mathrm{\ω}}^{2}M\widehat{\mathrm{\α}}+K\widehat{\mathrm{\α}})e^{j\mathrm{\ω}t}=0---\left(66\right)$

消去e^jωt后，有

$K\widehat{\mathrm{\α}}={\mathrm{\ω}}^{2}M\widehat{\mathrm{\α}}---\left(67\right)$

该方程有非零解的条件是

|K-ω²M|＝0 (68)

这就是本征方程。

F.求E步骤所获得的本征方程，获得特征值λ和对应的特征向量即振幅向量。

求解E步骤得到的广义本征方程(68)，得到一系列的特征值λ和对应的特征向量

G.对F步骤获得的特征值和对应特征向量进行后处理获得振动模态频率和对应振动振型。

对E步骤获得特征值λ进行处理，对应的振动模态频率为

$f_0=\frac{\sqrt{\mathrm{\λ}}}{2\mathrm{\π}}---\left(69\right)$

根据得到的本征方程(68)的特征向量结合基函数，由公式(31)可以得到求解域内的位移分布，这就是对应振动模态频率的振动振型。

Claims (3)

1.一种基于叠层基函数的三维力学模态仿真模拟方法，包括以下步骤：

A.将目标电子器件结构进行建模，引入位移边界条件或者应力边界条件建立对应的几何结构模型；

B.建立电子器件结构的线性弹性力学基本边值问题矩阵形式；

C.采用四面体网格剖分求解域；

D.选择标量叠层基函数，将位移在所有网格内用标量叠层基函数展开，并运用标准变分原理得到目标电子器件结构的有限元方程；

E.引入结构的惯性力，得到结构的自由振动有限元广义本征方程；

F.求E步骤所获得的本征方程，获得特征值λ和对应的特征向量即振幅向量；

G.对F步骤获得的特征值和对应特征向量进行后处理获得振动模态频率和对应振动振型。

2.如权利要求1所述基于叠层基函数的三维力学模态仿真模拟方法，其特征在于：所述步骤G中后处理具体为：运用处理特征值λ，处理对应特征向量u为位移矢量，f₀为振动模态频率，N为插值基函数的矩阵形式。

3.如权利要求1所述基于叠层基函数的三维力学模态仿真模拟方法，其特征在于：所述步骤D中，基函数为标量叠层基函数，基函数根据所求问题的复杂度和计算结果的精度要求按照如下叠层规则进行确定：

$W^p=\left\{{\mathrm{\λ}}_1^j{\mathrm{\λ}}_2^j{\mathrm{\λ}}_3^k{\mathrm{\λ}}_4^l\left|\begin{array}{c}i=0,1,...,p\\ j=0,1,...,p-i\\ k=0,1,...,p-i-j\\ l=0,1,...,p-i-j-k\end{array}\right.\right\}---\left(1\right)$

$Dim\left(W^p\right)=\frac{(p+1)(p+2)(p+3)}{6}---\left(2\right)$

式(1)(2)中，Wp表示所有基函数的集合，p表示选取基函数的阶数，Dim(W^p)表示基函数的个数，i,j,k,l表示上标，λ₁，λ₂，λ₃，λ₄为四面体网格最基本的标量基函数；

一般基函数会附着在点、边、面、体上则基函数集合表示成

$W^p=W_n^p\⊕W_e^p\⊕W_f^p\⊕W_v^p---\left(3\right)$

其中n,e,f,v表示点、边、面、体；

附着在点上的基函数为：

$\begin{array}{c}{W}_n^p=\{{\mathrm{\λ}}_1,{\mathrm{\λ}}_2,{\mathrm{\λ}}_3,{\mathrm{\λ}}_4\}\\ Dim\left(W_n^p\right)=4,(p=1)\end{array}---\left(4\right)$

附着在边上的基函数为：

$\begin{array}{c}{W}_e^p=\{e_{1,2}^p,e_{1,3}^p,e_{1,4}^p,e_{2,3}^p,e_{2,4}^p,e_{3,4}^p\}\\ Dim\left(W_e^p\right)=6(p-1),(p=2,3)\end{array}---\left(5\right)$

其中

附着在面上的基函数为：

$\begin{array}{c}{W}_f^p=\{f_{1,2,3}^p,f_{1,2,4}^p,f_{1,3,4}^p,f_{2,3,4}^p\}\\ Dim\left(W_f^p\right)=4\frac{(p-2)(p-1)}{2}\end{array}---\left(6\right)$

其中

附着在体上的基函数为：

$\begin{array}{c}{W}_v^p=\left\{v_{1,2,3,4}^p\right\}\\ Dim\left(W_v^p\right)=\frac{(p-3)(p-2)(p-1)}{6}\end{array}---\left(7\right)$

其中