@Bjorn Gustavsson maybe like this? could you help me to set the boundary condition and initial condition? since I dont know where to put that condition. Thank you in advance
Solving moisture diffusion equation
15 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Hi everyone, I am a little bit confused about this,
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/524524/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/524529/image.png)
Now Im stuck in the 'how to set the initial and boundary condition'. and the following figure is what I want to get and the trend should be like the following figure.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/524534/image.png)
and this is my temporary code,
clear; clc;
days=1000;
tcuring=14;
var(1)=0.001;
timeinit=0.001;
counter=2;
while (timeinit<days)
timeinit=timeinit*1.1;
var(counter)=timeinit;
counter=counter+1;
end
timeh=zeros(1,length(var)+1);
timeh(1)=0; %should be changed
for m=1:length(var)
timeh(m+1)=var(m)*24;
end
timeday=timeh/24;
%% Input
c3s=58.6;
c2s=17.2;
c3a=6.97;
c4af=9.9;
Eac3s=30000; %J/mol
Eac2s=56000; %J/mol
Eac3a=45000; %J/mol
Eac4af=32000; %J/mol
mc=0.4; %input mass of cement g/cm3
mw=0.157; %input mass of water g/cm3
ms=0.658; %input mass of sand g/cm3
mg=1.129; %input mass of gravel g/cm3
rh=1; %relative humidity
%year=100;
days=36500;%*year; %days
srfarea=330; %m2/g
refarea=330; %m2/g
wc=0.5; %water to cement ratio
To=293.15; %reference temperature
R=8.314; %J/mol/K
RH1=linspace(0.005,1,length(timeday));
T=20; %celcius
Tcal=T+273.15;
rhoeff=2.04; %density of HCP at 105C g/cm3
%Nucleation and Growth
Kc3s=1.7;%Maruyama's old paper
Nc3s=0.7;%Lothen
Hc3s=1.4;
Kc2s=0.1;
Nc2s=1;%Lothen
Hc2s=1.55;
Kc3a=1;
Nc3a=0.85;
Hc3a=1.5;
Kc4af=0.37;
Nc4af=0.7;
Hc4af=1.55;
%Diffusion
Kc3sd=1.7;
Kc2sd=0.1;
Kc3ad=0.4;
Kc4afd=0.1;
%Shell formation
Kc3ss=1.7;
Nc3ss=3.3;
Kc2ss=0.1;
Nc2ss=4.1;
Kc3as=2;
Nc3as=2.9;
Kc4afs=0.1;
Nc4afs=3.8;
ccem=0.78; %by lothenbach et al.(2008) and (maruyama et al 2006)
chyd=1.25; %J/gK
cwat=4.1852;
cagg=0.9; %ranging from 0.7-1
Qmax=411.15;%J/g
%% pre-allocating
toalfa=zeros(1,size(var,2)+1);
drc3s=zeros(1,size(var,2)+1);
drc2s=zeros(1,size(var,2)+1);
drc3a=zeros(1,size(var,2)+1);
drc4af=zeros(1,size(var,2)+1);
rc3sn=zeros(1,length(var)+1);
rc3ss=zeros(1,length(var)+1);
rc3sd=zeros(1,length(var)+1);
rc2sn=zeros(1,length(var)+1);
rc2ss=zeros(1,length(var)+1);
rc2sd=zeros(1,length(var)+1);
rc3an=zeros(1,length(var)+1);
rc3as=zeros(1,length(var)+1);
rc3ad=zeros(1,length(var)+1);
rc4afn=zeros(1,length(var)+1);
rc4afs=zeros(1,length(var)+1);
rc4afd=zeros(1,length(var)+1);
%% Parrot-Killoh model
for i=1:length(var) %boundary condition
if (i==1)
rc3sn(1)=0;
rc3ss(1)=0;
rc3sd(1)=0;
drc3s(1)=10^-15;
rc2sn(1)=0;
rc2ss(1)=0;
rc2sd(1)=0;
drc2s(1)=10^-15;
rc3an(1)=0;
rc3as(1)=0;
rc3ad(1)=0;
drc3a(1)=10^-15;
rc4afn(1)=0;
rc4afs(1)=0;
rc4afd(1)=0;
drc4af(1)=10^-15;
else
rc3sn(1)=0;
rc3ss(1)=0;
rc3sd(1)=0;
drc3s(1)=10^-15;
rc2sn(1)=0;
rc2ss(1)=0;
rc2sd(1)=0;
drc2s(1)=10^-15;
rc3an(1)=0;
rc3as(1)=0;
rc3ad(1)=0;
drc3a(1)=10^-15;
rc4afn(1)=0;
rc4afs(1)=0;
rc4afd(1)=0;
drc4af(1)=10^-15;
end
end
for j=1:length(var)
toalfa(j)=(drc3s(j)*c3s+drc2s(j)*c2s+drc3a(j)*c3a+drc4af(j)*c4af)/(c3s+c2s+c3a+c4af);
%C3S
if (wc*Hc3s>toalfa(j))
rc3sn(j+1)=Kc3s/Nc3s*(1-drc3s(j))*(-log(1-drc3s(j)))^(1-Nc3s)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc3sn(j+1)=Kc3s/Nc3s*(1-drc3s(j))*(-log(1-drc3s(j)))^(1-Nc3s)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*(1+3.333*Hc3s*wc-3.333*toalfa(j))^4*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
if (wc*Hc3s>toalfa(j))
rc3ss(j+1)=Kc3ss*(1-drc3s(j))^Nc3ss*(timeh(j+1)-timeh(j))/24*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc3ss(j+1)=Kc3ss*(1-drc3s(j))^Nc3ss*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc3s*wc-3.333*toalfa(j))^4*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
if (wc*Hc3s>toalfa(j))
rc3sd(j+1)=Kc3sd*(1-drc3s(j))^(2/3)/(1-(1-drc3s(j))^(1/3))*(timeh(j+1)-timeh(j))/24*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc3sd(j+1)=Kc3sd*(1-drc3s(j))^(2/3)/(1-(1-drc3s(j))^(1/3))*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc3s*wc-3.333*toalfa(j))^4*exp(Eac3s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
%C2S
if (wc*Hc2s>toalfa(j))
rc2sn(j+1)=Kc2s/Nc2s*(1-drc2s(j))*(-log(1-drc2s(j)))^(1-Nc2s)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc2sn(j+1)=Kc2s/Nc2s*(1-drc2s(j))*(-log(1-drc2s(j)))^(1-Nc2s)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*(1+3.333*Hc2s*wc-3.333*toalfa(j))^4*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
if (wc*Hc2s>toalfa(j))
rc2ss(j+1)=Kc2ss*(1-drc2s(j))^Nc2ss*(timeh(j+1)-timeh(j))/24*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc2ss(j+1)=Kc2ss*(1-drc2s(j))^Nc2ss*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc2s*wc-3.333*toalfa(j))^4*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
if (wc*Hc2s>toalfa(j))
rc2sd(j+1)=Kc2sd*(1-drc2s(j))^(2/3)/(1-(1-drc2s(j))^(1/3))*(timeh(j+1)-timeh(j))/24*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc2sd(j+1)=Kc2sd*(1-drc2s(j))^(2/3)/(1-(1-drc2s(j))^(1/3))*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc2s*wc-3.333*toalfa(j))^4*exp(Eac2s/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
%C3A
if (wc*Hc3a>toalfa(j))
rc3an(j+1)=Kc3a/Nc3a*(1-drc3a(j))*(-log(1-drc3a(j)))^(1-Nc3a)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc3an(j+1)=Kc3a/Nc3a*(1-drc3a(j))*(-log(1-drc3a(j)))^(1-Nc3a)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*(1+3.333*Hc3a*wc-3.333*toalfa(j))^4*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
if (wc*Hc3a>toalfa(j))
rc3as(j+1)=Kc3as*(1-drc3a(j))^Nc3as*(timeh(j+1)-timeh(j))/24*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc3as(j+1)=Kc3as*(1-drc3a(j))^Nc3as*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc3a*wc-3.333*toalfa(j))^4*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
if (wc*Hc3a>toalfa(j))
rc3ad(j+1)=Kc3ad*(1-drc3a(j))^(2/3)/(1-(1-drc3a(j))^(1/3))*(timeh(j+1)-timeh(j))/24*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc3ad(j+1)=Kc3ad*(1-drc3a(j))^(2/3)/(1-(1-drc3a(j))^(1/3))*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc3a*wc-3.333*toalfa(j))^4*exp(Eac3a/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
%C4AF
if (wc*Hc4af>toalfa(j))
rc4afn(j+1)=Kc4af/Nc4af*(1-drc4af(j))*(-log(1-drc4af(j)))^(1-Nc4af)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc4afn(j+1)=Kc4af/Nc4af*(1-drc4af(j))*(-log(1-drc4af(j)))^(1-Nc4af)*(timeh(j+1)-timeh(j))/24/refarea*srfarea*(1+3.333*Hc4af*wc-3.333*toalfa(j))^4*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
if (wc*Hc4af>toalfa(j))
rc4afs(j+1)=Kc4afs*(1-drc4af(j))^Nc4afs*(timeh(j+1)-timeh(j))/24*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc4afs(j+1)=Kc4afs*(1-drc4af(j))^Nc4afs*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc4af*wc-3.333*toalfa(j))^4*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
if (wc*Hc4af>toalfa(j))
rc4afd(j+1)=Kc4afd*(1-drc4af(j))^(2/3)/(1-(1-drc4af(j))^(1/3))*(timeh(j+1)-timeh(j))/24*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
else
rc4afd(j+1)=Kc4afd*(1-drc4af(j))^(2/3)/(1-(1-drc4af(j))^(1/3))*(timeh(j+1)-timeh(j))/24*(1+3.333*Hc4af*wc-3.333*toalfa(j))^4*exp(Eac4af/R*(1/To-1/Tcal))*((rh-0.55)/0.45)^4;
end
%C3S
if (rc3sn(j+1)<rc3ss(j+1))
if (rc3sn(j+1)<rc3sd(j+1))
drc3s(j+1)=drc3s(j)+rc3sn(j+1);
else
drc3s(j+1)=drc3s(j)+rc3sd(j+1);
end
else
if (rc3ss(j+1)<rc3sd(j+1))
drc3s(j+1)=drc3s(j)+rc3ss(j+1);
else
drc3s(j+1)=drc3s(j)+rc3sd(j+1);
end
end
%C2S
if (rc2sn(j+1)<rc2ss(j+1))
if (rc2sn(j+1)<rc2sd(j+1))
drc2s(j+1)=drc2s(j)+rc2sn(j+1);
else
drc2s(j+1)=drc2s(j)+rc2sd(j+1);
end
else
if (rc2ss(j+1)<rc2sd(j+1))
drc2s(j+1)=drc2s(j)+rc2ss(j+1);
else
drc2s(j+1)=drc2s(j)+rc2sd(j+1);
end
end
%C3A
if (rc3an(j+1)<rc3as(j+1))
if (rc3an(j+1)<rc3ad(j+1))
drc3a(j+1)=drc3a(j)+rc3an(j+1);
else
drc3a(j+1)=drc3a(j)+rc3ad(j+1);
end
else
if (rc3as(j+1)<rc3ad(j+1))
drc3a(j+1)=drc3a(j)+rc3as(j+1);
else
drc3a(j+1)=drc3a(j)+rc3ad(j+1);
end
end
%C4AF
if (rc4afn(j+1)<rc4afs(j+1))
if (rc4afn(j+1)<rc4afd(j+1))
drc4af(j+1)=drc4af(j)+rc4afn(j+1);
else
drc4af(j+1)=drc4af(j)+rc4afd(j+1);
end
else
if (rc4afs(j+1)<rc4afd(j+1))
drc4af(j+1)=drc4af(j)+rc4afs(j+1);
else
drc4af(j+1)=drc4af(j)+rc4afd(j+1);
end
end
%Rate of hydration
rc3s(j+1)=(drc3s(j+1)-drc3s(j))/((timeh(j+1)-timeh(j))/24);
rc2s(j+1)=(drc2s(j+1)-drc2s(j))/((timeh(j+1)-timeh(j))/24);
rc3a(j+1)=(drc3a(j+1)-drc3a(j))/((timeh(j+1)-timeh(j))/24);
rc4af(j+1)=(drc4af(j+1)-drc4af(j))/((timeh(j+1)-timeh(j))/24);
drc3s(2)=10^-15;
drc2s(2)=10^-15;
drc3a(2)=10^-15;
drc4af(2)=10^-15;
drc3s(2)=10^-15;
drc2s(2)=10^-15;
drc3a(2)=10^-15;
drc4af(2)=10^-15;
toalfa(length(var)+1)=(drc3s(length(var)+1)*c3s+drc2s(length(var)+1)*c2s+drc3a(length(var)+1)*c3a+drc4af(length(var)+1)*c4af)/(c3s+c2s+c3a+c4af);
end
for j=1:length(var)
rtoalfa(j+1)=(toalfa(j+1)-toalfa(j))/((timeh(j+1)-timeh(j))/24);
end
toalfa1=max(toalfa);
%% Heat related properties (assume the volume equals to 1 cm3)
wuc=(1-toalfa)*mc/(mw+mc+ms+mg);
wf0=mw/(mw+mc+ms+mg);
wf01=mw;
tetahpc=0.23; %based on arosa dabarera(2017)as a minimum water to binder ratio,The value for θhp,c was previously used in many studies for representingchemically bound water of cement (Powers, 1960; Neville, 1995; Lam, Wong, &Poon, 2000; Saengsoy, 2003)
wwhp=tetahpc*mc*toalfa/(mw+mc+ms+mg); %chemically bound water content
wwhp1=tetahpc*mc*toalfa; %g/cm3 blm dalam bentuk rate
tetagelc=0.0126+(0.0026/(-1.009+exp(0.1414*wc)));
wwgel=mc*toalfa*tetagelc/(mw+mc+ms+mg); %gel water content
wwgel1=mc*toalfa*tetagelc;
wfw=wf0-wwhp-wwgel; %weight of evaporable water (no unit)
wfw1=wf01-wwhp1-wwgel1;
wagg=(ms+mg)/(mw+mc+ms+mg);
whp=1-(wfw+wuc+wagg);
Ccon=wuc*ccem+whp*chyd+wfw*cwat+wagg*cagg; %J/gK
averageQt=toalfa*Qmax;%J/gK
averagedeltat=(averageQt*((c3s+c3a+c2s+c4af)/100)*mc)./((mc+ms+mg+mw)*(Ccon));
vct=0.9;%input(prompt);
if tcuring>5
if wc<=0.3
vmbet=(0.068-0.22/tcuring)*(0.85+0.45*0.3)*vct;%2.72*10^-4;
elseif wc>=0.7
vmbet=(0.068-0.22/tcuring)*(0.85+0.45*0.7)*vct;
else
vmbet=(0.068-0.22/tcuring)*(0.85+0.45*wc)*vct;
end
else
if wc<=0.3
vmbet=(0.068-0.22/5)*(0.85+0.45*0.3)*vct;%2.72*10^-4;
elseif wc>=0.7
vmbet=(0.068-0.22/5)*(0.85+0.45*0.7)*vct;
else
vmbet=(0.068-0.22/5)*(0.85+0.45*wc)*vct;
end
end
Cbet=exp(855/Tcal);
%prompt1='Cement type? type 1=1.1; type 2=1; type 3=1.15; type 4=1.5?';
type=1.1;%input(prompt1);
if tcuring>5
if wc<=0.3
n=(2.5+15/tcuring)*(0.33+2.2*0.3)*type;
elseif wc>=0.7
n=(2.5+15/tcuring)*(0.33+2.2*0.7)*type;
else
n=(2.5+15/tcuring)*(0.33+2.2*wc)*type;
end
else
if wc<=0.3
n=(2.5+15/5)*(0.33+2.2*0.3)*type;
elseif wc>=0.7
n=(2.5+15/5)*(0.33+2.2*0.7)*type;
else
n=(2.5+15/5)*(0.33+2.2*wc)*type;
end
end
kbet=(1-1/n)*(Cbet-1)/(Cbet-1);%0.83;
cm=100; % Prediction length
h=10; % the interval of distance (cm)
M=(cm)*(1/h); % Cut-off line of x-axis (days)
E=0.0001;
for mm=1:length(RH1)
if (RH1(mm)<=0.98)
wrht(mm)=(vmbet*Cbet*kbet*RH1(mm)/((1-kbet*RH1(mm))*(1+(Cbet-1)*kbet*RH1(mm))));
else
wrht(mm)=(vmbet*Cbet*kbet*RH1(mm)/((1-kbet*RH1(mm))*(1+(Cbet-1)*kbet*RH1(mm))));
end
end
for mm=1:length(RH1)
if (RH1(mm)<=0.98)
dwrht(mm)=(Cbet*RH1(mm)*kbet^2*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)^2) - (Cbet*kbet*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)) + (Cbet*RH1(mm)*kbet^2*vmbet*(Cbet - 1))/((RH1(mm)*kbet*(Cbet - 1) + 1)^2*(RH1(mm)*kbet - 1));
else
dwrht(mm)=(Cbet*RH1(mm)*kbet^2*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)^2) - (Cbet*kbet*vmbet)/((RH1(mm)*kbet*(Cbet - 1) + 1)*(RH1(mm)*kbet - 1)) + (Cbet*RH1(mm)*kbet^2*vmbet*(Cbet - 1))/((RH1(mm)*kbet*(Cbet - 1) + 1)^2*(RH1(mm)*kbet - 1));%50*w0293 - 50*w98ad;
end
end
ah=1.05-3.8*wc+3.56*(wc)^2;
bh=-14.4+50.4*wc-41.8*(wc)^2;
gamh=31.3-136*wc+162*(wc)^2;
Dh=ah+bh.*(1-2.^(-10.^(gamh*(RH1-1))));
RH=0.766; %relative humidity
RHH=zeros(M,length(var));
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% boundary condition %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for J=1:length(var)+1
RHH(:,1)=1;
if (J==1)
RHH(1,J)=0;
else
RHH(1,J)=0;
end
RHH(M+1,J)=0; %it's a boundary condition
end
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Crank & Nicolson method %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
r=1/(2*h^2);
for J=1:length(var)
for I=2:M
RHH(I,J+1)=RHH(I,J);
% %RHH(I+1,J)=RHH(I,J);
end
for L=1:500
LL=L; %number of times of calculation
for I=2:M
RHH(I,J+1)=(((r*((timeh(J+1)-timeh(J))/(24)))*(1/dwrht(J))*(Dh(J+1)*(RHH(I+1,J+1)-2*RHH(I,J+1)+RHH(I-1,J+1))+Dh(J)*(RHH(I+1,J)-2*RHH(I,J)+RHH(I-1,J))+0))+RHH(I,J));
end
S=2*Dh(J)*((timeh(J+1)-timeh(J))/(24))/h^2;
end
end
depth=0:h:cm;
% figure (1)
% plot(RH1,dwrht)
% grid on
% figure (2)
% plot(RH1,Dh)
% xlabel('RH')
% ylabel('Dh')
% grid on
%
figure (3)
plot(depth,RHH(:,58),depth,RHH(:,75),depth,RHH(:,123),depth,RHH(:,140),depth,RHH(:,147))
legend('5 hours','24 hours','100 days','500 days','1000 days')
grid on
xlabel('Depth (cm)')
ylabel('RH')
figure (4)
plot(timeday,RHH(2,:),timeday,RHH(3,:),timeday,RHH(4,:),timeday,RHH(5,:),timeday,RHH(6,:),timeday,RHH(7,:))
grid on
legend('10 cm','20 cm','30 cm','40 cm','50 cm')
xlabel('Days')
ylabel('RH')
Hope someone could help me to understand clearly how to perform this equation. Thankyou very much!
10 Kommentare
Bjorn Gustavsson
am 23 Feb. 2021
You "simply" have to make sure that you never take longer steps in time than:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/529289/image.png)
If your algorithm suggest taking a longer time-step you have to let the FTCS-scheme do that in sub-steps, and handle everything else accordingly. You don't need to save every step of the solution if that becomes too large, but if dt becomes too long you might get growing oscillations...
Antworten (0)
Siehe auch
Kategorien
Mehr zu Thermal Analysis finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!