5 基于 ADAMS的蜗杆传动机构多刚体动力学仿真方法.docx

5 基于 ADAMS的蜗杆传动机构多刚体动力学仿真方法.docx

《5 基于 ADAMS的蜗杆传动机构多刚体动力学仿真方法.docx》由会员分享，可在线阅读，更多相关《5 基于 ADAMS的蜗杆传动机构多刚体动力学仿真方法.docx（25页珍藏版）》请在冰豆网上搜索。

5基于ADAMS的蜗杆传动机构多刚体动力学仿真方法

5基于ADAMS的蜗杆传动机构多刚体动力学仿真方法

5.1多刚体动力学仿真

多刚体动力学仿真（MBS）是近十年发展起来的机械动力系统计算机数值仿真技术，它建立在多刚体动力学理论基础之上，应用多刚体动力学模型仿真机械系统运动过程中的运动学和动力学特性。

MBS提供了设计过程中设计方案的分析和优化，在机械设计领域获得越来越广泛的应用。

目前国内外较为流行的MBS软件主要有：

ADAMS、DADS、SIMPACK等。

5.1.1多刚体动力学的算法原理

设多刚体系统由n个刚体

组成，地球为零刚体

取定一个惯性参考基

和每个刚体的连体基

，

的原点

与质心

重合。

为了确定系统内每个刚体

相对惯性基的位形，可以用它的质心

的位置矢径

的三个分量

确定位置，连体基

的三个欧拉角

确定方位。

将这三个平动坐标和三个转动坐标写成

矢量列阵[25]

（4.1）

这种确定系统内每个刚体位形所采用的统一坐标称为笛卡尔广义坐标。

n个刚体组成的多刚体系统的位形由6n个笛卡尔广义坐标确定，可以写成如下

位置矢量列阵

（4.2）

在多刚体系统运动学中，刚体

的运动用质心

的

位置矢径列阵

和确定变换关系

的

方向余弦矩阵

来描述是方便的。

它们都可以用系统的笛卡尔广义坐标表示为

，

（4.3）

将式4.3的第一式对时间求一次和两次导数，得到刚体

质心的速度和加速度列阵

创建模型

创建几何模型

给模型添加约束和运动

给模型施加载荷

测试模型

定义测量量

对模型进行仿真

回放仿真动画

绘制仿真结果曲线

验证模型

输入测试数据

在绘制的曲线图上添加测试数据

仿真结果是否与试验结果一致？

模型的细化

添加摩擦力

改善施加的载荷函数

定义柔性体

定义控制

重新仿真分析

输入测试数据

在绘制的曲线图上添加测试数据

优化分析

进行主要影响因素的研究

完成试验设计分析

进行优化设计分析

定制用户环境

定制用户菜单

定制用户对话框

使用宏命令记录并重复复杂模型操作

图2-32用ADAMS软件进行虚拟样机设计的步骤

（4.4）

以上两式中定义了

（4.5）

则刚体

在惯性基

中的角速度矩阵与方向余弦矩阵之间有如下关系：

（4.6）

可以利用上式由方向余弦矩阵

计算角速度矢量的分量

。

在一般情况下角速度的阵列表达式可由已知的运动学方程得到

（4.7）

上式对时间求导，得到加速度阵列为

因此各刚体质心的速度、加速度和各刚体的角速度、角加速度的阵列均可用笛卡尔广义坐标表示出来。

上述算法原理一般以黑箱形式封装在类似于ADAMS这样的多刚体动力学软件之中。

5.1.2多刚体动力学建模要求

利用MBS技术进行机械系统动力学仿真的关键是建立与机械系统实际组成和运动状况一致的多刚体动力学模型。

多刚体动力学模型由刚体、刚体之间的连接、载荷、驱动源、输出几部分组成。

刚体是组成多刚体动力学模型最基本的构件，一个多刚体动力学模型至少由两个刚体构成，每个刚体由其几何特性参数（刚体的几何形状、铰接点、力的作用点、仿真结果的输出点等）、材料特性参数（材料密度、摩擦系数、弹性模量等）和质量特性参数（质心、质量和转动惯量矩阵）等来描述。

刚体之间的连接用来描述各刚体之间的相互连接关系，主要有以下几种形式[16]：

（1）铰接：

铰接是指刚体直接通过相互间的运动副进行连接，根据自由度的不同可分为滑动铰、转动铰、刚性铰、球铰等多种铰接形式。

（2）弹簧连接：

弹簧连接主要有线性弹簧、非线性弹簧、扭转弹簧等多种连接形式，它不仅可以仿真实际的弹簧，还可以用于仿真两个刚体之间的线性力、非线性力、扭转力矩等。

（3）阻尼连接：

通过阻尼连接，刚体在运动过程中的相互作用力与运动速度成正比。

（4）其它连接包括绳索连接、二力杆连接等。

载荷用来描述机构在运动过程中受到的外力，它既可以是时间的函数，也可以是机构中其它变量的函数。

驱动源用来描述驱动机械系统运动的原动力，它既可以是力和扭矩，也可以是位移、速度、加速度等，既可以是时间的函数，也可以是机构中其它变量的函数。

5.2动力学仿真分析软件ADAMS简介

ADAMS（AutomaticDynamicAnalysisofMechanicalSystem）软件是由美国机械动力公司（MechanicalDynamicsInc.）开发的多体机械系统动力学仿真分析软件，以刚性体为主要分析对象。

它使用交互式图形环境和零件库，约束库，力库，创建完全参数化的机械系统几何模型，其求解器采用多刚体系统动力学理论中的拉格朗日方程方法，建立系统动力学方程，对虚拟机械系统进行静力学，运动学和动力学分析，输出位移，速度，加速度和反作用力曲线。

ADAMS软件的仿真可用于预测机械系统的性能，运动范围，碰撞检测，峰值载荷以及计算有限元的输入载荷等。

ADAMS软件包括核心模块ADAMS/View和ADAMS/Solver，以及其他

扩展模块。

ADAMS/View（界面模块）是以用户为中心的交互式图形环境，它提供丰富的零件几何图形库、约束库和力库，将便捷的图形操作，菜单操作，鼠标点取操作与交互式图形建模、仿真计算、动画显示、优化设计、X-Y曲线图处理、结果分析和数据打印等功能集成在一起。

ADAMS/Solver（求解器）是ADAMS软件的仿真运算核心，它自动形成机械系统模型的动力学方程，提供静力学、运动学和动力学的解算结果。

ADAMS/View有各种建模和求解选项，以便精确有效地解决各种工程问题。

ADAMS/Controls（控制模块）可以通过简单的继电器、逻辑与非门、阻尼线圈等建立简单的控制机构，或者利用在通用控制系统软件（如：

MATLAB、MATRIX、EASY5）中建立的控制系统框图，建立包括控制系统、液压系统、气动系统和运动机械系统的仿真模型。

ADAMS/Linear（系统模态分析模块）可以在进行系统仿真时将系统非线性的运动学或动力学方程进行线性化处理，以便快速计算系统的固有频率（特征值）、特征向量和状态空间矩阵，更快更全面地了解系统的固有特性。

ADAMS/Flex（柔性分析模块）提供ADAMS软件与有限元分析软件之间的双向数据交换接口。

利用它与ANSYS、MSC/NASTRAN、ABAQUS、I-DEAS等软件的接口，可以方便地考虑零部件的弹性特征，建立多体动力学模型，以提高系统的仿真精度。

5.3ADAMS建模、仿真的步骤

用ADAMS进行建模、仿真和分析，一般要遵循以下步骤，各步骤简述如下：

（1）建造模型

a.创建零件。

有两种途径：

一是通过ADAMS/View的零件库来创建各种简单的单元零件；二是用ADAMS/Exchange引入复杂的其他CAD软件已建好的单元零件。

由于ADAMS/View的建模能力相当有限，所以大多数情况下选择第二种方法。

b.给模型施加约束和运动。

c.给模型施加各种作用力。

（2）测试模型

定义测量（definemeasure），对模型进行初步仿真，通过仿真结果检验模型中各个零件、约束及作用力是否正确。

（3）校验模型

导入实际实验测试数据，与虚拟仿真的结果进行比较。

（4）模型的细化（Refine）

经过初步仿真确定了模型的基本运动后，可以在模型中加入更复杂的单元，如在运动副上加入摩擦，用线性方程或一般方程定义控制系统，如加入柔性连接件等等，使模型与实际系统更加近似。

（5）定制用户自己的环境

用户可以定制菜单、对话框，或利用宏命令使许多重复工作可以自动进行。

[17]

5.4阿基米德蜗杆传动动力学模型的建立

5.4.1模型简化

对于蜗杆传动的几何模型，为了便于在输入ADAMS后进行动力学仿真分析，有必要对该模型进行简化。

简化主要体现在以下几个方面[18]：

（1）对系统动力学仿真分析影响不大的零件进行简化。

如将蜗杆减速器的外壳、联接螺栓等固定构件等简化为大地模型。

（2）输入、输出轴与虚支承通过销轴连接，二者转速相同，可将虚支承省略，仅保留输入、输出轴。

简化后的蜗杆传动模型只包括蜗杆、蜗轮、蜗杆转动轴承、输出轴、蜗轮转动轴承。

这样可使模型在输入ADAMS后仿真建模更为简明，添加约束方便，在保证仿真结果正确的前提下，减小了仿真计算的复杂性。

5.4.2模型数据的引入

为了在ADAMS中建立阿基米德蜗杆传动的刚体模型，我们采用模型引入的方式。

引入的模型数据是我们事先在UG中建立蜗杆传动的装配模型。

在ADAMS/View中输入和输出其它应用程序的CAD文件之前，还需要将数据文件转换成ADAMS/Exchange图形接口模块支持的IGES、STEP、DXF/DWG、SLP和Parasolid等几种格式图形文件中的某一种。

在输入CAD图形文件时，ADAMS/Exchage首先读入IGES、STEP、DXF/DWG、SLP和Parasolid等格式的图形文件，然后将这些图形转换为ADAMS/View中的各种几何单元，这些几何单元符合ADAMS/Solve中的GRAPHICS说明语句、多义线和壳体的要求。

可以根据需要，将输入的几何形体同任何ADAMS/View样机模型的构件联系起来。

[19]

我们把UG中将建立好的蜗杆传动装配模型输出为parasolid格式的模型文件，利用ADAMS的导入命令输入。

并且建议在UG中，添加装配体的蜗轮轴线以及蜗杆轴线与适当坐标轴平行的配合，这样可以保证导入ADAMS后的模型也与坐标轴一致，方便对其进行移动、添加约束等操作。

5.4.3设定模型零件的属性参数

阿基米德蜗杆传动是多齿啮合传动，各齿间的载荷分布很复杂。

其动力学性能除了受接触变形的影响外，还受制造误差、啮合间隙和轴受力变形等的影响。

为了便于研究，对蜗杆传动的刚体模型做如下假设：

[20]

（1）装配间隙为零，制造误差忽略不计；

（2）暂不考虑由于温度升高引起的蜗杆、蜗轮、轴的变形对动力学性能的影响；

（3）各部件均视为刚体。

对于输入的蜗杆传动虚拟样机模型，仿真分析时还需包括零部件质量、转动惯量等物理属性。

在几何模型数据转换过程中，如果数据文件包含体积信息，该零部件在导入ADAMS后可以自动计算其质量、转动惯量等属性。

反之，也可以通过手工定义零部件的材料类型，再由ADAMS根据各零部件的几何尺寸自动计算其质量及转动惯量。

当模型数据输入ADAMS后，原本在UG8.5中建立的各个零件的物理属性参数、装配约束关系等都会消失，所以在ADAMS中还要根据实际情况手工添加零件属性参数及约束关系。

根据实际要求，设定蜗轮零件属性为：

材料：

黄铜，密度＝

，弹性模量＝1.06E+11N/

，泊松比＝0.35。

其他零件属性为：

材料：

钢材，密度＝

，弹性模量＝2.07E+11N/

，泊松比＝0.29。

5.4.4添加模型零件之间的约束和载荷

根据平面二次包络环面蜗杆传动的运动规律，本系统加入的约束为：

（1）蜗杆轴承与大地之间添加固定约束（fixedJoint）；

（2）蜗杆轴承与蜗杆之间添加旋转副（RevoluteJoint）；

（3）蜗轮与输出轴之间用键联接，运动速度相同，添加固定约束（fixedJoint）可节省模型计算时间；

（4）输出轴与蜗轮转动轴承之间添加旋转副（RevoluteJoint）；

（5）蜗杆轴承与大地之间添加固定约束（fixedJoint）；

（6）蜗杆与蜗轮之间添加接触力（SolidtoSolidContact），这是一种基于力的接触模型。

图4.1在ADAMS中建立的蜗杆传动动力学模型

蜗杆的转速为1500r/min，蜗轮的最大扭矩是1845N·m。

所以给蜗杆施加一个固定转速的驱动，速度是9000deg/s。

再给蜗轮输出轴添加一个反方向的转

矩，大小是1845N·m。

ADAMS为接触作用的仿真提供了极大的方便：

通过修改接触副中接触刚度及其非线性指数、阻尼系数、最大阻尼时的击穿深度、接触面静态及动态摩擦系数，接触作用形式可以仿真至相当高的程度。

[21]

5.5阿基米德蜗杆传动动力学仿真结果分析

当我们在ADAMS/view中进行了动力学仿真后，需要分析整个仿真过程的数据变化，这时可以借助ADAMS/PostProcess模块。

ADAMS/PostProcess是ADAMS的后处理模块，绘制曲线的功能十分强大，利用ADAMS/PostProcess我们可以十分清晰的观察方针结果，也可以将所得到的结果转化成表格的形式。

5.5.1输入输出轴转速分析

定义蜗杆转速n输入＝30r/min＝180deg/s，z1＝1，z2＝40，按照此传动关系

需进行系统仿真验证：

蜗轮输出轴转速n输出＝0.75r/min＝4.5deg/s。

蜗杆传动仿真后，我们进入ADAMS/PostProcess模块，绘制蜗杆输入角速度和蜗轮输出角速度的特征曲线。

图4.2蜗杆输入角速度和蜗轮输出角速度的特征曲线

图中实线为蜗杆转速，虚线为蜗轮转速。

可以看出，蜗杆传动启动时从静止到转动的瞬间，动力学变化较大。

因而蜗轮在0～0.1s之间的角速度有一个明显得加速阶段，而在0.1s之后转速达到额定值，基本保持转速在4.5deg/s。

所以蜗杆仿真模型中蜗杆输入轴、蜗轮输出轴之间的转速及传动比关系符合设计要求及理论分析结果。

5.5.2输入输出轴角加速度分析

在ADAMS/PostProcess模块中绘制蜗杆输入角加速度和蜗轮输出角加速度的特征曲线。

图4.3蜗杆输入角加速度和蜗轮输出角加速度的特征曲线

图中实线为蜗杆输入角加速度，虚线为蜗轮输出角加速度。

从图可以看出，蜗杆角速度基本保持为零，而蜗轮在传动启动时从静止到转动的瞬间，动力学变化较大。

蜗轮在0～0.1s之间的角速度从零加速到4.5deg/s，在此时间段蜗轮有一个明显的角加速过程。

并且角加速度在0时刻为最大，达到28537deg/s2。

在0～0.1s加速度呈线性减小，在0.1s之后加速度基本保持0，说明在0.1s之后蜗轮的转速基本保持不变。

5.5.3蜗杆受力分析

在ADAMS/PostProcess模块中分析蜗杆受力的特征曲线。

我们分别绘制蜗杆的x轴（与蜗杆轴平行）的受力曲线。

图4.4蜗杆的x轴的受力曲线

蜗杆的x轴与蜗轮相切，为驱动蜗轮转动提供主要动力。

在启动推动蜗轮由静止到转动的瞬间，特别在0s时刻，蜗杆x轴方向的受力达到最大，经测量为415.0272N。

在0~0.1s时刻段，蜗杆x方向受力呈线性下降，在0.1s达到最小，此后基本保持在400N。

这与刚才研究的蜗轮的速度以及加速度的变化情况是相一致的。

我们再来绘制蜗杆的y向（垂至于蜗杆、蜗轮轴线）受力图。

图4.5蜗杆的y轴的受力曲线

我们可以看到在y方向，蜗杆的受力很小。

和上面的情况相似，0s时刻达到最大值3.7452N。

在0~0.1s时刻段，蜗杆y方向受力呈线性下降，在0.1s达到最小，此后基本保持为零的受力状态。

只是期间受力略由波动，但波动值也非常小，估计应该是由于啮合冲击造成的。

最后再来绘制蜗杆的z向（平行于蜗轮轴线）受力图。

图4.6蜗杆的z轴的受力曲线

我们可以看到在z方向，蜗杆的受力较小。

和上面的情况相似，0s时刻达到最大值39.895N。

在0~0.1s时刻段，蜗杆z方向受力呈线性下降，在0.1s达到最小，此后基本在35.714N受力状态。

5蜗杆传动机构柔性多体动力学仿真技术

5.1概述

在前面的章节中，我们已经对阿基米德蜗杆传动进行了多刚体动力学分析，得到了转速、角加速度等动力学参数。

为了更精确的分析在传动过程中，蜗杆蜗轮的轮齿具体受力，仅仅使用多刚体动力学分析是不够的。

因为在进行多刚体动力学建模时，把构件作为刚体，做运动分析时不考虑其弹性变形。

实际上，机构在较大的载荷作用下或者作加、减速运动时，机构就会产生较大的弹性变形。

因此要更真实地模拟机构的动态性能，必须将构件作为柔性体进行分析。

[22]

在机械系统中，柔性体将会对整个系统的运动产主重要影响，由此产生了基于刚体和柔体的多体动力学。

其主要目的是通过考虑柔性体的变形以及柔性体之问的高度耦合作，更好地确定各部件之间的载荷传递关系。

同样整个系统的运动情况也反过来决定了每个构件的受力状况和运动状态，从而决定了构件内部的应力应变分布。

ADAMS软件充分考虑了这个问题，在其众多的功能里，开发了进行柔性体仿真的功能。

ADAMS软件柔性化主要涉及两个问题：

柔性体振型叠加问题和柔性体的运动问题。

ADAMS软件柔性化主要使用了ADAMS/Flex模块，其原理是基于Craig-Bampton方法，又称为固定界面模态综合法，属于动态子结构法范畴，该方法由Craig和Bampton于1968年率先提出。

动态子结构法方法常用于求解多自由度系统的特征值问题，其基本思路是：

首先按照工程观点或结构的几何特征，遵循某些原则（主要是便于计算或实验）。

把完整的大型复杂结构人为抽象为若干个子结构。

对各个子结构进行模态分析，保留其主要模态信息，略去高阶模态已达到缩减自由度的目的。

根据各子结构交界面的位移协调条件，构成自由度大大缩减的总体系统方程。

求解此方程可以获得系统的固有频率和模态坐标下的主振型。

最后，进行坐标变换求得系统的位移、速度、加速度等动态响应。

[23]

本章采用多体动力学软件ADAMS和有限元分析软件ANSYS相结合，将有限元方法（FEM）和多体系统仿真分析（MSS）结合在一起建立阿基米德蜗杆传动的柔性多体动力学模型，全面、准确地获得蜗杆工作时的动态载荷，作后利用有限元法对蜗杆进行动态应力、应变分析。

5.2柔体动力学简介

多柔体系统动力学可以看作是多刚体系统动力学的自然延伸。

根据多柔体系统组成特点，一般以多刚体系统动力学的研究为基础，对系统中的柔性体进行不同的处理，在机械系统中常用的处理方法有离散法，模态分析法，形函数法和有限元法等。

将柔性体的分析结果与多刚体系统的研究方法相结合，最终得到系统的动力学方程。

[24]

柔性体的广义坐标描述采用相对描述法，见图。

柔性体上节点P的位置为3个矢量和：

式中

:

整体坐标系到柔性体局部坐标系位置矢量，

为节点

在局部坐标下未变形的位置矢量，

为点

的线性变形矢量。

上面的方程通常用整体坐标系下的矩阵形式表达：

其中

是局部坐标B到整体坐标的转换矩阵。

变形矢量

通常用模态叠加来表示：

其中

为节点

的线性自由度相对应的模态矩阵。

模态坐标

是柔性体的广义坐标，

是模态数。

因此柔性体的广义坐标可以表示为

将柔性体的广义坐标带入拉格朗日方程，得到柔性体的运动方程为：

式中，

、

、

为柔性体的广义坐标及其对时间的导数，

、

为柔性体的质量矩阵及其对时间的导数，

为质量矩阵对柔性体广义坐标的偏导数，它是一个

维张量，M为模态数。

5.3有限元分析原理概述

有限元分析（FEA，FiniteElementAnalysis）利用数学近似的方法对真实物理系统（几何和载荷工况）进行模拟。

还利用简单而又相互作用的元素，即单元，就可以用有限数量的未知量去逼近无限未知量的真实系统。

有限元分析是用较简单的问题代替复杂问题后再求解。

它将求解域看成是由许多称为有限元的小的互连子域组成，对每一单元假定一个合适的（较简单的）近似解，然后推导求解这个域总的满足条件（如结构的平衡条件），从而得到问题的解。

这个解不是准确解，而是近似解，因为实际问题被较简单的问题所代替。

由于大多数实际问题难以得到准确解，而有限元不仅计算精度高，而且能适应各种复杂形状，因而成为行之有效的工程分析手段。

有限元是那些集合在一起能够表示实际连续域的离散单元。

有限元的概念早在几个世纪前就已产生并得到了应用，例如用多边形（有限个直线单元）逼近圆来求得圆的周长，但作为一种方法而被提出，则是最近的事。

有限元法最初被称为矩阵近似方法，应用于航空器的结构强度计算，并由于其方便性、实用性和有效性而引起从事力学研究的科学家的浓厚兴趣。

经过短短数十年的努力，随着计算机技术的快速发展和普及，有限元方法迅速从结构工程强度分析计算扩展到几乎所有的科学技术领域，成为一种丰富多彩、应用广泛并且实用高效的数值分析方法。

5.4多体理论与有限元结合分析基本原理

5.4.1参考系

在有限元分析中，分析对象划分的网格节点，是相对与对象本身的某一惯性坐标。

为了描述方便，一般将惯性坐标连在物体的一端，并将坐标的一个轴线与物体的轴线重合。

将有限元分析的柔体加入多体系统中，需将柔体做相对的位移，即需将柔体的有限元方程乘以转换矩阵，从而实现由局部坐标向整个模型的惯性坐标的转换。

5.4.2多体系统中的有限元方程

采用拉格朗日方法建立的多体系统动力学方程，当代入其各项因子后，其简化形式形如有限元方程：

有限元分析的通用结构分析方程，一般是在有限元分析对象的局部坐标系下建立的。

方程中的变量是相对于局部坐标的坐标元素。

为了实现有限元分析方程与多体系统动力学分析方程的统一，可以通过多体理论中的坐标转换矩阵B，实现有限元分析方程向多体系统动力学的转换。

对于形如上式的有限元方程，其惯性坐标下的刚度、阻尼、质量矩阵，是由局部坐标系下的响应矩阵乘以局部坐标向惯性坐标的转换矩阵得到的。

即上述方程可以表示为[26]：

式中：

为静止状态下的单元质量矩阵；

为静止状态下的单元阻尼矩阵；

为静止状态下的单元刚度矩阵。

而力可以表示为位移、速度、加速度的函数

当单元的位移不大时，采用惯性坐标描述方程较为方便，因为转换矩阵不需更新。

但当单元位移较大时，转换矩阵必须更新，从而局部坐标向惯性坐标系转换的相关矩阵也必须更新，从而也就明显影响了系统的求解效率。

ADAMS软件中，所有的刚体与柔体都使用一个随刚体或柔体运动的浮动局部坐标系（LocalFloatingReferenceFrame），当刚体或柔体运动时，对于系统求解的每一步，从局部坐标向惯性坐标系转换的相关矩阵也必须更新，但惯性坐标系下的矩阵不需要重新组装，从而提高了系统的求解效率。

在ADAMS软件中，假定柔体的变形应力在材料的线性变化范围内，柔体的总位移便可以通过在局部坐标系中的相邻的单元变形叠加得到。

5.4.3求解方法

在多刚体模型的基础上，预求得柔性体与多体系统作用点的力、力矩与位移、速度、加速度等边界条件，柔性体以此边界条件求出变形与力、力矩，与多刚体模型的结果进行对比，如果误差较大，可进入下一次求解；此时以柔性体的变形或力、力矩为已知条件，求得系统对此作用点相应的力、力矩或位移，直到误差达到规定的范围。

总而言之，有限元（FEA）与多体机械系统仿真（MSS）拥有相同的系统动力学求解基础。

有限元分析使用惯性坐标系下的惯性矩阵能更好地适应小位移分析的需要。

而在机械系统仿真软件ADAMS/Solver中使用的局部坐标系下的局部单元矩阵，能及时有效地描述系统内不同部件的几何关系。

因此将多体理论与有限元