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提交計算后批量提取結果。