CN112199883B - 基于有效碰撞耗能计算的多颗粒阻尼器减振机理分析方法 - Google Patents

Abstract

本发明提出了一种基于有效碰撞耗能计算的多颗粒阻尼器减振机理分析方法，具体为建立了一种可定量反映颗粒阻尼器在工作状态下伴随发生的各类能量交换，及其内部颗粒团相对于空腔宏观运动状态的数值仿真方法，并给出了操作流程。包括有效碰撞模型的构建，颗粒阻尼器数值仿真模型的构建与求解，能量耗散及有效碰撞状态参数曲线提取方法等内容。基于提出的“有效碰撞”的计算模型，借助颗粒阻尼离散元程序(DEM)追踪颗粒与空腔在任意时间的运动状态与轨迹，给出不同响应加速度条件下颗粒团相对于空腔的宏观运动状态描述及各类能量耗散占能量交换总量的百分比，更为清晰地显示了各加速度条件下颗粒阻尼的减振机理。

Description

基于有效碰撞耗能计算的多颗粒阻尼器减振机理分析方法

技术领域

本发明属于被动减振技术数值仿真领域，具体为建立了一种可定量反映颗粒阻尼器在工作状态下伴随发生的各类能量交换，及其内部颗粒团相对于空腔宏观运动状态的数值仿真方法。

背景技术

由于极宽的温度适用范围，颗粒阻尼常被用于航空航天领域中极端工况下的振动抑制，并显示出很好的减振效果。已有研究表明，颗粒阻尼的减振效果可被绘制成随结构响应加速度量级变化的特征曲线，且具有高度非线性。为更好开展颗粒阻尼器设计，就必须对颗粒阻尼减振机理进行细致分析，找出影响减振效果的最大因素并加以利用。曾有研究者将阻尼器中所有离散颗粒等效为一个光滑的“集中质量”，利用颗粒与空腔碰撞次数的变化定性解释颗粒阻尼特征曲线的变化趋势。颗粒与空腔的碰撞次数可以体现两者之间发生动量交换的频次，但并不是只要颗粒与空腔发生碰撞，颗粒就会对空腔做负功，碰撞次数越多，也并不能完全表明颗粒对空腔所做的负功越多，阻尼器的减振效果越好。另外，从颗粒碰撞次数这种过于细观的角度去解释颗粒团体现出的宏观减振效果也还缺少一些可以反映颗粒宏观运动状态的参数作支撑，说服力不足。

本发明基于颗粒团相对于空腔的不同运动状态，建立了“有效碰撞”模型，同时借助颗粒阻尼离散元程序(DEM)追踪颗粒与空腔在任意时间的运动状态与轨迹，给出不同响应加速度条件下颗粒团相对于空腔的宏观运动状态描述及各类能量耗散占能量交换总量的百分比，更为清晰地显示了各加速度条件下颗粒阻尼的减振机理。在案例研究中，本发明利用DEM对已公开发布的试验进行了数值模拟，显示了较高的准确度，验证了模型的有效性。

发明内容

本发明的目的是提供一种基于颗粒阻尼离散元程序的数值分析方法，该方法既可反映颗粒阻尼器内部颗粒与主结构宏观作用状态，又可给出各类能量耗散占总能量交换百分比，很好地解释了颗粒阻尼呈现出高度非线性特点的内在能量交换机理。

一、有效碰撞模型的定义

(一)碰撞的分类

图1给出了颗粒与空腔发生斜碰前的示意。V_p为颗粒的绝对速度，V_pn、V_pt分别为颗粒沿空腔法向及切向的分量，空腔绝对速度V_c也沿法向和切向分为V_cn、V_ct。以空腔为研究对象，将颗粒与空腔沿法向的碰撞作用分为两类：

I类碰撞：碰撞产生的法向作用力与空腔的法向运动方向相反，此时颗粒对空腔做负功，可以起到减振的作用。发生这类碰撞时，颗粒相对于与空腔可能有两种运动情况：第一种，颗粒与空腔在发生碰撞前的绝对速度方向相反，即“迎面碰撞”；第二种，碰撞发生前颗粒与空腔的绝对速度方向相同，但该运动方向却与颗粒相对于空腔的运动方向相反，即空腔“追赶”颗粒。两种运动情况均可用V_cn·(V_cn-V_pn)≥0表示。

II类碰撞：碰撞产生的法向作用力与空腔的法向运动方向相同。此时颗粒对空腔做正功，无法起到减振作用。发生该类碰撞时，颗粒与空腔在碰撞发生前的绝对速度方向相同，且该运动方向与颗粒相对于空腔的运动方向相同，即颗粒“追赶”空腔。该种运动情况可用V_cn·(V_cn-V_pn)<0表示。

同理，按照碰撞产生的切向力(即摩擦力)对空腔所做的功也可将空腔与颗粒沿切向的碰撞分为I类和II类碰撞，发生I类碰撞时摩擦力沿切向对空腔做负功，颗粒与空腔的运动情况可用V_ct·(V_ct-V_pt)≥0表示，发生II类碰撞时摩擦力沿切线方向对空腔做正功，相应的运动情况可用V_ct·(V_ct-V_pt)<0表示。

(二)碰撞的能量交换模型

将单个周期内颗粒对空腔所做的功W_pel按照I类、II类进行细分(如图2所示)，进而写成公式(1)。其中W_pnel、W_pfel分别为颗粒通过法向碰撞和摩擦对空腔所做功，W_pnelh、W_pnelb分别为法向碰撞所做正功及负功、W_pfelh、W_pfelb分别为摩擦所做正功及负功。

W_pel＝W_pnel+W_pfel＝W_pnelh+W_pnelb+W_pfelh+W_pfelb (1)

通过对颗粒在单个周期内对空腔所做的功进行时域积分求解式(1)中的各项参数。假设空腔壁面数为n_b，空腔内颗粒数为n_p，在单个周期时间T内有个颗粒—空腔接触对，则颗粒对空腔做功的表达式可写为

其中计算起始时间为t₀，单个周期内时间步个数为N，单个时间步长为Δt，f_tbi、f_nbi分别为颗粒i与空腔b间的法向和切向作用力矢量，角标j为时间步编号，l为接触对编号。

对上述能量计算结果按照结构固有阻尼耗能进行归一。定义总耗能比R_el(Totalenergy dissipation ratio)为结构总耗能(结构固有阻尼耗能E_cdl与颗粒碰撞耗能E_pel之和)与结构固有阻尼耗能E_cdl的比值，根据能量转化关系，其表达式可写为：

其中W_cdl为固有阻尼做功，W_pel为颗粒碰撞对空腔做功。

同理给出法向碰撞耗能比R_nel、摩擦耗能比R_fel的定义，见式(4)、(5)。

其中E_penl，E_pfel分别为颗粒正碰与摩擦耗能，W_penl，W_pfel分别为颗粒正碰与摩擦对空腔做功。

(三)有效碰撞状态参数的定义

有效碰撞状态参数包括有效法向碰撞状态系数β_n及有效摩擦状态系数β_f，分别用来描述单个周期内颗粒通过法向碰撞和摩擦对空腔所做的负功在总能量交换中所占的比重，两系数的表达式如下：

将式(6)、(7)代入式(3)，则总耗能比R_el又可写为

二、颗粒阻尼器数值仿真模型

本发明基于离散单元法(Discrete Element Method，DEM)建立颗粒阻尼器数值仿真模型。模型采用的基本假设有：a)所有颗粒均为球形，且几何参数和物性参数均相同；b)单个计算时间步足够小，以至于任何时刻由于碰撞出现的微小扰动均不会随着时间的推进而放大。在该时间步内，每一个颗粒的位移、速度等状态信息均保持不变；c)发生碰撞时，颗粒之间会产生微小变形量(法向及切向叠合量)，但与颗粒的几何尺寸相比仍为小量。

图3给出了程序可描述的基本振动模型，包括一个颗粒空腔及内部的若干颗粒。令u和w分别表示空腔相对于平衡位置的水平和竖直方向位移，可得空腔的运动方程为：

式中“.”表示对时间的一阶导数，M，K，C分别表示主结构的等效质量、等效刚度及等效阻尼，F_cu表示颗粒沿水平方向对空腔的作用力。方程最右边的F_ext、ω、分别表示激振力幅值、激振频率及相位，其值可根据不同的激振方式设定。在实际应用中，图3所示的三维空腔用于表征安装颗粒阻尼器的主结构。

模型中的空腔只有1个平动自由度，颗粒包含6个自由度(3个方向平动，3个方向转动)。在某一瞬间，颗粒i可能同时与相邻的若干颗粒或者空腔壁面发生碰撞。为便于分析，把颗粒i所受的碰撞力分为两个部分(这里不考虑颗粒间的微观作用力，例如范德华力等)：一部分是颗粒与颗粒之间的碰撞导致的合力与合力矩，即(颗粒)内部作用力，另一部分则是颗粒与空腔壁面发生碰撞形成的合力及合力矩，以及颗粒所受的重力、离心力等，统称为(颗粒)外部作用力。图4(a)和图4(b)给出了颗粒与颗粒以及颗粒与空腔间的碰撞模型，表1则列出了描述碰撞的相关参数，其中r_i，r_j，P_i，P_j，V_i，V_j，ω_i，ω_j分别表示颗粒i及颗粒j的粒径、位移矢量、速度矢量及角速度矢量。空腔壁面b的位移矢量、速度矢量及壁面单位法向量矢量则分别表示为P_b，V_b及n_ib。

表1颗粒与颗粒及颗粒与空腔之间碰撞的相关参数

颗粒i的运动方程如下：

其中，m_i表示颗粒i的质量，I_i为颗粒绕质心的转动惯量，g为重力加速度。θ_i为颗粒的角速度矢量，f_nij、f_nib分别为颗粒与颗粒以及颗粒与空腔间的法向作用力矢量，f_tij、f_tib则为对应的切向作用力矢量。n_pair、b_pair分别表示在某个时间计算步内，与颗粒i发生接触的颗粒总数目及空腔总数目。从图4(a)和图4(b)中可以看到，当颗粒发生斜碰时，作用在接触点的切向力会产生力矩，使颗粒发生转动。

想要实现对颗粒运动的数值模拟，首先应具有可精确刻画颗粒碰撞过程的接触本构模型(即力—位移关系)。本发明采用的法向接触模型为基于Hertz接触理论提出的TS非线性粘弹性模型。以颗粒i及颗粒j间的接触为例，其法向接触模型的表达式为：

可以看到，法向作用力f_n为线弹性力f_ns和反映能量耗散的阻尼力f_nd的合力。式中的法向碰撞阻尼常数η可以表示为恢复系数R_n的函数，k_np则是两球体接触时的等效弹簧刚度，其表达式为：

式中E_i、E_j和υ_i、υ_j分别表示颗粒i、j的弹性模量和泊松比。当颗粒与空腔壁面发生碰撞时，对应的等效弹簧刚度可以写成

其中E_c和υ_c分别表示空腔的弹性模量和泊松比。

程序采用的切向接触模型是目前普遍采用的库伦摩擦模型。以颗粒i与颗粒j为例给出摩擦力模型的表达式(14)，颗粒与空腔之间的表达式与之相同。

f_tij＝-μ_d|f_nij|τ_ij (14)

其中μ_d为颗粒与颗粒以及颗粒与空腔之间的摩擦系数。

本发明优点及有益效果在于：基于上述“有效碰撞”模型，本发明可追踪颗粒与空腔在任意时间的运动状态与轨迹，给出不同响应环境下颗粒团相对于空腔的宏观运动状态描述，进而求得各类能量耗散占能量交换总量的百分比，更为清晰地显示了各加速度条件下颗粒阻尼的减振机理。

附图说明

图1为颗粒与空腔碰撞示意图。

图2为单个周期内颗粒对空腔所做功的分类。

图3为DEM仿真程序的基本振动模型。

图4(a)为颗粒间碰撞示意图。

图4(b)为颗粒与空腔之间的碰撞示意图。

图5为数值结果与试验值的比较。

图6为强迫振动条件下特征曲线与总耗能比曲线的比较。

图7(a)为强迫振动条件下的能量分析曲线。

图7(b)为强迫振动条件下的能量分析曲线。

图7(c)为强迫振动条件下的能量分析曲线。

图7(d)为强迫振动条件下的能量分析曲线。

图8为Γ＝5.35时，监测颗粒与空腔沿水平方向的速度时域响应曲线。

图9为Γ＝188.50时，监测颗粒与空腔沿水平方向的速度时域响应曲线。

具体实施方式

在颗粒阻尼器设计阶段，研究者通常可通过有限的试验获取较好的方案设计。但为进一步优化遴选方案，寻求当前设计减振效果的宏观和微观解释，进而获得设计方案的可靠优化路径，可以采用本发明的“有效碰撞”模型予以辅助，具体实施步骤如下：

步骤一：建立颗粒阻尼器数值仿真模型；

明确颗粒阻尼器的各项参数，其中空腔参数包括：空腔弹性模量E_c和泊松比υ_c，空腔尺寸参数长L、宽W、高H，空腔质量M、等效刚度K及等效阻尼C。颗粒参数包括：颗粒数量N，颗粒粒径d，颗粒弹性模量E，泊松比υ，颗粒密度ρ,库伦摩擦系数μ接触阻尼常数η。

步骤二：利用离散单元法程序求解颗粒阻尼器响应；

首先，确定空腔及其内部颗粒的位置。对于颗粒i，如果存在接触对且其法向叠合量δ_n>0，则利用公式(11)、(14)计算相应的法向及切向接触力；若法向叠合量δ_n≤0，则所有的接触力均设置为0；接着，对所有作用在颗粒i上的作用力(包括内部及外部作用力)求矢量和；然后，利用公式(9)、(10)求解颗粒i的平动及转动加速度，并按照上述过程对其余每一个颗粒的运动方程进行求解；最后，求解作用在颗粒团上的合力。根据牛顿第三定律得到颗粒对空腔的作用力，将其沿水平和竖直方向分解，得到F_cu，然后利用公式(8)更新外部作用力及空腔运动状态，进入下一个循环。

步骤三：提取能量耗散及有效碰撞状态参数曲线；

有效碰撞模型的定义

碰撞的分类

图1给出了颗粒与空腔发生斜碰前的示意。V_p为颗粒的绝对速度，V_pn、V_pt分别为颗粒沿空腔法向及切向的分量，空腔绝对速度V_c也沿法向和切向分为V_cn、V_ct。以空腔为研究对象，将颗粒与空腔沿法向的碰撞作用分为两类：

I类碰撞：碰撞产生的法向作用力与空腔的法向运动方向相反，此时颗粒对空腔做负功，可以起到减振的作用。发生这类碰撞时，颗粒相对于与空腔可能有两种运动情况：第一种，颗粒与空腔在发生碰撞前的绝对速度方向相反，即“迎面碰撞”；第二种，碰撞发生前颗粒与空腔的绝对速度方向相同，但该运动方向却与颗粒相对于空腔的运动方向相反，即空腔“追赶”颗粒。两种运动情况均可用V_cn·(V_cn-V_pn)≥0表示。

II类碰撞：碰撞产生的法向作用力与空腔的法向运动方向相同。此时颗粒对空腔做正功，无法起到减振作用。发生该类碰撞时，颗粒与空腔在碰撞发生前的绝对速度方向相同，且该运动方向与颗粒相对于空腔的运动方向相同，即颗粒“追赶”空腔。该种运动情况可用V_cn·(V_cn-V_pn)<0表示。

同理，按照碰撞产生的切向力(即摩擦力)对空腔所做的功也可将空腔与颗粒沿切向的碰撞分为I类和II类碰撞，发生I类碰撞时摩擦力沿切向对空腔做负功，颗粒与空腔的运动情况可用V_ct·(V_ct-V_pt)≥0表示，发生II类碰撞时摩擦力沿切线方向对空腔做正功，相应的运动情况可用V_ct·(V_ct-V_pt)<0表示。

碰撞的能量交换计算模型

将单个周期内颗粒对空腔所做的功W_pel按照I类、II类进行细分(如图2所示)，进而写成公式(1)。其中W_pnel、W_pfel分别为颗粒通过法向碰撞和摩擦对空腔所做功，W_pnelh、W_pnelb分别为法向碰撞所做正功及负功、W_pfelh、W_pfelb分别为摩擦所做正功及负功。

W_pel＝W_pnel+W_pfel＝W_pnelh+W_pnelb+W_pfelh+W_pfelb (1)

通过对颗粒在单个周期内对空腔所做的功进行时域积分求解式(1)中的各项参数。假设在单个周期内有N个时间步，每个时间步内有个颗粒—空腔接触对，则颗粒对空腔做功的表达式可写为

对上述能量计算结果按照结构固有阻尼耗能进行归一。定义总耗能比R_el(Totalenergy dissipation ratio)为结构总耗能E_loss与结构固有阻尼耗能E_cdl的比值，其表达式为：

同理给出法向碰撞耗能比R_nel、摩擦耗能比R_fel的定义，见式(4)、(5)。

有效碰撞状态参数的定义

有效碰撞状态参数包括有效法向碰撞状态系数β_n及有效摩擦状态系数β_f，分别用来描述单个周期内颗粒通过法向碰撞和摩擦对空腔所做的负功在总能量交换中所占的比重，两系数的表达式如下：

将式(6)、(7)代入式(3)，则总耗能比R_el又可写为

强迫振动条件下颗粒阻尼特征曲线的提取方法

颗粒阻尼依赖于主结构响应幅值的大小，本发明采用曲线作为颗粒阻尼的特征曲线。/>为归一化阻尼比，表达式为

其中，ζ为安装阻尼器后的结构阻尼比，ζ₀为空白对照的阻尼比。

Γ为无量纲响应加速度，表达式为

Γ＝aω² / g (10)

其中a为位移响应幅值，ω为激振频率，g为重力加速度。

当颗粒阻尼器发生强迫振动时，提取颗粒阻尼特征曲线的关键在于获得安装颗粒阻尼器后结构的主响应频率ω_d及其响应峰值A_max，这两项数据可以通过正弦扫频获得。当以固定位移a₀方式对空腔施加正弦扫频激励时，式(8)中激振频率ω为时变量，其表达式为其中ω₀为扫频起始频率，/>为扫频速率，t'为扫频时间。对颗粒空腔基础位移激励，ω为变量，力激励幅值F_ext及相位/>也均为变量，其表达式分别为：

采用正弦扫频的方式仅用一次数值仿真便可较为准确地获得安装颗粒阻尼器后结构的主响应频率及响应峰值。此时利用幅频响应曲线得到的峰值点可以获得无量纲加速度Γ，利用半功率带法得到归一化的结构阻尼比为获得多个无量纲加速度Γ条件下的还需通过多个量级的慢速正弦扫频得到颗粒阻尼特征曲线，然后利用求得的共振频率点对阻尼器施加单频简谐激励，方能得到稳定状态下颗粒阻尼的耗能比曲线。

案例

构建如下颗粒阻尼器模型：颗粒空腔质量M＝0.293kg，等效阻尼系数C＝0.116，等效刚度K＝1602.7N/m(可知结构固有模态频率为f_n＝11.771Hz，结构固有阻尼比ζ＝0.268％)。空腔弹性模量E_c＝281.85MPa，泊松比υ_c＝0.38。空腔长L＝58mm，宽H、高W均为38mm。颗粒空腔内填充200个直径d＝6mm的丙烯酸树脂颗粒，其质量率m_r＝9.2％。颗粒密度ρ＝1190kg/m³，颗粒弹性模量E_p＝331.37MPa，泊松比υ_c＝0.38，库伦摩擦系数μ＝0.52，接触阻尼常数η＝0.077。

对上述模型施加量级为a₀＝0.5mm、1.0mm的基础位移激励，并将计算结果与已有试验值进行了比较。由图5可知，建立的数值模型与试验值吻合较好，验证了利用DEM方法在计算强迫振动条件下的颗粒阻尼时，同样具有较高的准确性。

计算强迫振动状态下的总耗能比R_el，并将其与进行比较，如图6所示。可知总耗能比与归一化阻尼比两条曲线基本重合，因此可用能量分析的方法解释特征曲线的变化趋势。

图7(a)—图7(d)给出了对应的能量分析曲线。由图7(a)可知，强迫振动条件下颗粒阻尼特征曲线先增后减的趋势取决于颗粒耗能同固有阻尼耗能的相对关系，而在振动过程中，颗粒与空腔的法向碰撞和摩擦同样既包含对空腔做正功的成分，也包含做负功的成分(如图7(b)图所示)。当Γ较低时，摩擦也会起到耗能的作用。由图7(c)图可知，在Γ<10的范围内，摩擦耗能比R_fel大于0，其变化趋势与正碰耗能比R_nel一样，先增后减。在图7(d)图中可以看到，在强迫振动条件下，β_n及β_f随Γ先增后减，直到Γ>50后才逐渐趋于稳定。

由于图7(a)—图7(d)曲线上每一个有效碰撞状态参数的数据点均对应结构达到稳态后的结果，因此其对应的取值反映了在强迫振动条件下，与不同Γ对应的颗粒与空腔间不同的宏观相对运动状态。在颗粒团中随机选取了2个监测颗粒，用其速度曲线反映颗粒团整体的速度变化情况。

首先给出Γ＝5.35条件下，监测颗粒与空腔的时域曲线。此时β_n＝0.996，即几乎所有的法向碰撞均属于I类碰撞。图8中横坐标为周期数，纵坐标为水平方向的绝对速度值。可以看到，1#、2#监测颗粒所体现出的主要运动规律一致，此时可将颗粒团看做整体来描述其与空腔的相对运动状态。可以看到一个振动周期可被分成了6个阶段：

1～2：颗粒团与空腔壁板发生迎面碰撞(属于I类碰撞)，颗粒运动速度反向并达到极大值；2～3：空腔保持原有方向继续运动，持续与颗粒团作用(主要为I类碰撞)，颗粒团与空腔同向运动；3～4：空腔持续减速，颗粒于3点时刻被一侧壁板抛出，颗粒团基本保持匀速运动，保持时间为t₁；该段时间内当空腔速度大于0时颗粒与空腔之间的摩擦对空腔做正功，在空腔速度小于0后做负功；4～5：颗粒团与空腔另一侧壁板发生迎面碰撞(I类碰撞)，颗粒运动速度反向达到极小值；5～6：空腔继续沿反向运动，持续与颗粒团发生碰撞(主要为I类碰撞)，颗粒团与空腔同向运动；6～7：空腔持续加速，颗粒于5点时刻被壁板抛出。颗粒团保持基本匀速运动，保持时间为t₁’。该段时间内当空腔速度小于0时颗粒与空腔之间的摩擦对空腔做正功，在空腔速度大于0后做负功。

图9给出了Γ＝188.50时监测颗粒与空腔的速度时域响应曲线。此时β_n＝0.6025，即有近一半的法向碰撞在对空腔做正功。可以看到1#、2#监测颗粒的响应曲线基本完全重合，可将颗粒团作为整体进行运动分析：

1～2：颗粒团于1点位置被一侧壁板抛出后，基本保持匀速运动，保持时间为t₂；该段时间内颗粒与空腔之间的摩擦力始终对空腔做正功；2～3：颗粒与空腔运动方向相同，但颗粒速度更快(颗粒“追赶”空腔)，并与空腔另一侧壁板发生碰撞(II类碰撞)；3：颗粒团与空腔持续作用，直至空腔与颗粒团速度相同；3～4：颗粒团同空腔一同运动；4～5：颗粒团于4点时刻被该侧壁板抛出，基本保持匀速运动，保持时间为t₂’；该段时间内颗粒与空腔之间的摩擦力使用对空腔做正功；5～6：颗粒“追赶”空腔，并与空腔壁板发生碰撞(II类碰撞)；6：颗粒团与空腔持续作用，直至空腔与颗粒团速度相同。

由于颗粒团匀速运动的开始与结束分别对应颗粒自空腔的一端移动至另一端，因此比较颗粒团保持匀速运动的时间，就可以大致了解两种激励条件下颗粒相对于空腔运动速度的差异。由图8，图9可知，t₁>t₂，t₁’>t₂’，说明Γ较高时，颗粒相对于空腔的运动速度要高于Γ较低时的情况。通常我们认为，颗粒的动能越大，颗粒的活性越强，减振的效果就越强，而上述结果恰恰说明这种说法存在局限。Γ的大小不仅能够反映空腔可为颗粒传递动能的多少，更主要的是其会影响颗粒与空腔的碰撞方式，因此对于强迫振动，颗粒与空腔的碰撞方式是随Γ的变化而改变的。

Claims (1)

1.一种基于有效碰撞耗能计算的多颗粒阻尼器减振机理分析方法，其特征在于：

步骤一：建立颗粒阻尼器数值仿真模型

明确颗粒阻尼器的各项参数，其中空腔参数包括：空腔弹性模量E_c和泊松比υ_c，空腔尺寸参数长L、宽W、高H，空腔质量M、等效刚度K及等效阻尼C；颗粒参数包括：颗粒数量N，颗粒粒径d，颗粒弹性模量E，泊松比υ，颗粒密度ρ,库伦摩擦系数μ接触阻尼常数η；

步骤二：利用离散单元法程序求解颗粒阻尼器响应

首先，确定空腔及其内部颗粒的位置；接着对所有作用在颗粒上的作用力求矢量和；然后求解颗粒i的平动及转动加速度，并按照上述过程对其余每一个颗粒的运动方程进行求解；最后求解作用在颗粒团上的合力；根据牛顿第三定律得到颗粒对空腔的作用力，将其沿水平和竖直方向分解，得到颗粒作用与空腔的作用力，并更新外部作用力及空腔运动状态，进入下一个循环；

步骤三：提取能量耗散及有效碰撞状态参数曲线；

步骤三所述提取能量耗散及有效碰撞状态参数曲线，涉及有效碰撞模型的定义和强迫振动条件下特征曲线的提取两方面内容：

A:有效碰撞模型的定义:

碰撞的分类：

V_p为颗粒的绝对速度，V_pn、V_pt分别为颗粒沿空腔法向及切向的分量，空腔绝对速度V_c也沿法向和切向分为V_cn、V_ct；以空腔为研究对象，将颗粒与空腔沿法向的碰撞作用分为两类：

I类碰撞：碰撞产生的法向作用力与空腔的法向运动方向相反，此时颗粒对空腔做负功，起到减振的作用；发生这类碰撞时，颗粒相对于与空腔有两种运动情况：第一种，颗粒与空腔在发生碰撞前的绝对速度方向相反，即“迎面碰撞”；第二种，碰撞发生前颗粒与空腔的绝对速度方向相同，但该运动方向却与颗粒相对于空腔的运动方向相反，即空腔追赶颗粒；两种运动情况均用V_cn·(V_cn-V_pn)≥0表示；

II类碰撞：碰撞产生的法向作用力与空腔的法向运动方向相同；此时颗粒对空腔做正功，无法起到减振作用；发生该类碰撞时，颗粒与空腔在碰撞发生前的绝对速度方向相同，且该运动方向与颗粒相对于空腔的运动方向相同，即颗粒追赶空腔；该种运动情况用V_cn·(V_cn-V_pn)<0表示；

同理，按照碰撞产生的切向力对空腔所做的功也能将空腔与颗粒沿切向的碰撞分为I类和II类碰撞，发生I类碰撞时摩擦力沿切向对空腔做负功，颗粒与空腔的运动情况用V_ct·(V_ct-V_pt)≥0表示，发生II类碰撞时摩擦力沿切线方向对空腔做正功，相应的运动情况用V_ct·(V_ct-V_pt)<0表示；

碰撞的能量交换计算模型：

将单个周期内颗粒对空腔所做的功W_pel按照I类、II类进行细分，进而写成公式(1)；其中W_pnel、W_pfel分别为颗粒通过法向碰撞和摩擦对空腔所做功，W_pnelh、W_pnelb分别为法向碰撞所做正功及负功，W_pfelh、W_pfelb分别为摩擦所做正功及负功；

W_pel＝W_pnel+W_pfel＝W_pnelh+W_pnelb+W_pfelh+W_pfelb (1)

通过对颗粒在单个周期内对空腔所做的功进行时域积分求解式(1)中的各项参数；假设在单个周期内有N个时间步，每个时间步内有个颗粒—空腔接触对，则颗粒对空腔做功的表达式写为

对上述能量计算结果按照结构固有阻尼耗能进行归一；定义总耗能比R_el为结构总耗能E_loss与结构固有阻尼耗能E_cdl的比值，其表达式为：

同理给出法向碰撞耗能比R_nel、摩擦耗能比R_fel的定义，见式(4)、(5)；

有效碰撞状态参数的定义：

有效碰撞状态参数包括有效法向碰撞状态系数βn及有效摩擦状态系数βf，分别用来描述单个周期内颗粒通过法向碰撞和摩擦对空腔所做的负功在总能量交换中所占的比重，两系数的表达式如下：

将式(6)、(7)代入式(3)，则总耗能比R_el又写为

B：强迫振动条件下颗粒阻尼特征曲线的提取方法：

颗粒阻尼依赖于主结构响应幅值的大小，采用-Γ曲线作为颗粒阻尼的特征曲线；/>为归一化阻尼比，表达式为

其中，ζ为安装阻尼器后的结构阻尼比，ζ0为空白对照的阻尼比；

Γ为无量纲响应加速度，表达式为

Γ＝aω²/g (10)

其中a为位移响应幅值，ω为激振频率，g为重力加速度；

当颗粒阻尼器发生强迫振动时，提取颗粒阻尼特征曲线的关键在于获得安装颗粒阻尼器后结构的主响应频率ω_d及其响应峰值A_max，这两项数据通过正弦扫频获得；当以固定位移a₀方式对空腔施加正弦扫频激励时，式(8)中激振频率ω为时变量，其表达式为

其中ω0为扫频起始频率，为扫频速率，t'为扫频时间；对颗粒空腔基础位移激励，ω为变量，力激励幅值F_ext及相位/>也均为变量，其表达式分别为：

采用正弦扫频的方式仅用一次数值仿真便准确地获得安装颗粒阻尼器后结构的主响应频率及响应峰值；此时利用幅频响应曲线得到的峰值点获得无量纲加速度Γ，利用半功率带法得到归一化的结构阻尼比为获得多个无量纲加速度Γ条件下的/>还需通过多个量级的慢速正弦扫频得到颗粒阻尼特征曲线，然后利用求得的共振频率点对阻尼器施加单频简谐激励，得到稳定状态下颗粒阻尼的耗能比曲线。