how can I make this vode work?

3 views (last 30 days)
Richard Mogou
Richard Mogou on 8 Dec 2022
Answered: Image Analyst on 8 Dec 2022
clc;
clear all;
close all;
global phie e rhocp L dt tf N alpha1 dr drp alpha r
L=0.01;
dt=0.01; tf=50;
N=21; %M=round(tf/dt);
dr=L/(N-1); drp=2*dr;
alpha=0.5;
r=0.00;
r_out=0.01;
M=21; %nombre de pas entre r_in et r_out
dr=(L-r)/M;
%nr=M+1; %nombre total de noeuds radial
nphi=30;% nombre de pas d'angles
delta_phi=2*pi/nphi; %pas d'angle
t=500; %temps
nt=0.1; % pas dans l'espace temps
dt=t/nt;
T_in=28;
T_out=28;
err=10^(-3); Tavant=28; Tarriere=28;Ta=28;
%********************************** lameff en W/m.K ; I en Ampere ;
%X12=0.02887; c=1495; K=0.9898;
%rhocp en J/m.K X12=4.9054; c=8.4652; K=0.0082;
X12=2.887; c=14.95; K=0.009898;
% inputting the initial values******************** alfa est un coefficient entre 0 et 1 qui dépend de la tortuosité et de la porosité du matériau
%X12=4.7419; c=7.1141; K=0.0082; %X12,c et K sont les coefficients du modèle de GAB
alpha1=0.5;alpha=0.5; R=8.314; Mv=0.018; rhos=1000; rhoe=382; cps=1240; cpl=4181;% rhos est la masse volumique de bois de tali Ngono 2019 (changer la valeur en fonction de la teneur en eau initiale)
lambda=0.169; rhocp=8.08e5; p=1013.25; pva=1013.25; Lv=2504;% rhocp et lameff correspondent aux valeurs pour le Tali en transversale (Ngono et al 2019)
%********************************** rhos et cps sont respectivement la masse volumiq et la capacité thermik du materiau sec
%********************************** p est la Pression gazeuse totale %lameff la cond thermik effective du matériau
T=ones(N+1);Pvs=ones(N+1);HR=ones(N+1);Da=ones(N+1);pv=ones(N+1);xv=ones(N+1); X12=ones(N+1);c=ones(N+1);K=ones(N+1);
X=ones(N+1);dHR=zeros(N+1);dPvs=zeros(N+1);DVX=zeros(N+1);DVT=zeros(N+1);lambdae=ones(N+1);gamma=ones(N+1); HRa=ones(N+1);
T1=ones(M+1,N+1);X1=ones(M+1,N+1);DVX1=zeros(M+1,N+1);DVT1=zeros(M+1,N+1); lambdae1=ones(M+1,N+1); gamma1=ones(M+1,N+1);
%Dm=zeros(N+1);Dm1=zeros(M+1,N+1); Delta=zeros(N+1); Delta1=zeros(M+1,N+1);
%Ahr=ones(N+1);Bhr=ones(N+1);Chr=ones(N+1); c=ones(N+1);
re=3949.6159; %Rapport resistance electrique/section sonde (ohm/m2)
I=0.252;
phie=(re*I^2)/2;
rces=585.8002; %rces est la capacité thermique de la sonde
e=L; % epaisseur du matériau (m)
lam=lambda; %conductivité isolant (w/m.K)
roc=rhocp; %capacité isolant (J/m.K);
%lecture des paramètres
prompt={'epaisseur materiau (m)','Intensité (A)','conductivité thermique du materiau'};
def={num2str(e) num2str(I) num2str(lam)};
titre='Entrée des conditions expérimentales';
lineNo=[1 60];
answer=inputdlg(prompt,titre,lineNo,def);
Error using matlab.internal.lang.capability.Capability.require
Support for Java user interfaces is required, which is not available on this platform.

Error in inputdlg (line 55)
Capability.require(Capability.Swing);
stranswer1=deal(answer{1,:});
stranswer2=deal(answer{2,:});
stranswer3=deal(answer{3,:});
r0=str2double(stranswer1); %str2double plus rapide que str2num mais valide seulement pour les scalaires
r1=str2double(stranswer2);
I=str2double(stranswer3);
%fin lecture des paramètres
for i=2:1:N
T(i)=(Tavant+Tarriere)/2; Pvs(i)=exp(25.555-(5220/(T(i)+273))); %HR(i)=(HR(1)+HR(N+1))/2;
%teneur en eau initiale calculée par le modèle de GAB
X(i)=23;
X12(i)=0.0035*T(i)^2-0.1045*T(i)+4.7088; % X12(i)=-1.3786e-5*T(i)^2+0.0004558*T(i)+0.069;%
c(i)=-0.0001*T(i)^2-0.403*T(i)+19.312; % c(i)=0.0064*T(i)^2-0.58*T(i)+22; % % le modèle de X12 k et c en fonction de la temperature est tiré de la thèse de SAED RAJI
K(i)=3.7143e-5*T(i)^2-0.00077*T(i)+0.747; %K(i)=8e-6*T(i)^2-0.0007*T(i)+0.0214;
Ahr(i)=c(i)-1; Bhr(i)=(X12(i)*c(i))/X(i)-c(i)+2; Chr(i)=1;
%HR(i)=70;
HR(i)=(1/K(i))*(-Bhr(i)+sqrt(Bhr(i)^2-4*Ahr(i)*Chr(i))/(2*Ahr(i)));
pv(i)=HR(i)*Pvs(i); %pvs la Pression de vapeur d’eau saturante à la température T
Da(i)=2.17e-5*((T(i)+273)/273)^1.88; xv(i)=0.622*(pv(i)/(p)); %xv est la Fraction massique de vapeur d’eau
%X(i)=(X12(i)*c(i)*K(i)*HR(i))/(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i)); %modèle de GAB
X(i)=(X12(i)*c(i)*K(i)*HR(i))/((c(i)-2)*K(i)*HR(i)-1-(c(i)-1)*K(i)^2*HR(i)^2);
dHR(i)=(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))/(X12(i)*c(i)*K(i)+X(i)*K(i)*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))-X(i)*(1-K(i)*HR(i))*(c(i)*K(i)-K(i)));
dHR(i)=(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))/(X12(i)*c(i)*K(i)+X(i)*K(i)*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))-X(i)*(1-K(i)*HR(i))*(c(i)*K(i)-K(i)));
dPvs(i)=(5220/(T(i)+273)^2)*Pvs(i); %dHR(i)=abs(dHR(i));
DVX(i)=(alpha/rhos*(1-xv(i)))*((Da(i)*Mv)/(R*(T(i)+273)))*Pvs(i)*dHR(i);
DVT(i)=(alpha1/rhos*(1-xv(i)))*((Da(i)*Mv)/(R*(T(i)+273)))*HR(i)*dPvs(i); %rhocp(i)=rhos*(cps+X(i)*cpl);
lambdae(i)=(lambda+rhos*Lv*DVT(i)); gamma(i)=lambdae(i)/(rhocp); % lamapp est la conductivité thermique apparente; app la diffusivité apparente du materiau
%Dm(i)=DVX(i)+DVT(i);Delta(i)=DVT(i)/DVX(i);
%gamma(i)=lambdae(i)/(rhocp);
%X(i)=0.4061; %stockage des variables pour t=0
T1(1,i)=T(i)-Ta;
DVX1(1,i)=DVX(i); lambdae1(1,i)=lambdae(i);
DVT1(1,i)=(DVT(i)); gamma1(1,i)=gamma(i);
X1(1,i)=X(i); % Dm1(1,i)=Dm(i);
%Delta1(1,i)=Delta(i);
end
% Boundaries conditions
for i=1
T(i)=Tavant; X(i)=X(2);
end
T1(1,1)=T(1)-Ta; X1(1,1)=X(1);
for i=N+1
T(i)=Tarriere; X(i)=X(N);
end
T1(1,N+1)=T(N+1)-Ta; X1(1,N+1)=X(N+1);
%loop for time
% rhos=1044;Ta=28;
time=zeros(M,1);T_new = zeros(N+1,1); X_new = zeros(N+1,1);%New Temperature calcul.
DVX_new=zeros(N+1,1); DVT_new=zeros(N+1,1); lambdae_new=zeros(N+1,1); gamma_new=zeros(N+1,1);
%Dm_new=zeros(N+1,1); Delta_new=zeros(N+1,1)
for t=1:1:M % loop start for time
%New Temperature calcul.
for i=2:N % loop starts for space
%ct=((lambda+rhos*Lv*DVT(i+1))/rhocp)*(dt/(dr*dr)); dtt=(gamma(i))*(dt/(dr*dr));
%as=(dt/(dr*dr))*(rhos*Lv/rhocp)*DVX(i+1); bs=(dt/(dr*dr))*(rhos*Lv/rhocp)*DVX(i);
%T_new(i)=T(i,t)+ct*T(i+1,t)-(ct+dtt)*T(i,t)+dtt*T(i-1,t)-(rhos*Lv/rhocp)*(as+bs)*X(i,t)+(rhos*Lv/rhocp)*bs*X(i-1,t)+(rhos*Lv/rhocp)*as*X(i+1,t);
%T_new(i)=T(i,t)-(dtt+ct)*(2+(1/i))*T(i,t)+dtt*T(i-1,t)+(1+(1/i))*ct*T(i+1,t)-(as+bs)*(2+(1/i))*X(i,t)+as*(1+(1/i))*X(i+1,t)+bs*X(i-1,t);
T_new(i)=T(i,t)+a*((dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i+1,t)-a*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i,t)+a*(dt/(dr*dr))*T(i-1,t)+b*((dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i+1,t)-b*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i,t)+b*(dt/(dr*dr))*X(i-1,t);
end
%for the first node
T_new(1)=(phie/lambdae(2))*(dr)+(lambdae(2)/lambdae(2))*T_new(2);% +(Lv*rhos/lambdae(1))*(DVX(2)*X(2)-DVX(1)*X(1));
%for the last node
T_new(N+1)=Ta;
T(:,t+1)=T_new; % temperature update
%New water content calcul.
for i=3:N-1
at=DVX(i+1); bt=DVX(i-1); aw=DVT(i+1)*(dt/(drp*drp)); bw=DVT(i-1)*(dt/(drp*drp));
%X_new(i)=X(i,t)+as*X(i+1,t)-(as+bs)*X(i,t)+bs*X(i-1,t)+at*T(i+1,t)-(at+bt)*T(i,t)+bt*T(i-1,t);
%X_new(i)=X(i,t)+at*(1+(1/i))*X(i+1,t)-(at+bt)*(2+(1/i))*X(i,t)+bt*X(i-1,t)+aw*(1+(1/i))*T(i+1,t)-(as+bw)*(2+(1/i))*T(i,t)+bw*T(i-1,t);
X_new(i)=X(i,t)+aw*((dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i+1,t)-(aw+bw)*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i,t)+bw*(dt/(dr*dr))*T(i-1,t)+at*((dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i+1,t)-(at+bt)*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i,t)+bt*(dt/(dr*dr))*X(i-1,t);
end
%for the first node
% Pvs(1)=exp(25.555-(5220/(T_new(1)+273))); HRa(1)=100*pva/Pvs(1);
% X_new(1)=HRa(1)/T_new(1);
X_new(2)=X(2,t)+DVX(3)*(dt/(drp*drp))*(1+(1/3))*X(3,t)-DVX(2)*(dt/(drp*drp))*(2+(1/2))*X(2,t)+DVX(1)*(dt/(drp*drp))*X(1,t)+DVT(3)*(dt/(drp*drp))*(1+(1/3))*T(3,t)-DVT(2)*(dt/(drp*drp))*(2+(1/2))*T(2,t)+DVT(1)*(dt/(drp*drp))*T(1,t);
X_new(1)=X_new(2);
%for the last node
X_new(N)=X(N,t)+DVX(N+1)*(dt/(drp*drp))*(1+(1/N+1))*X(N+1,t)-DVX(N)*(dt/(drp*drp))*(2+(1/N))*X(N,t)+DVX(N-1)*(dt/(drp*drp))*X(N-1,t)+DVT(N+1)*(dt/(drp*drp))*(1+(1/N+1))*T(N+1,t)-DVT(N)*(dt/(drp*drp))*(2+(1/N))*T(N,t)+DVT(N-1)*(dt/(drp*drp))*T(N-1,t);
X_new(N+1)=X_new(N);
X(:,t+1)=X_new; %update of water content
% calculate the new diffusion coefficients
for i=1:1:N+1
Pvs(i)=exp(25.555-(5220/(T_new(i)+273)));%HR(i)=0.3;
X12(i)=-1.3786e-5*T_new(i)^2+0.0004558*T_new(i)+0.069;
c(i)=0.0064*T_new(i)^2-0.58*T_new(i)+22;
K(i)=3.7143e-5*T_new(i)^2-0.00077*T_new(i)+0.747;
Ahr(i)=c(i)-1; Bhr(i)=(X12(i)*c(i))/X_new(i)-c(i)+2; Chr(i)=1;
HR(i)=(1/K(i))*(-Bhr(i)+sqrt(Bhr(i)^2-4*Ahr(i)*Chr(i))/(2*Ahr(i)));
Da(i)=2.17e-5*((T_new(i)+273)/273)^1.88; pv(i)=HR(i)*Pvs(i); xv(i)=0.622*(pv(i)/(p));
%X(i)=(X12*c*K*HR(i))/((c-2)*K*HR(i)-1-(c-1)*K^2*HR(i)^2); %On utilise plus ceci car le regime devient dynamique
%dHR(i)=(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))/(X12(i)*c(i)*K(i)+X(i)*K(i)*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))-X(i)*(1-K(i)*HR(i))*(c(i)*K(i)-K(i)));
dHR(i)=(1-K(i)*HR(i))*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))/(X12(i)*c(i)*K(i)+X_new(i)*K(i)*(1+c(i)*K(i)*HR(i)-K(i)*HR(i))-X_new(i)*(1-K(i)*HR(i))*(c(i)*K(i)-K(i)));
dPvs(i)=(5220/(T_new(i)+273)^2)*Pvs(i); %dHR(i)=abs(dHR(i));%rhocp(i)=rhos*(cps+X(i)*cpl);
DVX_new(i)=(alpha/rhos*(1-xv(i)))*((Da(i)*Mv)/(R*(T_new(i)+273)))*Pvs(i)*dHR(i);
DVT_new(i)=(alpha1/rhos*(1-xv(i)))*((Da(i)*Mv)/(R*(T_new(i)+273)))*HR(i)*dPvs(i);
lambdae_new(i)=(lambda+rhos*Lv*DVT_new(i)); gamma_new(i)=lambdae_new(i)/(rhocp);
end
%update new coefficient
DVX=DVX_new; DVT=DVT_new;
lambdae=lambdae_new; gamma=gamma_new;
%Dm=Dm_new; Delta=Delta_new;
for i=1:1:N+1
T1(t+1,i)=T_new(i)-Ta; % data stocking
DVX1(t+1,i)=DVX(i); lambdae1(t+1,i)=lambdae(i);
DVT1(t+1,i)=(DVT(i)); gamma1(t+1,i)=gamma(i);
X1(t+1,i)=X_new(i); %Dm1(t+1,i)=Dm(i);
%Delta1(t+1,i)=Delta(i);
end
% if abs(v-T)<err
% break
% end
% end
end
X1; T1;
for t=1:1:M+1
time(t,1)=t*dt;
end
for i = 1:N+1
r(i) =(i)*dr;
end
figure('Name','coef de diff non isotherm','NumberTitle','off')
plot(time,DVX1(:,2),'linewidth',5)
%title('Coefficient de diffusion non isotherme'),
xlabel('temps(s)'),ylabel('DVX(m2/s)')
grid on
%subplot(2,2,1)
figure('Name','Coefficient de diffusion isotherme','NumberTitle','off')
plot(time,DVT1(:,2),'linewidth',5)
%title('Coefficient de diffusion isotherme'),
xlabel('time(s)'),ylabel('DVT(m2/s/K)')
grid on
%subplot(2,2,2)
figure('Name','Non isotherm diff coef vs water content','NumberTitle','off')
plot(X1(:,2),DVX1(:,2),'linewidth',5)
%title('Coefficient de diffusion non isotherme sur la face chauffée'),
xlabel('X(%)'),ylabel('DVX(m2/s)')
grid on
%subplot(2,2,3)
figure('Name','isotherm diff coef vs water content','NumberTitle','off')
plot(X1(:,2),DVT1(:,2),'linewidth',5)
%title('Coefficient de diffusion isotherme sur la face chauffée'),
xlabel('X(%)'),ylabel('DVT(m2/s/K)')
grid on
%subplot(2,2,4)
%subplot(1,2,1)
figure('Name','temperature vs lenght','NumberTitle','off')
plot(r,T1(100,:),'-',r,T1(500,:),'-',r,T1(800,:),'-',r,T1(1000,:),'-','linewidth',2.5)
%title('Temperature')%
xlabel('x(m)')
ylabel('T(°C)')
%subplot(1,2,2)
figure('Name','water vs lenght','NumberTitle','off')
plot(r,X1(100,:),'-',r,X1(500,:),'-',r,X1(800,:),'-',r,X1(1000,:),'-','linewidth',2.5)
%title('Teneur en eau')
xlabel('x(m)')
ylabel('X(%)')
%subplot(1,2,1)
figure('Name','temp in wall','NumberTitle','off')
mesh(T)
%title('Temperature dans la matériau')
xlabel('x(m)')
ylabel('temps(s)')
zlabel('T(°C)')
%subplot(1,2,2)
figure('Name','water in wall','NumberTitle','off')
mesh(X)
%title('Teneur en eau dans le matériau')
xlabel('x(m)')
ylabel('temps(s)')
zlabel('X(%)')

Answers (1)

Image Analyst
Image Analyst on 8 Dec 2022
In this line
T_new(i)=T(i,t)+a*((dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i+1,t)-a*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*T(i,t)+a*(dt/(dr*dr))*T(i-1,t)+b*((dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i+1,t)-b*(2*(dt/(dr*dr))+(1/r(i))*dt/(dr))*X(i,t)+b*(dt/(dr*dr))*X(i-1,t);
you don't yet have "a" defined. But it tries to use "a". You need to assign a value to "a" first before you attempt to use it.

Categories

Find more on Data Import from MATLAB in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by