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

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

## THE MATRIX FORM OF DIFFERENTIAL CONSTITUTIVE EQUATION FOR VISCOELASTIC MODEL1)

PENG Peihuo,*,2), HUANG Chaojun

*School of Science, Beijing University of Civil Engineering ＆ Architecture, Beijing 102612, China

China Construction Second Engineering Bureau Co., Ltd, Beijing 100160, China

 基金资助: 1)北京市属高校基本科研业务费专项资金项目(X20049)

Received: 2021-07-7   Revised: 2021-09-10

Abstract

This paper derived the general matrix form for the linear viscoelastic model for arbitrarily linked springs and dampers. Firstly, the underlying linear viscoelastic problem is modeled through the standard model composed of basic Maxwell units. Then the standard model is transformed into a directed graph, and the basic stress equation and strain equation are expressed in the form of independent path and closed enclosure. And then the general matrix form of a differential constitutive equation for an arbitrary viscoelastic model is derived. Finally, a universal paradigm suitable for computer programming is developed, and the relevant algorithm is realized with the Python language along with numerical results.

Keywords： viscoelasticity; mathematical model; differential form; constitutive equation; matrix form

PENG Peihuo, HUANG Chaojun. THE MATRIX FORM OF DIFFERENTIAL CONSTITUTIVE EQUATION FOR VISCOELASTIC MODEL1). Mechanics in Engineering, 2022, 44(2): 358-367 DOI:10.6052/1000-0879-21-267

Fig.1

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

$\dfrac{\dot{{\sigma }}}{E}+\dfrac{\sigma }{\eta }=\dot{{\varepsilon }}$

### 图2

Fig.2

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

$I_{r,j+1} {D}^{j}\varepsilon ={A}\times {D}^{j}\varepsilon_{i}$

${D}^{j}\varepsilon =\left[\varepsilon {\dot{{\varepsilon }}} \cdots {\dfrac{{\rm d}^{j}\varepsilon }{{\rm d}t^{j}}}\right]^{\rm T}$

${D}^{j}\varepsilon_{i} =\\\quad \bigg[ {\varepsilon_{1} } {\varepsilon_{2} } \cdots {\varepsilon_{i} } {\dot{{\varepsilon }}_{1} } {\dot{{\varepsilon }}_{2} }\cdots {\dot{{\varepsilon }}_{i} } \\\quad {\dfrac{{\rm d}^{j}\varepsilon_{1} }{{\rm d}t^{j}}} {\dfrac{{\rm d}^{j}\varepsilon_{2} }{{\rm d}t^{j}}} \cdots {\dfrac{{\rm d}^{j}\varepsilon_{i} }{{\rm d}t^{j}}}\bigg]^{\rm T}\\\quad$

$\boldsymbol{I}_{r, j+1}=\underbrace{\left[\begin{array}{cccc}\boldsymbol{d}_{r} & \mathbf{0} & \cdots & \mathbf{0} \\\mathbf{0} & \boldsymbol{d}_{r} & \cdots & \mathbf{0} \\\cdots & \cdots & \cdots & \cdots \\\mathbf{0} & \mathbf{0} & \cdots & \boldsymbol{d}_{r}\end{array}\right]}_{j+1}\}_{j+1}^{[}$

$\varepsilon =\varepsilon_{1} +\varepsilon_{3}\\$
$\varepsilon =\varepsilon_{1} +\varepsilon_{4} +\varepsilon_{5}\\$
$\varepsilon =\varepsilon_{2} +\varepsilon_{5}$

### 图3

Fig.3

$\boldsymbol{M}=\left[\begin{array}{ccccc}e_{1} & e_{2} & e_{3} & e_{4} & e_{5} \\1 & 1 & 0 & 0 & 0 \\-1 & 0 & 1 & 1 & 0 \\0 & -1 & 0 & -1 & 1 \\0 & 0 & -1 & 0 & -1\end{array}\right] \begin{array}{l}v_{1} \\v_{2} \\v_{3} \\v_{4}\end{array}$

$\left[\begin{array}{ccccc}1 & 1 & 0 & 0 & 0 \\-1 & 0 & 1 & 1 & 0 \\0 & -1 & 0 & -1 & 1 \\0 & 0 & -1 & 0 & -1\end{array}\right] \cdot \boldsymbol{x}=\left[\begin{array}{c}1 \\0 \\0 \\-1\end{array}\right] \begin{array}{l}v_{1} \\v_{2} \\v_{3} \\v_{4}\end{array}$

$\boldsymbol{x}_{1}=\left[\begin{array}{l}1 \\0 \\1 \\0 \\0\end{array}\right], \quad x_{2}=\left[\begin{array}{l}1 \\0 \\0 \\1 \\1\end{array}\right], \quad x_{3}=\left[\begin{array}{l}0 \\1 \\0 \\0 \\1\end{array}\right]$

$\boldsymbol{R}=\left[\begin{array}{lll}\boldsymbol{x}_{1} & \boldsymbol{x}_{2} & \boldsymbol{x}_{3}\end{array}\right]^{\mathrm{T}}=\left[\begin{array}{ccccc}1 & 0 & 1 & 0 & 0 \\1 & 0 & 0 & 1 & 1 \\0 & 1 & 0 & 0 & 1\end{array}\right]$

$A=\underbrace{\left[\begin{array}{cccc}R & 0 & \cdots & 0 \\0 & R & \cdots & 0 \\\cdots & \cdots & \cdots & \cdots \\0 & 0 & \cdots & R\end{array}\right]}_{j+1}\} j+1$

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

${I}_{n,m+1} {D}^{m}\sigma ={B}\times {D}^{m}\sigma_{i}$

$D^{m} \sigma=\left[\begin{array}{cccc}\sigma & \dot{\sigma} & \cdots & \frac{\mathrm{d}^{m} \sigma}{\mathrm{d} t^{m}}\end{array}\right]^{\mathrm{T}}$

$\begin{array}{l}\boldsymbol{D}^{m} \sigma_{i}= \\{\left[\begin{array}{llllllll}\sigma_{1} \sigma_{2} \cdots \sigma_{i} \dot{\sigma}_{1} \dot{\sigma}_{2} \cdots \dot{\sigma}_{i} \\\frac{\mathrm{d}^{m} \sigma_{1}}{\mathrm{~d} t^{m}} \frac{\mathrm{d}^{m} \sigma_{2}}{\mathrm{~d} t^{m}} \cdots \frac{\mathrm{d}^{m} \sigma_{i}}{\mathrm{~d} t^{m}}\end{array}\right]^{\mathrm{T}}}\end{array}$

$\sigma =\sigma_{1} +\sigma_{2}$
$\sigma =\sigma_{2} +\sigma_{3} +\sigma_{4}$
$\sigma =\sigma_{3} +\sigma_{5}$

$\left[\begin{array}{lllll}1 0 1 0 0 \\1 0 0 1 1 \\0 1 0 0 1\end{array}\right] \cdot \boldsymbol{x}=\left[\begin{array}{l}1 \\1 \\1\end{array}\right]$

$\boldsymbol{x}_{1}=\left[\begin{array}{l}1 \\1 \\0 \\0 \\0\end{array}\right], \quad \boldsymbol{x}_{2}=\left[\begin{array}{l}0 \\1 \\1 \\1 \\0\end{array}\right], \quad \boldsymbol{x}_{3}=\left[\begin{array}{l}0 \\0 \\1 \\0 \\1\end{array}\right]$

$\boldsymbol{E}=\left[\begin{array}{lll}\boldsymbol{x}_{1} \boldsymbol{x}_{2} \boldsymbol{x}_{3}\end{array}\right]^{\mathrm{T}}=\left[\begin{array}{lllll}1 1 0 0 0 \\0 1 1 1 0 \\0 0 1 0 1\end{array}\right]$

$B=\underbrace{\left[\begin{array}{cccc}\boldsymbol{E} \mathbf{0} \cdots 0 \\\mathbf{0} \boldsymbol{E} \cdots 0 \\\cdots \cdots \cdots \cdots \\0 0 \cdots \boldsymbol{E}\end{array}\right]}_{m+1}\} m+1$

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

${I}_{r,j+1} {D}^{j}{\varepsilon }'={H}\times {D}^{j+1}\sigma_{i}$

$\boldsymbol{D}^{j} \varepsilon^{\prime}=\left[\begin{array}{cccc}\dot{\varepsilon} & \ddot{\varepsilon} & \cdots & \frac{\mathrm{d}^{j+1} \varepsilon}{\mathrm{d} t^{j+1}}\end{array}\right]^{\mathrm{T}}$

$\boldsymbol{Y}=\underbrace{\left[\begin{array}{cccc}E_{1} & 0 & \cdots & 0 \\0 & E_{2} & \cdots & 0 \\\cdots & \cdots & \cdots & \cdots \\0 & 0 & \cdots & E_{i}\end{array}\right]}_{i}\} i$
$\boldsymbol{V}=\underbrace{\left[\begin{array}{cccc}\eta_{1} & 0 & \cdots & 0 \\0 & \eta_{2} & \cdots & 0 \\\cdots & \cdots & \cdots & \cdots \\0 & 0 & \cdots & \eta_{i}\end{array}\right]}_{i}\} i$

$\boldsymbol{H}=\underbrace{\left[\begin{array}{ccccc}\boldsymbol{R} \boldsymbol{V}^{-1} \boldsymbol{R} \boldsymbol{Y}^{-1} \ldots 0 0 \\0 \boldsymbol{R} \boldsymbol{V}^{-1} \cdots 0 0 \\\cdots \cdots \cdots \cdots \cdots \\0 \mathbf{0} \cdots \boldsymbol{R} \boldsymbol{V}^{-1} \boldsymbol{R} \boldsymbol{Y}^{-1}\end{array}\right]}_{j+2}\} j+1$

(4) 联立求解方程(14)与(24),消去中间变量${D}^{m}\sigma_{i}$,即可求得模型的应力$\sigma$与应变$\varepsilon$之间的微分关系式(需满足$m=j+1$,保证${D}^{m}\sigma_{i}={D}^{j+1}\sigma_{i} )$。具体求解过程如下。

${D}^{m}\sigma_{i} ={\rm null}({B})\times {k}+{P}$

${P}={F}\times {D}^{m}\sigma$

$\boldsymbol{F}=\underbrace{\left[\begin{array}{cccc}\boldsymbol{f} 0 \cdots 0 \\0 \boldsymbol{f} \cdots 0 \\\cdots \cdots \cdots \cdots \\0 0 \cdots \boldsymbol{f}\end{array}\right]}_{m+1}\} m+1$

${I}_{r,j+1} {D}^{j}{\varepsilon }'={H}\times {\rm null}({B})\times {k}+{H}\times {P}$

${C}={H}\times {\rm null}({B})$

${b}^{\rm T}\times {I}_{r,j+1} \times {D}^{j}{\varepsilon }'={b}^{\rm T}\times {H}\times {F}\times {D}^{m}\sigma$

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

$(j+1)r-(i-n)(m+1)\geqslant 1$

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

Fig.4

## 2 Python计算程序示例

### 图5

Fig.5   Program interface and results of 10-component model

### 图6

Fig.6   Calculation flow chart of constitutive equation of linear viscoelastic model

$\dfrac{1}{\eta_{3} }{\varepsilon }'+\dfrac{1}{E_{2} }{\varepsilon }"=\dfrac{1}{\eta_{3} \eta_{1} }\sigma +\\\qquad \lt(\dfrac{1}{E_{2} \eta_{3} }+\dfrac{1}{\eta_{1} E_{2} }+\dfrac{1}{\eta_{3} E_{1} }){\sigma }'+\dfrac{1}{E_{1} E_{2} }{\sigma }"\qquad$

Table 1  Calculation examples of two simple models

## 参考文献 原文顺序 文献年度倒序 文中引用次数倒序 被引期刊影响因子

Wineman A.

Nonlinear viscoelastic solids-a review

Mathematics and Mechanics of Solids, 2009, 14(3):300-366

Zheng Zaofeng, Cai Yuanqiang, Xu Changjie.

Axisymmetric dynamic response on layered saturated viscoelsatic soils under steady-state loads

Chinese Journal of Rock Mechanics and Engineering, 2005, 24(13):2380-2385 (in Chinese)

Yang Linde, Yan Jianping, Wang Yuezhao, et al.

Study on time-dependent properties and deformation prediction of surrounding rock

Chinese Journal of Rock Mechanics and Engineering, 2005, 24(2):212-216 (in Chinese)

Fan Jiashen.

Analytical solution for mode III dynamic rupture of standard linear viscoelastic solid with nonlinear damping

Applied Mathematics and Mechanics, 2000, 21(4):415-423 (in Chinese)

Ren S, Zhao G, Zhang S.

Elastic-viscoelastic composite sructures analysis with an improved burgers model

Journal of Vibration and Acoustics, 2018, 140(3):1-10

Huang C, Wang F, Gao T, et al.

A new viscoelastic mechanics model for the creep behaviour of fibre reinforced asphalt concrete

Frattura ed Integrità Strutturale, 2018, 12(45):108-120

Yin Zhiyong, Cai Shaoxi, Wang Zhengguo.

The 5-element viscoelasticity model of the blood and the physiological signification

Journal of Chongqing University $($Natural Sciences$)$, 1999(2):3-5 (in Chinese)

Hou Zhichao, Zhao Pei, Gao Jianghua.

Constitutive model with fractional derivatives for polyurethane foams under uniaxial compression

Sciencepaper Online, 2009, 4(10):765-770 (in Chinese)

Li Zhaosheng, Tan Yiqiu.

Research on asphalt mixture freezing-thawing damage performance

China Sciencepaper, 2014, 9(11):1279-1281 (in Chinese)

Zhan Xiaoli, Zhang Xiaoning, Wang Duanyi, et al.

Study on nonlinear viscoelastic constitutive equation of modified asphalt and its applications

Engineering Mechanics, 2009, 26(4):187-191 (in Chinese)

Pan Wenxiao, Tan Wenchang.

An unsteady flow of a viscoelastic fluid with the fractional Maxwell model between two parallel plates

Mechanics in Engineering, 2003, 25(1):19-22 (in Chinese)

Yang Tingqing. Theory of Viscoelastility. Wuhan: Huazhong University of Science and Technology Press, 1990 (in Chinese)

Zeng Xiangguo, Chen Huayan, Chen Jun. Theory of Viscoelastility. Chengdu: Sichuan University Press, 2016 (in Chinese)

Cai E. The Fundament of Viscoelastic Mechanics. Beijing: Beijing University of Aeronautics and Astronautics Press, 1989 (in Chinese)

/

 〈 〉