矿井提升系统数值仿真技术最新章节_刘义著

2.4 基于相对节点方程的钢丝绳建模方法

2.4.1 单元的划分

对摩擦式提升系统建模最关键的是建立钢丝绳的动力学方程。本节把钢丝绳看作由一系列连续三维梁单元组成的多柔性体系统，利用相对节点坐标方程，通过节点的相对节点位移建立钢丝绳的大变形结构的动力学平衡方程。同绝对节点法相似的是，相对节点坐标方程的节点位移变形也是沿着一定的路径累加得到的，并且在形成系统的运动方程之前，必须定义系统连续单元的信息。因此，同绝对节点坐标方程一样，首先把系统离散为有限单元组成的系统。本节利用欧拉-伯努力三维梁单元对钢丝绳的空间路径进行单元划分，得到一系列按照顺序分布的节点，如图2-10所示。

图2-10 单元划分及节点

显然，当考虑摩擦式提升机系统的尾绳时，划分单元后的钢丝绳路径形成一个闭环系统（即按照正向或者逆向路径形成的单元的节点首尾相连），当不考虑尾绳的影响时，形成一个开环系统。由于应用欧拉-伯努力三维梁所建立的动力学方程不考虑横截面剪切变形，扭转时认为横截面不发生翘曲。

2.4.2 相对节点方程

从空间结构中已经离散化的连续的一系列梁单元中任意取出两个单元 i -1和 i 用图2-11所示的坐标体系来说明做空间运动的柔性梁的相对节点的变形。单元 i -1由节点 i -1和 i 组成，单元 i 由节点 i 和节点 i +1组成。在空间结构中，在每个梁单元的节点上有6个自由度：沿着 X 、 Y 、 Z 3个方向的线位移和绕 X 、 Y 、 Z 3个轴的角位移。节点坐标的 X 轴方向为沿着梁单元的纵向的方向， Z 轴方向为沿着梁单元的横向方向， Y 轴方向由右手定则确定。变形前后的两个梁单元如图2-11所示。引入梁单元的节点参考坐标系，来描述节点的位置和方向矢量。图中的坐标系 XYZ 为总体惯性坐标系； x _k y _k z _k （ k 为梁单元的节点编号 i -1, i ）为节点 k 的节点参考坐标系； r _k 表示节点 k 在总体惯性坐标系的位置矢量。 x _（ _i _-1） _i y _（ _i _-1） _i z _（ _i _-1） _i 为固结在节点 i 上的关于节点 i -1的随体坐标系。随体坐标系与结构力学中有限元单元坐标系相同，是与节点上邻域的微元相固连的，而不是浮动的，即随体坐标系随着柔性体的变化而变化。在单元 i -1没有变形的情况下，节点随体参考坐标系 x _（ _i _-1） _i y _（ _i _-1） _i z _（ _i _-1） _i 和节点参考坐标系 x _i _-1 y _i _-1 z _i _-1 的方向是一致的。

图2-11 两个梁单元的变形示意图

柔性体在运动时，尤其是在发生小位移时，其空间的运动可以分解为刚性移动和刚性转动的组合运动来描述柔性体的变形。即在柔性体上任意选择一个基点，在这个基点上建立动参考系（或称为物体坐标系），这样就可以将柔性体的运动分解为随体坐标系的刚性运动和相对于动参考系的变形运动两部分来研究。然后将相对变形运动的各变量，通过空间旋转变换矩阵，转换到总体惯性参考坐标系下，从而建立柔性体系统的运动方程。

对于柔性体任意一点 p 的位置变化可以表示为

式中， r 为点 p 在总体惯性坐标系中关于 X 、 Y 、 Z 的位置矢量； r _o 为浮动坐标系原点在总体惯性坐标系中的位置矢量； A 为关于空间坐标变换的方向余弦矩阵； s _p 为柔性体在没变形时，点 p 在随体坐标系中的位置矢量； u _p 为相对变形矢量，包括位置变形及方向的变化。

基于以上的原理，对于柔性体单元 i -1，节点 i 在总体惯性坐标中的位置矢量可以利用节点 i -1的位置矢量以及与节点 i 的相对节点变形来表示。在这里，节点 i -1相当于浮动参考坐标系的原点；相应的坐标系 x _（ _i _-1） y _（ _i _-1） z _（ _i _-1）为柔性体（单元 i -1）的浮动参考坐标系；分别用 u ′ _（ _i _-1） _i 和 Θ ′ _（ _i _-1） _i 表示节点 i -1和 i 的位置和方位相对变化矢量。系统任意一节点在惯性坐标系的速度和虚位移分别定义为

相应的节点相对的速度和虚位移定义为

式中， A 为到惯性坐标系的方位转换矩阵。

则节点 i 在惯性坐标系 XYZ 的中的位置和方位分别用节点 i 的位置矢量以及与节点 i -1的相对变形表示为

A 为旋转转换矩阵。当单元 i -1发生变形时，关于节点 i 的参考坐标系 x ₍ _i _-1) _i y ₍ _i _-1) _i z ₍ _i _-1) _i 必然会产生一个相对于参考坐标系 x _i _-1 y _i _-1 z _i _-1 的转动角度。节点 i 在任一时刻的方位 A _i 可表示为相对于总体惯性坐标系的3个欧拉变换矩阵的乘积，即

式中， A _i _-1 表示节点坐标系 x _i y _i z _i 到总体惯性坐标系的变换矩阵；矩阵 D ₍ _i _-1) _i 为节点 i 的变形造成的方位变化形成的变换矩阵； C ₍ _i _-1) _i 表示由于随体坐标系置于柔性体上，而定义的由 x _i y _i z _i 到 x ₍ _i _-1) _i y ₍ _i _-1) _i z ₍ _i _-1) _i 转换的常矩阵。

相对角速度的关系为

节点 i 的参考坐标系 x ₍ _i _-1) y ₍ _i _-1) _i z ₍ _i _-1) _i 可通过3次旋转达到节点 i -1的节点坐标系 x _i _-1 y _i _-1 z _i _-1 的坐标轴位置。旋转顺序具有多种形式，但不能绕一个轴连续旋转两次，因为连续旋转两次等同于绕这个轴的一次旋转。在工程中经常使用经典欧拉转动顺序的主要是“3-1-3”旋转和“1-2-3”旋转。本文采用“1-2-3”旋转顺序。将依次旋转的3个角坐标分别记为 ψ′ ₍ _i _-1) _i ₁ 、 ψ′ ₍ _i _-1) _i ₂ 和 ψ′ ₍ _i _-1) _i ₃ ，方位相对变化矢量写成矩阵的形式

参考坐标系 x ₍ _i _-1) y ₍ _i _-1) _i z ₍ _i _-1) _i 3次转动后的位置矢量相对于转动前的位置矢量相当于3个方向余弦矩阵的乘积，相应的方位变化形成的变换矩阵 D ₍ _i _-1) _i 为3个转换矩阵的乘积，即

式中， D _τ ( ψ′ ₍ _i _-1) _iτ )表示参考坐标系 x ₍ _i _-1) _i y ₍ _i _-1) _i z ₍ _i _-1) _i 做 ψ′ ₍ _i _-1) _iτ ( τ =1,2,3)角度旋转时的变换矩阵。

根据三维刚体系统动力学可知

从而得到：

定义 r _k （ k = i , i -1）对应的虚位移为δ r _k ，对应的旋转虚位移为δ π ，利用虚位移原理，则式（2-103）存在如下的相应变分形式为

式（2-112）给出了 δ r ′ _i 与虚位移 δ r ′ _i _-1 的变换关系，式中的和分别表示由矢量分量组成相应的矩阵 s ′ _（ _i _-1） _i ₀ 和 u ′ _（ _i _-1） _i 的斜对称矩阵（或者称为反对称方阵）。

定义3个变换矩阵 A ₍ _i _-1) _i 、 A ₍ _i _-1) 、 A _i 的相互关系为

则

将式（2-114）代入式（2-104）得到

节点 i -1和节点 i 的相对旋转虚位移用δ π ′ _i 表示，式（2-103）相应变分形式经整理得到

式（2-116）表明了旋转虚位移δ π ′ _i 和δ π ′ _i _-1 之间的变换关系，其中

式（2-116）和式（2-115）联立，整理后写成矩阵的形式

式（2-118）体现了两个相邻单元节点的虚位移的关系。

这样沿着设定的路径顺序，反复利用式（2-118）就可以得到由 n 个单元组成的整个划分系统在绝对节点坐标下和相对节点坐标系下相对虚位移的关系

式（2-119）显示了虚位移在笛卡儿坐标系和相对节点坐标系的关系。

其中

同样对式（2-101）和式（2-102）求导联立得到

可知矩阵 B _（ _i _-1） _i ₁ 、 B _（ _i _-1） _i ₂ 就是关于相对速度到笛卡儿坐标系的转换矩阵。

得到系统在笛卡儿坐标系下的速度同相对速度的关系

式中

式中， nc 和 nr 分别表示笛卡儿坐标系和相对坐标系的个数。

对式（2-128）求微分得到节点的加速度

2.4.3 节点弹性力的计算

钢丝绳作为一个研究对象，从系统上讲，是一个典型柔性体的大变形，即几何非线性问题（不考虑材料的非线性）。这时结构的位移变化产生的二次内力不能够被忽略，这样结构的荷载、变形关系表现为非线性。对于这样的结构应该按照变形以后的位置来建立平衡方程。处理几何大变形的方法主要有总体拉格朗日法（T.L法）或者修正的拉格朗日法（U.L法）。U.L法不仅适用于小应变情况，由于在计算大应变时引入了分段线性化的几何关系以及非线性本构关系，使得它能够用于非弹性大应变分析，所以在工程中修正的拉格朗日法的适用范围更广，应用也更多一些。下面简单地说明应用修正的拉格朗日法建立空间梁单元的刚度矩阵的过程。对于本文的研究对象空间梁单元，根据Kirchhoff（克希霍夫）假定，等截面的梁单元应变能是由两部分组成的，即线性应变能和非线性应变能。同时注意到，对由大量节点组成的一个有限元计算模型，单元的应变能仅仅同单元内部节点的相对位移有关，而与结构的刚性运动无关。计算中应该注意的是，大部分缆索结构往往需要考虑初始预应力对刚度的影响。因此，在计算单元的应变能求结构的刚度矩阵时，必须考虑结构的几何非线性。

设系统中任意的一个空间梁单元 i -1上任一点在 t 时刻和 t +Δ t 时刻相对于节点参考坐标系 x _i _-1 y _i _-1 z _i _-1 的 x _i _-1 、 y _i _-1 、 z _i _-1 轴的变位增量和应变增量分别为Δ u 和Δ ε 。由于钢丝绳缆索类的结构仅端部受力，仅需考虑一个正应变和两个剪应变，由Green应变定义得

式（2-132）中的第1项是位移增量的线性函数，其他3项为位移增量的非线性函数。引入梁位移的型函数，则位移增量Δ u 与梁端的位移增量{Δ u ^e }（三个线位移和三个角位移）的关系可以用式（2-133）表示

则可以得到

即

令

则得到

在U.L.法的坐标系中，Kirchhoff应力 S 即为此时刻 t 的欧拉应力 σ 与应力增量Δ s 的和，即

不考虑材料的非线性，设材料的增量本构关系满足增量形式

式中，( D )为材料的弹性矩阵。

根据Kirchhoff假定，在 t +Δ t 时刻梁单元的虚功方程可以写作

式中，Δ s 表示Kirchhoff应力增量； p 、 q 分别表示作用于梁单元上的体力、面力。

得到

注意到式（2-143）中的｛Δ u ^e }为任意值，所以

式（2-144）的右侧表示单元上体力和面力的等效节点力。

把式（2-144）代入式（2-143）的左侧进行变换，得

对式（2-145）右侧第1项做进一步的变换得

令

把方程用矩阵的形式表示

其中，

式中，（ K _T ）为单元的切线刚度矩阵；（ K _L ）为单元小位移的线性刚度矩阵；（ K _σ ）为单元关于应力｛ σ }的对称矩阵，称为初应力矩阵或几何刚度；（ K _nL ）为大位移刚度矩阵。

对于非线性刚度矩阵，［ K _nL ］的求解分析在现有计算条件下十分困难，当结构划分的单元足够多时，单元的非线性刚度矩阵可以忽略不计，这对于大多数的工程问题是允许的。这样单元的刚度矩阵( K _T )就是由两部分组成：线性刚度矩阵和几何刚度矩阵。

即

在求得单元的刚度矩阵以后，单元 i -1在相对节点坐标下的应变能写成虚功的形式为

式中， q _（ _k _-1） _k 为节点 k -1与 k 的相对位移；表示由节点 k -1与 k 组成的梁单元 k -1的切线刚度矩阵。

单元在 i -1节点位置的矢量定义为

单元 i -1整个的位置矢量在相对节点坐标系统中可以写作

单元的弹性力可以由式（2-157）得到

系统的应变能可以写成单元应变能和虚功的形式

2.4.4 钢丝绳接触模型

钢丝绳的三维接触模型如图2-12所示。

图2-12 三维接触模型

设接触单元与摩擦轮的接触面的侵入值为 d ，相应地产生一个沿着矢量 n 与摩擦轮表面垂直的反力，此反力与侵入值 d 和成比例变化。根据前面梁单元的型函数可以知道，侵入值是一个沿单元长度方向上关于 x 的函数，同样在摩擦轮上产生的沿着切向与法向的接触力也是关于 x 的函数。其中法向的接触力可以写作

式中， n 为摩擦轮接触面的单位法向力； d 为允许侵入的最大值； k _p 、 c _p 分别表示弹性系数和阻尼系数。

侵入值定义为

式中， R 为摩擦轮的半径； r 为接触节点在总体坐标系中的矢量； r ₀ 为摩擦轮圆心在总体坐标系中的矢量。

法向矢量 n 可以定义为

根据三线蠕变速率摩擦模型，摩擦轮的切向摩擦力可以表示为

式中， μ ( v ( t ))为取决于相对速度 v ( t )的摩擦系数； t 表示垂直于法矢量 n 的单位矢量。

可以得到

接触力可以表示为

对于摩擦式提升机来说，假定摩擦轮的角速度为 ω ，在单元中线任意一点的相对切向速度为

摩擦轮的平衡方程为

式中， I 为摩擦轮的惯性矩；为摩擦轮的角加速度； I _i 为单元 i 的单元长度； T 为摩擦轮的阻力矩。

对于二维平面梁单元，其接触力的虚功为

式中， S ₀ 表示单元的型函数； Q _c 表示总接触力。

在系统中分散的接触力转化为总的接触力为

可知接触力为单元在长度方向上的积分。需要注意的是在接触情况下，对于同一个单元上不同的节点上的速度和受力是不同的。因为高斯积分方法用较少的积分点就能够达到较高的精度，从而大大节省了计算的时间，所以采用高斯求积法计算接触力，节点的接触力可以作为高斯积分点进行求解。由于接触力是非光滑曲线，因此必须保证足够多的积分点。由于采用了绝对节点坐标方程，在计算接触的过程中，相对很容易设定钢丝绳单元与摩擦轮发生接触的条件。钢丝绳单元除了搭载在摩擦轮上面的位置以外，都不需要考虑计算接触力。

2.4.5 系统运动方程求解

在相对独立的坐标系下，研究的对象为受约束的多柔性体系统，系统运动方程在总体惯性坐标系下的基本方程，又称欧拉-拉格朗日方程，可以写成下面微分代数方程的形式

把 δ Z = B δ q 代入得

式中， M 为质量矩阵。

式中， n 表示节点的个数； I 表示单位矩阵； J ′ 表示惯性矩； f ′ 表示外力； n ′ 表示外扭矩。对于梁单元的节点，其惯性矩为

式（2-171）中， Φ 为系统的运动学约束和驱动约束矩阵； λ ∈ R ^m 为拉格朗日乘子矢量； Q 为包含外力、接触力、应变引起的弹性力，以及由于速度变化引起的力的矢量矩阵

由于δ q 是任意的，所以可以得到

则受约束的运动动力方程可以用隐式的形式

式中， q ∈ R ⁿ ⁺¹ 为广义坐标矢量； λ ∈ R ^m 为约束拉格朗日乘子矢量； Φ ( q )∈ R ^m 为位置约束方程， Φ ( q )的雅克比矩阵 Φ _q ( q , t )=∂ Φ / ∂ q ∈ R ^m ⁽ ⁿ ⁺¹⁾ 为约束雅可比矩阵， t ∈ R 是时间。

引入速度矢量和加速度矢量 a

对式(2-181a)逐次取微分得到速度和加速度的约束矩阵

式（2-179）构成的方程组增广广义坐标 x 的维数为3（ n +1）+ m ，方程组维数为3（ n +1）+3 m ，这是一个超定问题，即所谓超定微分-代数方程组（ODAEs-overdetermined differential algebraic equations）。根据解的存在性定理中的假设，对式（2-181a）采用微分流形“投影”技术以消除超定性

式中， R ₁ 、 R ₂ ∈ R ^（ ⁿ ^+1）（ ⁿ ⁺¹ ^-m ^）。如果取得有效的 R ₁ 和 R ₂ ，那么式（2-180c）从 n +1维降到（ n +1 -m ）维，消除了系统冗余性。

做如上处理后，求解方程构成的方程组转化为求解方程组

式中， x 为增广广义坐标， x =（ q ^T ，， a ^T ， λ ^T ） ^T ，此时， x ∈ R ^3（ ⁿ ^+1）+ ^m ， H ₁ ∈ R ^3（ ⁿ ^+1）+ ^m ，消除了系统超定性。 U _i ∈ R ^（ ⁿ ^+1）（ ⁿ ⁺¹ ^-m ^）（ i =1，2），构成了速度和位移约束的基础。 U _i 的选择必须使得矩阵非奇异阵。因此 U _i 和的空间参数构成了整个空间 R ⁿ ⁺¹ 。