HGO是一个经典的描述生物组织力学行为的各向异性超弹本构，本文给出了实现HGO各向异性超弹UMAT的所有理论推导，免费分享给大家，码字不易，欢迎关注“九千CAE”公众号（b站、技术邻、仿真秀同昵称），欢迎动动小手一键三联（点赞+在看+分享）！点击“阅读原文”获取配套的视频课程和子程序源代码。
目录：
HGO超弹本构及应力应变关系
应变能密度函数
应力推导
摄动法计算切线刚度矩阵
ABAQUS UMAT中变形梯度的处理
局部坐标系下的变形梯度描述
ABAQUS UMAT中变形梯度
HGO超弹本构及应力应变关系
应变能密度函数
超弹模型假设一应变能密度函数，若该函数与取向有关，则相应的本构为各向异性超弹本构。HGO超弹本构的应变能密度函数如下
注意
应力推导
考虑
不变量关于变形梯度的偏导
由于应变能密度函数是表示为不变量的函数，所以先推导不变量关于变形梯度的偏导。
应变能密度关于不变量的偏导
若，则
注意，上式中每有一项，对应的求和项应划去。合并两种情形
若，则
若，则
合并两种情形有
注意，上式中没有求和符号
应力表达式
各向同性部分
各向异性部分
体积变形部分
摄动法计算切线刚度矩阵
由于HGO的切线刚度矩阵推导复杂，采用数值方法（摄动法）计算切线刚度矩阵。ABAQUS UMAT中采用Jaumann应力客观率
其中为旋转率张量，为变形率张量。由于，有
对变形梯度ij项进行摄动
其中为单位基矢。则
表明上述变形梯度的摄动只引起变形，不引起刚性转动。
同时
以 为例
有
故每摄动一次，可以确定刚度矩阵中的一列，共需摄动6次。最终
ABAQUS UMAT中变形梯度的处理
局部坐标系下的变形梯度描述
系统坐标系描述下，变形梯度张量定义为
其中和分别是参考构型和当前构型下某一物质点的位置向量，注意两者都是描述在系统坐标系下的。变形梯度极分解为
其中是右伸长张量，是左伸长张量，是旋转张量，描述刚体转动。设初始的单元局部坐标系为，变形后的共旋局部坐标为，其中。则定义相应的变形梯度
若将变形梯度完全表示在初始单元局部坐标系下
由于局部坐标是共旋的，有，则有，从而有
或
由于是正交矩阵，有，则
考虑极分解，有
ABAQUS UMAT中变形梯度
采用三维实体单元且设置局部坐标系时，ABAQUS在UMAT中传入的变形梯度为
因此，在编写相应的子程序时，需要先将分解获得，并将超弹应力表达式中的用代替。
