Matlab-FLUENT热物性离散化代码

clear;
clc;
%296.15:0.1:650;
T=500:0.01:2000;
P=8e3;%kPa
cpupper_limit=1.005;
cplower_limit=0.995;
muupper_limit=1.002;
mulower_limit=0.998;
lambdaupper_limit=1.002;
lambdalower_limit=0.998;
rhoupper_limit=1.002;
rholower_limit=0.998;
betaupper_limit=1.002;
betalower_limit=0.998;
hupper_limit=1.002;
hlower_limit=0.998;
prupper_limit=1.002;
prlower_limit=0.998;
rho_Tupper_limit=1.002;
rho_Tlower_limit=0.998;
n=length(T);
ss=length(T);
rho=zeros(size(T));
mu=zeros(size(T));
lambda=zeros(size(T));
cp=zeros(size(T));
rho_T=zeros(size(T));
beta=zeros(size(T));
h=zeros(size(T));
pr=zeros(size(T));
for i=1:n
rho(i)=refpropm('D','T',T(i),'P',P,'CO2');
mu(i)=refpropm('V','T',T(i),'P',P,'CO2');
lambda(i)=refpropm('L','T',T(i),'P',P,'CO2');
cp(i)=refpropm('Cp','T',T(i),'P',P,'CO2');
rho_T(i)=refpropm('W','T',T(i),'P',P,'CO2');
h(i)=refpropm('H','T',T(i),'P',P,'CO2');
beta(i)=-1/rho(i)*(rho_T(i));
pr(i)=mu(i)*cp(i)/lambda(i);
end
%% density drop
a=1;
t_rho=T(1);
for i=1:n
    for b=a+2:n
        radio=(rho(a)+rho(b))*(b-a)/(2*sum(rho(a:b))-rho(a)-rho(b));
        if (radio<rholower_limit)||(radio>rhoupper_limit)
            a=b-1;
            t_rho(end+1)=T(a);
            %disp(a);
            break;
        end
        if b==n
            a=n;
            t_rho(end+1)=T(a);
            break;  
        end
    end
    if a==n-2
        radio=(rho(a)+rho(n))*(n-a)/(2*sum(rho(a:n))-rho(a)-rho(n));
        if (radio<rholower_limit)||(radio>rhoupper_limit)
        t_rho(end+1)=T(n-1);
        t_rho(end+1)=T(n);
        break;
        else
          t_rho(end+1)=T(n);  
          break;
        end
    end
     if a==n-1
         t_rho(end+1)=T(n);
         break;
     end
end
%% rho_T drop
a=1;
t_rho_T=T(1);
for i=1:n
    for b=a+2:n
        radio=(rho_T(a)+rho_T(b))*(b-a)/(2*sum(rho_T(a:b))-rho_T(a)-rho_T(b));
        if (radio<rho_Tlower_limit)||(radio>rho_Tupper_limit)
            a=b-1;
            t_rho_T(end+1)=T(a);
            %disp(a);
            break;
        end
        if b==n
            a=n;
            t_rho_T(end+1)=T(a);
            break;  
        end
    end
    if a==n-2
        radio=(rho_T(a)+rho_T(n))*(n-a)/(2*sum(rho_T(a:n))-rho_T(a)-rho_T(n));
        if (radio<rho_Tlower_limit)||(radio>rho_Tupper_limit)
        t_rho_T(end+1)=T(n-1);
        t_rho_T(end+1)=T(n);
        break;
        else
          t_rho_T(end+1)=T(n);  
          break;
        end
    end
     if a==n-1
         t_rho_T(end+1)=T(n);
         break;
     end
end
%% viscosity drop
a=1;
t_mu=T(1);
for i=1:n
    for b=a+2:n
        radio=( mu(a)+ mu(b))*(b-a)/(2*sum( mu(a:b))- mu(a)- mu(b));
        if (radio<mulower_limit)||(radio>muupper_limit)
            a=b-1;
            t_mu(end+1)=T(a);
            %disp(a);
            break;
        end
        if b==n
            a=n;
            t_mu(end+1)=T(a);
            break;  
        end
    end
    if a==n-2
        radio=( mu(a)+ mu(n))*(n-a)/(2*sum( mu(a:n))- mu(a)- mu(n));
        if (radio<mulower_limit)||(radio>muupper_limit)
        t_mu(end+1)=T(n-1);
        t_mu(end+1)=T(n);
        break;
        else
          t_mu(end+1)=T(n);  
          break;
        end
    end
     if a==n-1
         t_mu(end+1)=T(n);
         break;
     end
end
%% lambda drop
a=1;
t_lambda=T(1);
for i=1:n
    for b=a+2:n
        radio=( lambda(a)+ lambda(b))*(b-a)/(2*sum( lambda(a:b))- lambda(a)- lambda(b));
        if (radio<lambdalower_limit)||(radio>lambdaupper_limit)
            a=b-1;
            t_lambda(end+1)=T(a);
            %disp(a);
            break;
        end
        if b==n
            a=n;
            t_lambda(end+1)=T(a);
            break;  
        end
    end
    if a==n-2
        radio=( lambda(a)+ lambda(n))*(n-a)/(2*sum( lambda(a:n))- lambda(a)- lambda(n));
        if (radio<lambdalower_limit)||(radio>lambdaupper_limit)
        t_lambda(end+1)=T(n-1);
        t_lambda(end+1)=T(n);
        break;
        else
          t_lambda(end+1)=T(n);  
          break;
        end
    end
     if a==n-1
         t_lambda(end+1)=T(n);
         break;
     end
end
%% Cp drop
a=1;
t_cp=T(1);
for i=1:n
    for b=a+2:n
        radio=( cp(a)+ cp(b))*(b-a)/(2*sum( cp(a:b))- cp(a)- cp(b));
        if (radio<cplower_limit)||(radio>cpupper_limit)
            a=b-1;
            t_cp(end+1)=T(a);
            %disp(a);
            break;
        end
        if b==n
            a=n;
            t_cp(end+1)=T(a);
            break;  
        end
    end
    if a==n-2
        radio=( cp(a)+ cp(n))*(n-a)/(2*sum( cp(a:n))- cp(a)- cp(n));
        if (radio<cplower_limit)||(radio>cpupper_limit)
        t_cp(end+1)=T(n-1);
        t_cp(end+1)=T(n);
        break;
        else
          t_cp(end+1)=T(n);  
          break;
        end
    end
     if a==n-1
         t_cp(end+1)=T(n);
         break;
     end
end
%% pr drop
a=1;
t_pr=T(1);
for i=1:n
    for b=a+2:n
        radio=( pr(a)+ pr(b))*(b-a)/(2*sum( pr(a:b))- pr(a)- pr(b));
        if (radio<prlower_limit)||(radio>prupper_limit)
            a=b-1;
            t_pr(end+1)=T(a);
            %disp(a);
            break;
        end
        if b==n
            a=n;
            t_pr(end+1)=T(a);
            break;  
        end
    end
    if a==n-2
        radio=( pr(a)+ pr(n))*(n-a)/(2*sum( pr(a:n))- pr(a)- pr(n));
        if (radio<prlower_limit)||(radio>prupper_limit)
        t_pr(end+1)=T(n-1);
        t_pr(end+1)=T(n);
        break;
        else
          t_pr(end+1)=T(n);  
          break;
        end
    end
     if a==n-1
         t_pr(end+1)=T(n);
         break;
     end
end
%% beta drop
a=1;
t_beta=T(1);
for i=1:n
    for b=a+2:n
        radio=( beta(a)+ beta(b))*(b-a)/(2*sum( beta(a:b))- beta(a)- beta(b));
        if (radio<betalower_limit)||(radio>betaupper_limit)
            a=b-1;
            t_beta(end+1)=T(a);
            %disp(a);
            break;
        end
        if b==n
            a=n;
            t_beta(end+1)=T(a);
            break;  
        end
    end
    if a==n-2
        radio=( beta(a)+ beta(n))*(n-a)/(2*sum( beta(a:n))- beta(a)- beta(n));
        if (radio<betalower_limit)||(radio>betaupper_limit)
        t_beta(end+1)=T(n-1);
        t_beta(end+1)=T(n);
        break;
        else
          t_beta(end+1)=T(n);  
          break;
        end
    end
     if a==n-1
         t_beta(end+1)=T(n);
         break;
     end
end
%% h drop
a=1;
t_h=T(1);
for i=1:n
    for b=a+2:n
        radio=( h(a)+ h(b))*(b-a)/(2*sum( h(a:b))- h(a)- h(b));
        if (radio<hlower_limit)||(radio>hupper_limit)
            a=b-1;
            t_h(end+1)=T(a);
            %disp(a);
            break;
        end
        if b==n
            a=n;
            t_h(end+1)=T(a);
            break;  
        end
    end
    if a==n-2
        radio=( h(a)+ h(n))*(n-a)/(2*sum( h(a:n))- h(a)- h(n));
        if (radio<hlower_limit)||(radio>hupper_limit)
        t_h(end+1)=T(n-1);
        t_h(end+1)=T(n);
        break;
        else
          t_h(end+1)=T(n);  
          break;
        end
    end
     if a==n-1
         t_h(end+1)=T(n);
         break;
     end
end
%% 
n=length(t_rho);
for i=1:n
rho1(i)=refpropm('D','T',t_rho(i),'P',P,'CO2');
end
n=length(t_rho_T);
for i=1:n
rho_T1(i)=refpropm('W','T',t_rho_T(i),'P',P,'CO2');
end
n=length(t_mu);
for i=1:n
mu1(i)=refpropm('V','T',t_mu(i),'P',P,'CO2');
end
n=length(t_lambda);
for i=1:n
lambda1(i)=refpropm('L','T',t_lambda(i),'P',P,'CO2');
end
n=length(t_cp);
for i=1:n
cp1(i)=refpropm('Cp','T',t_cp(i),'P',P,'CO2');
end
n=length(t_pr);
for i=1:n
pr1(i)=refpropm('^','T',t_pr(i),'P',P,'CO2');
end
n=length(t_h);
for i=1:n
h1(i)=refpropm('H','T',t_h(i),'P',P,'CO2');
end
n=length(t_beta);
for i=1:n
rho2(i)=refpropm('D','T',t_beta(i),'P',P,'CO2');
rho_T1(i)=refpropm('W','T',t_beta(i),'P',P,'CO2');
beta1(i)=-1/rho2(i)*(rho_T1(i));%%%%%%%%%%%%%%%%%%%%%%%%%
end
%% 
str1='FINEco2_P_';
str2='T0_';
filename=[str1,num2str(P),str2,num2str(T(1)),'_',num2str(T(ss)),'.scm'];
f=fopen(filename,'w');
fprintf(f,'(\n');
fprintf(f,'(%s fluid\n(chemical-formula . #f)\n',filename);
%***************************rho**********************************%
fprintf(f,'(density (polynomial piecewise-linear ');
n=length(t_rho);
for i=1:n
    fprintf(f,'(%.2f . %.2f) ',t_rho(i),rho1(i));
end
fprintf(f,'))\n');
%****************************cp**********************************%
fprintf(f,'(specific-heat (polynomial piecewise-linear ');
n=length(t_cp);
for i=1:n
    fprintf(f,'(%.2f . %.2f) ',t_cp(i),cp1(i));
end
fprintf(f,'))\n');
%**************************lambda*******************************%
fprintf(f,'(thermal-conductivity (polynomial piecewise-linear ');
n=length(t_lambda);
for i=1:n
    fprintf(f,'(%.2f . %f) ',t_lambda(i),lambda1(i));
end
fprintf(f,'))\n');
%*********************************b****************************%
fprintf(f,'(viscosity (polynomial piecewise-linear ');
n=length(t_mu);
for i=1:n
    fprintf(f,'(%.2f . %e) ',t_mu(i),mu1(i));
end
fprintf(f,'))\n');
fprintf(f,')\n)\n');
fclose(f);
clearvars a b i n;