,修晨曦,常江芳 武汉大学 土木建筑工程学院， 湖北 武汉 430072
收稿日期: 2014-09-22
基金项目: 国家自然科学基金资助项目(10802060,11172216,11472196);国家重点基础研究发展计划项目(2010CB731502); 湖北省自然科学基金资助项目(2013CFB287).
作者简介: 楚锡华(1977-),男,河南濮阳人,武汉大学副教授,博士生导师.
摘要: 在Volokh软化模型的基础上，考虑颗粒材料的变形和破坏机制，认为颗粒材料以剪切变形和破坏为主，提出了表示软化程度的软化因子概念，发展了两种能够部分表征颗粒材料力学行为特征的软化模式，并基于ABAQUS用户材料子程序(UMAT)开发相应程序.数值算例考察了上述模型模拟颗粒材料应变局部化现象和软化现象的能力.结果表明了所建议模型的有效性和可行性.
关键词: 颗粒材料 超弹性软化 应变软化 软化因子 应变局部化
Models of Hyperelasticity with Softening and Simulation of Strain Localization for Brittle Granular Materials
CHU Xi-hua
, XIU Chen-xi ,CHANG Jiang-fang School of Civil Engineering, Wuhan University, Wuhan 430072, China.
Corresponding author: CHU Xi-hua, associate professor, E-mail: Chuxh@whu.edu.cn
Abstract: Two modified softening models are developed based on the Volokh softening model to consider the deformation and failure of granular material which is governed by shear deformation and failure. Different treatments to the volume strain energy and the shear strain energy are used in the two modified models, and a concept of the softening factor to indicate the degree of softening is proposed. Based on user-defined material (UMAT) of ABAQUS, the simulation is implemented. Numerical examples illustrated the capability and performance of the presented models in modeling strain localization and strain softening of granular materials.
Key words: granular materials hyperelasticity with softening strain softening softening factor strain localization 颗粒材料在宏观上具有非常复杂的力学行为，如非线性、应变软化等.为描述颗粒材料的力学行为，基于连续介质模型已发展了大量弹塑性模型.这些模型需至少包含屈服面、流动法则和硬化(软化)法则3个基本要素.然而上述3个基本要素的确定通常具有一定难度.除弹塑性模型外，近年来，超弹性模型亦被用来描述岩土颗粒材料的力学行为[1, 2, 3].如Jiang和Liu基于热力学失稳和弹性模量的压力相关性推导了一个能够计算屈服和剪胀性的GE模型[2, 4, 5].Einav和Puzrin则将压力相关性的超弹性模型和相应塑性模型联合模拟土的力学行为[6].Houlsby等基于弹性模量为平均应力的幂函数推导了满足热力学要求的超弹性模型[3].考虑到颗粒材料力学特性，超弹性模型仅适用于应变非常小的颗粒材料体系，当需要描述颗粒材料的破坏行为时，通常将超弹性模型与塑性模型或损伤模型相结合.Volokh[7]提出的基于超弹性软化模型的概念则给出了另外一条可能的途径.超弹性软化模型的主要思想为[7]:①真实材料的应变能存在一个极限值；②真实材料的应变能为理想材料应变能的函数.随着应变增大(假定趋于无穷大)，真实材料的应变能将趋于极限值，其应力-应变关系表现出软化现象.必须指出Volokh超弹性模型是针对脆性材料或生物体材料提出的，且其模型假定真实材料和理想材料的应变相同，但应力不同.因此在超弹性软化模型中通常以应变描述应变能.
超弹性软化模型概念清晰，且形式简单，易于数值实现.本文试图将针对脆性材料和生物体材料的Volokh超弹性模型扩展至颗粒材料.考虑颗粒材料主要以剪切变形和破坏为主，文献[8]发展一个颗粒材料的超弹性软化模型，本文在此基础上，进一步建议体积应变能和剪切应变能区别处理，得到了两种修正的超弹性软化模型.并基于ABAQUS平台开发了相应的程序代码.数值算例考察了修正模型模拟应变局部化及应变软化现象的能力，验证了所建立模型的可行性与正确性.
1 Volokh超弹性软化模型为保持论文完整性，简单介绍Volokh超弹性软化模型的概念[7].以ψ(C)表示真实材料的应变能，其应力可表示为
其中：J=det(F)；C为Cauchy-Green变形张量，其与变形梯度F的关系为
假定存在一个与真实材料具有相同应变的理想材料，该理想材料的应变能以W表示，若真实材料性质与理想材料相同，则
此时，真实材料的应变能具有下列特征：
其中‖·‖表示张量的模.式(4)表明对理想材料(或与理想材料具有相同性质的真实材料)而言随着应变增加，应变能是无限制增加的.考虑到一般情况下，真实材料会在一定条件下发生破坏，意味着真实材料的应变能随应变增加到某一程度时具有一个临界值[7]，即
其中，Φ表示材料破坏时的极限应变能.式(5)所表述的材料破坏准则在本质上与传统的Mises准则或Tresca准则类似，均表示变形能存在临界值.满足式(5)的应变能具有多种表达式，一般可假定真实材料的应变能为理想材料应变能的函数，即
Volokh给出的一个表达式为
其中，
将式(7)代入式(1)，相应地有
对比式(1)可以看到，式(9)存在指数因子exp(-W/Φ)，随应变的增加，理想材料的应变能W不断增加，指数因子exp(-W/Φ)从1向0趋减，应力增加的趋势不断减小，并且应力在达到某一极限值后开始减小.这使得式(9)能够反应材料的应变软化现象.基于式(9)的模型称为超弹性软化模型.
对于线性各向同性超弹性材料，在小变形的情况下其应变能函数可表示为
式中，Δ和us2分别为应变第一主不变量和偏应变第二主不变量.
2 修正的超弹性软化模型考虑到岩土颗粒材料在破坏过程中体变部分与形变部分所起作用不同，因此在软化模型中采用不同的处理方式.
2.1 第一种软化方式：仅对形变部分软化考虑颗粒材料的变形和破坏特性，通常认为颗粒材料是以剪切变形和破坏为主，为此仅在应变能的剪切部分引入了软化函数，则有
进一步可认为颗粒材料的应变能由部分修正后的应变能(表征软化部分的材料)和部分理想材料的应变能组成，即
式中α为表示材料软化程度的因子，本文建议采用式(13)计算：
可见因子α在一定程度上可视为扰动状态因子[9, 10].式中α0为初始软化参数，β为模型参数，表征了软化对等效应变的依赖程度，ε为等效应变，采用
计算.对应变能函数ψM1求导可得到相应的应力，进一步对应变求导可得到相应的模量张量.
2.2 第二种软化方式：体变部分与形变部分采用不同的极限值对于颗粒材料而言，考虑虽然其更容易发生剪切破坏，但也可发生体变部分的破坏.为此本文假定极限体积应变能大于极限剪切应变能，分别考虑其软化，将第二种软化函数表示为
其中：Φ1为极限体积应变能；Φ2为极限剪切应变能.按照处理第一种软化方式相似的方法，进一步将应变能函数修正为
式中α计算见式(13).由式(15)可得到相应的应力和模量张量.
3 数值算例 算例1考虑图 1所示平面应变条件下的颗粒材料集合体，尺寸为0.6 m×0.8 m，上边界单元节点在水平方向被约束，下边界x,y 方向的自由度被约束，左右边界自由，颗粒集合体夹在两刚性板之间，承受由位移控制竖向加载.数值模拟所采用的参数为：体积模量K=60 MPa,剪切模量G=30 MPa，参数Φ=7.5 kJ/m3，参数Φ1=15.0 kJ/m3，参数Φ2=7.5 kJ/m3，α0=0.5，β=10.为考察式(12)与式(15)所建立模型对颗粒材料力学行为的模拟能力及参数α0和β的影响，图 2给出了Volokh软化模型与第一种和第二种修正的软化模型所预测的位移承载曲线.由图 2可以看出，随着位移加载，Volokh软化方式与两种软化方式所预测的颗粒集合体的承载力均先逐渐升高，当加载位移达到9 mm时，Volokh软化方式的承载力开始降低，材料进入软化阶段，当加载位移达到11 mm时，两种软化方式的承载力开始降低，材料进入软化阶段.图 3给出了加载位移为15 mm时两种软化方式所预测的等效应变分布图.可以看出，两种软化方式均预测出了明显的剪切带，第一种软化方式的高等效应变区域较第二种软化方式大.考察α0的变化对于材料软化的影响，图 4所示为第一种软化方式在加载位移为15 mm，α0分别取0.25，0.50和0.75且β取10时的位移承载曲线图.可以看到，α0在取0.50时软化明显.由于结果的相似性，未讨论第二种软化方式在不同α0取值时的位移承载曲线图.考察β的变化对于材料软化的影响，图 5给出第一种软化方式在加载位移为15 mm，β分别取5，10和15且α0取0.50时的位移承载曲线，当β取10时，材料软化明显.第二种软化方式结果类似，故未讨论第二种软化方式在不同β取值时的位移承载曲线图.从图 4可以看出α0取0.50时等效应变区域最大，剪切带最为明显.从图 5可以看到β取10时等效应变区域最大，且剪切带最为明显.图 6所示为不同加载时刻时第一种软化方式的等效应变分布图，可以看到，随着加载位移的增大，应变局部化逐渐增加，剪切带逐渐产生.由于第二种软化方式在不同α0与β取值时的等效应变分布与第一种软化方式类似，故未给出第二种软化方式的等效应变分布图.
图1(Figure 1)
 | 图1 平板压缩问题网格图 Fig. 1 Diagram of mesh for the rectangle plane subjectedto a prescribed vertical displacement |
图2(Figure 2)
 | 图2Volokh软化方式与两种软化方式所预测的位移承载曲线图(α0=0.5，β=10) Fig. 2 Load-displacement curves predicted by the Volokh softening model and the two softeningmodels (α0=0.5,β=10) |
图3(Figure 3)
 | (a)—第一种软化方式； (b)—第二种软化方式. 图3等效应变分布图(α0=0.5，β=10，加载位移=15 mm) Fig. 3 Equivalent strain distribution (α0=0.5,β=10load displacement=15 mm) |
图4(Figure 4)
 | (a)—位移-承载； (b)—α0=0.25； (c)—α0=0.50; (d)—α0=0.75. 图4第一种软化方式所预测的等效应变分布图 Fig. 4 Equivalent strain distribution predicted by the first softening model |
图5(Figure 5)
 | (a)—位移-承载; (b)—β=5； (c)—β=10; (d)—β=15.图5 第一种软化方式所预测的等效应变分布图 Fig. 5 Equivalent strain distribution predicted by the first softening model with β=5, 10, 15, α0=0.50under the load displacement of 15 mm |
图6(Figure 6)
 | (a)—加载位移=1.5 mm； (b)—加载位移=6.0 mm; (c)—加载位移=10.5 mm； (d)—加载位移=14.5 mm. 图6 第一种软化方式所预测的不同时刻的等效应变分布图(α0=0.5，β=10) Fig. 6 Equivalent strain distribution predicted by the first softening model at different step time (α0=0.5,β=10) |
算例2考虑图 7所示地基变形，尺寸为34 m×17 m.忽略自重，地基通过位于其上的刚性基础而受到荷载作用，地基和基础之间的接触面假定为理想粘接.所施荷载通过位于基础中心的有限元网格节点垂直作用于地基.地基的底部边界固定，左右两边界水平方向固定.模拟所用参数为体积模量K=60 MPa,剪切模量G=30 MPa，参数Φ=4 kJ/m3，参数Φ1=10 kJ/m3，参数Φ2=4 kJ/m3，α0=0.5，β=10.
图7(Figure 7)
 | 图7 地基平面网格图 Fig. 7 Diagram of mesh for the soil foundation |
图 8为加载位移为30 cm时理想模型与两种软化模型所预测的等效应变分布图.从图中看到，第一种软化方式的等效应变局部化区域最大，其应变局部化最明显.理想模型的局部化区域最小，应变局部化不明显.
图8(Figure 8)
 | (a)—第一种软化方式； (b)—第二种软化方式; (c)—理想模型.图8 等效应变分布图(α0=0.5，β=10，加载位移=30 cm) Fig. 8 Equivalent strain distribution (α0=0.5,β=10, load displacement=30 cm) |
4 结 论1) 本文将Volokh超弹性软化模型推广到脆性颗粒材料，所建议的模型具有模拟颗粒材料应变软化及应变局部化的能力.
2) 表征软化部分比例的控制因子α本质上类似于扰动状态概念中的扰动因子，其演化应伴随着材料从理想状态向软化状态调整.软化因子α采用随等效应变线性变化的演化方程，能够较好地模拟颗粒材料应变软化及应变局部化现象.
3) 第二种软化方式实质上包含了第一种软化方式.当第二种软化方式的参数Φ1取到无穷大时，意味着未对体变部分进行软化，这时第二种软化方式退化为第一种软化方式.
4) 相对于颗粒材料的其他本构的模型，超弹性软化模型主要以真实材料的应变能存在极限为基础，概念清晰，表达式简单，易于数值实现，希望能给颗粒材料本构行为的研究工作者一些借鉴.
参考文献
| [1] | Humrickhouse P W,Sharpe J P,Corradini M L.Comparison of hyperelastic models for granular materials[J].Physical Review E,2010,81(1):011303.( 1) |
| [2] | Jiang Y,Zheng H,Peng Z,et al.Expression for the granular elastic energy[J].Physical Review E,2012,85(5):051304.(1) |
| [3] | Houlsby G T,Amorosi A,Rojas E.Elastic moduli of soils dependent on pressure:a hyperelastic formulation[J].Géotechnique,2005,55(5):383-392.(2) |
| [4] | Jiang Y,Liu M.A brief review of “granular elasticity” [J].The European Physical Journal E,2007,22(3):255-260.(1) |
| [5] | Jiang Y,Liu M.Granular elasticity without the Coulomb condition[J].Physical Review Letters,2003,91(14):144301.(1) |
| [6] | Einav I,Puzrin A M.Pressure-dependent elasticity and energy conservation in elastoplastic models for soils[J].Journal of Geotechnical and Geoenvironmental Engineering,2004,130(1):81-92.(1) |
| [7] | Volokh K Y.Hyperelasticity with softening for modeling materials failure[J].Journal of the Mechanics and Physics of Solids,2007,55(10):2237-2264.(3) |
| [8] | Chang J F,Chu X H,Xu Y J.A softening hyperelastic model and simulation of the failure of granular materials[J].Geomechanics and Geoengineering,2014,7(4):335-353.(1) |
| [9] | 楚锡华,孔科,徐远杰.基于扰动状态概念与 EB 模型的粗粒料力学行为模拟[J].应用力学学报,2012,29(2):141-147.(Chu Xi-hua,Kong Ke,Xu Yuan-jie.Numerical simulation of mechanical behavior of coarse granular materials based on E-B model combining with disturbed state concept[J].Chinese Journal of Applied Mechanics,2012,29(2):141-147.)(1) |
| [10] | 付平,楚锡华,余村,等.基于扰动状态概念的颗粒材料应变局部化行为模拟[J].华南理工大学学报(自然科学版),2014,42(4):59-69.(Fu Ping,Chu Xi-hua,Yu Cun,et al.Simulation of strain localization of granular materials based on disturbed state concept[J].Journal of South China University of Technology (Natural Science Edition),2014,42(4):59-69.)(1) |
本文献在全文中的定位：
... 性模型外，近年来，超弹性模型亦被用来描述岩土颗粒材料的力学行为[1, 2, 3] ...
1
本文献在全文中的定位：
... 用来描述岩土颗粒材料的力学行为[1, 2, 3].如Jiang和Liu基于热力学失稳和弹性 ...
2
本文献在全文中的定位：
...[1, 2, 3].如Jiang和Liu基于热力学失稳和弹性模量的压力相关性推导了一个能够计 ...
... 基于弹性模量为平均应力的幂函数推导了满足热力学要求的超弹性模型[3].考虑到颗粒材料力学特性，超弹性模型仅适用于应变非常小的颗粒材 ...
1
本文献在全文中的定位：
... 一个能够计算屈服和剪胀性的GE模型[2, 4, 5].Einav和Puzrin则将压力相关性的超弹 ...
1
本文献在全文中的定位：
...[2, 4, 5].Einav和Puzrin则将压力相关性的超弹性模型和相应塑性模型联合模拟土的 ...
1
本文献在全文中的定位：
... 将压力相关性的超弹性模型和相应塑性模型联合模拟土的力学行为[6].Houlsby等基于弹性模量为平均应力的幂函数推导了满足热力学要求的超 ...
3
本文献在全文中的定位：
... 料的破坏行为时，通常将超弹性模型与塑性模型或损伤模型相结合.Volokh[7]提出的基于超弹性软化模型的概念则给出了另外一条可能的途径.超弹 ...
... 为保持论文完整性，简单介绍Volokh超弹性软化模型的概念[7].以ψ(C)表示真实材料的应变能，其应力可表示 ...
... 坏，意味着真实材料的应变能随应变增加到某一程度时具有一个临界值[7]，即 ...
1
本文献在全文中的定位：
... 性模型扩展至颗粒材料.考虑颗粒材料主要以剪切变形和破坏为主，文献[8]发展一个颗粒材料的超弹性软化模型，本文在此基础上，进一步建议体 ...
1
本文献在全文中的定位：
... 可见因子α在一定程度上可视为扰动状态因子[9, 10].式中α0为初始软化 ...
1
本文献在全文中的定位：
... 在一定程度上可视为扰动状态因子[9, 10].式中α0为初始软化参数，β为模型参数，表征了软 ...
