3AxialNew

Step 1Step 2Step 3

%step1: packing of model
clear;
%-------------initial parameters and build model--------------
A=obj_3AxisTester();%build tester
A.randomSeed=1;
A.waterPressure=0.1e6;%static water pressure
A.isClump=0;%particle is not clump
A.ballR=0.005;%average element radius
A.name='3AxialTest';
A.type='3Axial';%type='uniaxialC';%uniaxial compression
A.distriRate=0.2;%distribution coefficient
A.loadingType='stress';%displacement or stress
A.stressStep=-[50;100;200;400;800;1600;3200]*1e3;
A.maxStrain=0.15;
A.setType();
A.buildIntialModel();
A.setUIoutput();%set output of form
d=A.d;
d.mo.setGPU('auto');

sampleId=d.GROUP.sample;
%d.mo.aR(sampleId)=d.mo.aR(sampleId)*0.99;d.mo.setNearbyBall();%adjust element radius here

d.mo.setGPU('auto');%set auto GPU
A.gravitySediment();%sedimentation of elements
topPlatenId=d.GROUP.topPlaten;
d.mo.mGZ(topPlatenId)=d.mo.mGZ(topPlatenId)*10;
mfs.reduceGravity(d,5);
d.balance(50,20+d.SET.packNum);

d.status.dispEnergy();%display energy information
d.showFilter('SlideY',0.3,1);
d.show('StressZZ');

d.mo.setGPU('off');
d.clearData(1);%clear dependent data
d.recordCalHour('AxialStep1Finish');
save(['TempModel/' A.name '1.mat'],'A','d');
save(['TempModel/' A.name '1R' num2str(A.ballR) '-distri' num2str(A.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();

%step2: assign formal mechanical parameters
%-------------assign new properties--------------
clear;
load('TempModel/3AxialTest1.mat');
A.setUIoutput();%set the output message
balanceNum=100;%define the balance number in the following simulation
d=A.d;
d.mo.setGPU('off');%please always close the GPU when modifing the model
A.importModel(d);
d=A.getFinalModel();%get the model from the object
%-----------use low gravity to balance the model---------------
%Hertz Model
%d.mo.FnCommand='nFN1=obj.nKNe.*nIJXn;nR=obj.aR(1:m_Num)*nRow;nJR=obj.aR(obj.nBall);Req=nR.*nJR./(nR+nJR);nE=obj.aKN(1:m_Num)*nRow./(pi*nR);nJE=obj.aKN(obj.nBall)./(pi*nJR);Eeq=nE.*nJE./(nE+nJE);nFN2=-4/3*Eeq.*Req.^(1/2).*abs(nIJXn).^(3/2);f=nIJXn<0;nFN0=nFN1.*(~f)+nFN2.*f;';
d.mo.mGZ(:)=0;
d.status=modelStatus(d);
aMUp0=d.mo.aMUp;%record the friction coefficient
d.mo.aMUp(:)=0;%no friction between balls
d.mo.dT=d.mo.dT*4;%use greater step time
d.mo.setGPU('auto');%set auto GPU
d.balance(balanceNum,20+d.SET.packNum);%d.SET.packNum is the particle number along vertical direction
A.setTubeFixId();%set the fixed Id of tube
d.balance(balanceNum,20+d.SET.packNum);
d.mo.aMUp=aMUp0;%restore friction coefficient
d.mo.dT=d.mo.dT/4;%restore the step time

d.showFilter('SlideY',0.3,1);
d.show('StressZZ');

d.mo.setGPU('off');
d.clearData(1);%clear dependent data
d.recordCalHour('AxialStep2Finish');%record the time
save(['TempModel/' A.name '2.mat'],'A','d');
save(['TempModel/' A.name '2R' num2str(A.ballR) '-distri' num2str(A.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();

clear;
load('TempModel/3AxialTest2.mat');
A.setUIoutput();%set the output message
d=A.d;
d.calculateData();
d.mo.setGPU('off');%please always close the GPU when modifing the model
%----------set the parameters of numerical simulation
balanceRate=1;
balanceNum=20;%balance number in each step
loopNum=0;
totalCircle=5;%default value is 50
stepNum=100;
isSave=1;%save the .mat data

d.setStandarddT();%set the standard step time
d.mo.dT=d.mo.dT/balanceRate;
A.moveTop();%move top boundary to the top of the sample
%-----------set the loop-------------
platenArea=pi*A.sampleR^2;%the area of the platen
forceStep=A.stressStep*platenArea;%get the force step on the platen
botPlatenNum=length(d.GROUP.botPlaten);%element number of the platen
sampleH=mean(d.mo.aZ(d.GROUP.topPlaten))-mean(d.mo.aZ(d.GROUP.botPlaten));%original height of model
d.status=modelStatus(d);%initialize the status recorder
if strcmp(A.loadingType,'stress')
    totalCircle=length(forceStep);
end

botPlatenZZ=zeros(totalCircle+1,1);%initialize the position of bottom platen
topPlatenZZ=zeros(totalCircle+1,1);
botStressZZ=zeros(totalCircle+1,1);%initialize the position of bottom boundary
topStressZZ=zeros(totalCircle+1,1);
botPlatenZZ(1)=mean(d.mo.aZ(d.GROUP.botPlaten));
topPlatenZZ(1)=mean(d.mo.aZ(d.GROUP.topPlaten));
d.mo.mGZ(:)=0;
%-----------end set the loop-------------
%-----------the loop-------------
fName=['data/step/' A.name  num2str(A.ballR) '-' num2str(A.distriRate) 'loopNum'];
d.tic(totalCircle*stepNum);%record start time
save([fName '0.mat']);%return;
for circle=1:totalCircle
    if strcmp(A.loadingType,'stress')%when loading type is stress
        d.mo.mGZ(d.GROUP.topPlaten)=forceStep(circle)/botPlatenNum;%increase the stress
    end
    for step=1:stepNum
        for i=1:balanceNum*balanceRate
            d.mo.balance();
        end
        %d.status.SET.isWHT=1;%record the WHT in recordStatus
        d.recordStatus();
        d.toc();%show the note of time
    end
    botPlatenZZ(circle+1)=mean(gather(d.mo.aZ(d.GROUP.botPlaten)));
    topPlatenZZ(circle+1)=mean(gather(d.mo.aZ(d.GROUP.topPlaten)));
    botStressZZ(circle+1)=gather(d.status.bottomBFs(end,3))/platenArea;
    topStressZZ(circle+1)=gather(d.status.topBFs(end,3))/platenArea;
    if isSave==1
        d.clearData(1);
        save([fName num2str(circle) '.mat']);
        d.calculateData();
    end
end
A.data.topPlatenZZ=topPlatenZZ;%record the data
A.data.topStressZZ=topStressZZ;
A.data.botPlatenZZ=botPlatenZZ;
A.data.botStressZZ=botStressZZ;
A.data.sampleH=sampleH;

if  strcmp(A.loadingType,'stress')&&strcmp(A.type,'3Axial')
    loopNum=10;
end
for circle=1:loopNum%additional balance loops
    d.balance(100,5);%balance the model 500 times, record 5 times
    save(['data\step\' A.type '-R' num2str(A.ballR) '-' num2str(A.distriRate) 'circleLoop' num2str(circle) '.mat'],'A','d');
end
%-----------end the loop-------------

%------------return and save result--------------
d.setData();%set data for show()
data=A.data;
%calculate the strain and stress of sample
dPlatenZZ=data.topPlatenZZ-data.topPlatenZZ(1);
A.data.strain=-dPlatenZZ/data.sampleH;
A.data.stress=-data.botStressZZ;
figure;
plot(A.data.strain,A.data.stress,'-o');
xlabel('StrainZZ');
ylabel('StressZZ');

d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('AxialStep3Finish');
save(['TempModel/' A.name '3.mat'],'A','d');
save(['TempModel/' A.name '3R' num2str(A.ballR) '-distri' num2str(A.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();