【COMSOL】Marzas 材料模型 C 源文件代码解析

文件头

该材料模型输入为应变、材料参数、模型状态，输出为应力、雅可比矩阵，供 COMSOL 使用。

#include <math.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#ifdef _MSC_VER
#define EXPORT __declspec(dllexport)
#else
#define EXPORT
#endif
#define MIN(a,b) (((a)<(b))?(a):(b))
#define MAX(a,b) (((a)>(b))?(a):(b))

文件头里这些定义我只能看懂定义了MIN、和MAX

函数主体

EXPORT int eval(double e[6],       // Input: Green-Lagrange strain tensor components in Voigt order (xx,yy,zz,yz,zx,xy)double s[6],       // Output: Second Piola-Kirchhoff stress components in Voigt order (xx,yy,zz,yz,zx,xy)double D[6][6],    // Output: Jacobian of stress with respect to strain, 6-by-6 matrix in row-major orderint *nPar,         // Input: Number of material model parameters, scalar,double *par,       // Input: Parameters: par[0] = E0, par[1] = nu0, ...int *nStates,      // Input: Number of states, scalar               double *states) {  // States, nStates-vector

该计算应力结果的函数的输入有：
e[6]：格林-朗格朗日应变张量，以 Voigt 顺序存储，xx, yy, zz, yz, zx, xy。
s[6]: 皮奥拉·基尔霍夫应力，以 Voigt 顺序存储，xx, yy, zz, yz, zx, xy。
以上两者都是定义在参考构型上的。
D[6][6]: 应力-应变的雅可比矩阵。
nPar: 材料模型参数个数。
par: 材料模型参数数组。
nStates：状态数量。
states：状态数组。

在 COMSOL 软件中的外部材料选择类型“广义应力-应变关系”，就确定了输入量和输出量，还要求自定义材料属性明细，在里面，可以定义材料参数。在这个案例中，材料参数有10个：
E 0 , ν 0 E_0,\nu_0 E0,ν0：初始杨氏模量、泊松比
stchelp1、stchelp2、stchelp3：三个主拉伸伸长率
κ 0 \kappa_0 κ0：初始损伤阈值
A t , B t , A c , B c A_t,B_t,A_c,B_c At,Bt,Ac,Bc：模型中其它参数。（模型详情见上一篇）

  int i, j;double E0, nu0, stch[3], kappa0, At, Bt, Ac, Bc, kappa, ep[3], evol, lambLame, muLame;double sp[3], st[3], stsum, et[3], eef2, eef, ets, ecs, alphat, alphac, damt, damc, damage, E;

定义一些函数内变量，用于接受输入值，和中间计算过程。

  // Check inputsif (nPar[0] != 10)       // 10 inputs, E0, nu0, stretches and Mazars' parametersreturn 1;              // error code 1 = "Wrong number of parameters"if (nStates[0] != 2)     // 2 states to memorize damage variablesreturn 2;              // error code 2 = "Wrong number of states"

对输入参数数量进行检查。

  // Read input parameters from parameter vector, call convention://  {E0, nu0, comp1.solid.stchelp1, comp1.solid.stchelp2, comp1.solid.stchelp3, kappa0, At, Bt, Ac, Bc}E0 = par[0];    // undamaged Young's modulusif (E0 <= 0.0) return -1;nu0 = par[1];                  // undamaged Poisson's ratioif (nu0 >= 0.5 || nu0<= -1.0) return -1;        for (i = 0; i < 3; i++) {stch[i] = par[i+2];          // principal stretchesif (stch[i] <= 0.0) return -1;}kappa0 = par[5];               // Mazars' parameter: initial damage thresholdif (kappa0 <= 0.0) return -1;At = par[6];                   // Mazars' parameter Atif (At <= 0.0) return -1;Bt = par[7];                   // Mazars' parameter Btif (Bt <= 0.0) return -1;Ac = par[8];                   // Mazars' parameter Acif (Ac <= 0.0) return -1;Bc = par[9];                   // Mazars' parameter Bcif (Bc <= 0.0) return -1;kappa = MAX(states[0],kappa0); // kappa memorizes maximum equivalent tensile strain, initialized to kappa0

以上代码将输入参数值保存到定义的变量中，并对值的有效性进行检查。states的第一个值存储等效拉伸应变的历史最大值，第二个值存储损伤变量。

  // Define Lame parameterslambLame = E0 * nu0 / (1.0 + nu0) / (1.0 - 2.0 * nu0); // undamaged Lame parameter lambdamuLame = E0 / 2 / (1.0 + nu0);                         // undamaged Lame parameter mu

以上代码用初始杨氏模量和泊松比定义拉梅常数： λ , μ \lambda,\mu λ,μ.

  // Define principal engineering strains computed from principal stretchesevol = 0.0;for (i = 0; i < 3; i++) {ep[i] = stch[i] - 1.0;   // principal strainsevol += ep[i];           // volumetric strain}

以上代码由主拉伸率计算主工程应变，并相加计算体应变。

  stsum = 0.0, eef2 = 0.0;for (i = 0; i < 3; i++) {sp[i] = lambLame * evol + 2.0 * muLame * ep[i]; // principal stressesst[i] = MAX(sp[i], 0.0);                        // tensile stresses, positive only (Mazars' model)stsum += st[i];                                 // sum of tensile stresseset[i] = MAX(ep[i], 0.0);                        // tensile strains, positive only (Mazars' model)eef2 += et[i] * et[i];                          // equivalent tensile strain, squared.}eef = sqrt(eef2);                                 // equivalent tensile straineef2 += 1e-10;                                    // avoid divide by zero

sp[]存储由拉梅常数、体应变、主应变计算的主应力。
Mazars模型中定义的拉伸应力st[]，是对主应力取正值。然后计算等效拉伸应变。

ϵ e f = [ ϵ 1 ] 2 + [ ϵ 2 ] 2 + [ ϵ 3 ] 2 \epsilon_{ef}=\sqrt{[\epsilon_1]^2+[\epsilon_2]^2+[\epsilon_3]^2} ϵef=[ϵ1]2+[ϵ2]2+[ϵ3]2

  alphat = 0.0, alphac = 0.0;for (i = 0; i < 3; i++) {ets = (1.0 + nu0) / E0 * st[i] - nu0 / E0 * stsum; // define tensile strain from tensile stress and undamaged stiffness (Mazars' model)ecs = ep[i] - ets;                                 // define compressive strain (Mazars' model)alphat +=  ets * MAX(ep[i], 0.0) / eef2;           // tensile weight (Mazars' model)alphac +=  ecs * MAX(ep[i], 0.0) / eef2;           // compressive weight (Mazars' model)}

由拉伸应力通过初始材料参数计算拉伸应变，随后计算压缩应变。带入计算权重因子：

α t = ∑ i = 1 3 H i ϵ t i ϵ i / ϵ e f 2 \alpha_t=\sum_{i=1}^3H_i\epsilon_{ti}\epsilon_i/\epsilon_{ef}^2 αt=∑i=13Hiϵtiϵi/ϵef2

α c = ∑ i = 1 3 H i ϵ c i ϵ i / ϵ e f 2 \alpha_c=\sum_{i=1}^3H_i\epsilon_{ci}\epsilon_i/\epsilon_{ef}^2 αc=∑i=13Hiϵciϵi/ϵef2

  //Initialize damage variable with stored valuedamage = states[1];// Increase damage and update states only if the equivalent tensile strain increases its previous valueif(eef > kappa) { kappa = eef; // update kappa to the current value of equivalent strain// Define tensile and compressive damage variables (Mazars' model)damt = 1.0 - kappa0 * (1.0 - At) / kappa - At * exp(-Bt * (kappa - kappa0));damc = 1.0 - kappa0 * (1.0 - Ac) / kappa - Ac * exp(-Bc * (kappa - kappa0));// Define damage variabledamage = MAX(damage, alphat * damt + alphac * damc); // increase damage only if larger than its previous value stored in the state variabledamage = MIN(MAX(0.0, damage), 0.99); // Limit the damage variable between 0 and 0.99. Full damage means 0.99 stiffness reduction.states[0] = eef;     // use this state to memorize the historical maximum value of equivalent strainstates[1] = damage;  // update the state with the new value of the damage variable}

如果计算出的等效拉伸应变大于历史最大值，需要计算新的损伤因子 d d d，并且在 states 中更新等效拉伸应变历史最大值和新的损伤因子。

// Define damaged Young's modulusE = (1.0 - damage) * E0;   // Set up Jacobian matrix from damaged Young's modulus// 6x6 Elasticity matrix D based on E and nu0, initially filled by zerosfor (i = 0; i < 6; i++){for (j = 0; j < 6; j++) {D[i][j] = 0.0;}}  D[0][0] = D[1][1] = D[2][2] =  E * (1.0 - nu0) / (1.0 + nu0) / (1.0 - 2.0 * nu0);                 // upper diagonalD[3][3] = D[4][4] = D[5][5] =  E / (1.0 + nu0);                                                   // lower diagonal, note that this equals 2GD[0][1] = D[0][2] = D[1][0] = D[1][2] = D[2][0] = D[2][1] = E * nu0 / (1.0 + nu0) / (1.0 - 2.0 * nu0);// off diagonal

计算损伤后的模量与雅可比矩阵：

σ i j = 2 G ϵ i j + δ ( i , j ) λ ϵ j j \sigma_{ij}=2G\epsilon_{ij}+\delta(i,j)\lambda\epsilon_{jj} σij=2Gϵij+δ(i,j)λϵjj

  // Compute Hooke's law: s = D*efor (i = 0; i < 6; i++){s[i] = 0.0;for (j = 0; j < 6; j++) {s[i] += D[i][j] * e[j];}}  // Return value 0 if success, any other value trigger an exceptionreturn 0;

最后是计算应力。

总结

该模型使用当前应变状态去计算损伤因子，然后计算损伤后模量和输出应力。上一篇提到细化单元数量后，不收敛。尝试后发现，如果将单元类型设置为线性单元，就能够收敛。

【COMSOL】Marzas 材料模型 C 源文件代码解析相关推荐

对抗思想与强化学习的碰撞-SeqGAN模型原理和代码解析
GAN作为生成模型的一种新型训练方法,通过discriminative model来指导generative model的训练,并在真实数据中取得了很好的效果.尽管如此,当目标是一个待生成的非连续性序 ...
时间序列模型SCINet（代码解析）
前言 SCINet模型,精度仅次于NLinear的时间序列模型,在ETTh2数据集上单变量预测结果甚至比NLinear模型还要好. 在这里还是建议大家去读一读论文,论文写的很规范,很值得学习,论文地址 ...
Bert模型介绍及代码解析（pytorch）
Bert(预训练模型) 动机基于微调的NLP模型预训练的模型抽取了足够多的信息新的任务只需要增加一个简单的输出层注:bert相当于只有编码器的transformer 基于transformer ...
【COMSOL】在结构力学中使用自定义外部材料模型 · Mazars 损伤模型
前言虽然 COMSOL 提供了一些常用材料模型,但是,在软物质.断裂等研究中,需要使用材料库中不存在的材料模型形式.从5.2版本开始,支持用C语言自定义外部材料模型,编译成链接库使用. 这里通过一个 ...
关于Transformer你需要知道的都在这里------从论文到代码深入理解BERT类模型基石（包含极致详尽的代码解析！）
UPDATE 2.26.2020 为代码解析部分配上了Jay Ammar The Illustrated GPT-2 的图示,为想阅读源码的朋友缓解疼痛! 深入理解Transformer------从 ...
comsol分析时总位移代表什么_超弹性材料模型的压缩分析
为了表征超弹性材料,需要进行各种测试获取实验数据,包括承受单轴拉伸和压缩.双轴拉伸和压缩以及扭转测试.今天,我们向大家介绍如何使用通过单轴和双轴测试获得的拉伸和压缩测试数据,模拟由弹性泡沫材料制成的球 ...
【对比学习】CUT模型论文解读与NCE loss代码解析
标题:Contrastive Learning for Unpaired Image-to-Image Translation(基于对比学习的非配对图像转换) 作者:Taesung Park, Ale ...
DeepLearning | 图注意力网络Graph Attention Network（GAT）论文、模型、代码解析
本篇博客是对论文 Velikovi, Petar, Cucurull, Guillem, Casanova, Arantxa,et al. Graph Attention Networks, 2018 ...
代码解析之自行车模型在Apollo规划中的应用
大家好,我已经把CSDN上的博客迁移到了知乎上,欢迎大家在知乎关注我的专栏慢慢悠悠小马车(https://zhuanlan.zhihu.com/duangduangduang).希望大家可以多多交流, ...

【COMSOL】Marzas 材料模型 C 源文件代码解析

文件头

函数主体

总结

【COMSOL】Marzas 材料模型 C 源文件代码解析相关推荐

最新文章

热门文章