base_code.m

%% mod�lisation du swirl pot
clear all; close all
load('baseline2020results_damping.mat')
load('modelmoteur2019_2.mat')
% load('compare_data')
load('compare_data_wings_no_fan')
% load('compare_data_Lucka2019-04-13')
RPM=data{89}.Values
power=data{32}.Values
car_speed=data{91}.Values


RPM_to_power=[modelmoteur2019_2(:,1) modelmoteur2019_2(:,1)*2*pi/60.*modelmoteur2019_2(:,2)];
power=interp1(RPM_to_power(:,1),RPM_to_power(:,2),compare_rpm).*(compare_tps/100*0.9+0.1);
compare_time=compare_time-compare_time(1);
compare_speed=compare_speed/3.6;




m_vap=0;
m_dot=0.6;
P=101300 ;
T_in=compare_coolant_temp(1);
Q_out=10;
i=1;
% m_dot_out_s=m_dot;
V_specific=0.001043;
% T_out=20.5; % cas validation
T_out=18; %cas v
temps_max=200
alpha=0;%1.43e-7; %water conduction
H_speed=7.0178;

courant_desirer=0.8;

%Volume pour eau Composante m^3
dv=0.00005;
Volume_m=0.0007;
Volume_tube_1=0.0001;
Volume_swirl=0.0001;
Volume_tube_2=0.0001;
Volume_rad=0.0009;
Volume_tube_3=0.0001;
Volume_total=Volume_m+Volume_tube_1+Volume_swirl+Volume_tube_2+Volume_rad+Volume_tube_3;

sections=[Volume_m Volume_m+Volume_tube_1 Volume_m+Volume_tube_1+Volume_swirl Volume_m+Volume_tube_1+Volume_swirl+Volume_tube_2 Volume_m+Volume_tube_1+Volume_swirl+Volume_tube_2+Volume_rad Volume_m+Volume_tube_1+Volume_swirl+Volume_tube_2+Volume_rad+Volume_tube_3];


%Discretisation 
dt=0.02;

V=linspace(0,Volume_total,Volume_total/dv+1)
V(end)=[];
T=ones(length(V),1)*T_in;
X=zeros(length(V),1);
Volume_swirl=dv
T_2=T;
time=0
P3=P;
T3=[T]
V2=V;
V3=V';
X_2=X;
V_out=ones(length(V),1)*V_specific;
V_out2=V_out;
C_M=conduction_matrix(V,dv)*alpha;
I=eye(length(V));




figure(1)
figure(2)
dt_select=[0.04 0.03125 0.025 0.02 0.015625 0.01 ];

% hyper cube sampling 100
hyper_cube=lhsdesign(100,3)

H_speed_lhs=hyper_cube(:,2)*0.927225*2+7.845025-0.927225;
ratio_heat_power_lhs=hyper_cube(:,1)*0.2*2+0.8-0.2;
pompe_flow_lhs=hyper_cube(:,3)*0.1057*2+1.057-0.1057;


% H_speed_lhs=7.845025*1.08;
% ratio_heat_power_lhs=0.7;
% pompe_flow_lhs=1.057;

ratio_heat_power=1
for convergence=1:length(ratio_heat_power_lhs)
    m_vap=0;
    P=101300;
    V=linspace(0,Volume_total,Volume_total/dv+1)
    V(end)=[];
    T=ones(length(V),1)*T_in;
    X=zeros(length(V),1);
    Volume_swirl=dv
    T_2=T;
    time=0
    P3=P;
    T3=[T]
    V2=V;
    V3=V';
    X_2=X;
    V_out=ones(length(V),1)*V_specific;
    V_out2=V_out;
    C_M=conduction_matrix(V,dv)*alpha;
    I=eye(length(V));
    
    m_dot_max=pompe_flow_lhs(convergence)*0.0001*max(RPM.Data)
    
    dt=dv*courant_desirer/(V_specific*m_dot_max);
    dt=temps_max/round(temps_max/dt);
    
    T_pour_compare=T_in;
    
    H_speed=H_speed_lhs(convergence);
    %nombre de courant constant=
    for t=0:dt:temps_max-dt
%         Q_in=interp1(power.Time,power.Data,mod(t,max(RPM.Time)),'makima')*ratio_heat_power_lhs(convergence);
%         m_dot_in=pompe_flow_lhs(convergence)*interp1(RPM.Time,RPM.Data,mod(t,max(RPM.Time)),'makima')*0.0001;
%         vitesse=interp1(car_speed.Time,car_speed.Data,mod(t,max(RPM.Time)),'makima');
%         vitesse=20;
%         m_dot_in=1;
%         Q_in=8000;
%         
        Q_in=interp1(compare_time,power,mod(t,max(compare_time)),'makima')*ratio_heat_power_lhs(convergence);
        m_dot_in=pompe_flow_lhs(convergence)*interp1(compare_time,compare_rpm,mod(t,max(compare_time)),'makima')*0.0001;
        vitesse=interp1(compare_time,compare_speed,mod(t,max(compare_time)),'makima');
        
        cell_in_swirl2=find(V>=(Volume_tube_1+Volume_m-dv/10)& V<(Volume_tube_1+Volume_m+Volume_swirl));
        for i=1:length(cell_in_swirl2)
            cell_in_swirl=cell_in_swirl2(i);
            if i==1
            
        if isempty(cell_in_swirl)
%             cell_in_swirl=find(V>=(Volume_tube_1+Volume_m-dv/10)& V<=(Volume_tube_1+Volume_m+Volume_swirl));
            if isempty(cell_in_swirl)
                error('Time_step_to big')
            else
                [P2,T_2(cell_in_swirl),m_vap2,m_dot_out]=swirlpot(m_dot_in,X(cell_in_swirl),P,T(cell_in_swirl),Q_out,m_vap,dt);
                X_2(cell_in_swirl)=0;
                V2(cell_in_swirl)=V(cell_in_swirl)+V_specific*m_dot_in*dt;
            end
        else
            
            [P2,T_2(cell_in_swirl),m_vap2,m_dot_out]=swirlpot(m_dot_in,X(cell_in_swirl),P,T(cell_in_swirl),Q_out,m_vap,dt);
            X_2(cell_in_swirl)=0;
            V2(cell_in_swirl)=V(cell_in_swirl)+V_specific*m_dot_in*dt;
        end
            else
            V2(cell_in_swirl)=V(cell_in_swirl)+V_specific*m_dot_in*dt;
            end
        end
        m_dot_out=m_dot_in; % perte de masse negligable
        for i=1:length(T)
            if V(i)<Volume_m
                %Moteur
                
                [T_2(i),X_2(i)]=moteur(m_dot_in,P,T(i),X(i),Q_in/(Volume_m/(V_specific*m_dot_in*dt)));
                V2(i)=V(i)+V_specific*m_dot_in*dt;
                if V2(i)<Volume_m
                else
                    T_2(i)=(1-(V2(i)-Volume_m)/(V2(i)-V(i)))*(T_2(i)-T(i))+T(i);
                    X_2(i)=(1-(V2(i)-Volume_m)/(V2(i)-V(i)))*(X_2(i)-X(i))+X(i);
                end
                
            elseif V(i)>=Volume_m & V(i)<(Volume_tube_1+Volume_m)
                %Tube1
                T_2(i)=T(i);
                X_2(i)=X(i);
                V2(i)=V(i)+V_specific*m_dot_in*dt;
            elseif V(i)>=(Volume_tube_1+Volume_m+Volume_swirl) & V(i)<(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2)
                %Tube 2
                T_2(i)=T(i);
                X_2(i)=X(i);
                V2(i)=V(i)+V_specific*m_dot_in*dt;
                Vmax=(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2);
                if V2(i)<(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2)
                else
                    T_2(i)=rad(m_dot_in,P,T(i),T_out,V_specific,Volume_rad,dt*(V2(i)-Vmax)/(V2(i)-V(i)),vitesse,H_speed);
                end
                
            elseif V(i)>=(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2) & V(i)<(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2+Volume_rad)
                %rad
                T_2(i)=rad(m_dot_in,P,T(i),T_out,V_specific,Volume_rad,dt,vitesse,H_speed);
                X_2(i)=X(i);
                V2(i)=V(i)+V_specific*m_dot_in*dt;
                Vmax=(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2+Volume_rad);
                if V2(i)<(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2+Volume_rad)
                else
                    T_2(i)=rad(m_dot_in,P,T(i),T_out,V_specific,Volume_rad,dt*(1-(V2(i)-Vmax)/(V2(i)-V(i))),vitesse,H_speed);
                end
                
            elseif V(i)>=(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2+Volume_rad)
                %tube3
                T_2(i)=T(i);
                X_2(i)=X(i);
                V2(i)=V(i)+V_specific*m_dot_in*dt;
                
                if V2(i)<Volume_total
                else
                    [T_2(i),X_2(i)]=moteur(m_dot_in,P,T(i),X(i),Q_in/(Volume_m/(V_specific*m_dot_in*(dt*(V2(i)-Volume_total)/(V2(i)-V(i))))));
                end
                
            end
        end
        T3=[T3 T_2];
        T=T_2;
        X=X_2;
        V=V2;
        V3=[V3 V'];
        P3(length(P3)+1,1)=P;
        m_vap=m_vap2;
        P=P2;
        V_out=V_out2;
        time=[time t+dt];
        V(V>Volume_total)=V(V>Volume_total)-Volume_total;
        disp('time solved')
        disp(t+dt)
        hold off
        
        %Conduction in water
        
        T=(I/dt-C_M)\(T/dt);
        
        %check degazeur
        select=V>(Volume_tube_1+Volume_m+Volume_swirl+dv);
        check_degaz=find(X(select)~=0);
        if isempty(check_degaz)==1
        else
            disp('degaz_not_respected')
            error('degazage non respecte')
        end
        
%               plot(V(:)',T(:)','.')
%               pause(0.001)
            if mod(t,2)==0
%                 z3 = ones(size(V(:)));
%                 figure(1)
%             color_line3(V(:)',z3,X(:)',X(:)','LineWidth',4);
%         
%             colorbar;
%             figure(2)
%             plot(V(1),T(1));
%             hold on
%             z2 = zeros(size(V(:)));
%         
%             color_line3(V(:)',z2,T(:)',T(:)','LineWidth',4);
%         
%             plot([sections(1) sections(1)],[-1 1])
%             plot([sections(2) sections(2)],[-1 1])
%             plot([sections(3) sections(3)],[-1 1])
%             plot([sections(4) sections(4)],[-1 1])
%             plot([sections(5) sections(5)],[-1 1])
%             plot([sections(6) sections(6)],[-1 1])
%             xlim([0 Volume_total]);
%             colorbar;
%             ylim([-1 1]);
%             pause(0.01);
%               plot(V(:)',T(:)','.')
%               pause(0.001)
            end
            
        select_T=V<(Volume_tube_1+Volume_m+Volume_swirl+Volume_tube_2) & V >Volume_m;
        T_pour_compare=[T_pour_compare ; mean(T(select_T))];
    end
  
    
    T_final(convergence)=(mean(T));
    T_error(convergence)=abs(mean(T)-71.375630660228140);
    delta_t(convergence)=dt;
    delta_v(convergence)=dv;
    mean_T{convergence}=mean(T3)'
    dt_all{convergence}=0:dt:temps_max
    figure(3)
    hold on
    plot(time,mean(T3))
%     plot(time,T_pour_compare)
    plot(compare_time,compare_coolant_temp)
    xlabel('time (s)')
    ylabel('Temperature moyenne (C)')
%     title('Comparaison entre température mesurée et simulée pour calibration')
    title('Comparaison entre température mesurée et simulée pour validation')
    pause(0.1)
    E=0;
    norm_erreur=1;
    T_pour_compare_avreage=mean(T3);
    for i=1:length(time)
        T_data(i)=interp1(compare_time,compare_coolant_temp,time(i));
        E(i)=abs((T_data(i)-T_pour_compare_avreage(i)))^norm_erreur;
    end
    error_norm(convergence)=(sum(E')/length(E))^(1/norm_erreur);

end
T_pour_compare_avreage=mean(T3);

%%% mesure de l'erreur entre la la modélisation et la validation

%erreur d'experimentale =+- 0.5 degree C

norm_erreur=1;
for i=1:length(time)
    T_data(i)=interp1(compare_time,compare_coolant_temp,time(i));
    E(i)=abs((T_data(i)-T_pour_compare_avreage(i)))^norm_erreur;
end
error_norm=(sum(E')/length(E))^(1/norm_erreur);







%  T=delaunay(ratio_heat_power_lhs,H_speed_lhs)
%  trisurf(T,ratio_heat_power_lhs,H_speed_lhs,T_final')
%  shading interp
%  contour(ratio_heat_power_lhs,H_speed_lhs,T_final')
%  
% figure()
% plot(ratio_heat_power_lhs,T_final,'.')
% % figure()
% % plot(H_speed_lhs,T_final,'.')
% % figure()
% % plot(pompe_flow_lhs,T_final,'.')
% 
% linear_fit_ratio_heat_power=polyfit(ratio_heat_power_lhs,T_final',1);
% linear_fit_H_speed=polyfit(H_speed_lhs,T_final',1);
% linear_fit_pompe_flow=polyfit(pompe_flow_lhs,T_final',1);
% 
% 
% Correct_T_final_1=T_final'-(polyval(linear_fit_ratio_heat_power,ratio_heat_power_lhs)-polyval(linear_fit_ratio_heat_power,0.8));
% figure()
% plot(H_speed_lhs,Correct_T_final_1,'.')
% 
% 
% linear_fit_H_speed_2=polyfit(H_speed_lhs,Correct_T_final_1,1);
% Correct_T_final_2=Correct_T_final_1-(polyval(linear_fit_H_speed_2,H_speed_lhs)-polyval(linear_fit_H_speed_2,7.845025));
% figure()
% plot(pompe_flow_lhs,Correct_T_final_2,'.')
% Correct_T_final_3=T_final'-(polyval(linear_fit_H_speed_2,H_speed_lhs)-polyval(linear_fit_H_speed_2,7.845025));
% plot(ratio_heat_power_lhs,Correct_T_final_3,'.')
% 
% 
% 
% 
% 
% ordre_approx=log((T_final(3)-T_final(2))/(T_final(2)-T_final(1)))/log(0.5);
% 
% ordre_approx=log((T_final(4)-T_final(3))/(T_final(3)-T_final(2)))/log(0.5);