金属橡胶的本构关系研究★

王亚杰 赵亚哥白* 任小龙 黄新宇 王 璇 韩静怡

(东北林业大学，黑龙江哈尔滨 150040)

过去的许多试验表明，金属橡胶材料是一种具有良好的变形自恢复能力的非线性干摩擦阻尼材料。它具有含有良好三次非线性成分［1］的无记忆恢复力，除此还具有优秀的变形自恢复能力即记忆恢复力，金属橡胶耗能器的恢复能力具有非线性迟滞效应，其阻尼成分既有粘性阻尼又有干摩擦阻尼。国内外现在有宏观力学模型与微观力学模型两种建立金属橡胶本构关系比较著名的方法，在建立本构关系时需考虑金属橡胶材料诸多复杂性质。

1 建立宏观力学模型

通常采用理论法和实验法［2］两种方法来建立宏观数学本构模型，理论法常常被看做“白箱”问题，通过运用现有的定理规律，来推导判断系统中状态参数与作用的关系。而实验法是直接从实验数据出发，通过归纳总结来建立数学模型，这种情况又被称为“黑箱”问题。但是在实际中遇到的问题往往是介于两者之间的问题，必须将定理分析与实验归纳融合起来来寻找解决方法，遇到这种情况可以根据部分已知特性并通过特殊处理后来推导未知参数从而确定完整的状态方程，这种问题常常被称之为“灰箱”问题。目前常采用最小二乘法来确定未知参数，相比于极大似然法，辅助变量法等其他方法，最小二乘法具有概念简单，使用情况更加广泛等很多优点，用其求得的估计值具有最优统计特性，可以实现使得实验数据在平方误差意义上达到最小误差的效果。其实现过程为:假定一组变量含有m个元素x=(x1，x2…，xm)与某个变量z成某种确定的线性关系，如下:

方程中:θ=(θ1，θ2，…，θm)为一组参数并且是常数，式中 θi是未知的，该未知参数的取值可以通过在不同时刻观测变量z和x的关系来估计。假设已经在t1，t2，…，tn时刻对二者作了n次观测，其观测数据用 z(i)和 x1(i)，x2(i)，…，xm(i)，(i=1，2，…，n)两组参数来代表，如此便可以得到以下方程:

上述方程又可定义为回归函数，θm则为回归系数。其又可写成矩阵形式如下所示:

式中:

为了估计参数θ，前提是m≥n。如果m与n值相等，那么θ值则有唯一解，即是:

再设 J=(Z －Xθ)T(Z － Xθ)=ZTZ － θTXTZ － ZTXθ+θTXTXθ。

从而解得θ～=(XTX)－1XTZθ即为在最小二乘法基础上的估计值。

白鸿柏、黄协清［1］通过运用宏观力学分析方法，并使用最小二乘法对各参数进行识别取值来确定金属橡胶隔振器中的力—位移之三次非线性迟滞泛函本构关系，又得到三阶非线性多项式的弹性力和等效粘性阻尼力叠加组成了金属橡胶的恢复力的结论。姜洪源［3］曾提出建立金属橡胶构件的变形模型可以通过实验的理论方法来实现。

2 细观力学模型的研究

由于承受径向和轴向荷载的微螺旋弹簧冲压使得金属橡胶材料的细观结构特征成型，则金属橡胶的刚度由微螺旋弹簧的刚度决定。在分析中可忽略弹簧丝轴力和剪力的影响，只考虑弯矩影响。陈艳秋曾提出用小曲梁的串并连来描述材料总体刚度，可将杂乱金属丝看成多层微小弹簧，每一个弹簧都有一定刚度，它们共同组成金属橡胶总的刚度。彭威，白鸿柏［5］也曾通过引用新的材料参数铺层比例系数来建立微弹簧组合变形的细观本构模型。

3 建立硬化折线本构模型

为建立金属橡胶本构模型，可以运用实验法，该方法的简化分析模型不仅便于计算，简化了MR耗能减振结构的地震反应分析方法，而且能够表征材料的应变硬化特征，与试验曲线符合。根据实验法原理，在室温下对相对密度为0.27，加载频率为1 Hz，并处于最大应变幅值为20%的应力—应变曲线部分的金属橡胶材料进行拟合，其相应的应力应变曲线如图1所示。

构建模型可以利用MATLAB计算软件来实现:第一步将MR应力—应变曲线分成两部分，分别为加载段和卸载段;然后对前者采用局部依次拟合的方法，来得到图中的OA与AB曲线段，从而得到简化模型的强化应力、初始弹性模量等重要参数，最后确定BC段和CO段依据原则为能量相等，使简化模型与实际曲线包围相等的面积(见图2)。

图1 OMR-D“τ — γ”迟滞回线

图2 应变硬化折线本构模型拟合

下面将以此为例来编制计算程序，其详细步骤如下所述:

1)把名义应力、名义应变根据相应转换关系转换为真实应力与应变，将数据拷入txt文件中，以便MATLAB读取。

2)将原点和应力最大值作为拟合折线的两个点，近似读取中间加载点的坐标，写入系统，bx，by为选定刚度变化拐点坐标，cx，cy为最大横纵坐标。

3)输出面积差矩阵。

竖方向223个点表示应变点，i=1∶223，表示曲线最大横坐标值乘以10，纵坐标j=1∶20，表示有20个点。

4)计算步骤，假设已经确定了(15，04)为加载段拐点的坐标，并确定了(22，3)为曲线最大的点的坐标，这时卸载段的横坐标应该处于这两个坐标之间，若此时取x=19为卸载拐点横坐标，那么190那一行就是在矩阵里面对应的，找到这行最小面积差的那个点0.233 7，j=8，将其代入 for j=1∶20，y=(j/2 －1) × ( －1) ×0.1 进行处理计算，便可以得到y，此时，便确定了卸载点的纵坐标。

注意事项:

1)根据面积最小原理计算，若有两个符合要求的卸载点，则需选择横坐标位置符合真实规律的点。

2)卸载点对应的横坐标最好取中间的点(见图3)。

图3 本构拟合面积相等卸载点

最后根据上述数值可以调整ABAQUS关于金属橡胶的子程序参数，从而进一步建立本构模型。

4 结语

通过讨论细观力学模型并引入宏观力学模型，采用最小二乘法为参数估计的方法，研究变量y与若干元素的线性关系得到回归函数，定义误差矢量来求得回归系数的最小二乘的大致值。来建立金属橡胶构件的变形模型，在室温下相对密度为0.27，加载频率为1 Hz，并处于最大应变幅值为20%的应力—应变曲线部分的金属橡胶材料进行拟合，对实验结果按照包络面积相等原则，通过MATLAB计算软件来处理，从而得到应变硬化折线本构模型。下一步要研究如何建立金属橡胶阻尼器的剪力墙—连续梁结构分析模型，调整ABAQUS的子程序参数，使金属橡胶本构模型更加具体化，提高其可靠性。