Notebook

In [1]:

from TRIOMA.tools.Extractors.PAV import Component,Fluid,Membrane, Geometry
import numpy as np 
import TRIOMA.tools.correlations as corr
import TRIOMA.tools.materials as materials
import matplotlib.pyplot as plt
from scipy.interpolate import griddata
import matplotlib.colors as colors
from matplotlib.colors import LogNorm
import matplotlib.lines as mlines
from TRIOMA.tools.BreedingBlanket import BreedingBlanket

This analysis is a sweep on two parameters, in which the efficiency evaluated with the nodal model and the analytical efficiency are compared together, and the relative error in percentage is given. First we define the d_hyd_v, the first sweep vector. This is going to be sweeped against all others (otherwise the analysis would be very long if we sweep all vectors against each other). Then we define a vector of vectors to sweep. A str_v_vec is used in the "update_attribute" method to indicate which attribute is going to be updated, corresponding with the vector that is sweeped in the loop. The Fluid_v_bool vector is used in this code to check if the update attribute must be used for the Component class or the Component.fluid class. The same happens for the solid_v_bool. Color vectors will be used after for plotting contours. Eff_v_vec is used then to store the results.

In [2]:

##Define sweep vectors
color_vector = ["red", "blue", "green", "yellow", "purple", "orange"]
N_vec=50
d_hyd_v=np.logspace(np.log10(20.5E-3),np.log10(20E-2),N_vec)
T_vec=np.linspace(800,1000,N_vec)
U0_vec=np.linspace(1,5,N_vec)
# str_v_vec=['T','U0','Solubility','c_in',"thick","K_S","k_d","k_r"]
variables = {'T' : 'Temperature [K]', 'U0' : 'Velocity [m/s]', 'Solubility' : 'Solubility [mol/m^3]', 'c_in' : 'Concentration [mol/m^3]', 'thick' : 'Thickness [m]', 'K_S' : 'Partition Coefficient', 'k_d' : 'Desorption Rate [1/s]', 'k_r' : 'Reaction Rate [1/s]'}
c_in_vec=np.logspace(-6,-1,N_vec)
D_vec=np.logspace(-10,-8,N_vec)
thick_vec=np.logspace(-3,-1,N_vec)
K_S_vec=np.logspace(-3,-1,N_vec)
k_d_vec=np.logspace(2,6,N_vec)
k_r_vec=np.logspace(2,6,N_vec)
solubility_vec=np.logspace(-3,-1,N_vec)
v_vec=np.array([T_vec, U0_vec, solubility_vec, c_in_vec,  thick_vec, K_S_vec])
eff_v_vec=np.array([])
fluid_v_bool=np.array([True , True , True ,False,False,False,False,False,False])
solid_v_bool=np.array([False, False, False,False,True ,True ,True ,True ,True ])

Here variables for the Breeding Blanket class are defined.

In [3]:

#Define other HX constraints

T_hot_prim=900
T_hot_sec=838
T_cold_prim=800
T_cold_sec=581
T_sec_ave=(T_hot_sec+T_cold_sec)/2
rho_sec=2263.628-0.636*T_sec_ave
mu_sec=0.075439-2.77E-4*(T_sec_ave-273.15)+3.49E-7*(T_sec_ave-273.15)**2-1.474E-10*(T_sec_ave-273.15)**3
k_sec=0.45
cp_sec=1396.044+0.172*(T_sec_ave)
N_HX=3
Q=1E9
m_in_sec=Q/(cp_sec*(T_hot_sec-T_cold_sec))

Empty arrays to store results with the append method are defined, and the double sweep takes place. In each iteration, the same component is defined and then one attribute is changed according to the sweep. Then color map plots are displayed, with additional isovariable contours.

In [4]:

from turtle import color

from matplotlib import contour
from matplotlib.pyplot import sca
from pyparsing import line


for j,vec in enumerate(v_vec):
        eff_v=np.array([])
        var_str=list(variables.keys())[j]
        T=800
        res_vec=np.array([])
        d_hyd_v_res=np.array([])
        var_vec=np.array([])
        L_vec=np.array([])
        n_pipes_v=np.array([])
        H_v=np.array([])
        W_v=np.array([])
        power_v=np.array([])
        for var in vec:
                for i,d_hyd in enumerate(d_hyd_v):
                        if var_str=='T':
                                T=var
                        mat=materials.Flibe(T)
                        flibe=Fluid(T=T, Solubility=mat.Solubility, MS=False,D=mat.D, d_Hyd=d_hyd ,mu=mat.mu,rho=mat.rho,U0=1,k=mat.k,
                         cp=mat.cp)
                        BB=BreedingBlanket(Q=Q,TBR=1.08,T_out=T_hot_prim,T_in=T_cold_prim,fluid=flibe, c_in=0 )
                        BB.get_cout()
                        c0=BB.c_out
                        Steel = Membrane( T=T,
                                # D=7E-10,
                                D=1E-4,
                                thick=2.1E-3,
                                K_S=4.41E-3,
                                k_d=1E6,
                                k_r=1E6,k=21)
                        geom_PAV=Geometry(thick=2.1E-3, D=d_hyd,L=1)
                        PAV = Component(c_in=c0, geometry=geom_PAV, fluid=flibe, membrane=Steel)
                        if fluid_v_bool[j]==True:
                                PAV.fluid.update_attribute(var_str,var)
                        elif solid_v_bool[j]==True:
                                PAV.membrane.update_attribute(var_str,var)
                        else:
                                PAV.update_attribute(var_str,var)
                        n_pipes=BB.m_coolant/(PAV.fluid.rho*PAV.fluid.U0*PAV.fluid.d_Hyd**2/4)
                        PAV.geometry.n_pipes=n_pipes
                        Q_HX=Q/n_pipes/N_HX
                        HX=Component(c_in=1, eff=0.08, fluid=flibe, membrane=Steel, geometry=Geometry(thick=2.1E-3, D=1E-4, L=1))
                        PAV.get_adimensionals()
                        d_h_sec=2.5E-2 # from HX forecasts
                        V_sec=m_in_sec/(rho_sec*d_h_sec**2*np.pi/4)/N_HX/n_pipes
                        Re_sec=corr.Re(rho=rho_sec,u=V_sec,mu=mu_sec,L=d_h_sec)
                        Pr_sec=corr.Pr(c_p=cp_sec,mu=mu_sec,k=k_sec)
                        h_coeff_sec=corr.get_h_from_Nu(corr.Nu_DittusBoelter(Re_sec, Pr_sec), k_sec,d_h_sec)
                        R_sec=1/h_coeff_sec
                        U = PAV.get_global_HX_coeff(0)
                        L= corr.get_length_HX(corr.get_deltaTML(T_hot_prim, T_cold_prim, T_cold_sec, T_hot_sec), PAV.fluid.d_Hyd, PAV.U, Q_HX)
                        PAV.geometry.update_attribute('L',L)
                        PAV.get_efficiency(plotvar=False)
                        PAV.analytical_efficiency()
                        out_flux=(PAV.c_in*(1-PAV.eff)*PAV.fluid.U0*PAV.fluid.d_Hyd**2/4)
                        eff_v=np.append(eff_v, abs(PAV.eff-PAV.eff_an)/PAV.eff_an)
                        d_hyd_v_res=np.append(d_hyd_v_res, d_hyd)
                        var_vec=np.append(var_vec, var)
                        # res_vec=np.append(res_vec, PAV.eff)
                        res_vec=np.append(res_vec, abs(PAV.eff-PAV.eff_an)/PAV.eff_an)
                        L_vec=np.append(L_vec, L)
                        n_pipes_v=np.append(n_pipes_v, n_pipes)
                        H_v=np.append(H_v, PAV.H)
                        W_v=np.append(W_v, PAV.W)
                        PAV.get_pumping_power()  
                        power_v=np.append(power_v, PAV.pumping_power)
                eff_v_vec=np.append(eff_v_vec, eff_v)
        plt.figure(j)
        plt.title(list(variables.values())[j])
        plt.yscale('log')
        plt.xscale('log')
        x = np.logspace(np.log10(min(d_hyd_v_res[:])), np.log10(max(d_hyd_v_res[:])), num=100)
        y = np.logspace(np.log10(min(var_vec[:])), np.log10(max(var_vec[:])), num=100)
        X, Y = np.meshgrid(x, y)
        Z = griddata((d_hyd_v_res, var_vec), (res_vec)*100, (X, Y), method='nearest')
        dZdX, dZdY = np.gradient(Z, x, y, edge_order=2)
        ZL=griddata((d_hyd_v_res, var_vec), L_vec, (X, Y), method='cubic')
        Zn_pipes=griddata((d_hyd_v_res, var_vec), n_pipes_v, (X, Y), method='cubic')
        ZH=griddata((d_hyd_v_res, var_vec), H_v, (X, Y), method='cubic')
        ZW=griddata((d_hyd_v_res, var_vec), W_v, (X, Y), method='cubic')
        Zratio=griddata((d_hyd_v_res, var_vec), H_v/W_v, (X, Y), method='cubic')
        Zpower=griddata((d_hyd_v_res, var_vec), power_v, (X, Y), method='cubic')
        contour1=plt.contour(X, Y, Z, levels=[1E-4,1E-3,1E-2,1E-1,1,10],norm=LogNorm(),colors=color_vector[0])
        contour2=plt.contour(X, Y, ZL, levels=range(1,30,1),colors=color_vector[1],linestyles='dotted')
        contour3=plt.contour(X, Y, Zn_pipes, levels=range(500,3000,500),colors=color_vector[2],linestyles='dashed')
        contour4=plt.contour(X, Y, ZH, levels=[0.01,0.2,100],colors=color_vector[3],linestyles='dotted')
        contour5=plt.contour(X, Y, ZW, levels=[0.1,0.2,10,20],colors=color_vector[4],linestyles='dotted')
        contour6=plt.contour(X, Y, Zratio, levels=[0.0001,0.2,1000],colors=color_vector[5],linestyles='dotted')
        contour7=plt.contour(X, Y, Zpower, levels=[0.5E4,1E4,1.5E4,2E4],colors='black',linestyles='dotted')
        scatter=plt.scatter(X,Y, c=Z ,norm=colors.LogNorm()) #,norm=colors.LogNorm() 
        plt.clabel(contour1, inline=True, fontsize=8)
        plt.clabel(contour2, inline=True, fontsize=8)
        plt.clabel(contour3, inline=True, fontsize=8)
        plt.clabel(contour4, inline=True, fontsize=8)
        plt.clabel(contour5, inline=True, fontsize=8)
        plt.clabel(contour6, inline=True, fontsize=8)
        plt.clabel(contour7, inline=True, fontsize=8)
        if scatter:
                if np.nanmin(scatter.get_array()) < np.nanmax(scatter.get_array()):
                        cbar = plt.colorbar(scatter)  # Show color scale for scatter
                        cbar.set_label('Relative error')
                else:
                        print("Cannot create colorbar: data does not have a valid range of values.")
        line1 = mlines.Line2D([], [], color=color_vector[0], markersize=15, label='Error %')
        line2 = mlines.Line2D([], [], color=color_vector[1], markersize=15, label='Length HX(m)')
        line3 = mlines.Line2D([], [], color=color_vector[2], markersize=15, label='number of pipes')
        line4 = mlines.Line2D([], [], color=color_vector[3], markersize=15, label='Fluid/Surface')
        line5 = mlines.Line2D([], [], color=color_vector[4], markersize=15, label='Diffusion/Surface')
        line6= mlines.Line2D([], [], color=color_vector[5], markersize=15, label='Fluid/Diffusion')
        line7= mlines.Line2D([], [], color='black', markersize=15, label='Pumping Power')
        plt.legend(handles=[ line1,line2, line3], loc='upper right')
        plt.legend(handles=[ line1,line2, line3,line4,line5,line6,line7], bbox_to_anchor=(1.3, 1), loc='upper left')
        plt.show()

Plots indicate very low relative error between the methods. Discontinuites in the plot indicate different if branches in the axial method (eg: fluid/diffusion>1000 in the solubility sweep), or different techniques to evaluate Sherwood (based on Reynolds number)