CN115081146A - 一种颗粒阻尼动力吸振器耦合数值计算方法 - Google Patents

Abstract

本发明属于减振降噪领域，具体地说，是一种颗粒阻尼动力吸振器耦合数值计算方法，使用MBD和DEM算法分别对每个机械部分和每个颗粒的运动方程进行数值积分析，用直接积分法将二阶运动微分方程转化为一阶微分方程，然后用龙格‑库塔方法求解初值，再用耦合多体动力学和离散元模型计算，本发明采用了多体动力学和离散元耦合的方法，从多体动力学和离散元的原理出发，所建立的颗粒阻尼动力吸振器耦合方法，结果准确，能够充分的反映出结构的运动状态，为颗粒阻尼动力吸振器计算研究提供有效的技术支持。

Description

一种颗粒阻尼动力吸振器耦合数值计算方法

技术领域

本发明属于减振降噪领域，具体地说，是一种颗粒阻尼动力吸振器耦合数值计算方法。

背景技术

目前，基于调谐质量的强烈振动可以改善颗粒阻尼的阻尼效果，而阻尼效果的改善也将促进动力吸振器的振动控制性能，因此动力吸振器和颗粒阻尼的结合具有一定的研究意义。

耦合的颗粒阻尼动力吸振器的总阻尼效果是颗粒阻尼和动力吸振器两系统相互促进的，而不是动力吸振器和颗粒阻尼的线性叠加。因此，有必要讨论颗粒阻尼动力吸振器的阻尼特性。目前，现有技术中还没有对颗粒阻尼动力吸振耦合的复合阻尼器的仿真学进行研究过，也没有提出相关的算法。

发明内容

为了解决颗粒阻尼动力吸振器的耦合仿真问题，本发明提出一种颗粒阻尼动力吸振器耦合数值计算方法，该计算方法是基于多体动力学和离散元的仿真方法，可以高效准确的对颗粒阻尼动力吸振器进行计算仿真研究，对颗粒阻尼动力吸振器的数值计算具有重要的意义，具有良好的市场应用前景，弥补市场空白。

本发明采用的具体技术方案如下：

一种颗粒阻尼动力吸振器耦合数值计算方法，包括以下步骤：

步骤一：给出机械零件和阻尼颗粒的几何参数和材料特性，并设置了初始条件和边界条件；MBD计算t_MBD和DEM计算t_DEM的物理时间在计算开始时设置为零；

步骤二：由于MBD计算的时间步长Δt_MBD远大于DEM计算的时间步长Δt_DEM，因此DEM计算首先开始，并遵循DEM循环，DEM循环流程如下：

流程a、通过接触算法计算颗粒之间以及颗粒与机械零件之间的接触位移；

流程b、利用方程式确定接触力：

流程c、计算每个颗粒的线加速度和角加速度；

流程d、更新每个颗粒的位置和速度；

流程e、如果DEM计算t_DEM的物理时间小于下一个MBD步骤的物理时间，计算仍然遵循DEM循环，从流程a重复到流程d，否则，调用MBD方案，计算转到步骤三。

步骤三：DEM将机械零件上的等效颗粒力传递给MBD。在DEM循环之间的每个MBD时间步长执行MBD方案。MBD循环流程如下：

流程a、计算广义力。

流程b、通过化简方程求解MBD的控制方程。

流程c、利用RK4方法更新机械零件的位置和速度。

流程d、一旦MBD计算在每个时间步Δt_MBD结束，回到步骤二，计算相应的平移和旋转的机械零件发送回DME循环，然后确定阻尼颗粒的运动行为和等效颗粒部分再次机械零件。

本发明的有益效果：本发明相对于实验来说能明显的缩短产品开发周期、减少实验次数、降低产品的开发成本；动力吸振器和颗粒阻尼器作为低频振动控制的重要措施，动力吸振器能在不移动机械设备、不改变基座结构、不修改船体结构的情况下，添加颗粒阻尼系统吸收主系统振动能量，抑制结构振动，因此研究颗粒阻尼和动力吸振器耦合是抑制低频噪声的一个突破点。

附图说明

图1是本发明颗粒阻尼动力吸振器示意图。

图2a是颗粒与颗粒碰撞原理图。

图2b是颗粒与壁面碰撞原理图。

图3a是颗粒与颗粒碰撞原理图。

图3b是颗粒与壁面碰撞原理图。

图4是颗粒阻尼动力吸振器耦合流程图。

图5是多体动力学和离散元耦合仿真和实验对比图。

具体实施方式

为了加深对本发明的理解，下面将结合附图和实施例对本发明做进一步详细描述，该实施例仅用于解释本发明，并不对本发明的保护范围构成限定。

本发明实施例是对颗粒阻尼动力吸振器进行多体动力学和离散元耦合数值计算，对其方法进行分析，通过本实施例的分析，对本发明做详细的说明，具体实施方式分为以下步骤：

MBD和DEM使用算法，分别对每个机械部分和每个颗粒的运动方程进行数值积分析。用直接积分法将二阶运动微分方程转化为一阶微分方程。然后用龙格-库塔方法求解初值，图4为耦合多体动力学和离散元(MBD-DEM)模型的计算流程图，具体流程如下：

步骤一：给出了机械部件和阻尼颗粒的几何参数和材料特性，并设置了初始条件和边界条件；如表1，MBD计算t_MBD、DEM计算t_DEM的物理时间，在计算开始时设置为零。

步骤二：由于MBD计算的时间步长(Δt_MBD)远大于DEM计算的时间步长(Δt_DEM)，因此DEM计算首先开始，并遵循DEM循环。DEM循环如下：

a)通过接触算法计算颗粒之间以及颗粒与机械零件之间的接触位移(重叠)；

b)利用方程式(1)-(13)确定接触力，具体步骤如下：

1、颗粒在离散元(DEM)中，每个颗粒的运动通常由颗粒i表示，遵守牛顿第二运动定律，方程(1)描述质心G的平移运动，方程(2)描述整个质点的旋转运动。

其中m⁽ⁱ⁾是颗粒i的质量，I⁽ⁱ⁾是相应的关于质心G的质量惯性矩，

是颗粒i的位置矢量，

是平移加速度，

是角速度，

是角加速度。

是我在接触点的接触力作用在颗粒j,

是重力加速度,

上的接触点颗粒的位置向量i在接触j。

是颗粒上的联系人数量我和象征“×”表示向量叉乘。

2、在颗粒动力学中，颗粒之间以及颗粒与机械部件之间的接触是一种直接且简单的方法。图2和图3表示出颗粒之间以及颗粒与机械部件之间的接触示意图。颗粒间接触的判据为

否则为无接触。同样，颗粒-壁面接触的接触判据由

给出，否则为无接触。在上述不等式中，

为颗粒j的位置向量，

为接触点在壁上的位置向量，r_i，r_j分别表示颗粒i和j的半径。

3、考虑两个接触的颗粒，即颗粒i和j，m_i和m_j分别表示颗粒i、j的质量；r_i和r_j表示这些颗粒各自的半径。分别定义m^*和r^*为等效质量和半径，表示为m^*＝m_im_j/(m_i+m_j)和r^*＝r_ir_j/(r_i+r_j)。两个颗粒具有相同的材料性质，剪切模量为G，泊松比为v。当颗粒i与颗粒j以非零的相对速度斜向碰撞时，颗粒i与颗粒j之间存在接触力。作用在颗粒

上的接触力包括法向分量

和切向分量

表示为：

4、各分量用弹簧力和阻尼力之和表示，表示为：

F_cn＝F_cn，s+F_cn，d (4)

F_ct＝F_ct，s+F_ct，d (5)

5、在法向上，式(4)中的F_cn，s，F_cn，d分别为法向弹簧力和法向阻尼力。

式中，δ_n为法向接触位移，

为相对法向速度(j颗粒相对于i颗粒)，ξ_n为法向阻尼比，k_n为法向接触刚度。在切向方向上，式5中的F_ct，s，F_ct，d分别为切向弹簧力和切向阻尼力，表示为：

F_ct，s＝F_{ct，s，(n-1)}+ΔF_ct，s (8)

6、ΔF_ct，s和F_{ct，s，(n-1)}分别为当前和之前时间步长的切向弹簧力，ΔF_ct，s为切向弹簧力的增量。ξ_t为切向阻尼比，k_t为切向接触刚度，

为相对切向速度(颗粒j相对颗粒i)。增量ΔF_ct，s进一步由ΔF_ct，s＝-k_sΔδ_t给出，Δδ_t为增量切向位移，表示为

其中Δt为时间步长。

7、为简单起见令ξ_n＝ξ_t＝ξ，ξ可以表示为

其中e为恢复系数。k_n和k_t可进一步表示为：

切向力受库仑摩擦定律限制，其表达式为：

|F_ct|≤μF_cn (11)

其中μ是摩擦系数，该接触力模型所需的参数包括剪切模量(G)、泊松比(υ)、颗粒间

的摩擦系数(μ_b)、颗粒与机械部件之间的摩擦系数(μ_w)、颗粒间的恢复系数(e_b)、颗粒与机械部件之间的恢复系数(e_w)。

8、根据牛顿力学方程，含颗粒阻尼动力吸振器方程可以表示为下式：

在方程式(12)中，FCV(t)分别表示颗粒对盒的垂直合力，如式(13)所示。

式中，F_cv，j为颗粒盒第j个接触点接触力的垂直分量，N_contact为颗粒盒上的总接触数，法向力F_cn一般可以表示为

其中K为接触刚度；χ为滞回阻尼因子；a和b是幂指数，右边的第一项和第二项分别表示弹性力分量和阻尼力分量。

c)利用式(1)、(2)计算每个颗粒的线加速度和角加速度；

d)通过显式积分更新每个颗粒的位置和速度。平动速度和角速度用中间间隔的速度值表示为方程式(14)和(15)，该位置如式(16)所示。

e)如果DEM(t_DEM)的物理时间小于下一个MBD步骤的物理时间(t_MBD+Δt_MBD)，计算仍然遵循DEM循环，从步骤2(a)重复到步骤2(d)。否则，调用MBD方案，计算转到步骤三。

步骤三：DEM将机械零件上的等效颗粒力传递给MBD。在DEM循环之间的每个MBD时间步长执行MBD方案，MBD循环如下：

a)计算广义力；

b)通过化简方程求解MBD的控制方程；

c)利用(RK4)方法更新机械零件的位置和速度；

d)一旦MBD计算在每个时间步Δt_MBD结束，回到步骤(2)，计算相应的平移和旋转的机械零件发送回DME循环，然后确定阻尼颗粒的运动行为和等效颗粒部分再次机械零件；

数值计算具体实施例数据如下表1所示，振动台对基座施加0.12g的加速度，用上述方法，对颗粒阻尼动力吸振器系统进行数值计算，分别对无动力吸振器、有动力吸振器、颗粒阻尼动力吸振器进行仿真和实验的对比，实验与仿真对比如图5所示，可以看出有动力吸振器和无动力吸振器相比加速度的峰值下降了37.5％，加颗粒之后，加速度峰值下降明显，含颗粒动力吸振器和有动力吸振器相比两个峰的加速度峰值分别下降了78％和89％，实验和仿真也吻合较好。

表1动力吸振和颗粒阻尼耦合研究模型参数

以上显示和描述了本发明的基本原理、主要特征及优点。本行业的技术人员应该了解，本发明不受上述实施例的限制，上述实施例和说明书中描述的只是说明本发明的原理，在不脱离本发明精神和范围的前提下，本发明还会有各种变化和改进，这些变化和改进都落入要求保护的本发明范围内。本发明要求保护范围由所附的权利要求书及其等效物界定。

Claims (8)

1.一种颗粒阻尼动力吸振器耦合数值计算方法，其特征在于，使用MBD和DEM算法分别对每个机械部分和每个颗粒的运动方程进行数值积分析，用直接积分法将二阶运动微分方程转化为一阶微分方程，然后用龙格-库塔方法求解初值，再用耦合多体动力学和离散元模型计算。

2.根据权利要求1所述的颗粒阻尼动力吸振器耦合数值计算方法，其特征在于，具体包括以下步骤：

步骤一、给出了机械部件和阻尼颗粒的几何参数和材料特性，并设置了初始条件和边界条件，MBD计算t_MBD、DEM计算t_DEM的物理时间，在计算开始时设置为零；

步骤二、由于MBD计算的时间步长Δt_MBD大于DEM计算的时间步长Δt_DEM，因此DEM计算首先开始，并遵循DEM循环；

步骤三、DEM将机械零件上的等效颗粒力传递给MBD，在DEM循环之间的每个MBD时间步长执行MBD循环方案。

3.根据权利要求2所述的颗粒阻尼动力吸振器耦合数值计算方法，其特征在于，所述步骤二中，DEM循环流程如下：

流程a、通过接触算法计算颗粒之间以及颗粒与机械零件之间的接触位移；

流程b、利用方程式确定接触力：

流程c、计算每个颗粒的线加速度和角加速度；

流程d、更新每个颗粒的位置和速度；

流程e、如果DEM计算t_DEM的物理时间小于下一个MBD步骤的物理时间，计算仍然遵循DEM循环，从流程a重复到流程d，否则，调用MBD方案，计算转到步骤三。

4.根据权利要求3所述的颗粒阻尼动力吸振器耦合数值计算方法，其特征在于，所述步骤二中DEM循环的流程b利用方程式确定接触力，具体方法如下：

颗粒在离散元DEM中，每个颗粒的运动由颗粒i表示，遵守牛顿第二运动定律：方程(1)描述质心G的平移运动，方程(2)描述整个质点的旋转运动：

其中m⁽ⁱ⁾是颗粒i的质量，I⁽ⁱ⁾是相应的关于质心G的质量惯性矩，

是颗粒i的位置矢量，

是平移加速度，

是角速度，

是角加速度，

是i在接触点的接触力作用在颗粒j，

是重力加速度，

上的接触点颗粒的位置向量i在接触j；

颗粒间接触的判据为

否则为无接触，同样，颗粒-壁面接触的接触判据由

给出，否则为无接触，在上述不等式中，

为颗粒j的位置向量，

为接触点在壁上的位置向量，r_i，r_j分别表示颗粒i和j的半径；

考虑两个接触的颗粒，即颗粒i和j，m_i和m_j分别表示颗粒i、j的质量；r_i和r_j表示这些颗粒各自的半径；分别定义m^*和r^*为等效质量和半径，表示为m^*＝m_im_j/(m_i+m_j)和r^*＝r_ir_j/(r_i+r_j)，两个颗粒具有相同的材料性质，剪切模量为G，泊松比为v，当颗粒i与颗粒j以非零的相对速度斜向碰撞时，颗粒i与颗粒j之间存在接触力，作用在颗粒

上的接触力包括法向分量

和切向分量

表示为：

各分量用弹簧力和阻尼力之和表示，表示为：

F_cn＝F_cn，s+F_cn，d (4)

F_ct＝F_ct，s+F_ct，d (5)

在法向上，式(4)中的F_cn，s，F_cn，d分别为法向弹簧力和法向阻尼力；

式中，δ_n为法向接触位移，

为相对法向速度，ξ_n为法向阻尼比，k_n为法向接触刚度，在切向方向上，式5中的F_ct，s，F_ct，d分别为切向弹簧力和切向阻尼力，表示为：

F_ct，s＝F_{ct，s，(n-1)}+ΔF_ct，s (8)

ΔF_ct，s和F_{ct，s，(n-1)}分别为当前和之前时间步长的切向弹簧力，ΔF_ct，s为切向弹簧力的增量，ξ_t为切向阻尼比，k_t为切向接触刚度，

为相对切向速度，增量ΔF_ct，s进一步由ΔF_ct，s＝-k_sΔδ_t给出，Δδ_t为增量切向位移，表示为

其中Δt为时间步长；

令ξ_n＝ξ_t＝ξ，ξ表示为

其中e为恢复系数，k_n和k_t进一步表示为：

切向力受库仑摩擦定律限制，其表达式为：

|F_ct|≤μF_cn (11)

其中μ是摩擦系数，该接触力模型所需的参数包括剪切模量(G)、泊松比(v)、颗粒间的摩擦系数(μ_b)、颗粒与机械部件之间的摩擦系数(μ_w)、颗粒间的恢复系数(e_b)、颗粒与机械部件之间的恢复系数(e_w)。

5.根据权利要求4所述的颗粒阻尼动力吸振器耦合数值计算方法，其特征在于，所述步骤二中，根据牛顿力学方程，含颗粒阻尼动力吸振器方程表示为下式：

在方程式(12)中，F_CV(t)分别表示颗粒对盒的垂直合力，如式(13)所示：

式中，F_cv，j为颗粒盒第j个接触点接触力的垂直分量，N_contact为颗粒盒上的总接触数，法向力F_cn表示为

其中K为接触刚度；χ为滞回阻尼因子；a和b是幂指数，右边的第一项和第二项分别表示弹性力分量和阻尼力分量。

6.根据权利要求5所述的颗粒阻尼动力吸振器耦合数值计算方法，其特征在于，所述步骤二中DEM循环的流程c中利用式(1)、(2)计算每个颗粒的线加速度和角加速度。

7.根据权利要求6所述的颗粒阻尼动力吸振器耦合数值计算方法，其特征在于，所述步骤二中DEM循环的流程d通过显式积分更新每个颗粒的位置和速度，平动速度和角速度用中间间隔的速度值表示为方程式(14)和(15)，该位置如式(16)所示：

8.根据权利要求7所述的颗粒阻尼动力吸振器耦合数值计算方法，其特征在于，所述步骤三中MBD循环方案的流程具体如下：

流程a、计算广义力；

流程b、通过化简方程求解MBD的控制方程；

流程c、利用RK4方法更新机械零件的位置和速度；

流程d、一旦MBD计算在每个时间步Δt_MBD结束，回到步骤二，计算相应的平移和旋转的机械零件发送回DME循环，然后确定阻尼颗粒的运动行为和等效颗粒部分再次机械零件。

Images