clc close all hidden tic tend = 3600 * lunar_day_duration ; %end of simulated lunar day (seconds) Instabilities in the simulation may appear if tend is set above 100 hours (360000s) tendnight = 3600 * lunar_night_duration; %end of simulated lunar night (seconds) %% MATERIALS PROPERTIES %NATIVE REGOLITH PROPERTIES % k_nat=0.01; thermal conductivity (W/Km) Uncomment to use constant properties k_nat = @(~,state) 9.3e-3*(1 + 0.073*((state.u.*state.u.*state.u)/(350*350*350))); %temperature-dependent thermal conductivity (W/km) density_nat = 1800; %density (kg/m3) cp_nat = 840; %specific heat (J/(kg*K)) %THERMAL MASS PROPERTIES % k_TM = 2.1; thermal conductivity (W/Km) Uncomment to use constant properties k_TM = @(~,state) 6e-7*state.u.*state.u - 0.0028*state.u + 3.3753; %temperature-dependent thermal conductivity density_TM = 3000; %density (kg/m3) % cp_TM = 800; %specific heat (J/(kg*K)) Uncomment to use constant properties cp_TM = @(~,state) -0.0005*state.u.*state.u + 1.4332*state.u + 371.5; %temperature-dependent specific heat (J/(kg*K)) emissivity_TM = 0.96; %top surface emissivity(-) TM_absorptivity = 0.85; %top surface absorptivity (-) %FLUFF LAYER PROPERTIES density_fluff = 1300; %Density of the fluff layer of regolith (kg/m3) cp_fluff = 854; %specific heat of the fluff layer of regolith (J/(kg.K)) % k_fluff=2.287e-3; thermal conductivity (W/Km) Uncomment to use constant properties k_fluff = @(~,state) 9.22e-4*(1 + 1.48*((state.u.*state.u.*state.u)/(350*350*350))); %temperature-dependent thermal conductivity (W/km) %% CREATE THE MODEL numberOfPDE=1; %number of partial differential equqtions to be solved model=createpde(numberOfPDE); %create the model object %% CREATE THE GEOMETRY NATREG = [3; 4; -fluff_width; TM_diam+fluff_width; TM_diam+fluff_width; -fluff_width; -depth_nat_reg; -depth_nat_reg; TM_height-fluff_thickness; TM_height-fluff_thickness]; TM = [2; 4; 0; TM_diam; TM_diam; 0; 0; 0; TM_height; TM_height]; FLUFF = [2; 4; -fluff_width; TM_diam+fluff_width; TM_diam+fluff_width; -fluff_width; TM_height-fluff_thickness; TM_height-fluff_thickness; TM_height; TM_height]; TEG1 = [3; 4; TM_diam; TM_diam+0.01; TM_diam+0.01; TM_diam; TEG_y_bot; TEG_y_bot; TEG_y_top; TEG_y_top]; TEG2 = [3; 4; 0; 0-0.01; 0-0.01; 0; TEG_y_bot; TEG_y_bot; TEG_y_top; TEG_y_top]; HSP1 = [3; 4; TM_diam-HSP_thickness; TM_diam; TM_diam; TM_diam-HSP_thickness; HSP_y_bot; HSP_y_bot; HSP_y_top; HSP_y_top]; HSP2 = [3; 4; 0; 0+HSP_thickness; 0+HSP_thickness; 0; HSP_y_bot; HSP_y_bot; HSP_y_top; HSP_y_top]; TB = [3; 4 ; (TM_diam/2)-(thermal_beam_thickness/2); (TM_diam/2)+(thermal_beam_thickness/2); (TM_diam/2)+(thermal_beam_thickness/2); (TM_diam/2)-(thermal_beam_thickness/2); TM_height-thermal_beam_y_top-thermal_beam_length; TM_height-thermal_beam_y_top-thermal_beam_length; TM_height-thermal_beam_y_top; TM_height-thermal_beam_y_top]; gd = [NATREG TM FLUFF TEG1 TEG2 HSP1 HSP2 TB]; sf = 'NATREG - TEG1 -TEG2 + TM + FLUFF + HSP1 + HSP2 + TB'; ns = char('NATREG','TM','FLUFF','TEG1','TEG2','HSP1','HSP2','TB'); ns = ns'; dl = decsg(gd,sf,ns); geometryFromEdges(model,dl); generateMesh(model,'Hmax',mesh_size); %DISPLAY THE GEOMETRY % f1 = figure('Position',[30 500 600 500]) % pdegplot(model,'EdgeLabels','on','FaceLabels','on') % xlim([-0.3 0.5]) % ylim([-0.3 TM_height+0.1]) % axis equal; %% INITIAL CONDITIONS u0 = 254; %initial temperature of the system (K) Thot = 254; %initial temperature of the hot side (K) Tcold = 250; %Temperature at the cold side of TEG (K) (lower than 254K to avoid instabilities a the start of the simulation) Pr = 0; %initial power rejected by the TEG module (W) setInitialConditions(model,u0); %set initial temperatures to u0 value %% BOUNDARY CONDITIONS %Boundary conditions applyBoundaryCondition(model,'dirichlet','Edge',28,'u',254); %Temperature obtained over a certain penetration depth on the moon %Flux definition heat_flux_TM = RSI * m_factor * fresnel_eff * TM_absorptivity; %flux redirected and concentrated (W/m2) heat_flux_fluff = DSI; %flux from Direct Solar Irradiance (W/m2) %BOUNDARIES CONDITIONS applyBoundaryCondition(model,'neumann','Edge',8,'g',heat_flux_TM); %The top surface of the TM recieves the heat flux applyBoundaryCondition(model,'neumann','Edge',[7 9],'g',heat_flux_fluff); %The lunar soil which surround the TM recieves direct solar irradiance %% COEFFICIENTS OF THE PDE EQUATION % Native regolith m1 = 0; d1 = density_nat * cp_nat; c1 = k_nat; a1 = 0; f1 = 0; specifyCoefficients(model,'m',m1,'d',d1,'c',c1,'a',a1,'f',f1,'face',1); %Face 1 represents the native regolith % TM m2 = 0; d2 = @(~,state) density_TM*(-0.0005*state.u.*state.u + 1.4332*state.u + 371.5); %d2 = density * cp_TM c2 = k_TM; a2 = 0; f2 = 0; specifyCoefficients(model,'m',m2,'d',d2,'c',c2,'a',a2,'f',f2,'face',[2 4]); % Fluff-layer m3 = 0; d3 = density_fluff * cp_fluff; c3 = k_fluff; a3 = 0; f3 = 0; specifyCoefficients(model,'m',m3,'d',d3,'c',c3,'a',a3,'f',f3,'face',[3 5]); %Hot Sink Plate (HSP) m4 = 0; d4 = HSP_density * HSP_cp; c4 = HSP_conductivity; a4 = 0; f4 = 0; specifyCoefficients(model,'m',m4,'d',d4,'c',c4,'a',a4,'f',f4,'face',[6 7]); % %Thermal beam m5 = 0; d5 = thermal_beam_density * thermal_beam_cp; c5 = thermal_beam_conducitivy; a5 = 0; f5 = 0; specifyCoefficients(model,'m',m5,'d',d5,'c',c5,'a',a5,'f',f5,'face',8); %% SIMULATION %Declare matrices that will be used % Smatrix=[]; %to store nodal solution (temperature) at each iteration % SOL_DAY=[]; %store nodal solutions (temperatures) for the entire simulation % F=[]; %store a figure as frame % F2=[]; %stores all frames in a matrix % TTOP=[]; %stores the temperatures (K) at the top surface of the TM for the entire simulation % TBOT=[]; %stores the temperatures (K) at the bottom surface of the TM for the entire simulation % TMID=[]; %stores the temperatures (K) at the middle horizontal of the TM for the entire simulation % TTM=[]; %stores the a fctitious mean value of the TM temperature (K) for the entire simulation % DeltaT=[]; %stores the temperature difference between the hot and cold side % effTEG=[]; %stores the efficiency of the TEG for the entire simulation % PowerOut=[]; %stores the power output from one TEG element for the entire simulation (W) % PowerAbs=[]; %stores the power absorbed by one TEG element for the entire simulation (W) % THOT=[]; %stores the values of temperatures at the hot side of the TEG element for the entire simulation (K) % TCOLD = []; %stores the values of temperatures at the cold side of the TEG element for the entire simulation (K) % TRAD = []; %stores the values of the radiator temperature assumed to be homogeneous (K) % PowerRej=[]; %stores the power rejected by one TEG element for the entire simulation (W) %% Matrix Pre-allocation nNodes = size(model.Mesh.Nodes,2); ncolDay = lunar_day_duration*3600/step; ncolNight = lunar_night_duration*3600/step; ncolTot = ncolDay + ncolNight; F=[]; F2=[]; Smatrix = zeros(nNodes,2) ./0; SOL_DAY = zeros( nNodes,ncolDay)./0; SOL_NIGHT = zeros( nNodes,ncolNight)./0; vide = zeros(1,ncolTot)./0; TTOP = vide; TBOT = vide; TMID = vide; TTM = vide; DeltaT = vide; effTEG = vide; PowerOut = vide; PowerAbs = vide; PowerRej = vide; THOT = vide; TCOLD = vide; TRAD = vide; %% DEFINE LINES AND GEOMETRY WERE WE WANT TO KNOW THE TEMPERATURE %Line for the top surface temperature Xtop=linspace(0,TM_diam,5); %x coordinates of 5 points Ytop=[TM_height TM_height TM_height TM_height TM_height]; %y coordinates of 5 points %line for the temperature at the top of the fluff layer X2=linspace(-fluff_width,-0.1,5); %x coordinates of 5 points spaced wrt to the top temperature of the TM Y2=[TM_height TM_height TM_height TM_height TM_height]; %x coordinates of 5 points spaced wrt to the top temperature of the TM %line for the temperature at the bottom of the TM Xbot=linspace(0,TM_diam,5); %x coordinates of 5 points Ybot=[0 0 0 0 0]; %y coordinates of 5 points %line for the temperature at the middle of the TM Xmid=linspace(0,TM_diam,5); %x coordinates of 5 points Ymid=[TM_height/2 TM_height/2 TM_height/2 TM_height/2 TM_height/2]; %x coordinates of 5 points %line for the temperature at the hot sink of the TEG (TM-TEG interface) Xteg=[TM_diam TM_diam TM_diam TM_diam TM_diam]; %x coordinates of 5 points Yteg=linspace(TEG_y_bot,TEG_y_top,5); %y coordinates of 5 points %query points to compute the TM block mean temperature Xtm=[0 TM_diam/4 TM_diam/2 TM_diam/(4/3) TM_diam 0 TM_diam/4 TM_diam/2 TM_diam/(4/3) TM_diam 0 TM_diam/4 TM_diam/2 TM_diam/(4/3) TM_diam 0 TM_diam/4 TM_diam/2 TM_diam/(4/3) TM_diam 0 TM_diam/4 TM_diam/2 TM_diam/(4/3) TM_diam ]; Ytm=[0 0 0 0 0 TM_height/4 TM_height/4 TM_height/4 TM_height/4 TM_height/4 TM_height/2 TM_height/2 TM_height/2 TM_height/2 TM_height/2 TM_height/(4/3) TM_height/(4/3) TM_height/(4/3) TM_height/(4/3) TM_height/(4/3) TM_height TM_height TM_height TM_height TM_height]; format compact for i=step:step:tend %iteration instructions t1=i-step; %initial time (0 in the first iteration, 0+step at the second, etc...) t2=i %final time after iteration time=t1:step:t2; %time vector for the current iteration %final time after iteration (example: [600 900] %if step = 300s, next iteration will be[900 1200] results = solvepde(model,time); %solve the previously specified model over the specified time Smatrix = results.NodalSolution; %obtain the nodal solution (at t1 and t2) SOL_DAY(:,t2/step) = Smatrix(:,end); %t2 to the matrix which stores all solutions setInitialConditions(model,results); %Update the initial conditions based on the results. New IC are needed for the next iteration. Tintrpbot = interpolateSolution(results,Xbot,Ybot,1:length(time)); %Interpolate the temperature over the previously defined vector bot_temp=mean(Tintrpbot); %compute the average temperature on that location at t1 qnd t2 TM_temp_bot=bot_temp(:,end); %keep the temperature at t2 TBOT(t2/step) = TM_temp_bot; %stores those temperatures for the entire simulation Tintrpmid = interpolateSolution(results,Xmid,Ymid,1:length(time)); %Idem mid_temp=mean(Tintrpmid); TM_temp_mid=mid_temp(:,end); TMID(t2/step)=TM_temp_mid; Tintrptop = interpolateSolution(results,Xtop,Ytop,1:length(time)); %Idem rad_temp=mean(Tintrptop); TM_temp_top=rad_temp(:,end); TTOP(t2/step)=TM_temp_top; %The temperature at the top of the thermal mass should be determined to %compute the losses due to radiation heat_flux_loss=emissivity_TM*StefanBoltz*(TM_temp_top*TM_temp_top*TM_temp_top*TM_temp_top); heat_flux_net=heat_flux_TM-heat_flux_loss; %we reduce the incoming heat flux by radiative losses applyBoundaryCondition(model,'neumann','Edge',8,'g',heat_flux_net); %Boundaries conditions are updated based on the new net heat flux Tintrp2 = interpolateSolution(results,X2,Y2,1:length(time)); %Idem for the fluff layer incoming flux and radiative losses rad_temp2=mean(Tintrp2); mean_rad_temp2=rad_temp2(:,end); heat_flux_loss2=emissivity_TM*StefanBoltz*(mean_rad_temp2*mean_rad_temp2*mean_rad_temp2*mean_rad_temp2); heat_flux_net2=heat_flux_fluff-heat_flux_loss2; applyBoundaryCondition(model,'neumann','Edge',[7 9],'g',heat_flux_net2); %Boundaries conditions on the fluff layer are updated for the next iteration Tintrptm = interpolateSolution(results,Xtm,Ytm,1:length(time)); %Interpolation of the temperature for the query points previously defined tm_temp=mean(Tintrptm); %compute mean temperature of the thermal at t1 and t2 TM_temp=tm_temp(:,end); %we only keep temperature at t2 TTM(t2/step)=TM_temp; %stores the mean TM temperature in the matrix if boolFig fig = figure('InnerPosition',[30 500 600 500]); pdeplot(model,'XYData',SOL_DAY(:,t2/step),'Contour','off','ColorMap','hot'); %display the temperature profile at the end of the simulated day period. xlim([-0.3 0.5]) ylim([-0.3 TM_height+0.1]) caxis([250 1100]) %temperature range definition to adjust color scale axis equal; F=getframe; %capture the last figure display as a frame. It is the temperature profile at time t2. F2=[F2,F]; %stores the frame for each iteration if the F2 matrix. close(fig); end %Interpolate the hot temperature of the hot sie of the TEG Tintrteg = interpolateSolution(results,Xteg,Yteg,1:length(time)); temp_teg=mean(Tintrteg); %mean temperature at the hot plate Thot=temp_teg(:,end); THOT(t2/step)=Thot; %store hot temperature over the simulation of the day period %Simulink Cold Sink + radiator simulation %IMPORTANT: this section of the script was used when the simulation %used the simulink model of the cold sink coupled with a radiator. % The following lines are no more used for the advanbced model. % Due to instabilities problem induced, the cold sink temperature has been %modelled by a simple fonction which mimique the behavior previously %observed with the simulink model. If soneone wants to integrate a %simulink model, keep in mind that he could reuse and adapt the below %written script. % simName='Cold_sink_advanced_model'; %calls the simulink model % % configSet = getActiveConfigSet(simName); % configSetNames = get_param(configSet,'ObjectParameters'); % % paramNameValStruct.SaveOutput = 'on'; %indicates that % simulation output shall be saved % paramNameValStruct.OutputSaveName = 'yout'; %name of the output % object % simout = sim(simName,paramNameValStruct); %simulation execution % yout = simout.get('yout'); %get the ouput values out of the % simulation % % DataSim_cold = yout{1}.Values.Data; %get the first output data % Tcold = DataSim_cold(end,:); %retain only the last value of the % simulation (in that case, it was the temperature of the cold sink) % DataSim_rad = yout{2}.Values.Data; %get the second ouput data % Trad = DataSim_rad(end,:);%retain only the last value of the % simulation (in that case, it was the temperature of the radiator) Psun=RAD_area*RAD_abs*DSI; P_rad_space = RAD_area*RAD_emissivity*StefanBoltz*(Tcold*Tcold*Tcold*Tcold - Tspace*Tspace*Tspace*Tspace); Delta_Tcold = (1/(RAD_mass*RAD_specific_heat))*( Pr*2 + Psun - P_rad_space)*step; Tcold = Tcold + Delta_Tcold; TCOLD(t2/step) = Tcold; % THERMOELECTRIC FUNCTION FOR SIMULATION [ TEG_DeltaT,Voc,Uoc,TEG_eff,Riteg,Iload,Uload,tcond,TEG_Pout,Prej,TEG_Pabs]=tegmoduleSiGe(Thot,Tcold,n); DeltaT(t2/step) = TEG_DeltaT; %Adds the last DeltaT value to the relevant warehouse matrix effTEG(t2/step) = TEG_eff; %Adds the last TEG_eff value to the relevant warehouse matrix PowerOut(t2/step) = TEG_Pout; %Adds the last TEG_Pout value to the relevant warehouse matrix PowerAbs(t2/step) = TEG_Pabs; %Adds the last TEG_Pabs value to the relevant warehouse matrix TEG_abs_flux = PowerAbs(:,t2/step)/TEG_side; %Conpute the absorbed flux based on the TEG side (considering a 1m length perpendicular to the screen). applyBoundaryCondition(model,'neumann','Edge',[22 36],'g',-TEG_abs_flux); %generates a boundary condition %on the hot plate which is the absorbed flux by the TEG Pr = PowerAbs(:,(t2/step)) - PowerOut(:,(t2/step)); %Power rejected at the cold side (W) to be evacuated towards radiator PowerRej(t2/step) = Pr; %Stores the rejected power in the relevant warehouse matrix end transient_results = results; %stores the results in a transient result object. The next section of the script %is the night, and the initial conditions will depend on the results of the day period. %% NIGHT SIMULATION emissivity_TM_night=emissivity_TM/emissivity_reduction_night; numberOfPDE=1; %IDEM model=createpde(numberOfPDE); NATREG = [3; 4; -fluff_width; TM_diam+fluff_width; TM_diam+fluff_width; -fluff_width; -depth_nat_reg; -depth_nat_reg; TM_height-fluff_thickness; TM_height-fluff_thickness]; TM = [2; 4; 0; TM_diam; TM_diam; 0; 0; 0; TM_height; TM_height]; FLUFF = [2; 4; -fluff_width; TM_diam+fluff_width; TM_diam+fluff_width; -fluff_width; TM_height-fluff_thickness; TM_height-fluff_thickness; TM_height; TM_height]; TEG1 = [3; 4; TM_diam; TM_diam+0.01; TM_diam+0.01; TM_diam; TEG_y_bot; TEG_y_bot; TEG_y_top; TEG_y_top]; TEG2 = [3; 4; 0; 0-0.01; 0-0.01; 0; TEG_y_bot; TEG_y_bot; TEG_y_top; TEG_y_top]; HSP1 = [3; 4; TM_diam-HSP_thickness; TM_diam; TM_diam; TM_diam-HSP_thickness; HSP_y_bot; HSP_y_bot; HSP_y_top; HSP_y_top]; HSP2 = [3; 4; 0; 0+HSP_thickness; 0+HSP_thickness; 0; HSP_y_bot; HSP_y_bot; HSP_y_top; HSP_y_top]; TB = [3; 4 ; (TM_diam/2)-(thermal_beam_thickness/2); (TM_diam/2)+(thermal_beam_thickness/2); (TM_diam/2)+(thermal_beam_thickness/2); (TM_diam/2)-(thermal_beam_thickness/2); TM_height-thermal_beam_y_top-thermal_beam_length; TM_height-thermal_beam_y_top-thermal_beam_length; TM_height-thermal_beam_y_top; TM_height-thermal_beam_y_top]; gd = [NATREG TM FLUFF TEG1 TEG2 HSP1 HSP2 TB]; sf = 'NATREG - TEG1 -TEG2 + TM + FLUFF + HSP1 + HSP2 + TB'; ns = char('NATREG','TM','FLUFF','TEG1','TEG2','HSP1','HSP2','TB'); ns = ns'; dl = decsg(gd,sf,ns); geometryFromEdges(model,dl); generateMesh(model,'Hmax',mesh_size); % figure(3) % pdegplot(model,'EdgeLabels','on','FaceLabels','on') % xlim([-0.3 0.5]) % ylim([-0.3 TM_height+0.1]) % axis equal; % % figure(4) % pdeplot(model) % xlim([-0.3 0.5]) % ylim([-0.3 TM_height+0.1]) % axis equal; %Set initial conditions setInitialConditions(model,transient_results); %Boundary conditions applyBoundaryCondition(model,'dirichlet','Edge',28,'u',254); heat_flux_night=0; %During the night, there are no incoming flux from the sun heat_flux_night2=0; %during the night, there are no incoming flux from the sun applyBoundaryCondition(model,'neumann','Edge',8,'g',heat_flux_night); applyBoundaryCondition(model,'neumann','Edge',[7 9],'g',heat_flux_night2); % COEFFICIENTS % Native regolith m1=0; d1=density_nat*cp_nat; c1=k_nat; a1=0; f1=0; specifyCoefficients(model,'m',m1,'d',d1,'c',c1,'a',a1,'f',f1,'face',1); % TM m2=0; d2 = @(~,state) density_TM*(-0.0005*state.u.*state.u + 1.4332*state.u + 371.5); k_TM = @(~,state) 6e-7*state.u.*state.u - 0.0028*state.u + 3.3753; a2=0; f2=0; specifyCoefficients(model,'m',m2,'d',d2,'c',k_TM,'a',a2,'f',f2,'face',[2 4]); % Fluff m3=0; d3=density_fluff*cp_fluff; c3=k_fluff; a3=0; f3=0; specifyCoefficients(model,'m',m3,'d',d3,'c',c3,'a',a3,'f',f3,'face',[3 5]); %Hot Sink Plate (HSP) m4=0; d4=HSP_density*HSP_cp; c4=HSP_conductivity; a4=0; f4=0; specifyCoefficients(model,'m',m4,'d',d4,'c',c4,'a',a4,'f',f4,'face',[6 7]); %Thermal beam m5=0; d5=thermal_beam_density*thermal_beam_cp; c5=thermal_beam_conducitivy; a5=0; f5=0; specifyCoefficients(model,'m',m5,'d',d5,'c',c5,'a',a5,'f',f5,'face',8); % NIGHT SIMULATION for i=step:step:tendnight; t1=i-step; t2=i time=t1:step:t2; ind = lunar_day_duration*(3600/step) + t2/step; results = solvepde(model,time); Smatrix = results.NodalSolution; SOL_NIGHT(:,t2/step) = Smatrix(:,end); setInitialConditions(model,results); Tintrpbot = interpolateSolution(results,Xbot,Ybot,1:length(time)); bot_temp=mean(Tintrpbot); TM_temp_bot=bot_temp(:,end); TBOT(ind) = TM_temp_bot; Tintrpmid = interpolateSolution(results,Xmid,Ymid,1:length(time)); mid_temp=mean(Tintrpmid); TM_temp_mid=mid_temp(:,end); TMID(ind) = TM_temp_mid; Tintrptop = interpolateSolution(results,Xtop,Ytop,1:length(time)); rad_temp=mean(Tintrptop); TM_temp_top=rad_temp(:,end); TTOP(ind) = TM_temp_top; %Although there is no solar incoming flux, radiative losses have to be %taken into account. We take into account the reduced emissivity of the %thermal mass. heat_flux_loss = emissivity_TM_night*StefanBoltz*(TM_temp_top*TM_temp_top*TM_temp_top*TM_temp_top); heat_flux_night_net = heat_flux_night - heat_flux_loss; applyBoundaryCondition(model,'neumann','Edge',8,'g',heat_flux_night_net); Tintrp2 = interpolateSolution(results,X2,Y2,1:length(time)); rad_temp2=mean(Tintrp2); mean_rad_temp2=rad_temp2(:,end); heat_flux_loss2=emissivity_TM*StefanBoltz*(mean_rad_temp2*mean_rad_temp2*mean_rad_temp2*mean_rad_temp2); heat_flux_night_net2 = heat_flux_night2-heat_flux_loss2; applyBoundaryCondition(model,'neumann','Edge',[7 9],'g',heat_flux_night_net2); Tintrptm = interpolateSolution(results,Xtm,Ytm,1:length(time)); tm_temp=mean(Tintrptm); TM_temp=tm_temp(:,end); TTM(ind) = TM_temp; if boolFig fig3 = figure('InnerPosition',[30 500 600 500]); %display in figure 5 the last temperature profile of the night pdeplot(model,'XYData',SOL_NIGHT(:,t2/step),'Contour','off','ColorMap','hot'); xlim([-0.3 0.5]); ylim([-0.3 TM_height+0.1]); caxis([250 1100]) %temperature range to set the colormap axis equal; F=getframe; F2=[F2,F]; close(fig3); end %Interpolate the hot temperature of the hot side of the TEG Tintrteg = interpolateSolution(results,Xteg,Yteg,1:length(time)); temp_teg=mean(Tintrteg); Thot=temp_teg(:,end); THOT(ind) = Thot; %Simulink Cold Sink + radiator simulation Psun=0; P_rad_space = RAD_area*RAD_emissivity*StefanBoltz*(Tcold*Tcold*Tcold*Tcold - Tspace*Tspace*Tspace*Tspace); Delta_Tcold = (1/(RAD_mass*RAD_specific_heat))*( Pr*2 + Psun - P_rad_space)*step; Tcold = Tcold + Delta_Tcold; TCOLD(ind) = Tcold; %stores the values of the cold temperature (constant during the day but not during the night) %IF WE WANT TO USE A SIMULINK %IDEM % simName='Cold_sink_advanced_model'; %calls the simulink model % % configSet = getActiveConfigSet(simName); % configSetNames = get_param(configSet,'ObjectParameters'); % % paramNameValStruct.SaveOutput = 'on'; %indicates that % simulation output shall be saved % paramNameValStruct.OutputSaveName = 'yout'; %name of the output % object % simout = sim(simName,paramNameValStruct); %simulation execution % yout = simout.get('yout'); %get the ouput values out of the % simulation % % DataSim_cold = yout{1}.Values.Data; %get the first output data % Tcold = DataSim_cold(end,:); %retain only the last value of the % simulation (in that case, it was the temperature of the cold sink) % DataSim_rad = yout{2}.Values.Data; %get the second ouput data % Trad = DataSim_rad(end,:);%retain only the last value of the % simulation (in that case, it was the temperature of the radiator) % TRAD = [TRAD,Trad]; %stores the temperature of the radiator %TEG FUNCTION [ TEG_DeltaT,Voc,Uoc,TEG_eff,Riteg,Iload,Uload,tcond,TEG_Pout,Prej,TEG_Pabs]=tegmoduleSiGe(Thot,Tcold,n); DeltaT(ind) = TEG_DeltaT; effTEG(ind) = TEG_eff; PowerOut(ind) = TEG_Pout; PowerAbs(ind) = TEG_Pabs; TEG_abs_flux = PowerAbs(:,ind)/TEG_side; applyBoundaryCondition(model,'neumann','Edge',[22 36],'g',-TEG_abs_flux); Pr = PowerAbs(:,ind) - PowerOut(:,ind); %Power rejected at the cold side (W) PowerRej(ind) = Pr; %Stores rejected power values end %this result matrix stores all warehouse matrices to facilitate %copy-past fron MATLAB workspace to excel spreadsheet if needed. MATRIX_RESULTS =[ TTOP; TMID; TBOT; TTM; THOT; TCOLD; DeltaT; PowerOut; PowerAbs; PowerRej; effTEG]; %Return the last values computed by the TEG function. This values are given %for one TEG element. Be careful, there are 2 TEG element presently %implemented in the model: one on each side of the TM. [ DeltaT_end,Voc_end,Uoc_end,effteg_end,Riteg_end,Iload_end,Uload_end,tcond_end,Pelec_end,Pr_end,Pa_end ] = tegmoduleSiGe( Thot,Tcold,n ); %PLOT FIGURES FOR RESULTS time_axis = [step : step : tend+tendnight]; %create a time vector including day and night time time_axis = time_axis./3600; %convert seconds into hours for readibility purpose a = zeros(1,length(time_axis))+PowReq; %creates a matrix with power requirement b = PowerOut.*2*(1-performance_margin); %multiply the power output of one TEG element by the number of TEG element (one on each side of the TM) TES_number_sizing = (a./b)*(1+design_margin); %Divide the power requirement by the total power output for one TES System. %It returns the number of TES System with this design which shall be implemented to satisfy the power requirement. fig = figure(10) %creates a figure to display all plot in one scree (compliance to MR6) subplot(2,3,1); plot(time_axis,TTOP,time_axis,TMID,time_axis,TBOT,time_axis,TTM); title('Temperature profile of the TM (K)'); legend('TTOP','TMID','TBOT','TTM'); xlabel('Time (hrs)'); ylabel('Temperature (K)'); grid on grid minor subplot(2,3,2); plot(time_axis,THOT,time_axis,TCOLD,time_axis,DeltaT); title('Hot and cold temperatures (K)'); legend('THOT','TCOLD','DeltaT'); xlabel('Time (hrs)'); ylabel('Temperature (K)'); grid on grid minor subplot(2,3,3); plot(time_axis,PowerOut); title('Power Output from one TEG element'); xlabel('Time (hrs)'); ylabel('Power (W)'); grid on grid minor subplot(2,3,4); plot(time_axis,PowerAbs,time_axis,PowerRej); title('Heat flow through each TEG element (W)'); legend('Power absorbed','Power rejected'); xlabel('Time (hrs)'); ylabel('Power (W)'); grid on grid minor subplot(2,3,5); plot(time_axis,effTEG); title('Efficiency of the TEG elements'); xlabel('Time (hrs)'); ylabel('Efficiency (-)'); grid on grid minor subplot(2,3,6); plot(time_axis,TES_number_sizing) title('Number of TESS needed including margins'); xlabel('Time (hrs)'); ylabel('Number of TESS (-)'); ylim([0 20000]); grid on grid minor toc