跳转至

ParticleMigration

user_ParticleMigration1.m
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
clear;clc
fs.randSeed(2);
bBox=[0,0,0;1,0.5,0.5];
WLH=bBox(2,:)-bBox(1,:);
boxV=prod(WLH);

grainV=0.25*0.7;%packing rate: 0.7
grainR=0.5*mfs.getGradationDiameter([0.01,0.02,0.02;0.05,0.10,0.98],grainV);
%grainR=0.5*mfs.getGradationDiameter([0.005,0.010,0.015;0.025,0.050,0.985],grainV);

fineId=find((2*grainR)<0.04);

tic;sampleObj=pfs3d.genAssembly(grainR,bBox);toc
%fs.showObj(sampleObj);
%--------------------------
B=obj_Box();
B.name='ParticleMigration';
B.ballR=0.02;%mean(grainR.^3).^(1/3);
B.isSample=0;
B.sampleW=WLH(1);
B.sampleL=WLH(2);
B.sampleH=WLH(3);

B.buildInitialModel();
d=B.d;
%--------------------------
sampleId=d.addElement(1,sampleObj);
d.addGroup('sample',sampleId);
d.delElement(d.GROUP.topPlaten);
d.moveGroup('topB',0,0,-2*B.ballR);
d.show('aR')

moNum=2;
rRanges=fs.getAutoRanges(d.mo.aR,moNum);
aMoId=fs.getMoId(d.mo.aR,rRanges);
aMoId(end)=0;
d.mo.setChildModel(aMoId,rRanges);
d.addGroup('Coarse',d.mo.moCell{1}.parent_mId);
d.addGroup('Fine',d.mo.moCell{2}.parent_mId);
%---------------------------
%d.mo.mGZ(:)=0;
n=20;
for ii=1:n
   %d.mo.afterBalance='';
   d.connectGroup();
   d.deleteConnection('boundary');
   %d.mo.mVX(:)=0;
   %d.mo.mVY(:)=0;
   %d.mo.mVZ(:)=0;
   fs.limitFrame(d,d.GROUP.sample,0,1,0,0.5,0,0.5);
   d.balance('Standard',1/n);
end
d.show('-aR','EnergyCurve','ForceCurve','HeatCurve');
%histogram(d.mo.aZ(d.GROUP.Fine),linspace(bBox(1,3),bBox(2,3),21));

%---------save the data
d.mo.setGPU('off'); 
d.clearData(1);
d.recordCalHour('Step1Finish');
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();
user_ParticleMigration2.m
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
clear
load('TempModel/ParticleMigration1.mat');
B.setUIoutput();%set output of message
d=B.d;
d.calculateData();
d.mo.setGPU('off'); 
d.getModel();%get xyz from d.mo
d.resetStatus();
%---------------------
d.mo.aMUp(:)=0;
d.breakGroup();

coarseId=d.GROUP.Coarse;
fineId=d.GROUP.Fine;
%TR=getMesh(d,coarseId);
d.addFixId('X',coarseId);
d.addFixId('Y',coarseId);
d.addFixId('Z',coarseId);

% rigB=d.GROUP.rigB;
% rigBXYZ=[d.mo.aX(rigB),d.mo.aY(rigB),d.mo.aZ(rigB)];
% %F=vecnorm(rigBXYZ(:,[2,3])-[0.4,0.4],2,2)<0.1;
% fh=@(YZ,yz0,r0)vecnorm(YZ-yz0,2,2)<r0;
% F=fh(rigBXYZ(:,[2,3]),[0.2,0.2],0.05)|fh(rigBXYZ(:,[2,3]),[0.2,0.6],0.05);
% F=F|fh(rigBXYZ(:,[2,3]),[0.6,0.2],0.05)|fh(rigBXYZ(:,[2,3]),[0.6,0.6],0.05);
% F=F|fh(rigBXYZ(:,[2,3]),[0.4,0.4],0.05);
% 
% d.mo.aKN(rigB(F))=0;
% d.mo.aKS(rigB(F))=0;
% d.mo.setKNKS();
%d.showB=2;d.show('aKN');
d.moveGroup('rigB',0.1,0,0);

%p=build2pore(d,'Coarse');
%p.show;

w1=0.70;w2=0.20;
coarseWC=w2+(w1-w2)*(1-rescale(d.mo.aX(coarseId)));%

showFilter=abs(d.mo.aX-0.2)<0.05|abs(d.mo.aX-0.4)<0.05|abs(d.mo.aX-0.6)<0.05|abs(d.mo.aX-0.8)<0.05;
showFilter(end)=false;
showFilter(d.GROUP.Fine)=true;

saveDir=['data\',B.name,'_multifield'];
mkdir(saveDir);
save([saveDir,'\loopNum',num2str(0),'.mat']);
totalCircle=50;
steps=100;
d.tic(totalCircle);
for ii=1:totalCircle
    for jj=1:steps
        %1.场的演化
        %
        %2.建立粗颗粒场 & 梯度计算
        aPos=[d.mo.aX,d.mo.aY,d.mo.aZ];
        if ~exist('TR','var')
            TR=pfs3d.getMesh(aPos(coarseId,:),d.mo.aR(coarseId));
        end
        triGrad=pfs3d.calGrad(TR.ConnectivityList,aPos(coarseId,:),coarseWC);
        %3.映射细颗粒 & 受力计算
        fineLoc=TR.pointLocation(aPos(fineId,:));
        f=isnan(fineLoc);%out of boundingbox
        fineLoc(f)=1;%avoid index error
        fGrad=triGrad(fineLoc,:);
        fGrad(f,:)=0;
        kF=-2e4*median(d.mo.aR(coarseId))*9.8*d.mo.mM(fineId);%based on gravity
        d.mo.mGX(fineId)=fGrad(:,1).*kF;
        d.mo.mGY(fineId)=fGrad(:,2).*kF;
        d.mo.mGZ(fineId)=fGrad(:,3).*kF;
        %4.颗粒运动
        d.balance(1);
        %d.recordStatus();
    end
    %5.额外数据记录
    d.figureNumber=1;
    d.showFilter('Filter',showFilter);
    d.show('Displacement');
    frames(ii)=getframe();
    save([saveDir,'\loopNum',num2str(ii),'.mat']);
    d.toc();
end
fs.movie2gif('1.gif',frames,0.2);
%movie('1.gif');

% n=20;
% bins=linspace(0,1,n+1);
% N=histcounts(d.mo.aX(d.GROUP.Fine),bins);
% N=N/length(d.GROUP.Fine);
% plot(movmean(bins,2,'Endpoints','discard'),N,'-^')
% ylim([0.4,1.6]*1/n);
% yline(1/n,'--')

%---------save the data
d.mo.setGPU('off'); 
d.clearData(1);
d.recordCalHour('Step2Finish');
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']);
d.calculateData();