[TRNSYS-users] W-editor Runtine Error
Denis GARASA
denis.garasa at etu.univ-nantes.fr
Wed Dec 7 12:18:05 PST 2011
Hello,
We're student and we're working on a project on TRNSys. We're are buiding
a program on W-editor that can be used on TRNSys's W-editor boxes.
We already have a working program on Matlab. We've adapted this program to
feet the W-editor language, but we always get a Runtime Error, without
more precision.
The program always us to calculate the evolution of the temperature with a
time dependence into a wall.
Is W-editor able to run a time dependent program?
When we'll get the program to work on W-editor, we want to put all entries
in TRNSys inputs and the solution matrix into an output.
So I attach to this mail the Matlab and the W-editor programs.
Thanks in advance
Best regards
Denis GARASA
-------------- next part --------------
%programme de bilan enthalpique sur un composant discrétisé en 4 noeuds
%interne (dont deux de surface).
%Les conditions limites sont définies dans les milieux 1 et 2;
%
%Le schéma equivalent est le suivant
%Apports radiatifs (IR1) | | IR2
%Echanges radiatifs (Tr1 et hr1) | T1-----T2-----T3-----T4 | Tr2 et hr2
%Ecanges convectifs (Te1 et hc1) | | Te2 et hc2
%Définition des conditions limites
%-------------milieu 1
IR1 = 0; % apports radiatifs (W/m²)
Tr1 = 30+273; % Température equivalente pour les échanges radiatifs (°C)
Te1 = 30+273; % Température équivalente pour les échange convectifs (°C)
hc1 = 20; % Coefficient d'échange convectif coté mimieu 1 (W/m²/K)
%-------------milieu 2
IR2 = 0; % apports radiatifs (W/m²)
Tr2 = 10+273; % Température equivalente pour les échanges radiatifs (°C)
Te2 = 10+273; % Température équivalente pour les échange convectifs (°C)
hc2 = 20; % Coefficient d'échange convectif coté mimieu 2 (W/m²/K)
%Définition des caractéristiques radiatives
eps1 = 0.9; % Emissivité sur la surface coté milieu 1
eps2 = 0.9; % Emissivité sur la surface coté milieu 2
alpha1 = 0.5 ; % Absorptivité du matériau coté milieu 1
alpha2 = 0.5 ; % Absorptivité du matériau coté milieu 2
tho1 = 0.5 ; % Transmitivité du matériau coté milieu 1
tho2 = 0.5 ; % Transmitivité du matériau coté milieu 2
%Définition des caractiristiques du matériau
lambda = 0.201; % Conductivité thermique (W/m/K)
rho = 400; % densité (kg/m3)
Cp = 1400; % Capacité calorifique (J/kg/K)
ep = 5E-2; % Epaisseur du matériau (m)
%Définition des conditions initiales
Tinit = 0; % Température de condition initale (°C)
%definition des cond et capacité en chaque noeud
cond = lambda / ep * 3;
capa = rho * Cp * ep / 3;
capa2 = capa / 2;
condneg = -1 * cond; %matrice négative des conductances
boltz = 5.67e-8;
% Temps final physique
tempsfinal = 1 * 3600;
% Données numériques
npas = 1;
NNpas = 1E+3;
dt = tempsfinal / npas;
conv = matrix(4, 1, 0);
% Initialisation du tableau de sauvegarde
Tsauv = matrix( 4, 1, 1) * 273;
Tnv = matrix( 4, npas, 1) * Tinit + 273;
% Choix de la progression
% Geométrique = 0;
% Arithmétique != 0;
rgeo = 2; % Raison de la progression géométrique
rarithm = 1; % Raison de la progression arithmétique
progression = 1;
%programme de convergence des températures
convergence = 1e4;
while convergence > 1e-3
%while npas < NNpas
% Sélection de la progression
if ( progression )
% Incrémentation du nombre de pas de temps
npas = npas + rarithm;
else
npas = npas * rgeo;
end
% Calcul de la valeur du pas de temps
dt = tempsfinal / ( npas - 1 );
% Initialisation
Tnv = matrix( 4, npas, 1) * Tinit + 273;
for n = 2 : 1 : npas
%calcul des coefficients d'échange radiatifs
hr1 = eps1 * boltz * ( Te1^2 + Tnv(1,n-1)^2 ) * ( Te1 + Tnv(1,n-1) );
hr2 = eps2*boltz*(Te2^2+Tnv(4,n-1)^2)*(Te2+Tnv(4,n-1));
%definition de la matrice des capacités
MAT(1,1) = capa2 / dt + hr1 + hc1 + cond;
MAT(1,2) = condneg;
MAT(1,3) = 0;
MAT(1,4) = 0;
MAT(2,1) = condneg;
MAT(2,2) = capa / dt + 2 * cond;
MAT(2,3) = condneg;
MAT(2,4) = 0;
MAT(3,1) = 0;
MAT(3,2) = condneg;
MAT(3,3) = capa / dt + 2 * cond;
MAT(3,4) = condneg;
MAT(4,1) = 0;
MAT(4,2) = 0;
MAT(4,3) = condneg;
MAT(4,4) = capa2/ dt + cond + hr2 + hc2;
% Définition du veteur résultat
SOL(1) = capa2 / dt * Tnv(1,n-1) + hc1 * Te1 + hr1 * Tr1 + alpha1 * IR1;
SOL(2) = capa/dt * Tnv(2,n-1);
SOL(3) = capa/dt * Tnv(3,n-1);
SOL(4) = capa2/dt*Tnv(4,n-1)+hc2*Te2+hr2*Tr2+alpha2*IR2;
% Résolution du système
IMAT = inv(MAT);
SOLt = [SOL(1),SOL(2),SOL(3),SOL(4)];
for k = 1 : 1 : 4;
Tnv(k,n) = IMAT(k,1) * SOLt(1) + IMAT(k,2) * SOLt(2) + IMAT(k,3) * SOLt(3) + IMAT(k,4) * SOLt(4);
end
end
% Calcul de la norme infinie
for k = 1 : 1 : 4;
conv(k) = Tnv(k,npas) - Tsauv(k);
end
% Calcul du carré de la norme euclidienne d'ordre 2
convergence = 0;
for n = 1 : 1 : 4
convergence = convergence + conv(n)^2;
end
% Sauvegarde du champ de température
for k = 1 : 1 : 4;
Tsauv(k) = Tnv(k,npas);
end
end
-------------- next part --------------
%programme de bilan enthalpique sur un composant discrétisé en 4 noeuds
%interne (dont deux de surface).
%Les conditions limites sont définies dans les milieux 1 et 2;
%
%Le schéma equivalent est le suivant
clear all;
clc;
tic;
%Apports radiatifs (IR1) | | IR2
%Echanges radiatifs (Tr1 et hr1) | T1-----T2-----T3-----T4 | Tr2 et hr2
%Ecanges convectifs (Te1 et hc1) | | Te2 et hc2
%Définition des conditions limites
%-------------milieu 1
IR1 = 0; % apports radiatifs (W/m²)
Tr1 = 30+273; % Température equivalente pour les échanges radiatifs (°C)
Te1 = 30+273; % Température équivalente pour les échange convectifs (°C)
hc1 = 20; % Coefficient d'échange convectif coté mimieu 1 (W/m²/K)
%-------------milieu 2
IR2 = 0; % apports radiatifs (W/m²)
Tr2 = 10+273; % Température equivalente pour les échanges radiatifs (°C)
Te2 = 10+273; % Température équivalente pour les échange convectifs (°C)
hc2 = 20; % Coefficient d'échange convectif coté mimieu 2 (W/m²/K)
%Définition des caractéristiques radiatives
eps1 = .9; % Emissivité sur la surface coté milieu 1
eps2 = .9; % Emissivité sur la surface coté milieu 2
alpha1 = .5 ; % Absorptivité du matériau coté milieu 1
alpha2 = .5 ; % Absorptivité du matériau coté milieu 2
tho1 = .5 ; % Transmitivité du matériau coté milieu 1
tho2 = .5 ; % Transmitivité du matériau coté milieu 2
%Définition des caractiristiques du matériau
lambda = 0.201; % Conductivité thermique (W/m/K)
rho = 400; % densité (kg/m3)
Cp = 1400; % Capacité calorifique (J/kg/K)
ep = 5E-2; % Epaisseur du matériau (m)
%Définition des conditions initiales
Tinit = 0; % Température de condition initale (°C)
%definition des cond et capacité en chaque noeud
cond = lambda / ep * 3;
capa = rho * Cp * ep / 3;
capa2 = capa / 2;
condneg = -1 * cond; %matrice négative des conductances
boltz = 5.67e-8;
% Temps final physique
tempsfinal = 1 * 3600;
% Données numériques
npas = 1;
NNpas = 1E+3;
dt = tempsfinal / npas;
conv = zeros( 4, 1 );
% Initialisation du tableau de sauvegarde
T = ones( 4, 1) * inf;
Tnv = ones( 4, npas ) * Tinit + 273;
% Choix de la progression
% Geométrique = 0;
% Arithmétique != 0;
rgeo = 2; % Raison de la progression géométrique
rarithm = 1; % Raison de la progression arithmétique
progression = 1;
%programme de convergence des températures
convergence = inf;
while ( ( convergence > 1E-3 ) && ( npas < NNpas ) )
% Sélection de la progression
if ( progression )
% Incrémentation du nombre de pas de temps
npas = npas + rarithm;
else
npas = npas * rgeo;
end
% Calcul de la valeur du pas de temps
dt = tempsfinal / ( npas - 1 );
% Initialisation
Tnv = ones( 4, npas ) * Tinit + 273;
for n = 2 : 1 : npas
%calcul des coefficients d'échange radiatifs
hr1 = eps1 * boltz * ( Te1^2 + Tnv(1,n-1)^2 ) * ( Te1 + Tnv(1,n-1) );
hr2 = eps2 * boltz * ( Te2^2 + Tnv(4,n-1)^2 ) * ( Te2 + Tnv(4,n-1) );
%definition de la matrice des capacités
MAT(1,1) = capa2 / dt + hr1 + hc1 + cond;
MAT(1,2) = condneg;
MAT(1,3) = 0;
MAT(1,4) = 0;
MAT(2,1) = condneg;
MAT(2,2) = capa / dt + 2 * cond;
MAT(2,3) = condneg;
MAT(2,4) = 0;
MAT(3,1) = 0;
MAT(3,2) = condneg;
MAT(3,3) = capa / dt + 2 * cond;
MAT(3,4) = condneg;
MAT(4,1) = 0;
MAT(4,2) = 0;
MAT(4,3) = condneg;
MAT(4,4) = capa2/ dt + cond + hr2 + hc2;
% Définition du veteur résultat
SOL(1) = capa2 / dt * Tnv(1,n-1) + hc1 * Te1 + hr1 * Tr1 + alpha1 * IR1;
SOL(2) = capa / dt * Tnv(2,n-1);
SOL(3) = capa / dt * Tnv(3,n-1);
SOL(4) = capa2 / dt* Tnv(4,n-1) + hc2 * Te2 + hr2 * Tr2 + alpha2 * IR2;
% Résolution du système
IMAT = inv(MAT);
SOLt = [SOL(1),SOL(2),SOL(3),SOL(4)];
for k = 1 : 1 : 4;
Tnv(k,n) = IMAT(k,1) * SOLt(1) + IMAT(k,2) * SOLt(2) + IMAT(k,3) * SOLt(3) + IMAT(k,4) * SOLt(4);
end
end
% Calcul de la norme infinie
for k = 1 : 1 : 4;
conv(k) = Tnv(k,npas) - T(k);
end
% Calcul du carré de la norme euclidienne d'ordre 2
convergence = 0.;
for n = 1 : 1 : 4
convergence = convergence + conv(n)^2;
end
% Sauvegarde du champ de température
for k = 1 : 1 : 4;
T(k) = Tnv(k,npas);
end
end
duree = toc;
fprintf( '\n\n Temps de calcul : %.2E s', duree );
% Abscisse de traçage des graphiques
temps = linspace( 0, tempsfinal, npas );
% Ordonnées initiales
VTe1 = ones(1,length(Tnv(4,:))).*Te1-273;
VTe2 = ones(1,length(Tnv(4,:))).*Te2-273;
Tnv = Tnv - 273;
figure(1);
clf;
hold on;
plot( temps, Tnv(1,:),'g', temps,Tnv(2,:),'y',temps,Tnv(3,:),'b',temps,Tnv(4,:),'r')
plot( temps, VTe1, 'k','linewidth', 2 );
plot( temps, VTe2, 'k', 'linewidth', 2 );
grid on;
box on;
xlabel('Temps (h)');
ylabel('Température (°C)');
legend('Tsurf1','Tint12','Tint21','Tsurf2','Tmilieu1','Tmilieu2', 'location', 'best' );
axis( [ 0 tempsfinal min(Tinit,5) max(Te1,Te2)-273+5 ] );
fprintf( '\n\n\t\tFIN\n\n' );
More information about the TRNSYS-users
mailing list