diff --git a/.github/workflows/formatting.yml b/.github/workflows/formatting.yml index 480c3bf..2f02d52 100644 --- a/.github/workflows/formatting.yml +++ b/.github/workflows/formatting.yml @@ -1,5 +1,5 @@ name: Check code formatting -on: [push, pull_request] +on: [push, pull_request, workflow_dispatch] jobs: check-formatting: runs-on: ubuntu-latest diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 6b17bca..6ec1043 100644 --- a/.github/workflows/test.yml +++ b/.github/workflows/test.yml @@ -1,5 +1,5 @@ name: tests -on: [push, pull_request] +on: [push, pull_request, workflow_dispatch] jobs: test-linux: runs-on: ubuntu-latest diff --git a/requirements.txt b/requirements.txt index d599151..731d36b 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,5 +1,5 @@ numpy -scipy==1.9.2 +scipy>=1.9.2 h5py vtk pytest diff --git a/srlife/damage.py b/srlife/damage.py index 4b9053e..56043e5 100644 --- a/srlife/damage.py +++ b/srlife/damage.py @@ -1033,7 +1033,7 @@ def calculate_volume_flaw_flattened_eq_stress( else: sigma_e[sigma_e < 0] = 0 g = ( - np.trapz( + np.trapezoid( (sigma_e / (sigma_e_max + self.tolerance)) ** Navg[..., None, None], time, axis=0, @@ -1196,7 +1196,7 @@ def calculate_surface_flaw_flattened_eq_stress( else: sigma_e[sigma_e < 0] = 0 g = ( - np.trapz( + np.trapezoid( (sigma_e / (sigma_e_max + self.tolerance)) ** Navg[..., None, None, None], time, @@ -1474,7 +1474,7 @@ def calculate_volume_flaw_element_log_reliability( # Calculating ratio of cyclic stress to max cyclic stress (in one cycle) # for time-independent and time-dependent cases g = ( - np.trapz( + np.trapezoid( (pstress / (pstress_max + 1.0e-14)) ** Navg[..., None], time, axis=0, @@ -1556,7 +1556,7 @@ def calculate_surface_time_dep_pstress( # Calculating ratio of cyclic stress to max cyclic stress (in one cycle) # for time-independent and time-dependent cases g = ( - np.trapz( + np.trapezoid( (surf_pstress / (surf_pstress_max + 1.0e-14)) ** Navg[..., None, None], time, axis=0, @@ -1743,7 +1743,7 @@ def calculate_volume_flaw_avg_normal_stress( sigma_n_0 = sigma_n else: g = ( - np.trapz( + np.trapezoid( (sigma_n / (sigma_n_max + self.tolerance)) ** Navg[..., None, None], time, axis=0, @@ -1890,7 +1890,7 @@ def calculate_surface_flaw_time_dep_normal_stress( sigma_n_0 = sigma_n else: g = ( - np.trapz( + np.trapezoid( (sigma_n / (sigma_n_max + self.tolerance)) ** Navg[..., None, None, None], time, diff --git a/srlife/data/damage/SiC.xml b/srlife/data/damage/SiC.xml index 64e6f87..f75fc6a 100644 --- a/srlife/data/damage/SiC.xml +++ b/srlife/data/damage/SiC.xml @@ -46,12 +46,12 @@ - 298.15 1073.15 1273.15 1473.15 1673.15 1773.15 - 0.0 0.0 0.0 0.0 0.0 0.0 + 298.15 811.15 1089.15 1366.15 1566.15 + 0.0 0.0 0.0 0.0 0.0 - 298.15 1073.15 1273.15 1473.15 1673.15 1773.15 - 0.0 0.0 0.0 0.0 0.0 0.0 + 298.15 811.15 1089.15 1366.15 1566.15 + 0.0 0.0 0.0 0.0 0.0 298.15 811.15 1089.15 1366.15 1566.15 diff --git a/srlife/interface.py b/srlife/interface.py new file mode 100644 index 0000000..3ac0a95 --- /dev/null +++ b/srlife/interface.py @@ -0,0 +1,1061 @@ +""" + This module has all of the helper functions needed to run an srlife analysis. + It is meant to allow a simple script interface to define a problem and run the + analysis. + + One such example is defined in test/test_interface.py +""" + +import numpy as np +from scipy.interpolate import RegularGridInterpolator + +# srlife specific functions +from srlife import ( + receiver, + solverparams, +) + + +# UNIT CONVERSION FUNCTIONS +def convert_mm_to_m(mm_in): + """ + Unit conversion from mm to m + + Args: mm_in (double): qty to convert + """ + return mm_in / 1000 + + +def convert_m_to_mm(m_in): + """ + Unit conversion from m to mm + + Args: m_in (double): qty to convert + """ + return m_in * 1000 + + +def convert_kWm2_to_Wmm2(kWm2_in): + """ + Unit conversion from kW/m^2 to W/mm^2 + + Args: kWm2_in (double): qty to convert + """ + return kWm2_in * 1000e-6 + + +def convert_Pa_to_MPa(Pa_in): + """ + Unit conversion from Pa to MPa + + Args: Pa_in (double): qty to convert + """ + return Pa_in * 1e-6 + + +def convert_MPa_to_Pa(MPa_in): + """ + Unit conversion from MPa to Pa + + Args: MPa_in (double): qty to convert + """ + return MPa_in * 1e6 + + +def convert_C_to_K(C_in): + """ + Unit conversion from C to K + + Args: C_in (double): qty to convert + """ + return C_in + 273.15 + + +def read_month_day_hour_flux_file(flux_file_name, data_shape, z_offset_in): + """ + Read and process the flux data from a solarPILOT file + + Args: + flux_file_name (string): path and name of flux file from solarpilot + data_shape (list[int]): list with n_rows and n_cols of flux data + z_offset_in (double): z offset between srlife model and solarpilot model + this will shift z data that is read up or down + + Returns: + z_data (np.array(double)): vector with z positions from flux sampling + theta_data (np.array(double)): vector with theta positions from flux sampling + flux_data (np.array(double)): matrix with flux data at sampled positions + """ + # read flux data + # massage into format + # output + row, col = data_shape + z_data = np.zeros([row, 1]) + theta_data = np.zeros([1, col]) + flux_data = np.zeros([row, col]) + try: + theta_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=6, usecols=range(1, col + 1) + )[0] + z_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=7, usecols=range(0, 1) + ) + z_data += z_offset_in + flux_data = np.loadtxt( + flux_file_name, delimiter=",", skiprows=7, usecols=range(1, col + 1) + ) + # Data is on cylinder, so 180 will repeat + theta_data = np.append(theta_data, -180) + flux_data = np.hstack((flux_data, flux_data[:, 0][:, None])) + except FileNotFoundError: + # if we cannot find the file + # print a message and just return all zeros + print( + f"File could not be found!!! Returning all zeros\n filename: {flux_file_name}" + ) + # convert height to mm + z_data = convert_m_to_mm(z_data) + return z_data, theta_data, flux_data + + +def get_flux_interpolators_from_data_files( + times, start_time, end_time, month, day, flux_data_dir, flux_data_shape, z_offset +): + """ + Read flux files and create interpolator functions for each month, day, hour + + Args: + times (list(int)): list of index of hours used for analysis + Could probably remove this if (end_time - start_time) == len(times) + start_time (int): first hour to read flux data + end_time (int): last hour to read flux data + month (int): month to read flux data + day (int): day to read flux data + flux_data_dir (string): path to flux data files + flux_data_shape (list(int)): n_rows and n_cols of flux data in flux files + z_offset_in (double): z offset between srlife model and solarpilot model + this will shift z data that is read up or down + + Returns: + flux_interpolators_by_hour (list(RegularGridInterpolators)): list with the interpolating + functions for each hour in times. + """ + print("Making flux interpolators from data") + # Time zero is NOT analyzed, so we do not make a fluxFn for that + flux_interpolators_by_hour = [0] * (len(times) - 1) + for i_count, i_time in enumerate(range(start_time, end_time + 1)): + flux_filename = f"{flux_data_dir}/data_{month}_d{day}_hr{int(i_time)}.csv" + z_pts, theta_pts, flux_pts = read_month_day_hour_flux_file( + flux_filename, flux_data_shape, z_offset + ) + interpFun = RegularGridInterpolator( + (z_pts, theta_pts), + flux_pts, + method="linear", + bounds_error=False, + fill_value=0.0, + ) + flux_interpolators_by_hour[i_count] = interpFun + return flux_interpolators_by_hour + + +def apply_tube_flux_bcs( + tube, + tube_absorbance, + tube_rec_theta, + flux_interpolators_by_hour, +): + """ + Actually apply the heat flux to the tube across width and length + + Args: + tube (srlife.Tube): srlife object for this tube + tube_z (list[double]): list of z coordinates of tube nodes + tube_theta (list[double]): list of theta coordinates of tube nodes + tube_absorbance (double): absorbance of tubes + tube_rec_theta (double): receiver theta value for this tube + flux_interpolators_by_hour (list[RegularGridInterpolator]): + list of interpolating functions for heat flux on tubes + """ + num_steps = len(flux_interpolators_by_hour) + tube_theta = np.linspace(0, 360, tube.nt + 1)[:-1] + tube_z = np.linspace(0, tube.h, tube.nz) + # add extra time step for time zero + tube_flux = np.zeros([num_steps + 1, len(tube_theta), len(tube_z)]) + for i_hour in range(num_steps): + # for each hour of flux data + # sample tube at tube_rec_theta and z + for i_z, z in enumerate(tube_z): + # for each zPt + pt_z_theta = [z, tube_rec_theta] + rec_flux = convert_kWm2_to_Wmm2( + flux_interpolators_by_hour[i_hour](pt_z_theta) + ) + # apply across crown so only sunSide gets flux + for i_theta, theta in enumerate(tube_theta): + # for each thetaPt around tube + if 0 < theta < 180: + tube_flux[i_hour + 1, i_theta, i_z] = rec_flux[0] * np.cos( + (theta - 90) * np.pi / 180 + ) + else: + tube_flux[i_hour + 1, i_theta, i_z] = 0.0 + # make bc object for tube with heat flux data + # NOTE: this is done here and not in cerate_tube because if we initialize + # bc data to np.zeros, flux will always be zero unless we reload from a file + tube.set_bc(apply_heatflux_bc(tube, tube_flux * tube_absorbance, "outer"), "outer") + # while doing this, also apply convective BCs to inner face + tube.set_bc(apply_convective_bc(tube, "inner"), "inner") + + +def apply_convective_bc(tube, loc): + """ + Defines and applies a convective BC to a surface of a tube + specificed by loc + + Args: + tube (srlife.Tube): srlife object for this tube + loc (string): "inner" or "outer" to specify which surface BC is applied to + + """ + nz = tube.nz + if loc == "outer": + r = tube.r + else: + r = tube.r - tube.t + return receiver.ConvectiveBC( + r, tube.h, nz, tube.times, np.zeros((len(tube.times), nz)) + ) + + +def apply_heatflux_bc(tube, data, loc): + """ + Defines and applied a heat flux BC to a surface of a tube object + specified by loc + + Args: + tube (srlife.Tube): srlife object for this tube + data (list[double]): 3d list of heat flux bcs [hour, theta, z] + loc (string): "inner" or "outer" to specify which surface BC is applied to + + """ + if loc == "outer": + r = tube.r + else: + r = tube.r - tube.t + nt = tube.nt + nz = tube.nz + return receiver.HeatFluxBC(r, tube.h, nt, nz, tube.times, data) + + +def calc_and_write_tube_flux_bcs( + rec, + ass_tube_per_panel, + num_panels, + flux_interpolators_by_hour, + tube_absorbance, +): + """ + Applies flux boundary conditions to analyzed tubes. Flux value is interpolated based on + centerline tube location. It is applied along the tube face as a cosine distribution + with max at sun-facing centerline to zero at 90 deg in either direction. + + Args: + rec (Receiver): the receiver object that the tubes belong to + ass_tube_per_panel (int): the assumed number of tubes per panel, or analyzed number of tubes + num_panels (int): number of panels in the receiver #BPMToDo: remove this and get from rec + flux_interpolators_by_hour (list(RegularGridInterpolators)): list of flux interpolators + on receiver by hour of analysis + tube_z (list(double)): z nodes points along tube + tube_theta (list(double)): theta nodes points along tube + tube_absorbance (double): absorbance of tubes + write_tube_flux_to_csv (bool): flag to dump flux data per tube to a csv file to debug + + Returns: None + """ + print("Assigning flux BCs to tubes") + for panel_key, panel in rec.panels.items(): + # in each panel of rec + for tube_key, tube in panel.tubes.items(): + print(f"Panel: {panel_key} Tube: {tube_key}") + # in each tube of panel + tube_rec_theta = 180 - (360 * int(panel_key) / num_panels) + if ass_tube_per_panel == 1: + # if only one tube, place at center of panel + tube_rec_theta += 360 / (2 * num_panels) + else: + # if more than one tube, evenly space them including edges + tube_rec_theta -= ( + 360 / num_panels * int(tube_key) / (ass_tube_per_panel - 1) + ) + if tube_rec_theta == -180: + # wrap data around tube + tube_rec_theta = 180 + apply_tube_flux_bcs( + tube, + tube_absorbance, + tube_rec_theta, + flux_interpolators_by_hour, + ) + + +def set_rec_flow_paths(rec, panel_flow_path, mass_flow_per_path, T_in_per_path): + """ + Set flowpaths within a receiver object + + Args: + rec (Receiver): the receiver object to set flowpaths for + panel_flow_path (list[string]): panel names indicating the flowpath. Can have a list + of lists if more that one flowpath + mass_flow_per_path (list[double]): list with entries of initial mass flow rate for + each path + T_in_per_path (list[double]): list with entries of inlet temp for each path + """ + print("Setting receiver flowpaths") + T_in_per_path = convert_C_to_K(T_in_per_path) + # add flowpaths to receiver + for i_path, flow_path_panels in enumerate(panel_flow_path): + # for each path in given list + print(f"Path {i_path}") + flow_path = [0] * len(flow_path_panels) + for i_panel, panel in enumerate(flow_path_panels): + # for each panel in path + print(f"Panel: {panel}") + flow_path[i_panel] = panel + panel = rec.panels[panel] + for tube in panel.tubes.values(): + # for each tube in panel + times = tube.times + # initialize mass flow for path + mass_flow = mass_flow_per_path[i_path] + # convert mass flow from kg/s to kg/hr + mass_flow *= 3600 + # set inlet temp + T_in_flow_path = T_in_per_path[i_path] * np.ones_like(times) + rec.add_flowpath(flow_path, times, mass_flow, T_in_flow_path) + + +def check_outlet_temps_and_step_mass_flow( + rec, path, path_key, T_out_target, pct_err_outlet_temp, breaker +): + """ + Check outlet temp of all tubes from last panel in flowpath. + Average their value and use it to step mass flow rate. + + Args: + rec (srlife.Receiver): Receiver object who owns the flowpath we are + optimizing + path (OrderedDict): dictionary of flowpath information + path_key (string): name of flowpath + T_out_target (double): target outlet temperature + pct_err_outler_temp (double): acceptable error % for outlet temp (for convergence) + breaker (list[bool]): check whether all flowpaths have converged + """ + # All tubes in last panel go to manifold, so estimate + # temp as avg of all tubes output + panels = path.get("panels") + T_in_path = path["inlet_temp"][1:] + last_panel = rec.panels[panels[-1]] + num_tubes = last_panel.ntubes + num_time_steps = next(iter(last_panel.tubes.values())).ntime + T_out_path = np.zeros([num_tubes, num_time_steps - 1]) + for i_tube, tube in enumerate(last_panel.tubes.values()): + # for tubes in the last panel of the flow path + # results for tube are (time, nz) shape + # want all time except 0 and the last zPos + T_out_path[i_tube] = tube.axial_results.get("fluid_temperature")[1:, -1] + T_out_path = T_out_path.mean(axis=0) + print(f"T_out = {T_out_path - 273.15}") + if ( + any(T_out_path > T_out_target * (1 + pct_err_outlet_temp / 100)) is False + and any(T_out_path < T_out_target * (1 - pct_err_outlet_temp / 100)) is False + ): + # we have converged this flowpath + breaker[int(path_key)] = True + else: + mass_flow = path.get("mass_flow") + mass_flow[1:] = ( + mass_flow[1:] * (T_out_path - T_in_path) / (T_out_target - T_in_path) + ) + mass_flow[0] = mass_flow[1] + rec.flowpaths[path_key]["mass_flow"] = mass_flow + print(f"Mass flow: {mass_flow/3600}") + + +def optimize_mass_flow_rate_per_path( + rec, rec_filename, T_out_target, pct_err_outlet_temp, solver, save_heat_to_vtu +): + """ + Iteratively calculate the ideal mass flow rate for each path such that the + outlet temperature matches a target temperature. + + Args: + rec (Receiver): srlife Receiver object we are optimizing flow for + rec_filename (string): file name for the receiver object to optimize + T_out_target (double): target output temperature + pct_err_outlet_temp (double): acceptable percentage of error for the outlet temp + solver (managers.SolutionManager): srlife solver object + save_heat_to_vtu (bool): whether or not to save vtu files of results + """ + print("Optimizing mass flow rate to match outlet temp") + T_out_target = convert_C_to_K(T_out_target) + # make sure receiver is up to date + solver.receiver = rec + # Mass flow optimization + # limit of error on outlet temp + max_mass_flow_iter = 5 + for i_opt in range(max_mass_flow_iter): + print(f"Iteration = {i_opt}") + breaker = [False] * len(rec.flowpaths) + if i_opt > 0: + for path_key, path in rec.flowpaths.items(): + print(f"FlowPath: {path_key}") + # check last panel tube outlet temps + check_outlet_temps_and_step_mass_flow( + rec, path, path_key, T_out_target, pct_err_outlet_temp, breaker + ) + for panel_key, panel in rec.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"ht-tube-{panel_key}-{tube_key}") + + if all(breaker): + print(f"Converged!!! solved in {i_opt} iterations") + if save_heat_to_vtu: + for panel_key, panel in rec.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"ht-tube-{panel_key}-{tube_key}") + break + # solve heat transfer problem + solver.solve_heat_transfer() + filename = rec_filename + "_ht_iter_" + str(i_opt) + ".hdf5" + rec.save(filename) + + +def calc_fluid_velocity(tube_mass_flow, fluid_rho, tube_Dh): + """ + Calculate the velocity of the fluid in a tube based on mass flow rate + + Args: + tube_mass_flow (double): mass flow rate in tube + fluid_rho (double): fluid density + tube_Dh (double): hydraulic diameter of tube + """ + return tube_mass_flow / fluid_rho * (4.0 / (np.pi * tube_Dh**2)) + + +def calc_reynolds_number(fluid_rho, fluid_vel, fluid_mu, tube_Dh): + """ + Calculate the Reynolds number for tube flow + + Args: + fluid_rho (double): fluid density + fluid_vel (double): fluid velocity + fluid_mu (double): fluid dynamic viscosity + tube_Dh (double): tube hydraulic diameter + + """ + return fluid_rho * fluid_vel * tube_Dh / fluid_mu + + +def calc_fd(Re, tube_eta, tube_Dh): + """ + Calculate Darcy friction factor from: + Zigrange and Sylvester 1985, A review of explicit friction factor equations. + J of Energy Resources Technology) + Equation 13 + + Args: + Re (double): Reynolds number + tube_eta (double): tube roughness + tube_Dh (double): tube hydraulic diameter + Returns: + fd (double): darcy friciton factor + """ + if Re < 4000: + # laminar flow + fd = 64 / Re + else: + fd = ( + 1 + / ( + -2 + * np.log10( + (tube_eta / tube_Dh) / 3.7 + - 5.02 / Re * np.log10((tube_eta / tube_Dh) / 3.7 + 13 / Re) + ) + ) + ** 2 + ) + return fd + + +def calc_friction_p_loss(fluid_rho, fluid_mu, tube_mass_flow, tube_dict, delta_l): + """ + Calculate pressure loss from friction in tubes + + Args: + fluid_rho (double): fluid density in tubes + fluid_mu (double): fluid dynamic visco in tubes + tube_mass_flow (double): mass flow rate in tube + tube_dict (dict): dictionary of tube specs + delta_l (double): length of tube + + Returns: + friction_p_loss (double): pressure loss from friction + """ + # tube roughness + tube_eta = tube_dict["eta"] + tube_Dh = convert_mm_to_m(tube_dict["od"] - 2 * tube_dict["t"]) + fluid_vel = calc_fluid_velocity(tube_mass_flow, fluid_rho, tube_Dh) + Re = calc_reynolds_number(fluid_rho, fluid_vel, fluid_mu, tube_Dh) + fd = calc_fd(Re, tube_eta, tube_Dh) + return fd * (0.5 * fluid_rho) * (fluid_vel**2) / tube_Dh * delta_l + + +def calc_manifold_bend_loss( + max_mass_flow, manifold_dict, manifold_rho, manifold_mu, num_bends +): + """ + Calculate pressure loss from flow through bends in manifold tubes + + Args: + max_mass_flow (double): mass flow in manifold + manifold_dict (dict): dictionary with manifold tube specs + manifold_rho (double): fluid density in manifold + manifold_mu (double): dynamic viscosity in manifold + num_bends (int): total number of bends in manifold tubes + """ + manifold_eta = manifold_dict["eta"] + manifold_pipe_t = convert_mm_to_m(manifold_dict["t"]) + manifold_pipe_id = convert_mm_to_m(manifold_dict["od"] - 2 * manifold_pipe_t) + manifold_pipe_bend_radius = convert_mm_to_m(manifold_dict["bend_radius"]) + manifold_eta = manifold_dict["eta"] + manifold_Dh = manifold_pipe_id + if manifold_eta / manifold_Dh > 0.0001: + K0 = 0.42 + elif manifold_eta / manifold_Dh == 0: + K0 = 0.21 + else: + K0 = 0.21 * (1 + 1000 * manifold_eta / manifold_Dh) + + manifold_vel = calc_fluid_velocity(max_mass_flow, manifold_rho, manifold_Dh) + Re = calc_reynolds_number(manifold_rho, manifold_vel, manifold_mu, manifold_Dh) + manifold_fd = calc_fd(Re, manifold_eta, manifold_Dh) + K = ( + K0 ** ((manifold_pipe_id / manifold_pipe_bend_radius) ** 0.5) + + 0.5 * np.pi * (manifold_pipe_bend_radius / manifold_pipe_id) * manifold_fd + ) + return num_bends * K * manifold_rho * manifold_vel**2 + + +def calc_path_total_p_loss( + rec, + path, + path_key, + delta_T, + tube_dict, + manifold_dict, + fluid, + rec_circumference, + num_panels, +): + """ + Calculate the total pressure loss in a flowpath from flow friction, + head loss, and flow through pipe bends. + + Args: + rec (Receiver): receiver object to analyze + path (dict): flowpath to get pressure loss through + path_key (string): name of flowpath + delta_T (double): temp step through flowpath + tube_dict (dict): dictionary of usefule tube info + manifold_dict (dict): dict of useful info for manifold tubes + fluid (thermalfluid.ThermalFluidMaterial): working fluid in receiver + rec_circumference (double): circumference of receiver + num_panels (int): number of panels in receiver + + Returns: + total_p_loss (double): total pressure loss in flowpath + """ + g = 9.81 # m/s^2 + # NOTE: this function works mostly with meters + act_tube_per_panel = tube_dict["ass_tube_per_panel"] * tube_dict["tube_mult"] + num_bends = num_panels * manifold_dict["bends_per_panel"] + + inlet_temp = path["inlet_temp"][0] + # get outlet avg outlet temp from each tube in last panel + # then average down again to single value + outlet_temp = np.mean( + [ + tube.axial_results.get("fluid_temperature")[:, -1][1:-1].mean() + for tube_key, tube in rec.panels[path["panels"][-1]].tubes.items() + ] + ) + max_mass_flow = path["mass_flow"].max() / 3600 # kg/s + tube_mass_flow = max_mass_flow / act_tube_per_panel + flow_path_length = convert_mm_to_m(tube_dict["h"]) * len(path["panels"]) + temp_steps = np.arange(inlet_temp, outlet_temp, delta_T) + rho_with_T = fluid.rho(temp_steps) * 1e9 # kg/m^3 + mu_with_T = fluid.mu(temp_steps) * 1000 / 3600 # Pa-s + + # Friction Loss + # pos spacing of temp change on tubes + delta_l = flow_path_length / np.size(temp_steps) + friction_p_loss = sum( + calc_friction_p_loss(rho, mu, tube_mass_flow, tube_dict, delta_l) + for rho, mu in zip(rho_with_T, mu_with_T) + ) + # Head Loss + avg_temp = 0.5 * (inlet_temp + outlet_temp) + head_p_loss = fluid.rho(avg_temp) * g * convert_mm_to_m(tube_dict["h"]) + # Manifold Loss + manifold_friction_p_loss = calc_friction_p_loss( + fluid.rho(avg_temp) * 1e9, + fluid.mu(avg_temp) * 1000 / 3600, + max_mass_flow, + manifold_dict, + (len(path["panels"]) - 1) * (rec_circumference / num_panels), + ) + manifold_bend_loss = calc_manifold_bend_loss( + max_mass_flow, + manifold_dict, + fluid.rho(avg_temp) * 1e9, + fluid.mu(avg_temp) * 1000 / 3600, + num_bends, + ) + return friction_p_loss + head_p_loss + manifold_bend_loss + manifold_friction_p_loss + + +def calc_p_loss_from_flows_temps(rec, tube_dict, manifold_dict, fluid, outlet_p): + """ + Calculate the pressure loss across the receiver based on the results + of a flow simulation. This captures friction loss, head loss, and bend losses + in the tubes and manifold pipes of the receiver + + Args: + rec (Receiver): the receiver object to calculate pressure loss + tube_dict (dict): a dictionary containing info about rec tubes + manifold_dict (dict): a dictionary containing info about manifold tubes + fluid (thermalfluid.ThermalFluidMaterial): thermal material of working fluid + outlet_p (double): pressure at flowpath outlets + """ + print("Calculating pressure loss from flow at temps") + # receiver details + num_panels = len(rec.panels) + # back out rec_diam + rec_diam = ( + tube_dict["ass_tube_per_panel"] + * tube_dict["tube_mult"] + * num_panels + * (tube_dict["od"] + tube_dict["spacing"]) + / np.pi + ) + rec_circumference = np.pi * convert_mm_to_m(rec_diam) + print(f"Rec diam = {convert_mm_to_m(rec_diam)}") + # This is number of tempsteps along tube + delta_T = 5 + flow_path_p_loss = np.zeros(len(rec.flowpaths)) + for path_key, path in rec.flowpaths.items(): + # for each flow path in rec + print(f"Flow path: {path_key}") + total_p_loss = calc_path_total_p_loss( + rec, + path, + path_key, + delta_T, + tube_dict, + manifold_dict, + fluid, + rec_circumference, + num_panels, + ) + flow_path_p_loss[int(path_key)] = total_p_loss + print(f"Total pressure loss= {(total_p_loss/1e6)} MPa") + print(f"Inlet pressure= {(outlet_p + total_p_loss/1e6)} MPa") + return flow_path_p_loss + + +def update_tube_pressure_bcs(rec, inlet_p_per_path, outlet_p): + """ + Set tube pressure BCs based on pressure loss calc after thm solve + + Args: + rec (Receiver): receiver object to set the pressure for + inlet_p_per_path (np.array): array of inlet pressures for each flowpath + outlet_p (double): outlet pressure, same for all paths + """ + # set tube pressures and temperatures + for path_key, path in rec.flowpaths.items(): + # for each flowpath in the model + tube_pressures = np.linspace( + inlet_p_per_path[int(path_key)], outlet_p, len(path["panels"]) + 1 + )[:-1] + for i_panel, panel_key in enumerate(path["panels"]): + # for each panel in flowpath + for tube in rec.panels[panel_key].tubes.values(): + # for each tube in panel + # get tube times + times = tube.times + pressure = np.ones_like(times) + pressure[0] = 0.0 + tube_pressure = tube_pressures[i_panel] * pressure + tube_pressure_bc = receiver.PressureBC(times, tube_pressure) + tube.set_pressure_bc(tube_pressure_bc) + + +def cycle_tube_pressure_bcs(tubes_dict, num_cycles, cyclic_times): + """ + Set pressure bcs on internal face of tube from flow + + Args: + tubes_dict (dict): rec.panels[panel].tubes dictionary tubes from a panel + tube_pressure (list[double]): pressure values for tubes in a given panel with height + pressure (list[double]): properly sized array of pressure data to define bc + times (list[int]): list of analysis times + + """ + for tube in tubes_dict.values(): + # NOTE: as soon as I cycle results + # ghost temp, fluid temp, fluid velocity results are no longer good + # so lets delete themo + try: + tube.quadrature_results.pop("ghost_temperature") + tube.axial_results.pop("fluid_temperature") + tube.axial_results.pop("fluid_velocity") + except KeyError: + # We do not need to do anything here, results are already gone + pass + # for each tube in panel + press_bc = tube.pressure_bc + pressure = press_bc.data + tube_pressure = np.tile(pressure[1:], num_cycles) + tube_pressure = np.append(0, tube_pressure) + tube_pressure_bc = receiver.PressureBC(cyclic_times, tube_pressure) + tube.set_pressure_bc(tube_pressure_bc) + + +def set_and_downsample_tube_temp_bcs( + tubes_dict, + set_init_T_to_inlet_T, + inlet_T, + cyclic_times, + num_cycles, + analysis_type, + loc, +): + """ + Handle specific details for temperature BCs in a tube. + Set the initial temperature + Find and set the appropriate temp for 2d analysis + + Args: + tubes_dict (dict): rec.panels[panel].tubes dictionary tubes from a panel + set_init_T_to_inlet_T (bool): whether to initialize tube temps as the + inlet temperature + inlet_T (double): inlet temperature for tubes + cyclic_times (list[int]): list of times across all cycles of analysis + num_cycles (int): number of cycles in analysis + analysis_type (string): "2d" or "3d" to handle turn on or off downsampling + of temperature results for simpler analysis + loc (string): "max_T" or "max_avg_T" to set where temp is sampled in 2d analysis + """ + for tube in tubes_dict.values(): + # for each tube in panel + T = tube.results["temperature"] + if set_init_T_to_inlet_T: + T[0] = inlet_T + tube.T0 = inlet_T + _, _, _, r4 = np.shape(T) + # this will tile stack all temp results num_cycles times + # It does not account for cycle heuristic + T_0 = np.array([T[0]]) + T = np.tile(T[1:], (num_cycles, 1, 1, 1)) + T = np.append(T_0, T, axis=0) + tube.results["temperature"] = T + tube.set_times(cyclic_times) + T_3d = tube.results["temperature"] + if analysis_type == "2d": + if loc == "max_T": + T_max = 0 + T_max_h_index = -1 + for i_temp in range(1, r4 - 1): + if np.max(T_3d[:, :, :, i_temp]) > T_max: + T_max = np.max(T_3d[:, :, :, i_temp]) + T_max_h_index = i_temp + tube.make_2D(tube.h / (r4 - 1) * T_max_h_index) + elif loc == "max_avg_T": + T_3d_avg = np.average(T_3d, axis=0) + T_max = 0 + T_max_h_index = -1 + for i_temp in range(1, r4 - 1): + if np.max(T_3d_avg[:, :, i_temp]) > T_max: + T_max = np.max(T_3d_avg[:, :, i_temp]) + T_max_h_index = i_temp + tube.make_2D(tube.h / (r4 - 1) * T_max_h_index) + tube.results["temperature"] = T_3d[:, :, :, T_max_h_index] + + +def process_single_panel_analysis(rec, struct_output_dict, single_panel_analysis_id): + """ + Remove all panels that will not be analyzed + + Args: + rec (Receiver): object we are analyzing + struct_output_dict (dict): dictionary with output information + save_struct_to_vtu (bool): whether or not to write files to vtu + st_fname (string): filename for structure + tube_fname (string): filename for tube + single_panel_analysis_is (string): name of panel to analyze + """ + # file name mods + st_fname = ( + f"sing_panel_{single_panel_analysis_id}_{struct_output_dict['st_filename']}" + ) + struct_output_dict["st_filename"] = st_fname + tube_fname = ( + f"sing_panel_{single_panel_analysis_id}_{struct_output_dict['tube_filename']}" + ) + struct_output_dict["tube_filename"] = tube_fname + # update model to use just one panel + single_panel_model = rec.panels[single_panel_analysis_id] + rec.panels.clear() + rec.flowpaths.clear() + rec.add_panel(single_panel_model) + + +def save_structural_results(rec_struct, st_fname, tube_fname): + """ + Save results from analysis to files + + Args: + rec_struct (Receiver): receiver object we analyzed + st_fname (string): hdf5 filename for structure + tube_fname (string): base filename for tube vtus + """ + # post process + rec_struct.save(st_fname + ".hdf5") + for panel_key, panel in rec_struct.panels.items(): + for tube_key, tube in panel.tubes.items(): + tube.write_vtk(f"{tube_fname}_{panel_key}_{tube_key}") + + +def run_struct_analysis( + rec_filename, + set_init_T_to_inlet_T, + loc, + analysis_type, + is_single_panel_analysis, + single_panel_analysis_id, + num_cycles, + solver, + struct_output_dict, +): + """ + Runs the structural analysis on the receiver with a given set of thermal + hydraulics results. + + Args: + rec_filename (string): filename of the receiver to analyze + set_init_T_to_inlet_T (bool): whether to initialize tube temps as the + inlet temperature + loc (string): "max_T" or "max_avg_T" used to select location to downsample + analysis domain if lower dimensional analysis is done + analysis_type (string): "2d" or "3d" used to decide whether full analysis + or downsampled analysis is used + is_single_panel_analysis (bool): whether a single panel is analyzed or the + full receiver + single_panel_analysis_id (string): name of single panel to analyze if + is_single_panel_analysis == True + num_cycles (double): number of times to repeat load cycles for the analysis + solver (managers.SolutionManager): system solver for the receiver + struct_output_dict (dict): dictionary with info about output for st files + save_to_vtu (bool): save files to vtu output + st_filename (string): hdf5 filename for structural output + tube_filename (string): filename for vtk output of tube analysis + + ToDos: + rename loc to better name, make it and analysis_type bools? + """ + # passing filename here because we may change receiver significantly + # i.e. removing many panels + rec_struct = receiver.Receiver.load(rec_filename + ".hdf5") + rec_struct.days *= num_cycles + inlet_T = rec_struct.flowpaths["0"]["inlet_temp"][0] + + # set tube pressures and temperatures + for path in rec_struct.flowpaths.values(): + # for each flowpath in the model + for panel_key in path["panels"]: + # for each panel in flowpath + # times we actually analyze receiver for thm + first_tube = next(iter(rec_struct.panels[panel_key].tubes.values())) + analysis_times = first_tube.times[1:] + # create cyclic time for life + cyclic_times = np.tile(analysis_times, num_cycles) + for i_cyc in range(num_cycles): + cyclic_times[ + i_cyc * len(analysis_times) : (i_cyc + 1) * len(analysis_times) + ] += (i_cyc * rec_struct.period) + # initial steps for pressure and time + cyclic_times = np.append([0], cyclic_times) + pressure = np.ones_like(cyclic_times) + pressure[0] = 0.0 + + cycle_tube_pressure_bcs( + rec_struct.panels[panel_key].tubes, + num_cycles, + cyclic_times, + ) + set_and_downsample_tube_temp_bcs( + rec_struct.panels[panel_key].tubes, + set_init_T_to_inlet_T, + inlet_T, + cyclic_times, + num_cycles, + analysis_type, + loc, + ) + # remove unnecessary panels if needed + if is_single_panel_analysis: + process_single_panel_analysis( + rec_struct, struct_output_dict, single_panel_analysis_id + ) + solver.receiver = rec_struct + # run structural problem + solver.solve_structural() + if struct_output_dict["save_to_vtu"]: + save_structural_results( + rec_struct, + struct_output_dict["st_filename"], + struct_output_dict["tube_filename"], + ) + rec_struct.save(rec_filename + "_struct.hdf5") + + +def create_receiver(tube_dict, num_days, times, period, panel_k, num_panels, results): + """ + Create an srlife.Receiver object based on user defined inputs + + Args: + tube_dict (dict): Python dictionary that holds useful tube definitions + "od", "t", "h", "nr", "nt", "nz", "spacing", "T0", "tube_k", + "eta", "tube_mult", "ass_tube_per_panel" + num_days (int): number of days to analyze + times (list[int]): list of index of hours to analyze + period (int): number of hours to analyze in a day + panel_k (double): panel stiffness used in SystemSolver object + num_panels (int): number of panels + results (list): set of results for tube objects + """ + # array of 1's with leading and trailing zeros multiplied by input pressure + rec = receiver.Receiver(period, num_days, panel_k) + for _ in range(num_panels): + rec.add_panel(create_panel(tube_dict, times, results)) + return rec + + +def create_panel(tube_dict, times, results): + """ + Create an srlife.Panel object to add to a receiver + + Args: + tube_dict (dictionary): Python dictionary that holds useful tube definitions + times (list[int]): list of index of hours to analyze + results (list): list of results for tube objects + + Returns: + panel (srlife.Panel): panel object + """ + panel = receiver.Panel(tube_dict["tube_k"]) + for _ in range(tube_dict["ass_tube_per_panel"]): + panel.add_tube(create_tube(tube_dict, times, results)) + return panel + + +def create_tube(tube_dict, times, results): + """ + Create an srlife.Tube object to add to panel + + Args: + tube_dict (dictionary): Python dictionary that holds useful tube definitions + times (list[int]): list of index of hours to analyze + results (list): list of results for tube object + + Returns: + tube (srlife.Tube): tube object + """ + tube = receiver.Tube( + 0.5 * tube_dict["od"], + tube_dict["t"], + tube_dict["h"], + tube_dict["nr"], + tube_dict["nt"], + tube_dict["nz"], + tube_dict["T0"], + ) + tube.set_times(times) + tube.multiplier_val = tube_dict["tube_mult"] + tube.T0 = tube_dict["T0"] + # initialize pressure data + pressure_data = np.ones(len(times)) * 1.0 + pressure_data[0] = 0.0 + pressure_data[-1] = 0.0 + pressure = receiver.PressureBC(times, pressure_data) + tube.set_pressure_bc(pressure) + for res in results: + tube.add_results( + res, + np.zeros((len(times), tube_dict["nr"], tube_dict["nt"], tube_dict["nz"])), + ) + return tube + + +def sample_parameters(num_threads, verbose, rtol, atol): + """ + Create a set of srlife.SolverParams to be used by srlife solvers + + Args: + num_threads (int): number of threads for analysis + verbose (bool): whether solvers should print verbose info + rtol (double): relative tol for solver convergence + atol (double): abs tol for solver convergence + Returns: None + """ + params = solverparams.ParameterSet() + + params["nthreads"] = num_threads + params["progress_bars"] = True + # If true store results on disk (slower, but less memory) + params["page_results"] = False + + params["thermal"]["miter"] = 200 + params["thermal"]["verbose"] = verbose + params["thermal"]["steady"] = True + params["thermal"]["substep"] = 2 # 10 + + params["thermal"]["solid"]["rtol"] = rtol + params["thermal"]["solid"]["atol"] = atol + params["thermal"]["solid"]["miter"] = 200 + params["thermal"]["solid"]["verbose"] = verbose + params["thermal"]["solid"]["substep"] = 2 # 10 + + params["thermal"]["fluid"]["rtol"] = rtol + params["thermal"]["fluid"]["atol"] = atol + params["thermal"]["fluid"]["miter"] = 200 + params["thermal"]["fluid"]["verbose"] = verbose + params["thermal"]["fluid"]["substep"] = 2 # 100 + + params["structural"]["rtol"] = rtol + params["structural"]["atol"] = atol + params["structural"]["miter"] = 50 + params["structural"]["verbose"] = verbose + + params["system"]["rtol"] = rtol + params["system"]["atol"] = atol + params["system"]["miter"] = 10 + params["system"]["verbose"] = verbose + + # If true store results on disk (slower, but less memory) + params["page_results"] = False + + return params diff --git a/srlife/receiver.py b/srlife/receiver.py index fc43b01..7dee2c2 100644 --- a/srlife/receiver.py +++ b/srlife/receiver.py @@ -1509,7 +1509,7 @@ def __init__(self, radius, height, nz, fluid_T, film): self.fluid_T = fluid_T self.film = film - if fluid_T.shape != (nz,) or film.shape != (nz): + if fluid_T.shape != (nz,) or film.shape != (nz,): raise ValueError( "Film coefficient and fluid temperature data must have size (nz,)" ) diff --git a/srlife/thermal.py b/srlife/thermal.py index 14f21c6..3c4243a 100644 --- a/srlife/thermal.py +++ b/srlife/thermal.py @@ -86,6 +86,7 @@ def __init__(self, pset=solverparams.ParameterSet()): self.solid_params = pset.get_default("solid", solverparams.ParameterSet()) self.thermo_params = pset.get_default("fluid", solverparams.ParameterSet()) + self.conv_fluid_nodal_T = np.array([]) def solve_receiver( self, @@ -311,7 +312,9 @@ def solve_fluid(self, i, time, dt): ) # Solve for the fluid temperatures - nodal_temps = model.solve(time) + nodal_temps = model.solve(time, self.conv_fluid_nodal_T) + # save these values to use as initial guess in next solve + self.conv_fluid_nodal_T = nodal_temps # Update the stored fluid temperature and flow velocities flow_rates, tube_temperatures = model.recover_tube_results( diff --git a/srlife/thermohydraulics/flowpath.py b/srlife/thermohydraulics/flowpath.py index 40a17b0..1b11bf0 100644 --- a/srlife/thermohydraulics/flowpath.py +++ b/srlife/thermohydraulics/flowpath.py @@ -275,6 +275,11 @@ def __init__( self.miter = 100 self.verbose = verbose + # Save previous solution for use as initial guess + # to start, this is empty + self.T_prev_sol = np.array([]) + self.has_prev_conv_T = False + def add_panel(self, weights, ri, h, metal_temp, material): """ Construct and add the standard panel -> manifold link @@ -342,7 +347,7 @@ def _setup(self): self.dof_map.append(list(range(self.nvals, self.nvals + obj.size))) self.nvals += obj.size - def solve(self, t): + def solve(self, t, T_init_guess=np.array([])): """ Solve for the current fluid temperatures @@ -353,7 +358,10 @@ def solve(self, t): self.t = t # Significant decision... - T = np.zeros((self.nvals)) + if T_init_guess.size != 0: + T = T_init_guess + else: + T = np.zeros((self.nvals)) # Initial residual R, J = self.RJ(T) diff --git a/test/Rec_Gold.hdf5 b/test/Rec_Gold.hdf5 new file mode 100644 index 0000000..23ef1c7 Binary files /dev/null and b/test/Rec_Gold.hdf5 differ diff --git a/test/data_6_d20_hr6.csv b/test/data_6_d20_hr6.csv new file mode 100644 index 0000000..314e01e --- /dev/null +++ b/test/data_6_d20_hr6.csv @@ -0,0 +1,32 @@ +Flux intensity plot (units kW/m2) +Max. flux,845.368290,kW/m2 +Min. flux,5.230420,kW/m2 +Ave. flux,317.202681,kW/m2 + +,Receiver circumferential position (0=North; +East) [deg] --> +Receiver vertical position [m],180.0,165.6,151.2,136.8,122.4,108.0,93.6,79.2,64.8,50.4,36.0,21.6,7.2,-7.2,-21.6,-36.0,-50.4,-64.8,-79.2,-93.6,-108.0,-122.4,-136.8,-151.2,-165.6, +11.5,5.865428,11.283991,13.880661,16.842514,22.815152,32.645430,43.186667,58.565392,66.347012,72.230932,85.402586,97.840197,109.429092,121.305551,141.017264,153.904936,151.934404,162.649725,138.853229,117.681587,74.587489,43.714894,27.364659,20.800167,6.491756, +11.0,17.658028,29.897828,34.924294,38.709276,45.014505,56.841670,73.787590,98.083153,110.704619,123.382102,152.472569,186.575658,219.479074,254.866113,308.945140,333.658735,333.745597,368.698913,315.307864,273.352880,179.640739,110.583814,70.516157,54.587778,27.069226, +10.5,36.674288,53.725361,61.776878,66.756774,72.478764,85.964130,111.261308,144.960179,163.938502,186.376298,235.371592,297.652387,357.103579,420.466061,513.749207,548.871754,543.427509,563.915575,478.943700,407.648136,272.732278,164.844613,104.091638,79.118222,59.918282, +10.1,58.385827,75.340911,86.230176,94.034261,99.839781,115.841957,150.200820,191.523915,218.579634,253.160567,321.114763,408.338732,488.118858,568.696977,680.669811,724.857991,707.606809,678.833907,585.772046,491.729620,344.654885,204.458647,127.979380,92.955093,89.964197, +9.6,72.547265,89.155024,104.661356,116.229239,123.349889,143.254884,184.851984,230.298822,267.376752,315.788281,398.582772,501.378057,590.364240,668.949451,769.061696,824.309470,794.997917,698.053072,611.227515,502.933148,360.085565,211.882003,130.134056,94.298304,104.189166, +9.1,73.377691,92.286521,116.005702,132.032528,142.310761,166.880555,211.843447,257.966925,306.702048,369.585016,462.757405,570.680349,657.203287,715.064570,778.664687,840.188562,803.299265,678.705748,599.267354,483.599009,346.786118,204.969999,125.356205,92.589900,98.652784, +8.6,68.126975,88.699430,118.705590,140.683250,156.752226,185.942657,231.035091,276.232795,336.840817,412.816882,511.412373,614.727913,685.250575,718.887212,751.434871,810.145711,772.694578,659.723400,585.398800,463.408221,333.813909,199.803179,124.830549,91.094809,91.101483, +8.1,63.940598,85.630442,116.109772,143.114612,166.228001,199.362980,244.063480,289.067585,359.612452,444.521903,542.487606,635.177725,682.745331,704.719091,734.388058,785.658757,755.393345,657.349880,579.319679,449.218987,329.253533,188.562395,114.133972,84.836870,89.174000, +7.7,65.938220,87.002191,114.267323,142.468584,170.958670,206.581963,252.037516,299.190928,375.814454,463.690297,555.999563,637.906126,671.859390,694.371924,738.305565,794.258378,769.387549,671.946760,588.653729,467.603572,344.652994,194.437526,116.424957,90.995471,93.698518, +7.2,72.749707,90.784625,114.938310,142.247289,172.041149,208.226388,255.404607,307.232328,385.451697,470.630232,555.793486,631.760997,669.048445,694.993634,748.389640,819.531479,788.565498,685.087510,601.891633,496.902157,360.697599,208.416296,128.624520,104.597274,98.665569, +6.7,75.284866,92.306701,116.572879,142.513004,170.562028,206.046851,255.191811,312.701706,389.458540,468.575243,547.455033,622.168154,672.309036,702.227319,750.151506,827.207277,787.296141,686.564217,600.968183,498.867845,354.642957,214.084027,132.532529,105.996993,99.333239, +6.2,70.965130,90.513041,117.968740,141.772216,167.594225,202.900224,254.045631,316.943235,391.111456,463.535258,537.986543,612.409394,673.077829,709.334630,740.888804,807.798666,764.351904,675.427188,582.730503,469.032813,332.308408,202.373653,120.507249,91.521556,94.799816, +5.7,65.369504,87.770128,118.570521,140.517851,165.446535,201.413718,254.932867,321.221018,393.579647,461.522886,533.271927,608.270106,672.165140,712.171492,736.292361,788.051437,755.380651,664.216106,562.734177,445.375915,320.775324,191.000050,112.753438,80.974581,90.969674, +5.3,65.242545,87.630401,117.939154,140.059561,165.567066,202.408342,258.663648,324.952294,397.262377,465.003605,535.895485,612.184177,673.197282,709.606967,745.951442,798.345822,773.845989,670.894856,571.483626,466.254142,338.391464,198.446811,125.640660,89.312089,95.177955, +4.8,70.196015,90.199699,116.108580,140.232651,167.285845,205.241671,262.334833,324.229833,398.012697,471.886049,544.303563,620.755932,671.838768,701.648864,754.459887,818.378749,785.713573,689.098141,592.352752,496.439935,352.445625,210.427393,141.827345,99.332781,104.781586, +4.3,73.220712,91.635071,114.157709,140.770533,169.512385,207.393869,260.669920,313.829163,389.590378,474.782080,554.314978,629.917527,667.109611,692.568176,748.590553,814.474930,771.477045,691.755474,595.052072,497.307045,346.584656,207.445635,140.437135,96.522706,107.539481, +3.8,70.409957,90.916623,114.300921,142.161480,170.575333,205.693454,249.149957,291.363464,367.493169,464.898938,558.319548,638.446104,669.411450,691.472540,736.063023,792.200397,744.409295,675.053970,581.309402,471.844117,329.698648,201.432119,121.422787,84.903556,98.904763, +3.3,66.780681,89.950648,115.486836,142.499273,167.874759,197.647559,227.030377,258.131103,331.951006,436.722451,545.794823,639.946223,682.582719,702.179375,734.972898,783.185103,743.191115,655.591600,564.161077,455.173329,322.085878,201.880120,117.572959,83.055553,93.472703, +2.9,69.395461,90.794585,116.332867,138.255296,159.014698,182.116590,197.061057,218.908629,286.921045,390.521892,509.117976,620.450519,688.392962,715.721760,754.817587,808.516245,776.752821,660.333729,571.716366,467.782506,332.096122,208.290729,130.425582,90.604223,102.996816, +2.4,76.096444,93.912987,115.017278,127.997929,143.141059,159.085041,162.804739,178.222814,237.049681,330.724159,447.850780,569.836090,660.615014,710.692050,779.728021,845.368290,807.663824,685.291615,590.610631,485.669678,335.146460,201.751498,131.093149,90.737506,114.298249, +1.9,76.145988,92.622307,105.890143,111.298611,120.551093,129.631886,127.064487,138.460137,185.790386,263.167141,367.977154,488.420249,590.661192,662.443853,765.339580,833.682320,789.255551,705.651209,606.208101,503.061993,344.327190,207.680376,129.159426,88.261077,107.457543, +1.4,63.130895,80.433153,85.132597,87.938759,92.506930,96.425126,92.207886,101.001451,136.037633,194.051011,278.519378,384.219591,485.234091,559.577868,674.698926,729.640393,693.747841,679.421832,588.045764,492.662784,339.485696,217.604536,119.520737,85.949469,77.003034, +0.9,41.958724,59.259016,56.733468,59.984219,62.544511,63.997906,60.940287,67.682370,91.338342,130.400483,190.666579,270.068952,353.606763,411.888101,509.713656,544.559336,528.545285,554.225394,484.146205,404.930489,279.459182,187.582891,94.763647,74.438300,43.539273, +0.5,20.874981,34.346603,29.849904,33.548026,36.352097,37.111611,35.968006,40.788031,55.310048,78.728883,116.325708,165.596506,218.557759,252.014828,310.135687,325.431913,323.981447,358.292285,317.543371,267.717581,185.578683,129.770336,62.883900,51.794553,19.053695, +-0.0,6.838274,13.739533,11.878137,14.792917,17.984004,18.461297,18.648988,21.709624,29.880444,42.480334,62.727140,88.124764,112.684268,125.552273,146.843451,150.444623,149.552792,156.682002,138.905323,112.543353,75.687404,50.695634,23.663315,19.962505,5.230420, diff --git a/test/data_6_d20_hr7.csv b/test/data_6_d20_hr7.csv new file mode 100644 index 0000000..1beb8a5 --- /dev/null +++ b/test/data_6_d20_hr7.csv @@ -0,0 +1,32 @@ +Flux intensity plot (units kW/m2) +Max. flux,983.427384,kW/m2 +Min. flux,8.922926,kW/m2 +Ave. flux,396.461019,kW/m2 + +,Receiver circumferential position (0=North; +East) [deg] --> +Receiver vertical position [m],180.0,165.6,151.2,136.8,122.4,108.0,93.6,79.2,64.8,50.4,36.0,21.6,7.2,-7.2,-21.6,-36.0,-50.4,-64.8,-79.2,-93.6,-108.0,-122.4,-136.8,-151.2,-165.6, +11.5,8.922926,15.235459,15.685769,21.794866,31.948658,45.217232,54.469739,77.368883,83.637944,99.849242,119.121781,128.287710,151.895427,170.283334,164.679445,207.063765,211.601128,202.007188,170.970138,148.781458,95.484908,56.229191,34.401579,15.069059,11.714641, +11.0,29.528849,41.431577,39.220303,46.098350,61.591009,82.280134,98.675091,135.926759,149.747750,178.183071,220.106817,249.931875,311.494423,367.480888,359.621292,466.268497,489.703524,466.166995,392.716792,350.138107,232.988908,147.912204,96.486663,47.187907,39.618773, +10.5,60.938384,76.785420,71.634903,77.705316,99.077692,126.279086,152.031612,204.699260,230.785744,274.426236,343.744364,404.745647,505.332961,604.958286,603.398440,732.057587,746.701333,698.254957,587.610869,510.333073,340.939173,219.483638,155.357787,97.907739,83.279784, +10.1,87.588346,109.179247,103.311393,110.488520,137.335142,168.011649,204.633129,270.291947,311.761649,372.172976,464.858497,558.507991,677.500809,796.079043,820.935327,902.933964,898.445635,835.036405,712.882125,607.175918,410.498907,262.293785,183.012993,154.015067,122.510062, +9.6,96.287036,127.714499,127.011989,139.574885,169.231515,200.762637,248.367286,321.850411,379.415697,456.763456,563.392400,679.621901,791.182740,895.914871,958.689898,956.132834,916.681507,847.301445,731.708069,603.305992,409.601966,254.966704,172.228053,175.158471,131.236359, +9.1,90.143800,129.039677,140.242220,160.077581,190.294213,224.047020,281.724333,357.231639,429.676482,521.138343,629.886225,750.983334,835.693442,904.219990,983.427384,939.968280,883.380317,812.272659,711.540924,575.027630,396.005912,242.918301,157.120404,163.854282,117.177716, +8.6,83.518119,122.080194,142.697478,167.944854,200.294942,241.236328,306.475531,379.959550,464.564908,564.062626,662.694268,771.383697,832.136784,870.824521,936.427630,903.316042,848.285491,778.114304,694.105988,563.496718,397.992773,255.916526,177.988589,149.415893,111.110618, +8.1,84.844018,118.463194,138.658459,165.599334,202.250879,253.897090,322.203681,392.994969,485.475046,586.946228,669.439170,758.922648,816.250320,851.797831,891.460919,882.148651,841.094116,768.139258,676.133189,549.975923,387.048359,255.158521,190.101168,160.930912,126.711105, +7.7,88.606675,121.368241,134.598874,160.105816,201.085459,261.182535,326.877128,397.053303,491.820843,593.060860,664.836097,741.967334,808.589537,860.190711,899.164891,896.083232,868.887095,802.126067,693.041056,562.966183,398.573746,251.822365,183.396783,176.402568,137.802650, +7.2,88.710764,125.769425,133.170083,157.922124,200.522821,261.621018,321.133330,393.126654,485.161978,587.240577,661.752582,737.851675,811.417253,876.622836,935.762797,923.708858,892.948306,836.065051,723.803863,591.259520,419.344867,242.410757,166.677255,168.130270,128.310612, +6.7,86.065741,124.753428,134.301412,160.309443,200.867271,255.742796,310.372835,384.517840,471.875633,574.340288,662.382007,745.739866,820.111667,881.109089,947.579867,929.113098,888.759499,832.776927,739.211230,587.313809,433.472993,244.231068,177.563916,154.157871,116.069420, +6.2,86.851176,119.416419,137.085627,164.416711,200.750613,247.559469,301.901607,376.635029,460.551694,559.939780,662.615454,756.177377,828.487583,870.598640,922.294510,905.099617,854.867439,796.788237,718.561180,558.587341,412.199248,263.844891,199.314350,161.916323,122.734418, +5.7,89.944198,117.985859,139.770124,166.781516,199.735027,242.690870,301.664359,374.567335,457.908739,550.657663,660.314397,762.249796,832.090315,862.507878,898.352879,878.265216,838.176828,783.117760,690.109151,526.153235,387.453526,260.598841,194.389202,163.884132,130.553481, +5.3,89.997151,122.085845,140.100728,166.342130,198.567216,245.294549,311.857081,381.799322,466.912137,551.666036,656.961141,760.431916,830.252376,873.655748,915.045446,887.889543,865.355925,808.681142,694.228783,549.774008,393.866204,256.784323,169.279408,147.322551,121.243026, +4.8,89.016938,125.803005,137.797867,163.431607,198.236588,253.544994,328.770256,396.809319,484.309183,563.091723,654.056414,749.780361,821.193348,888.449898,944.178123,920.186307,893.232929,835.925933,720.447840,586.505016,411.281675,240.192916,164.269971,136.802472,106.931973, +4.3,89.585605,125.020490,135.508295,161.239002,199.411261,261.879500,343.723023,410.684704,501.928364,580.298203,655.398990,738.531490,807.999827,885.633338,941.081684,926.871260,891.109574,832.541985,723.872008,593.237626,413.601192,238.403508,190.497411,159.829692,113.365312, +3.8,89.184180,119.888233,136.709522,163.201599,201.690289,264.297254,345.662312,411.045793,507.533261,594.457950,664.975344,741.203071,801.621025,870.108223,910.592801,904.258621,858.150510,807.080057,699.575834,579.324819,400.026988,257.210962,196.531796,172.245066,126.789928, +3.3,88.066721,115.749899,141.537761,168.288468,202.118211,257.171147,328.525994,390.148364,490.999016,592.495575,676.400264,761.210249,810.063869,859.240873,897.831413,875.389555,835.998638,792.386018,673.780323,561.338860,393.633684,259.267434,179.446006,153.165240,123.084747, +2.9,89.637651,118.346971,146.327505,170.610268,196.608815,239.559905,294.101691,350.298785,449.723285,562.537575,669.778693,779.674951,830.921057,874.765262,930.612874,880.033816,849.270790,812.020584,679.540443,570.250074,411.013310,258.931136,157.142770,131.240173,109.023860, +2.4,95.625073,126.155695,143.126580,162.824216,181.968365,212.401973,249.688600,300.479925,389.101990,501.243356,625.314710,762.514919,840.506703,908.864416,970.330468,914.828105,882.660814,839.002695,703.170130,575.029215,411.689479,245.738549,174.539984,150.499486,109.352259, +1.9,98.987652,125.427274,123.163849,139.814752,157.063236,177.196743,201.199731,246.809825,316.670459,414.383581,537.723478,685.399883,796.632680,900.007702,949.631834,937.812431,917.069825,860.379968,731.356639,597.231557,420.878011,260.531246,197.470836,170.412804,119.088469, +1.4,88.497300,105.834032,88.435603,104.431550,122.911758,135.594408,151.351234,190.013396,238.658210,313.131211,418.000464,550.699542,676.507731,793.646209,826.210678,894.664366,902.993490,832.705332,720.638009,599.666548,414.077457,276.754368,179.212555,153.595975,114.637919, +0.9,63.003641,73.883265,51.880555,66.051805,84.542109,92.352880,103.037852,131.823296,162.457595,211.949464,287.643763,384.811562,496.902131,595.715128,619.813791,727.392819,755.892001,682.046833,597.090584,501.748379,337.279137,237.624803,118.780378,91.264270,84.510060, +0.5,32.955955,40.059096,24.578039,34.593556,49.616879,54.264599,61.490290,79.176869,97.472717,126.205968,172.268150,228.552909,303.243562,360.032313,377.496057,466.057780,497.897431,449.949907,397.557133,339.144048,226.251746,160.985653,61.623834,35.757786,44.002549, +-0.0,11.499806,14.719394,9.367245,14.635576,24.200211,26.711428,31.423467,39.928058,50.755003,65.652015,89.589039,115.675556,150.847196,170.674352,179.282804,209.431456,220.527056,196.411236,171.605219,144.060590,90.502142,60.620801,20.447250,9.534841,13.907637, diff --git a/test/test_interface.py b/test/test_interface.py new file mode 100644 index 0000000..4b186a6 --- /dev/null +++ b/test/test_interface.py @@ -0,0 +1,555 @@ +import unittest +import tempfile + +import numpy as np +from scipy.interpolate import RegularGridInterpolator +import multiprocess +import subprocess +import h5py + +from srlife import ( + receiver, + thermal, + structural, + system, + damage, + library, + managers, + interface, +) + +header = b"\n\n\n\n\n\n" + +fake_data = b"""-,180,135,90,45,0,-45,-90,-135, +0,1,1,1,1,1,1,1,1 +1,1,1,1,1,1,1,1,1 +2,1,1,1,1,1,1,1,1 +3,1,1,1,1,1,1,1,1 +""" + +nr_unit = 2 +nt_unit = 4 +nz_unit = 4 +tube_dict_unit = { + "od": 21.3, + "t": 1.25, + "h": 3, + "nr": nr_unit, + "nt": nt_unit, + "nz": nz_unit, + "spacing": 1, + "T0": 300, + "tube_k": "rigid", + "eta": 0.1, + "tube_mult": 58.5, + "ass_tube_per_panel": 2, +} + + +def calc_fd(Re, tube_eta, tube_Dh): + """ + Calculate Darcy friction factor from: + Zigrange and Sylvester 1985, A review of explicit friction factor equations. + J of Energy Resources Technology) + Equation 13 + + Args: + Re (double): Reynolds number + tube_eta (double): tube roughness + tube_Dh (double): tube hydraulic diameter + Returns: + fd (double): darcy friciton factor + """ + if Re < 4000: + # laminar flow + fd = 64 / Re + else: + fd = ( + 1 + / ( + -2 + * np.log10( + (tube_eta / tube_Dh) / 3.7 + - 5.02 / Re * np.log10((tube_eta / tube_Dh) / 3.7 + 13 / Re) + ) + ) + ** 2 + ) + return fd + + +class InterfaceUnitTests(unittest.TestCase): + """ + Run unit tests on some funcitons of interface.py + """ + + def test_file_not_found(self): + flux_filename = "no_file_existing.csv" + r = 10 + c = 17 + data_shape = [r, c] + z_offset_in = 0.0 + z_data, theta_data, flux_data = interface.read_month_day_hour_flux_file( + flux_filename, data_shape, z_offset_in + ) + print(theta_data) + self.assertTrue(np.array_equal(z_data, np.zeros([r, 1]))) + self.assertTrue(np.array_equal(theta_data, np.zeros([1, c]))) + self.assertTrue(np.array_equal(flux_data, np.zeros([r, c]))) + + def test_fake_file(self): + # Create a mock CSV file content to test the parsing + with tempfile.NamedTemporaryFile(delete=True) as f: + f.write(header) + f.write(fake_data) + f.seek(0) + z_data, theta_data, flux_data = interface.read_month_day_hour_flux_file( + f.name, [4, 8], 0 + ) + # spot check data + print(z_data) + self.assertEqual(z_data[2], 2000) + self.assertEqual(theta_data[-1], -180) + self.assertEqual(flux_data[2, 5], 1) + + def test_apply_tube_flux_bcs(self): + # test applying BCs via interpolator + times = np.array([0, 1]) + tube = interface.create_tube(tube_dict_unit, times, []) + tube_absorbance = 1.0 + tube_rec_theta = 12 + z_pts = np.linspace(10, 20, nz_unit) + theta_pts = np.linspace(15, -15, nt_unit) + data = np.array([[(2 * ix + 4 * iy) for ix in z_pts] for iy in theta_pts]) + interpFun = RegularGridInterpolator( + (z_pts, theta_pts), + data, + method="linear", + bounds_error=False, + fill_value=0.0, + ) + tube_rec_theta = 4 + interface.apply_tube_flux_bcs( + tube, tube_absorbance, tube_rec_theta, [interpFun] + ) + # spot check tube BCs + self.assertAlmostEqual(tube.outer_bc.data[0, 0, 2], interpFun([0.0, 3.0])) + self.assertEqual(tube.outer_bc.data[0, 3, 3], 0.0) + + def test_calc_fluid_velocity(self): + mfr = 10 + rho = 5 + Dh = 3 + vel = mfr / rho * (4.0 / (np.pi * Dh**2)) + self.assertAlmostEqual(interface.calc_fluid_velocity(mfr, rho, Dh), vel) + + def test_calc_reynolds_number(self): + mfr = 28 + rho = 5 + mu = 5 + Dh = 3 + vel = mfr / rho * (4.0 / (np.pi * Dh**2)) + Re = rho * vel * Dh / mu + self.assertAlmostEqual(interface.calc_reynolds_number(rho, vel, mu, Dh), Re) + + def test_calc_fd(self): + Re = [2000, 10000] + eta = 0.2 + Dh = 0.012 + for iRe in Re: + fd = interface.calc_fd(iRe, eta, Dh) + fd_test = calc_fd(iRe, eta, Dh) + self.assertAlmostEqual(fd, fd_test, places=12) + + def test_calc_friction_p_loss(self): + mfr = 0.001 + rho = 50 + mu = 0.004 + delta_l = 1 + p_loss = interface.calc_friction_p_loss(rho, mu, mfr, tube_dict_unit, delta_l) + test_val = 26.09267992590423 + self.assertAlmostEqual(p_loss, test_val) + + def test_calc_manifold_bend_loss(self): + mfr = 0.001 + rho = 1550 + mu = 0.004 + manifold_tube_dict = { + "od": 21.3, + "t": 1.25, + "bend_radius": 18, + "bends_per_panel": 4, + "eta": 0.01, + } + num_bends = 4 + + p_loss = interface.calc_manifold_bend_loss( + mfr, manifold_tube_dict, rho, mu, num_bends + ) + test_val = 0.12351409567895541 + self.assertAlmostEqual(p_loss, test_val) + + def test_create_receiver(self): + num_days = 1 + times = np.array([0, 1, 2]) + period = 2 + panel_k = "disconnect" + num_panels = 12 + results = [] + rec = interface.create_receiver( + tube_dict_unit, num_days, times, period, panel_k, num_panels, results + ) + self.assertEqual(num_panels, rec.npanels) + for panel in rec.panels.values(): + self.assertEqual(tube_dict_unit["ass_tube_per_panel"], panel.ntubes) + + def test_cycle_tube_pressure_bcs(self): + times = np.array([0, 1]) + tube_pressure = 12 + pressure = np.ones_like(times) * tube_pressure + pressure[0] = 0.0 + tube = interface.create_tube(tube_dict_unit, times, []) + tube_pressure_bc = receiver.PressureBC(times, pressure) + tube.set_pressure_bc(tube_pressure_bc) + + tubes_dict = {"0": tube} + num_cycles = 3 + cyclic_times = np.array([0, 1, 2, 3]) + interface.cycle_tube_pressure_bcs(tubes_dict, num_cycles, cyclic_times) + gold = np.array([0, 12, 12, 12]) + self.assertTrue(np.array_equal(gold, tube.pressure_bc.data)) + + def test_set_and_downsample_tube_temp_bcs(self): + times = np.array([0, 1, 2]) + cyclic_times = np.array([0, 1, 2]) + num_cycles = 1 + tube = interface.create_tube(tube_dict_unit, times, []) + tubes_dict = {"0": tube} + data = 245 * np.ones([tube.ntime, tube.nr, tube.nt, tube.nz]) + tube.add_results("temperature", data) + # Test 1: inlet T + interface.set_and_downsample_tube_temp_bcs( + tubes_dict, + True, + 245, + cyclic_times, + num_cycles, + "3d", + "", + ) + self.assertEqual(245, tube.T0) + analysis_type = "2d" + locs = ["max_T", "max_avg_T"] + for i_loc, loc in enumerate(locs): + tube = interface.create_tube(tube_dict_unit, times, []) + tubes_dict = {"0": tube} + data = (i_loc + 245) * np.ones([tube.ntime, tube.nr, tube.nt, tube.nz]) + tube.add_results("temperature", data) + interface.set_and_downsample_tube_temp_bcs( + tubes_dict, + False, + 245, + cyclic_times, + num_cycles, + analysis_type, + loc, + ) + self.assertTrue( + np.array_equal(tube.results["temperature"], data[:, :, :, 0]) + ) + + +class InterfaceRegressionTest(unittest.TestCase): + def test_simple_model(self): + multiprocess.set_start_method("spawn", force=True) + # Receiver filename to write to + rec_filename = "./test/Rec_Test" + # file extension used by srlife + hdf5_ext = ".hdf5" + + # GEOMETRY + # rec geom + num_panels = 2 # number of panels in receiver + rec_diam = 10000 # mm + rec_height = 12000 # mm + + # rec tube geom + tube_od = 21.3 # mm (Custom, Schedule 5) + tube_t = 1.24 # mm (Custom, Schedule 5) + tube_h = 12000 # mm (also rec height) + tube_spacing = 1 # mm + tube_eta = 0.01e-3 # tube roughness + tube_absorbance = 0.98 # absorbance of tube + z_offset = 0.25 # z offset for flux data + + # manifold tube geom + manifold_tube_od = 502.15 # mm + manifold_tube_t = 45.24 # mm + manifold_bend_radius = 2 * manifold_tube_od # used in p loss calc + manifold_tube_eta = tube_eta # roughness + manifold_bends_per_panel = 4 + # END GEOMETRY + + # ASSUMPTIONS + ass_tube_per_panel = 2 # number of tubes analyzed at each panel + outlet_p = 0.5 # MPa + + # tube assumption work + act_tube_per_panel = int( + np.pi * rec_diam / (num_panels * (tube_od + tube_spacing)) + ) + tube_multiplier = act_tube_per_panel / ass_tube_per_panel + print(f"Actual tubes per panel = {act_tube_per_panel}") + print(f"Tube Multiplier = {tube_multiplier}") + + # Structure assumptions + panel_k = "disconnect" # spring stiffness panel connect + tube_k = "rigid" # stiffness of tube connect + set_init_T_to_inlet_T = True + loc = "max_avg_T" # "max_T" + analysis_type = "2d" # "3d" + is_single_panel_analysis = True + single_panel_analysis_id = "0" + # END ASSUMPTIONS + + # BC/IC INFO + # flux data input and shape + flux_data_dir = "./test" + flux_data_shape = [25, 25] + tube_initial_temp = 300.0 # K + # END BC/IC INFO + + # ANALYSIS TIMES and CYCLE INFO + num_days = 1 + start_time = 6 + end_time = 7 + time_step = 1 # hour, can do nonInt + # stimes is index of hours to analyze, starting with first hour of flux BCs + stimes = np.arange(1, (end_time + 1) - start_time + time_step, time_step) + period = len(stimes) + # array of sample times of across number of days + times = np.tile(stimes, num_days) + # add time zero for initial conditions + times = np.append(np.array([0]), times) + # month and date from flux data + month = 6 + day = 20 + print("Times analyzed: \n", times) + num_cycles = 2 + # END ANALYSIS TIMES AND CYCLE INFO + + # TUBE DISC INFO + nz = int(interface.convert_mm_to_m(rec_height) * 0.5 + 1) + nt = 8 + nr = 2 + # END TUBE DISC INFO + + # SOLVERS # + # setup solver parameters + num_threads = 1 + # NOTE: using large tolerances to assure quick convergence + rtol = 1.0e-2 + atol = 1.0e-3 + params = interface.sample_parameters(num_threads, False, rtol, atol) + # thermal + thermal_solver = thermal.ThermohydraulicsThermalSolver(params["thermal"]) + # Structural solver + structural_solver = structural.PythonTubeSolver(params["structural"]) + # receiver system solver + system_solver = system.SpringSystemSolver(params["system"]) + # material damage model + damage_model = damage.TimeFractionInteractionDamage(params["damage"]) + # END SOLVERS # + + # MATERIALS + # fluid model + mat_fluid = library.load_thermal_fluid("32MgCl2-68KCl", "base") + # material models + st_mat_model = "elastic_creep" # "elastic_model" "base" + mat_thermal, mat_deformation, mat_damage = library.load_material( + "740H", "base", st_mat_model, "base" + ) + # END MATERIALS + + # structure output + save_struct_files_to_vtu = False + st_filename = "structure_" + st_mat_model + tube_filename = "tube_" + st_mat_model + + # FLOW SYSTEM SETUP + T_out_target = 720 # target outlet temperature, unit: C + # tolerance for mass flow optimization + # NOTE: making very large for test for single iter convergence + pct_err_outlet_temp = 100.25 + panel_flow_path = [["1"], ["0"]] + mass_flow_per_path = np.ones((len(panel_flow_path), len(times))) * 650 # kg/s + T_in_per_path = np.array([500, 500]) # Celcius + use_cycle_reset_heuristic = False + save_heat_to_vtu = False + # END FLOW SYSTEM SETUP + + # START ANALYSIS WORK + # rec tube data structure + tube_dict = { + "od": tube_od, + "t": tube_t, + "h": tube_h, + "nr": nr, + "nt": nt, + "nz": nz, + "spacing": tube_spacing, + "T0": tube_initial_temp, + "tube_k": tube_k, + "eta": tube_eta, + "tube_mult": tube_multiplier, + "ass_tube_per_panel": ass_tube_per_panel, + } + # manifold tube data structure + manifold_tube_dict = { + "od": manifold_tube_od, + "t": manifold_tube_t, + "bend_radius": manifold_bend_radius, + "bends_per_panel": manifold_bends_per_panel, + "eta": manifold_tube_eta, + } + + # STEP 1: create a receiver + rec = interface.create_receiver( + tube_dict, num_days, times, period, panel_k, num_panels, results=[] + ) + + # STEP 2: read and assign flux BCs to receiver tubes + # read flux files, make interpolator functions by hour and save + # NOTE: not doing multiple days right now + # NOTE: can we assume zData and thetaData are constant in time??? + flux_interpolators_by_hour = interface.get_flux_interpolators_from_data_files( + times, + start_time, + end_time, + month, + day, + flux_data_dir, + flux_data_shape, + z_offset, + ) + # Take flux interpolators and write flux BCs for each tube + # flux is sampled at tube centerline and applied smoothly + # across sunSideSurface via cosine + interface.calc_and_write_tube_flux_bcs( + rec, + ass_tube_per_panel, + num_panels, + flux_interpolators_by_hour, + tube_absorbance, + ) + + rec.save(rec_filename + hdf5_ext) + + # STEP 3: add flow paths to receiver + # add flow paths to receiver + interface.set_rec_flow_paths( + rec, panel_flow_path, mass_flow_per_path, T_in_per_path + ) + rec.save(rec_filename + hdf5_ext) + + # make solver + solver = managers.SolutionManager( + rec, + thermal_solver, + mat_thermal, + mat_fluid, + structural_solver, + mat_deformation, + mat_damage, + system_solver, + damage_model, + pset=params, + ) + + # STEP 4: optimize mass flow rate to match target outlet temp + # now optimize mass flow rate for each path + # this takes a long time, can skip if it has already been done + run_mass_flow_opt = True + # Heuristics + if use_cycle_reset_heuristic: + solver.add_heuristic(managers.CycleResetHeuristic()) + if run_mass_flow_opt: + interface.optimize_mass_flow_rate_per_path( + rec, + rec_filename, + T_out_target, + pct_err_outlet_temp, + solver, + save_heat_to_vtu, + ) + # Saving optimized results to base filename + rec.save(rec_filename + hdf5_ext) + else: + rec_filename = "./Rec_Test_ht_iter_2" + rec = receiver.Receiver.load(rec_filename + hdf5_ext) + + mass_flow = np.array([path["mass_flow"] for path in rec.flowpaths.values()]) + print(f"Mass flow: {mass_flow/3600}") + + # STEP 5: calc pressure loss from flow + # calculate pressure loss from + # pipe friction + # head loss + # manifold loss + flow_path_p_loss = interface.calc_p_loss_from_flows_temps( + rec, tube_dict, manifold_tube_dict, mat_fluid, outlet_p + ) + inlet_p = outlet_p + interface.convert_Pa_to_MPa(flow_path_p_loss) + interface.update_tube_pressure_bcs(rec, inlet_p, outlet_p) + rec.save(rec_filename + hdf5_ext) + + # STEP 6: solve structure and life receiver + struct_output_dict = { + "save_to_vtu": save_struct_files_to_vtu, + "st_filename": st_filename, + "tube_filename": tube_filename, + } + interface.run_struct_analysis( + rec_filename, + set_init_T_to_inlet_T, + loc, + analysis_type, + is_single_panel_analysis, + single_panel_analysis_id, + num_cycles, + solver, + struct_output_dict, + ) + # END ANALYSIS WORK + # do some checks on hdf5 files + file_gold = h5py.File("./test/Rec_Gold.hdf5", "r+") + file_test = h5py.File(rec_filename + hdf5_ext, "r+") + # flux bcs for a panel/tube + data_name = "/panels/1/tubes/0/outer_bc/data" + print("Flux BC") + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # fluid temps for two tubes + print("Fluid Temp") + data_name = "/panels/1/tubes/1/axial_results/fluid_temperature" + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # temps for a tube + print("Temp") + data_name = "/panels/1/tubes/0/results/temperature" + self.assertTrue( + np.allclose( + file_test[data_name], file_gold[data_name], rtol=rtol, atol=atol + ) + ) + # NOTE: in interface.run_struct_analysis we make this a + # single panel model, and don;t save it in this run + # so structure results are not available + # all else matches and run_struct_analysis finished, + # so call it cleared