Abaqus中材料参数随机场实现
1、前言
考虑材料参数空间变异性的岩土工程对象的数值分析是岩土工程研究中重要分支。当前，考虑材料参数空间变异性（即参数随机场）的分析手段中，除了极少数非主流的数值分析软件可以实现一键式随机场分析以外，大多数复杂的随机场实现都存在较高的门槛，且难以实现复杂的岩土对象相互作用分析。在主流岩土工程分析软件中，如，flac和abaqus中，前者需要使用fish编程，且基本要借助第三方软件计算随机场才能实现；后者可以利用python进行前处理，实现随机场的过程相对较为简单，且适合进行复杂的岩土工程对象的数值分析，因此基于abaqus的数值分析更具有普适性。
在abaqus中实现随机场的过程是先利用其他语言（如matlab,python）生成随机场结果文件，然后在abaqusCAE中将结果文件中材料参数分别赋值给每个单元；然后批量生成inp文件，最后批量计算inp文件，其中，生成inp文件可以在CAE中进行，也可以在matlab或者python里面直接编辑生成inp文件，批量计算可以在command或者cae中进行。
2、随机场文件生成
随机场的生成主要参考的文献是《考虑自相关函数影响的边坡可靠度分析》，这篇文献后面列出了生成随机场的matlab代码，其核心的算法是采用乔列斯基分解，5000个单元以内时候，matlab的计算速度是很快的。以下是以函数形式调用的随机场生成算法。
function [field]=midpoint_RF (Coord, mu,cov,dh,dv,Nsim,ACF)
%考虑自相关函数影响的边坡可靠度分析,李典庆
sigma=mu.*cov;
% rxy=0;
% rou=[1 rxy
% rxy 1];
rou=eye(size(mu,1));
L1=chol(rou); %L1==rou
% 读取随机场单元中心点横坐标x和纵坐标y
mLem=length(Coord);
% ACF=1 指数型(SNX); ACF=2 高斯型(SQX); ACF=3 二阶自回归型(CSX); ACF=4 指数余弦型(SMK)
for k=1:size(mu,1)
pxy=zeros(mLem);
for i=1:mLem
for j=1:mLem
dx=abs(Coord(i,1)-Coord(j,1));
dy=abs(Coord(i,2)-Coord(j,2));
switch ACF
case 1 % 指数型自相关函数(SNX)
pxy(i,j)=exp(-2*(dx/dh(k)+dy/dv(k)));
case 2 % 高斯型自相关函数(SQX)
pxy(i,j)=exp(-pi*((dx/dh(k))^2+(dy/dv(k))^2));
case 3 % 二阶自回归型自相关函数(CSX)
pxy(i,j)=exp(-4*(dx/dh(k)+dy/dv(k)))*(1+4*dx/dh(k))*(1+4*dy/dv(k));
case 4 % 指数余弦型自相关函数(SMK)
pxy(i,j)=exp(-(dx/dh(k)+dy/dv(k)))*cos(dx/dh(k))*cos(dy/dv(k));
case 5 % 三角型自相关函数(BIN)
if dx<dh(k) & dy<dv(k)
pxy(i,j)=(1-dx/dh(k))*(1-dy/dv(k));
else
pxy(i,j)=0;
end
end
end
end
PXY(:,:,k)=pxy;
end
for k=1:size(mu,1)
L2(:,:,k)=chol(PXY(:,:,k))';
end
% 拉丁超立方抽样产生独立标准正态随机样本（固定样本种子）
randn('state',0);rand('state',0);
% UU=lhsnorm(zeros(2*mLem,1),eye(2*mLem), Nsim)';
UU=lhsnorm(zeros(size(mu,1)*mLem,1),eye(size(mu,1)*mLem), Nsim)';
% 计算ln-xi的标准差和ln-xi的均值
sLn=sqrt(log(1+(sigma./mu).^2));mLn=log(mu)-sLn.^2/2;
for imod=1:Nsim
U=[];
for k=1:size(mu,1)
U=[U UU(((k-1)*mLem+1):k*mLem,imod)];
end
U_=U*L1;
H0=[];
for k=1:size(mu,1)
H0=[H0 L2(:,:,k)*U_(:,k)];
field(:,imod,k)=exp(mLn(k)+sLn(k)*H0(:,k));
end
end
end
这个函数需要设置参数有
mu=[2e7 ]'; cov=[0.3 ]';%均值及变异系数 dh=[20 ]; dv=[5 ]; %自相关距离 Nsim=10;%随机场数量 ACF=1;%自相关函数类型
此外还需要所有的需要赋值随机参数的单元及节点坐标，这两个文件可以在inp文件中直接提取，其源头可以通过在abaqus中建立的确定性分析模型得到。程序运行后即可生成指定个数的参数随机场，每一个参数的虽有随机场数据保存为一个txt文件，每一个txt文件中的每一列为一个参数的一个随机场数据。
上述过程也可以在python中实现。关于python和matlab的对比，我个人觉得matlab在计算速度、变量值查看方面比python优秀，python载入库是我觉得比较繁琐的事情，因此，复杂的计算和文件操作都用matlab进行。
3、随机参数的赋值
在abaqus中，前处理都是利用python进行的，且在CAE中每一个操作，都会在abaqus的工作目录中的abaqus.rpy中生成相对应的代码。不熟悉abaqus某一步操作代码时候，都可以在cae中进行一步操作，然后去看abaqus.rpy中生成的最新的一句代码。另外，还可以使用file->macro manager录制相应cae中操作的宏代码，录制完成后会在工作目录下面生成abaqusMacros.py文件，里面是相应操作的python代码。Abaqus中运行python代码有两种方式，第一种是在abaqus下面的输入框中输入代码，如下图；还有一种方法是在abaqusCAE中，file->run script。
在随机场的模拟中，模型中每一个单元都有相应的材料参数，每个单元赋予材料参数值的python代码如下：
for ele in elelist: ele=int(ele) mdb.models['Model-1'].Material(name='TempM-'+str(ele)) mdb.models['Model-1'].materials['TempM-'+str(ele)].Elastic(table=((para1[ele-1][inp],mu[1]),)) mdb.models['Model-1'].HomogeneousSolidSection(name='TempS-'+str(ele),material='TempM-'+str(ele),thickness=None) p.Set(elements=myele[(ele-1):(ele)],name='Tempset-'+str(ele)) region=p.sets['Tempset-'+str(ele)] p.SectionAssignment(region=region,sectionName='TempS-'+str(ele),offset=0.0, offsetType=MIDDLE_SURFACE,offsetField='', thicknessAssignment=FROM_SECTION)
Abaqus是可以进行inp提交计算的。赋予材料参数后，需要进一步得到inp文件，以便后续批量提交inp进行计算。Abaqus中，inp的生成代码如下：
mdb.Job(name='Job-'+str(inp), model='Model-1', description='', type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, echoPrint=OFF, modelPrint=OFF, contactPrint=OFF, historyPrint=OFF, userSubroutine='', scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=1, numGPUs=0) mdb.jobs['Job-'+str(inp)].writeInput(consistencyChecking=OFF)
4、inp文件批量生成
Abaqus中run script的方法获得inp文件的方法耗时很长，特别是单元数量超过1000时，利用abaqusCAE生成inp文件的做法基本是不可取的。我们可以在abaqus中生成一个inp，然后在修改inp文件中材料参数部分得到一个新的inp。
Importdata是matlab中读取数据的函数，数据类型可以使txt、csv等，其中当类型为txt文件时，可以设置文件头的行数deadline，headline之前的文本或者数据会读取为cell，后面的数据会存为矩阵，其中cell是将txt文件中每一行存为cell中一个元素。因此，只要设置headline大于inp文件的最大行数，就可以利用将inp全部读取为cell。
例如，利用下面的语句，
inpfile=importdata('Job-0.inp',',',25000);
得到的inp读取结果为：
上述变量inpfiel中，关于材料参数的部分第一行是关键词Material，第二行是力学参数的类型，如是弹性还是塑性参数，紧接着第三行是参数值，如果还有其他类型的参数，如摩尔库伦参数，会在弹性参数后面叠加参数类型关键词和材料参数值。
不同随机场模型的inp，由于只有随机场不一样，也就是每个单元的材料参数值不一样，因此只要修改上述inp文件中材料参数值的部分就可以得到一个新的inp。读取inp然后修改inp的关键代码如下：
numberofINP=10;
maxlineINP=25000;%数值要大于最大行数
materialline=[14994 21360];%材料参数设置起始行和结束行
%% 载入数据
inpfile=importdata('Job-0.inp',',',maxlineINP);
param1=importdata('param1.txt');
%% 生成inp
for i=0:(numberofINP-1)
name_INP={'Job-',num2str(i),'.inp'};
INPname=[name_INP{1} name_INP{2} name_INP{3}];
%% 修改材料参数
for j=materialline(1):3:(materialline(2)-3)
line=(j-materialline(1))/3+1;
value={num2str(param1(line,i+1)),', ','0.13'};
assignmaterial=[value{1} value{2} value{3}];
inpfile(j,1)=cellstr(assignmaterial);
end
最后将修改好的inpfile写入到一个新的inp文件中即生成了一个新的inp文件，这个过程相比在abaqusCAE中生成inp文件快百倍。同样的，批量生成inp使用其他语言一样可以达成，如python。
5、inp文件批量提交计算
在abaqus中，有两种方式提交模型计算，一种是从model的方式，另一种是从inp文件提交。利用python提交inp文件进行计算的代码如下：
for i in range(numinp): mdb.JobFromInputFile(name='Job-'+str(i),inputFileName='E:\\temp\\Job-'+str(i)+'.inp', type=ANALYSIS, atTime=None, waitMinutes=0, waitHours=0, queue=None, memory=90, memoryUnits=PERCENTAGE, getMemoryFromAnalysis=True, explicitPrecision=SINGLE, nodalOutputPrecision=SINGLE, userSubroutine='', scratch='', resultsFormat=ODB, multiprocessingMode=DEFAULT, numCpus=1, numGPUs=0) mdb.jobs['Job-'+str(i)].submit(consistencyChecking=OFF)
提交计算后就是提取计算结果，这些可以先在cae中操作一遍，以得到提取结果的代码，然后将代码另存，以便在inp提交计算后批量提取结果。
查看更多评论 >