黏弹性模型微分型本构方程的矩阵形式1)

彭培火 黄朝军

*(北京建筑大学理学院,北京 102612)

†(中国建筑第二工程局有限公司,北京 100160)

黏弹性材料广泛应用于各类工程中,如塑料、橡胶等工业材料,土壤、沥青等地质材料,肌肉、血液等生物材料,不胜枚举。由于黏弹性材料兼具弹性与黏性两种特性,其力学性能依赖于环境温度、加载时间、加载速率和幅值等条件,因此建立合适的模型,准确地描述材料的黏弹性行为,对工程应用及理论研究都具有十分重要的意义。材料黏弹性特性的描述和本构方程的表达是黏弹性力学的重要研究内容之一。为了描述黏弹性材料的应力−应变规律,常利用弹性元件和黏性元件以及它们的各种组合构成的模型来进行模拟。如Maxwell 模型(图1(a))、Kelvin 模型(图1(b))、Poynting–Thomson 模型(图1(c))、Burgers 模型(1(d))、广义Maxwell 模型、广义Kelvin 模型等等。多年来,人们开发了各种各样的模型来表征材料的黏弹性特性,对黏弹性材料的本构方程的数学描述进行了大量的研究[1]。利用这些模型,诸多学者在各个工程领域开展了大量的研究。例如,将土体视为黏弹性材料,郑灶锋等[2]利用Merchant 模型研究了稳态载荷频率和黏滞性对地基振动的影响,研究电站地下厂房洞室围岩的变形规律时,杨林德等[3]基于Kelvin–Voigt 模型来模拟和计算围岩的黏弹性变形特征。范家参[4]采用Poynting–Thomson 模型对断裂动力学问题进行了研究。Ren 等[5]利用改进的Burgers 模型,探究了黏弹性复合材料结构的振动特征。Huang 等[6]利用改进的Burgers 模型来描述沥青混凝土的黏弹性本构行为。为了更准确地表征血液的流变特性,更好地与高切变率下的实验数据相吻合,人们也常常采用由更多元件组成的黏弹性模型来进行相关研究[7]。候之超等[8]利用分数微分方程描述泡沫材料的黏弹性特征,研究了单向恒速率加载下模型的力学响应。李兆生等[9]利用广义Maxwell 模型构建损伤力学本构模型,研究了沥青材料的冻融损伤过程。詹小丽等[10]采用分数阶导数Maxwell 模型对沥青的动态黏弹性力学性能进行了研究。潘文潇等[11]研究了两平板之间的广义Maxwell 黏弹性流体非定常流动的解析解。不同的模型具有不同的优缺点,如Maxwell 模型可以模拟应力松弛现象,但不能表示蠕变过程；Kelvin模型能够模拟蠕变,却不能表示应力松弛[12-13]。一般来说,随着弹簧和阻尼器数量的增多,模型能够模拟的黏弹性特征也会越来越贴近实际材料[14]。但是,对于多个弹簧和阻尼器以任意连接方式组成的黏弹性模型,其微分型本构方程如何快速计算出,是各种基于黏弹性材料结构计算的基础。图论算法可方便快捷地计算路径等问题,而有向图的路径和闭包围又与黏弹性模型的应力方程和应变方程存在联系,以计算由10 个弹性元件和黏性元件构成的较为复杂的模型的微分本构方程为例,推导得到了任意线性黏弹性模型的微分型本构方程的一般矩阵形式,总结了一套适合计算机编程的固定范式。并利用Python语言编写了图形界面用户程序,可快速地构建任意复杂的线性黏弹性模型和计算其应力−应变微分关系式。该计算方法对准确地构建实际工程材料的黏弹性计算模型、精确地拟合实际工程材料的黏弹性特性及应力−应变规律具有一定的参考意义。

1 本构方程的矩阵形式表述

一般的黏弹性力学教材中都给出了各种基本模型的微分型本构方程[14],如：Maxwell 模型、Kelvin模型、Poynting–Thomson 模型、Burgers 模型等等。由于弹性元件和黏性元件数量不多,它们的本构方程推导过程并不难。为了实现计算机编程计算本构方程,以下利用图论算法和矩阵运算,通过推导由10 个弹性元件和黏性元件构成的较为复杂模型的应力−应变微分关系,总结和建立任意线性黏弹性模型的微分型本构方程的一般形式。

如图2(a) 所示的10 元件黏弹性模型,可以看成由Maxwell 单元连接组成,如图2(b) 中所示。单个Maxwell 单元的应力−应变微分关系为

图2 Fig.2

基于公式(1),只要写出10 元件黏弹性模型的应变ε、模型的应力σ与各个Maxwell 单元的应变εi(i= 1,2,··· ,5),应力σi(i= 1,2,··· ,5) 之间的关系式,将它们适当求导,并写成矩阵形式,模型的微分型本构方程即可以通过各系数矩阵的运算得到。为了计算便捷和程式化,借助图论中有向图的相关工具,可迅速将求解过程编写成程序。现将整个求解过程描述如下。

(1)根据黏弹性模型中的弹簧和阻尼器的连接方式,写出模型的总应变ε与各个Maxwell 单元的应变εi之间的关系式,并求j阶导数。将模型的总应变ε及其导数与各Maxwell 单元的应变εi及其导数的关系式以矩阵形式记为

式中,记Djε表示由模型的总应变ε及其1～j阶导数组成的列向量,有j+ 1 个元素,形如式(3)所示

记Djεi表示各个Maxwell 单元的应变εi及其1～j阶导数组成的列向量,有i×(j+1) 个元素,i为Maxwell 单元的数量,形如式(4) 所示

系数矩阵A有(j+1)×r行,i×(j+1) 列,为行满秩矩阵,即R(A) = (j+1)r。Ir,j+1是形如式 (5) 所示的矩阵。因模型的总应变ε与各个Maxwell 单元的应变εi之间的关系式共有r个独立方程,故乘以该Ir,j+1矩阵来将其合写成一个矩阵方程。

式中,dr为共有r个元素且值全为1 的向量。将模型转化为有向图,如图3 所示。模型的总应变ε与各个基础元件的应变εi之间的r个独立应变关系式可以根据由Maxwell 单元组成的模型中由起点到终点的独立路径的数量进行确定。有向图中独立路径的数量即为r,如图3(a) 中所示,r= 3。可以写出3 个总应变ε与各个Maxwell 单元的应变εi之间的关系式,即

图3 Fig.3

根据图论中求解有向图路径的方法,先写出图的关联矩阵为

列出方程Mx=b,然后求其解(x的元素只能是1 或者0,向量b表示从起点v1到终点v4)。

方程(10) 的满足条件的解中,有3 个线性无关的解向量为

以式(11) 中的每一个解表示一条从起点v1到终点v4的路径,解的第几个元素为1 表示该路径经过第几条有向边。如：向量x2的第1,4,5 行元素为1,其余元素为0,表示图3(a)中路径 2○是一条经过e1,e4,e5从起点v1到终点v4的路径。由x1,x2,x3构成的矩阵即为独立路径矩阵,记为R。

对比式(2)～式(8),式(12),可以得出式(2)中的矩阵A可由独立路径矩阵R构造而成。

(2)根据黏弹性模型中的弹簧和阻尼器的连接方式,写出模型的总应力σ与各个Maxwell 单元的应力σi之间的关系式,并求m阶导数,将模型的总应力σ及其导数与各Maxwell 单元的应力σi及其导数的关系式以矩阵形式记为

式中,模型的总应力σ与各个基础元件的应力σi之间的关系式共有n个独立方程,利用乘以In,m+1来将其合写成一个矩阵方程。记Dmσ表示由模型的总应力σ及其1～m阶导数组成的列向量,有m+1个元素,形如式(15) 所示。

记Dmσi表示各个Maxwell 单元的应力σi及其1～m阶导数组成的列向量,有i×(m+1)个元素,i为Maxwell 单元的数量,形如式(16) 所示。

系数矩阵B有(m+1)×n行,i×(m+1) 列,为行满秩矩阵,即R(B)=(m+1)n。模型的总应力σ与各个基础元件的应力σi之间的n个独立关系式可以根据有向图中包含起点v1的闭包围的数量进行确定(其中闭包围须满足条件：与闭包围相交的有向边的方向应指向外)。如图3(b)中所示,n=3。可以写出3 个总应力σ与各个Maxwell 单元的应力σi之间的应力关系式,即

因为有向图的任一闭包围必与每一条路径相交,故闭包围矩阵的任一行与独立路径矩阵R的每一行求点积必为1,于是闭包围矩阵可通过求方程Rx=c的解获得(非齐次项c为全1 向量,其行数与R的行数相同,解x的元素只能是1 或者0)。

方程(20) 的满足条件的解中,有3 个线性无关的解向量为

以上式(21)中的每一个解表示有向图中一条包含起点v1的闭包围,解的第几个元素为1 表示该闭包围与第几条有向边相交。如：向量x2的第2,3,4 行元素为1,表示图3(b)中闭包围 2○与有向边e2,e3,e4相交。由以上解x1,x2,x3构成的矩阵即为闭包围矩阵,记为E。

对比式(14)～式(19),式(22),可以得出式(14)中的矩阵B可由闭包围矩阵E构造而成。

(3) 将Maxwell 单元的应力−应变微分关系式(1) 求适当阶数的导数后代入式(2),可得到模型的总应变ε与各个Maxwell 单元的应力σi之间的关系式,以矩阵形式记为

式中,记Djε′表示对向量Djε再求一阶导数,故Djε′是由模型的总应变ε的1～j+1 阶导数组成的列向量,有j+1 个元素,形如式(25) 所示。

由黏弹性模型中各弹簧的弹性模量Ei及各阻尼器的黏度系数ηi构造Y矩阵和V矩阵分别为

不难得出,式(24) 中的系数矩阵H为

(4) 联立求解方程(14) 与(24),消去中间变量Dmσi,即可求得模型的应力σ与应变ε之间的微分关系式(需满足m=j+1,保证Dmσi=Dj+1σi)。具体求解过程如下。

将式(14) 看成是非齐次线性方程组,将向量Dmσi看成待求量,In,m+1Dmσ为非齐次量,则Dmσi可由Dmσ表示为

式中,null(B) 表示矩阵B的零空间矩阵,P为式(14) 的任意一特解,可按式(30) 方法构造

其中矩阵F为

式中,列向量f可取独立路径矩阵R的任一行通过转置得到,F矩阵有i×(m+1) 行,m+1 列。将式(29) 代入式(24),可得

记

求出CT的零空间矩阵 null(CT)。 式 (33)中,null(B) 为i ×(m+ 1) 行,(i −n)×(m+ 1)列矩阵,H为r ×(j+ 1) 行,i ×(j+ 2) 列矩阵(因为m=j+ 1,H ×null(B) 符合矩阵乘法规则)。C为r ×(j+1) 行,(i −n)×(m+1) 列矩阵,故CT的零空间矩阵null(CT) 有r ×(j+ 1)行,r×(j+1)−(i −n)×(m+1) 列。记null(CT)的某一个列向量为b,将其转置并左乘式(32),则可以得到黏弹性模型的应力σ与应变ε之间的微分关系式的矩阵形式表述为

(5)上述计算过程中,关于r个应变方程及n个应力方程求多少阶导数最合适的问题,即m和j的取值为多少合适。从上述的应变方程和应力方程以及推导过程可知,应力的求导阶数比应变应该高一阶,故应有m=j+1。为了使CT的零空间矩阵null(CT) 至少含有一个向量b,则应有

式中,i为黏弹性模型中Maxwell 单元的数量。

(6) 为了将上述计算方法程式化, 需要将黏弹性模型中的所有弹性元件和黏性元件都配对成Maxwell 基本单元。当由弹簧和阻尼器构成的连接图中存在有非Maxwell 基本单元时,先把非Maxwell单元的部分扩展成Maxwell 单元。例如,图1(c) 所示的Poynting–Thomson 模型,可以通过添加相应的弹簧和阻尼器进行扩展成如图4(a) 所示,则该模型中的非Maxwell 单元的部分都扩展成了Maxwell 单元,如图4(b) 所示。先计算出图4 所示的黏弹性模型的微分本构方程,然后令η2为∞,即与原模型等价。

图4 Fig.4

2 Python 计算程序示例

基于以上推导的计算方法,利用Python 编写了图形用户界面程序(如图5 所示),可以实现自由搭建模型和计算本构方程等功能。如将软件界面右侧的弹簧或阻尼器拖入左侧图框中,即可构造以任意的连接方式组成的线性黏弹性模型。程序在鼠标拖拽的同时会自动识别各个弹性元件或黏性元件的连接结构,若自动识别的连接结构不符合要求,还可以通过<编辑>菜单进行调整。黏弹性模型构建完毕后即可转换得到与之相对应的有向图,并得到关联矩阵M。在<计算>菜单中可设置相关参数并进行计算,如采用符号计算或是数值计算,零空间矩阵的求解方法等等。程序按照图6 所示的流程图可快速计算出黏弹性模型的微分形式的本构方程。并在此基础上,由本构方程计算松弛模量、蠕变柔度、复模量等描述黏弹性材料的基本物理量。计算结果输出在程序界面下方的图框内,也可输出中间计算结果(如独立路径矩阵R、闭包围矩阵E等等)。通过多个模型的测试,利用该程序均能快速正确地计算出黏弹性模型的本构方程。其中表1 为2 个简单黏弹性模型的算例,依据图6 所示的计算流程,给出了中间过程及最终结果。程序还可继续扩展,开发具有其他功能的代码。

图5 程序界面及10 元件模型的计算结果Fig.5 Program interface and results of 10-component model

图6 线性黏弹性模型本构方程计算流程图Fig.6 Calculation flow chart of constitutive equation of linear viscoelastic model

表1 中右侧所示模型,若令η2和E3为无穷大,表示阻尼η2处和弹簧E3处为刚性连接,即与图1(d)所示的Burgers 模型等价。按照式(34) 计算得到的本构方程中,η2和E3都出现在分母的位置,令η2和E3为无穷大,则倒数为零,可以很方便计算得到Burgers 模型的本构方程为

表1 两个简单模型的计算示例(续)Table 1 Calculation examples of two simple models (continued)

3 结语

诸多工程材料如塑料、橡胶、沥青等具有随时间变化的力学响应,如恒载荷下的蠕变、恒定变形下的应力松弛和卸载时的延迟应变恢复等现象。材料的这些黏弹性特征对它们承载时的力学性能有着重要的影响。为了从数学上对黏弹性材料的应力−应变规律进行描述,常利用弹性元件和黏性元件以及它们的组合来建立模型,模拟材料的黏弹性行为。一般来说,弹簧和阻尼器数量越多,模型能够模拟的黏弹性特征也会越来越贴近实际材料。但是,元件越多,模型的本构方程的计算也越复杂。对于多个弹簧和阻尼器以任意连接方式组成的黏弹性模型,本文推导出了其微分型本构方程的一般矩阵形式的计算方法。首先通过添加弹性参数或黏性参数为无穷大的虚加元件,将存在非Maxwell 基本单元的模型扩展为由Maxwell 基本单元构成的标准模型。然后将标准模型转化为与之对应的有向图,根据图论的有关知识,基于有向图的特征,由独立路径的形式可以写出基本应变方程,由闭包围的形式可以写出基本应力方程,并根据式(35) 对其求适当阶数的导数,按照图6 所示的计算流程可以得到独立路径矩阵、闭包围矩阵及其他中间结果,通过各系数矩阵的运算,即可根据式(34)直接得到模型的微分型本构方程的矩阵形式表述。本文总结了一套适合计算机编程的固定范式,编写了Python 程序,可方便搭建黏弹性模型及快速计算其本构关系,并在此基础上可继续开发具有其他功能的代码。