首页/文章/ 详情
banner

黏弹性:广义Maxwell模型本构理论推导

3月前浏览6135

本文摘要:(由ai生成)

本文深入探讨了二阶及以上的广义Maxwell黏弹性模型理论及其实现方法。文章介绍了高阶微分形式的广义Maxwell本构,包括一维到三维模型的推广,并提出了简化模型控制方程的方法。文章还将理论应用于UMAT列式,提供了一种更通用、易扩展的黏弹性格式。该科研笔记旨在帮助有需要的读者,同时作者欢迎读者反馈并承诺提供UMAT本构代码。

0. 上回书说到

2024年1月30日,小喵发了一篇文章:

黏弹性:从一维模型到三维UMAT本构实现

里面介绍了基础的黏弹性模型,包括Kelvin和Maxwell形式的标准线性模型(Zener模型),如何从一维推导向三维,以及如何将微分形式的本构关系与Prony级数形式的公式联系起来。

这篇文章在此基础上更进一步,讨论一下二阶及以上的广义Maxwell模型的理论,以及如何实现。

 

备注:这篇文章是小喵出于自身科研需要而总结。我自己也是学习者,一边推导一边写成了这篇推送。限于能力和水平,无法100%保证全部推导都是正确的,只能说我已尽我所能保证推导和讲解过程的可读性和准确性。

这篇文章,本想把UMAT代码都写完后再一起发出来。但鉴于后续时间安排,恐怕就又要延后一星期甚至更久。写到这里Markdown源文件已经有超过18,000个字符,考虑到篇幅因素,还是决定把公式推导部分提前单独发出来。

正文目录:

  • 0. 上回书说到

  • 1. 高阶微分形式的广义Maxwell本构

  • 2. 引入局部应变的简化形式

  • 3. 广义Maxwell黏弹性UMAT列式

  • 4. 结尾和预告

观前提醒,如果读到某个公式卡壳,千万不要死磕。第1节的公式形式复杂,到第2节会变简单。

1. 高阶微分形式的广义Maxwell本构

让我们回顾一下,从一根Maxwell黏弹性元件开始。

 

弹簧和阻尼器串联,两个元件应力相等,总应变等于二者叠加。因此Maxwell元件的控制方程可以用应变率的形式写为:

那么,对于一根刚度为     的弹簧和一根     的Maxwell元件并联组成的,Maxwell形式标准线性模型(这里把弹簧刚度的下标修改成0,是为了和后面广义Maxwell模型一致),它的弹簧和Maxwell两个元件并联,应变相等,应力为二者叠加。

 

对于弹簧元件,有:

对于Maxwell元件,和前面公式    一样:

这里我们使用一个比较通用的推导方式,方便向高阶推广:从应力的零阶导数项开始,试着让方程中不要出现    这种局部元件应力,只有整体应力    

为了实现这个目的,我们首先写出:

将公式    做简单的变形,结合公式    ,可以写出:

很显然,公式    还不够干净。它的右侧还多出来一个     项。

为了把它平衡掉,我们回头来看公式     。把公式     左右两边同时对时间求导,形式还是不变的:(求导符号我也懒得改正体了,太麻烦了,毁灭吧……摆烂.jpg)

这样就有了     。我们用它来平衡掉公式     中应力对时间的一阶导数项,就得到:

引入弛豫时间     ,简单整理一下,就大功告成:

式     就是一阶广义Maxwell黏弹性模型的控制微分方程,描述应力和应变(以及它们的时间导数)之间的关系。很简单,很好理解,对不对?甚至对于刚才我们推导的一阶广义Maxwell模型来说,仅靠观察也差不多能得到这个结果,随便怎么推都能推出公式    来。

——————

那现在让我们考虑二阶广义Maxwell模型。模型里的参数分别为     

对于弹簧元件,显然有(这次我们直接准备到二阶导数):

对于两个Maxwell元件,和公式    一样。但这次我们也给它多准备一阶的时间导数:

还是和刚才一样,首先写出     

果然,这次的公式     里面想要再直接用弹簧元件的一阶导数是不能摆平了。因为这次两个局部应力导数前面的系数不同。

那么为了伺候这两位爷,我们构建出下面的形式:

在    里面,把    那两项刨除掉,剩下的几个导数项,把    里面对时间求过一阶导的公式拿出来,简单变变形:

综上,我们把公式     和    放在一起,合并同类项:(这公式已经有点长了)

诶,现在我们可以发现,在公式    中不平衡的    两个时间导数项,在公式    里面系数变得一致了。这样,只需要在公式    中补充上    的项,即可完全消去局部应力,得到整体应力-应变之间的微分关系。

所以最终,二阶Maxwell黏弹性模型的控制方程,我们这里就可以参考之前黏弹性的推导,引入弛豫时间    

最终形成公式    。类比一下,和前面的公式    形式很像,对不对。

到这一步,我们可以大致用文字总结一下。

对于包含1个弹簧和N个Maxwell元件的广义Maxwell黏弹性模型,为了推导出其总应力-总应变之间的关系,需要列出每个局部元件的控制方程,并将每一行方程对时间求导。然后从应力对时间的零阶导数开始,逐层将局部元件的应力加总。这一过程总会残留下应力对时间的更高阶导数项。重复此过程,一直到应力对时间的第N阶导数项才可以被完全配平。

(好像用文字怎么描述都说不太清楚。不过如果你认真跟着前面的推导读下来,相信你一定可以理解这段话。)

推广到N阶,对于包含共计1个弹簧和N个Maxwell元件的广义Maxwell黏弹性模型,其控制方程的一般形式为:

——————————

前方高能。这个形式来自维 基 百 科,可难写了……

——————————

稍微展开一点,另外一种通用形式:

我已经懒得再仔细解释上面那两坨公式了。总之,按照前面的推导逻辑,N阶广义Maxwell模型的控制方程里确实需要包含应力和应变对时间的最高N阶导数。

推导结束,可喜可贺……是吗?

如果有限元软件里的广义Maxwell黏弹性模型要拿上面的控制方程去写代码,那这代码的繁琐程度恐怕就不是一般人类能看懂的了。

2. 引入局部应变的简化形式

还有没有更简洁一点的形式呢? 就,类似之前我们看到过的Prony级数,对于广义Maxwell模型的松弛模量,多一阶只需要加一项(而不是提升一阶导数):

——那肯定是可以的。

仔细想一想前面的推导,我们为什么要这么麻烦地去不断提高应力和应变对时间的求导阶次?

答案是为了得到总应力和总应变之间的关系,在方程里抹掉所有局部元件的信息。

那么,如果我们在方程里保留每一个局部的信息,会怎么样呢?

实际上,这就相当于在控制方程里引入更多额外的变量,提高方程组的维度,从而降低方程组的阶数。

来,一起来见证魔法。

仍然以最简单的一阶广义Maxwell模型为例。

 

这一次,我们在列出     的时候,不再追求把两部分应力合并在一起,就这么分开放着。本构关系是应力与应变之间的关系。那么既然要把     和     分开,也就要多定义一个独立的应变。很显然——在这个模型里,我们可以把Maxwell元件中阻尼器的应变单独写出来。

将阻尼参数为     的阻尼器应变定义为

那么   那根弹簧的应变就应该是   

这样,弹簧部分本构仍然是

Maxwell部分的本构:

这次在式    中,我们要用     和     来表示     ,而不要出现高阶导数     。这也不难,只需要把式    对时间的导数降低一阶:

就可以写出一阶广义Maxwell模型的控制方程了:

很显然,模型中的     实际对应Maxwell元件的应力卸载到零时,黏弹性模型的长期模量;而     对应的则是黏弹性模型的瞬时模量。

投共一念起,刹那天地宽。)引入黏壶应变     后,我们发现式     这样的黏弹性控制方程形式一下子简洁了许多。同理二阶广义Maxwell黏弹性模型的控制方程就可以写为:

直接就可以把N阶控制方程的通用形式都写出来了。

当然,增加了N个变量,还需要相应地补充N条方程,说清楚     和     之间的关系。这并不难。对式     做一个变形:

式     构成了完整的N阶广义Maxwell黏弹性模型的控制方程。

绕了一大圈,原来这么简单就解决了啊。那你早说呀,我前面还吭哧吭哧推那些公式干嘛呢。

3. 广义Maxwell黏弹性UMAT列式

现在,有了前面的公式     ,再来审视一下黏弹性本构模型的UMAT。让我们把它修改成更通用、更方便扩展的广义Maxwell黏弹性格式。

从式     出发,这个本构直接可以推广到三维。三维的广义Maxwell模型本构关系可以用张量记法写为:

其中,上标显然不表示乘方,主要是下标被张量标记占了只能写成上标。     是和第n个Maxwell元件关联的阻尼器内应变,    是对应的弹簧刚度。第一个     是瞬时模量,    和前面的    一样都是四阶弹性张量,具有相同的形式。每一个     都服从以下关于时间的微分方程:

根据公式     ,其实就可以直接写出相应的具体有限元列式。为了使方程形式简洁,使用     两个变量作为弹性常数。将初始(瞬时)Lamé常数记为     ,一阶广义Maxwell元件对应的刚度为     ,相应的弛豫时间为     。则式     可以用张量的形式展开为:

 

这里的符号记法与前文略有出入:前面描述广义Maxwell模型的时候,我们将并联的单个弹簧单元的刚度写成     ,下标为0以和其他Maxwell元件形成区分。但实际上我们知道那根弹簧的刚度    其实是整个模型的长期模量。

在后面的推导中,我们为了和Abaqus文档里的黏弹性公式保持一致,将初始的瞬时模量下标写为0,而将长期模量的下标写为     。例如,瞬时Lamé常数为    ,对应的长期模量数值为     (读者如果要参考本公 众  号文章来组织自己论文中的公式,请自行斟酌标记的前后逻辑一致性)

写成矩阵形式(哎呀这段LaTeX可难写了):

 

要注意,在Abaqus文档页

Abaqus > Introduction & Spatial Modeling > Conventions > Convention used for stress and strain components

里面,介绍了Abaqus以向量形式存储张量分量的顺序。上面使用的顺序是Abaqus/Standard的存储顺序。Abaqus/Explicit使用了不同的顺序:    这样。

另外,同一页文档还写明,Abaqus始终使用工程剪应变,    

所以后面的     剪切分量,我这里也自动乘了2.

同样使用中心差分列式(此处下标表示该变量所处的时刻),

式     可以写为:

化简后得:

式     的形式就可以简单地扩展到2阶乃至更高阶的Maxwell黏弹性模型: