TunnelHeat

Step 1Step 2Step 3

%build the geometrical model
clear;
fs.randSeed(1);%build random model
B=obj_Box;%build a box object
B.name='BoxTunnelHeat';
B.GPUstatus='auto';
ballR=0.1;%width, length, height, radius
distriRate=0.2;%define distribution of ball radius, 
isClump=0;
%--------------initial model------------
B.isUI=0;%when run the code in UI_command, isUI=1
B.ballR=ballR;
B.isClump=isClump;
B.distriRate=distriRate;
B.sampleW=5;
B.sampleL=5;
B.sampleH=5;
%B.BexpandRate=4;
%B.PexpandRate=4;
B.type='topPlaten';
%B.type='TriaxialCompression';
B.setType();
B.buildInitialModel();%B.show();

d=B.d;
%d.show('aR');return;
%--------------end initial model------------
B.gravitySediment();
B.compactSample(1);%input is compaction time
%------------return and save result--------------
mfs.reduceGravity(d,5);
d.balance('Standard');
d.status.dispEnergy();%display the energy of the model
d.clearData(1);%clear dependent data
d.recordCalHour('BoxTunnelHeat1Finish');
save(['TempModel/' B.name '1.mat'],'B','d');
save(['TempModel/' B.name '1R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();
d.show('aR');

clear
load('TempModel/BoxTunnelHeat1.mat');
%------------initialize model-------------------
B.setUIoutput();
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%get xyz from d.mo
%------------end initialize model-------------------

%---------------delele elements on the top
mZ=d.mo.aZ(1:d.mNum);%get the Z of elements
topLayerFilter=mZ>max(mZ)*0.9;
d.delElement(find(topLayerFilter));%delete elements according to id

%-------------set new material----------------
matTxt=load('Mats\StrongRock.txt');
Mats{1,1}=material('StrongRock',matTxt,B.ballR);
Mats{1,1}.Id=1;
d.Mats=Mats;
d.groupMat2Model({'sample'},1);
%-------------end set new material----------------

%-----------------make pipe--------------------
mX=d.mo.aX(1:d.mNum);mY=d.mo.aY(1:d.mNum);mZ=d.mo.aZ(1:d.mNum);mR=d.mo.aR(1:d.mNum);
innerR=0.2;
layerNum=2;
minBallR=min(mR)*0.9;
Rrate=1;
ringObj=f.run('fun/makeRing.m',innerR,layerNum,minBallR,Rrate);
tubeL=4;
dZ=minBallR*sqrt(3)/2;
ringNum=(tubeL-minBallR)/dZ+1;
ringZ=(1:ringNum)*dZ;
dAngle=180/(length(ringObj.R)/layerNum);
tubeObj=ringObj;
for i=1:ringNum-1
    newObj=mfs.rotate(ringObj,'XY',dAngle*i);
    newObj=mfs.move(newObj,0,0,dZ*i);
    tubeObj=mfs.combineObj(tubeObj,newObj);
end
outerR=(max(tubeObj.X)-min(tubeObj.X))/2;
%-----------------end make tube--------------------

%--------------------make tunnel-----------------------
sampleId=d.getGroupId('sample');
sX=d.mo.aX(sampleId);sY=d.mo.aY(sampleId);sZ=d.mo.aZ(sampleId);
dipD=0;dipA=0;radius=outerR;height=5;
columnFilter=f.run('fun/getColumnFilter.m',sX,sY,sZ,dipD,dipA,radius+B.ballR,height);
zFilter=sZ>1;
tunnelId=find(columnFilter&zFilter);
d.delElement(tunnelId);
%--------------------end make tunnel-----------------------

%-------------------insert the ringTube-------------------
tubeId=d.addElement(1,tubeObj);%add a slope boundary
d.addGroup('ringTube',tubeId);%add a new group
d.setClump('ringTube');%set the pile clump
d.moveGroup('ringTube',(max(mX)+min(mX))/2,(max(mY)+min(mY))/2,1);
d.minusGroup('sample','ringTube',0.4);%remove overlap elements from sample
innerTubeId=find(d.mo.aR==minBallR);
d.addGroup('innerTube',innerTubeId);%add a new group
%-------------------end insert the ringTube-------------------

d.addFixId('X', d.GROUP.ringTube);%fix the coordinates of elements of the pile
d.addFixId('Y', d.GROUP.ringTube);
d.addFixId('Z', d.GROUP.ringTube);
d.balanceBondedModel0();%bond the elements and balance the model
d.breakGroup();%break all connections
d.balance('Standard');%balance the model before saving

d.show();
%--------------------save data-----------------------
d.mo.setGPU('off');
d.clearData(1);
d.recordCalHour('BoxTunnelHeat2Finish');
save(['TempModel/' B.name '2.mat'],'B','d');
save(['TempModel/' B.name '2R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
%--------------------end save data-----------------------

d.calculateData();
d.show('ZDisplacement');

%note: this is an example to show thermomechanical coupling
%it is not completed, specified heat must considered in real applications
%Chun Liu, Nanjing University
%-------------------user_mxSlope3.m;
clear;
load('TempModel/BoxTunnelHeat2.mat');
B.setUIoutput();
%------------initialize model-------------------
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();%d.setModel();%reset the initial status of the model
d.resetStatus();
d.mo.isHeat=1;%calculate heat in the model
d.setStandarddT();
d.mo.isCrack=1;
%------------end initialize model-------------------
%d.show();return;

%--------------------set block force--------------------
innerTubeId=d.GROUP.innerTube;
d.mo.SET.aT=ones(d.aNum,1)*15;%initial temperature is 15 degrees
initialT=d.mo.SET.aT(1:d.mNum);%record the initial temperatures
innerTubeT=25;
d.mo.SET.aT(d.mNum+1:d.aNum)=-1000;%boundary is insulated
mdR0=zeros(d.mNum,1);%deviation of radius of active elements
%--------------------end set block force--------------------

totalCircle=10;
stepNum=20;
fName=['data/step/' B.name  num2str(B.ballR) '-' num2str(B.distriRate) 'loopNum'];
save([fName '0.mat']);%return;
gpuStatus=d.mo.setGPU('auto');
d.tic(totalCircle);%record the initial time of loop
for i=1:totalCircle
    d.mo.setGPU(gpuStatus);
    for j=1:stepNum
        %-----1. calculate the temperature difference, nTempDiff
        d.mo.SET.aT(innerTubeId)=innerTubeT;%assign temperture to innerTube group
        nRow=ones(1,size(d.mo.nBall,2));
        NBall=gather(d.mo.nBall);
        nTempDiff=d.mo.SET.aT(NBall)-d.mo.SET.aT(1:d.mNum)*nRow;%difference of water content bewteen elements and neighbors

        %-----2. heat conduction
        nThermalC=0.02;%Coefficient of thermal conductivity
        cbFilter=gather(d.mo.cFilter|d.mo.bFilter);%define the contact filter
        nBoundaryFilter=(d.mo.SET.aT(NBall)==-1000);%temperature of neighboring elements
        inslatedFilter=(~cbFilter|nBoundaryFilter);
        nTempFlow=nTempDiff.*nThermalC;%temperature various
        nTempFlow(inslatedFilter)=0;        
        mTempFlow=sum(nTempFlow,2);
        d.mo.SET.aT(1:d.mNum)=d.mo.SET.aT(1:d.mNum)+mTempFlow;
        %-----3. update element radius and properties according to temperature, i.e. aT
        expandRate=0.001;%expanded by 0.1% when temperature increased by 1 degree
        mdR1=gather(d.mo.aR(1:d.mNum).*(d.mo.SET.aT(1:d.mNum)-initialT)*expandRate);%describe how radius varies with T
        d.mo.aR(1:d.mNum)=d.aR(1:d.mNum)+mdR1;
        maxDis=gather(max(max(abs(d.mo.dis_mXYZ),[],2)+mdR1-mdR0));%transform data to CPU
        if (maxDis>0.5*d.mo.dSide)
            fs.disp(['balanceTime' num2str(d.mo.totalT) '->expand->setNearbyBall']);
            d.mo.setNearbyBall();
            mdR0(:)=mdR1;
        end
        %here, you may modify other properties of elements when temperature changed
        d.mo.balance();%calculation
        d.recordStatus();
    end
    %clear and save data
    d.clearData(1);
    d.mo.setGPU('off');
    save([fName num2str(i) '.mat']);
    d.calculateData();
    d.toc();%show the note of time
end

d.recordCalHour('BoxTunnelHeat3Finish');
d.mo.setGPU('off');
save(['TempModel/' B.name '3.mat'],'B','d');
save(['TempModel/' B.name '3R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.show('aR');