diff --git a/LICENSE.txt b/LICENSE.txt
index ab35250..e0baf6a 100644
--- a/LICENSE.txt
+++ b/LICENSE.txt
@@ -1,6 +1,6 @@
The MIT License (MIT)
-Copyright (c) 2021-2022 Leif-Thore Deck, David Ochsenbein
+Copyright (c) 2021-2024 Leif-Thore Deck, Andraž Košir, David Ochsenbein
Permission is hereby granted, free of charge, to any person obtaining a copy
of this software and associated documentation files (the "Software"), to deal
diff --git a/README.md b/README.md
index cb7c776..ec737cd 100644
--- a/README.md
+++ b/README.md
@@ -3,7 +3,7 @@
  
      
-SNOW (**S**tochastic **N**ucleation **O**f **W**ater) is an open source modeling framework to simulate the freezing process in a large number of vials, tailored towards pharmaceutical freeze-drying processes. SNOW was developed as part of a research collaboration between ETH Zurich's Separation Processes Laboratory and Janssen (pharmaceutical companies of J&J). It is brought to you by Leif-Thore Deck and Dr. David Ochsenbein.
+SNOW (**S**tochastic **N**ucleation **O**f **W**ater) is an open source modeling framework to simulate the freezing process in a large number of vials, tailored towards pharmaceutical freeze-drying processes. SNOW was developed as part of a research collaboration between ETH Zurich's Separation Processes Laboratory and Janssen (pharmaceutical companies of J&J). It is brought to you by Leif-Thore Deck, Andraž Košir and Dr. David Ochsenbein.
## Description
@@ -11,9 +11,15 @@ SNOW is a model capable of simulating the entire freezing process for a batch of
In addition to the simulation of the thermal evolution of vials during freezing, SNOW tracks a number of characteristic quantities of the freezing process for each vial. These quantities are the nucleation time, the nucleation temperature and the solidification time. Both nucleation and solidification are correlated with attributes of the frozen product such as the mean crystal size and the product activity of an API. Variability in these quantities among vials may thus be used as surrogate for batch heterogeneity. SNOW consequently may be used in process design and optimization to identify suitable cooling protocols with sufficiently small batch heterogeneity.
-A more detailed description of the modeling framework is presented in a recent publication, which is publicly accessible under https://doi.org/10.1016/j.ijpharm.2021.121276 While version 1.0 is tailored towards simulations of the freezing stage in freeze-drying, version 1.1. was developed for pallet freezing; pallet freezing is for example used in the manufacturing of the Janssen COVID-19 vaccine, which served as case study for the model. Version 1.1 and the case study on the COVID-vaccine are discussed in detail in a recent pre-print at ChemRxiv: https://doi.org/10.26434/chemrxiv-2022-gnwhf
+A more detailed description of the modeling framework is presented in a recent publication, which is publicly accessible under https://doi.org/10.1016/j.ijpharm.2021.121276 While version 1.0 is tailored towards simulations of the freezing stage in freeze-drying, version 1.1. was developed for pallet freezing; pallet freezing is for example used in the manufacturing of the Janssen COVID-19 vaccine, which served as case study for the model. Version 1.1 and the case study on the COVID-vaccine are discussed in detail in a recent scientific publication: https://doi.org/10.1016/j.ijpharm.2022.122051
-### Additional features currently supported (as of version 1.1)
+An extension of the modeling framework for the simulation of the freezing process with spatial resolution within individual containers has been implemented in version 2.0 and is discussed in detail within the scope of a recent scientific publication: https://doi.org/10.1016/j.cej.2024.148660.
+
+### Additional features supported by version 2.0
+- Spatial simulation of freezing within individual containers with different dimensionalities (0D, 1D, 2D)
+- Simulation of three process configurations (shelf-ramped freezing, vacuum-induced surface freezing, jacket-ramped freezing)
+
+### Additional features supported by version 1.1
- Simulation of systems comprising vials arranged in three spatial dimensions
- Arbitrary initial temperature of the vials
- Improvements in the numerical implementation (Second method to compare the initial amount of ice formed, faster data saving)
@@ -35,6 +41,6 @@ A more detailed description of the modeling framework is presented in a recent p
Please report bugs as [Github issues](https://github.com/SPL-ethz/snow/issues/new?assignees=ltdeck&labels=bug) or via [Email](mailto:deckl@ethz.ch), preferably
including a screenshot that illustrates the problem.
-Copyright (c) 2021-2022 Leif-Thore Deck, David Ochsenbein
+Copyright (c) 2021-2024 Leif-Thore Deck, Andraž Košir, David Ochsenbein
The snow package logo has been designed using resources from Flaticon.com.
diff --git a/docs/.buildinfo b/docs/.buildinfo
index c984184..9356573 100644
--- a/docs/.buildinfo
+++ b/docs/.buildinfo
@@ -1,4 +1,4 @@
# Sphinx build info version 1
# This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done.
-config: 2ccd3234bc1616e1d4b456b969de6cf8
+config: c1e8aa707df1a9e0485c5a57df96b0be
tags: 645f666f9bcd5a90fca523b33c5a78b7
diff --git a/docs/.doctrees/api/ethz_snow.doctree b/docs/.doctrees/api/ethz_snow.doctree
index 7fcf564..4cbfdec 100644
Binary files a/docs/.doctrees/api/ethz_snow.doctree and b/docs/.doctrees/api/ethz_snow.doctree differ
diff --git a/docs/.doctrees/api/modules.doctree b/docs/.doctrees/api/modules.doctree
index 9651d44..77c74cf 100644
Binary files a/docs/.doctrees/api/modules.doctree and b/docs/.doctrees/api/modules.doctree differ
diff --git a/docs/.doctrees/authors.doctree b/docs/.doctrees/authors.doctree
index ae3a8e7..46246cb 100644
Binary files a/docs/.doctrees/authors.doctree and b/docs/.doctrees/authors.doctree differ
diff --git a/docs/.doctrees/changelog.doctree b/docs/.doctrees/changelog.doctree
index 8bb639f..66f7011 100644
Binary files a/docs/.doctrees/changelog.doctree and b/docs/.doctrees/changelog.doctree differ
diff --git a/docs/.doctrees/development.doctree b/docs/.doctrees/development.doctree
index 2f30c8f..b58a3fc 100644
Binary files a/docs/.doctrees/development.doctree and b/docs/.doctrees/development.doctree differ
diff --git a/docs/.doctrees/index.doctree b/docs/.doctrees/index.doctree
index a3c8e2b..a28c9d9 100644
Binary files a/docs/.doctrees/index.doctree and b/docs/.doctrees/index.doctree differ
diff --git a/docs/.doctrees/license.doctree b/docs/.doctrees/license.doctree
index 6bca634..6a9803d 100644
Binary files a/docs/.doctrees/license.doctree and b/docs/.doctrees/license.doctree differ
diff --git a/docs/.doctrees/readme.doctree b/docs/.doctrees/readme.doctree
index 728e9e3..f0802d9 100644
Binary files a/docs/.doctrees/readme.doctree and b/docs/.doctrees/readme.doctree differ
diff --git a/docs/.doctrees/tutorial.doctree b/docs/.doctrees/tutorial.doctree
index 87e5459..93edc58 100644
Binary files a/docs/.doctrees/tutorial.doctree and b/docs/.doctrees/tutorial.doctree differ
diff --git a/docs/_modules/ethz_snow/constants.html b/docs/_modules/ethz_snow/constants.html
index 0012e6f..481a87f 100644
--- a/docs/_modules/ethz_snow/constants.html
+++ b/docs/_modules/ethz_snow/constants.html
@@ -8,7 +8,7 @@
- ethz_snow.constants — ethz_snow 1.1.0 documentation
+ ethz_snow.constants — ethz_snow 1.1.0.post1.dev34+ga1dec07.d20230707 documentation
@@ -37,7 +37,7 @@
-
+
@@ -67,7 +67,7 @@
T_eq_l=T_eq-depression
- # bundle things into a dict now
- # don't do it earlier for readability
+ # latent heat of fusion
+ Dh=float(config["water"]["Dh"])
+
+ # needed for solving nucleation
+ k_f=float(config["solution"]["k_f"])
+ M_s=float(config["solution"]["M_s"])
+
+ # general constants
+ sigma_B=float(config["general"]["sigma_B"])
+ k_B=float(config["general"]["k_B"])
+
+ # mass of solute and water in the solution
+ mass_solute=mass*solid_fraction
+ mass_water=mass*(1-solid_fraction)
+
+ # freezing configuration
+ configuration=str(config["configuration"])
+
constVars=[
+ "dimensionality","T_eq","T_eq_l","kb",
@@ -351,6 +427,7 @@
[docs]classSnowing:
+ """A class to handle a single SNOWing (Stochastic Nucleation of Water
+ with INternal Gradients) simulation.
+
+ More information regarding the equations and their derivation can be found in
+ XXX, Deck et al. (2023).
+
+ Parameters:
+ k (dict): Heat transfer coefficients.
+ opcond (OperatingConditions): Operating conditions of the run.
+ temperature (str): Model dimensionality.
+ configuration (str): Freezing configuration.
+ plotting (bool): Whether evolutions are plotted or not.
+ Nrep (int): Number of repetitions.
+ configPath (Optional[str]): The path of the (optional) custom config yaml.
+ """
+
+ def__init__(
+ self,
+ k:dict={"int":0,"ext":0,"s0":50,"s_sigma_rel":0},
+ opcond:OperatingConditions=OperatingConditions(),
+ Nrep:int=1,
+ configPath:Optional[str]=None,
+ ):
+ """Construct a Snowing object.
+
+ Args:
+ k (dict, optional): Heat transfer coefficients.
+ Defaults to {"int": 0, "ext": 0, "s0": 50, "s_sigma_rel": 0}.
+ opcond (OperatingConditions, optional):
+ Operating conditions of the run. Defaults to OperatingConditions().
+ Nrep (int, optional): Number of repetitions. Defaults to 1.
+ configPath (Optional[str], optional): The path of the (optional)
+ custom config yaml. Defaults to None.
+
+ Raises:
+ TypeError: If k is not a dict.
+ ValueError: If a aspecific key is missing in dict k.
+ TypeError: If opcond is not of type operatingConditions.
+
+ Examples:
+ The following command will create and instance of the Snowing class with
+ default operating conditions and heat transfer parameters:
+
+ >>> S = Snowing()
+
+ To customize the operating conditions, first define the operating conditions
+ and the heat transfer parameters:
+
+ >>> d = {"int": 0, "ext": 0, "s0": 50, "s_sigma_rel": 0}
+ >>> c = {"rate": 0.5 / 60, "start": 20, "end": -50}
+ >>> h = [dict(duration=60*60, temp=-10)]
+ >>> op = OperatingConditions(t_tot=5*3600, cooling=c, holding=h)
+
+ Then an instance of the Snowing class can be created:
+
+ >>> S = Snowing(k=d, opcond=op)
+
+ See tutorial for more information on how to use the spatial model functionalities.
+ """
+
+ ifnotisinstance(k,dict):
+ raiseTypeError(f"Input k must be of type dict. Was given {type(k)}.")
+ elifnotall([keyink.keys()forkeyinHEATFLOW_REQUIREDKEYS]):
+ raiseValueError(
+ (
+ f"A required key was missing from dictionary k, specifically "
+ +f"{set(HEATFLOW_REQUIREDKEYS)-set(k.keys())}."
+ )
+ )
+ else:
+ self.k=k
+
+ # for saving the results
+ self._domain=None
+ self._prop=None
+ self._time=None
+ self._shelfTemp=None
+ self._temp=None
+ self._iceMassFraction=None
+
+ # dictionary of results
+ self._stats={}
+
+ # number of repetitions
+ self.Nrep=Nrep
+
+ # dictionary of results for multiple simulations
+ ifself.Nrep>1:
+ self._statsMultiple={}
+
+ # path to the parameters file
+ self.configPath=configPath
+
+ # simulation progress
+ self._simulationStatus=0
+
+ # operating conditions
+ ifnotisinstance(opcond,OperatingConditions):
+ raiseTypeError(
+ "Input opcond must be an instance of class OperatingConditions."
+ )
+ else:
+ self.opcond=opcond
+
+ @property
+ defconfigPath(self)->Optional[str]:
+ """Get configPath property."""
+ returnself._configPath
+
+ @configPath.setter
+ defconfigPath(self,value):
+ """Set configPath property."""
+ self.const=calculateDerived(value)
+ self._configPath=value
+
+ @property
+ defsimulationStatus(self)->int:
+ """Return simulation status of instance.
+
+ Returns:
+ int: Simulation status. 0 = not run, 1 = run.
+ """
+ if(self._simulationStatus==1)or(self._tempisnotNone):
+ self._simulationStatus=1
+
+ returnself._simulationStatus
+
+ @property
+ deftime(self)->np.ndarray:
+ """Return time array in the simulation of instance.
+
+ Returns:
+ np.ndarray: Time array in hours.
+ """
+ assert(self._simulationStatus==1)or(
+ self._timeisnotNone
+ ),"Simulation not yet carried out or completed."
+ returnself._time
+
+ @property
+ deftemp(self)->np.ndarray:
+ """Return temperature evolution in the simulation of instance.
+
+ Returns:
+ np.ndarray: Temperature evolution in °C.
+ """
+ assert(self._simulationStatus==1)or(
+ self._tempisnotNone
+ ),"Simulation not yet carried out or completed."
+ returnself._temp
+
+ @property
+ defshelfTemp(self)->np.ndarray:
+ """Return shelf temperature evolution in the simulation of instance.
+
+ Returns:
+ np.ndarray: Shelf temperature evolution in °C.
+ """
+ assert(self._simulationStatus==1)or(
+ self._shelfTempisnotNone
+ ),"Simulation not yet carried out or completed."
+ returnself._shelfTemp
+
+ @property
+ deficeMassFraction(self)->np.ndarray:
+ """Return ice mass fraction evolution in simulation of instance.
+
+ Returns:
+ np.ndarray: Ice mass fraction evolution (no units).
+ """
+ assert(self._simulationStatus==1)or(
+ self._iceMassFractionisnotNone
+ ),"Simulation not yet carried out or completed."
+ returnself._iceMassFraction
+
+ @property
+ defresults(self)->dict:
+ """Return minimum, mean, kinetic mean and maximum nucleation temperature
+ in simulation of instance.
+
+ Returns:
+ dict: 4 values of nucleation temperature (minimum, kinetic mean,
+ mean and maximum). Irrelevant for homogeneous model simulation.
+ """
+ assert(
+ self._simulationStatus==1
+ ),"Simulation not yet carried out or completed."
+ ifself.Nrep==1:
+ self._stats_df=pd.DataFrame.from_dict(self._stats,orient="index").T
+ returnself._stats_df
+ elifself.Nrep>1:
+ self._statsMultiple_df=pd.DataFrame.from_dict(
+ self._statsMultiple,orient="index"
+ )
+ self._statsMultiple_df.columns=list(self._keys)
+ returnself._statsMultiple_df
+
+ # simulate the entire process
+
+
+ # -------------------------------------------------------------------------- #
+ # 0D implementation of the freezing model
+ # -------------------------------------------------------------------------- #
+
+ # homogeneous model
+ def_run_0D(self,seed=0):
+ # write into local vars for convenience: geometry
+ A=self.const["A"]
+ V=self.const["V"]
+
+ # density
+ rho_l=self.const["rho_l"]
+
+ # mass of individual components
+ mass=self.const["mass"]
+ mass_water=self.const["mass_water"]
+ mass_solute=self.const["mass_solute"]
+
+ # solute mass fraction
+ w_s=mass_solute/mass
+
+ # specific heat capacities
+ cp_w=self.const["cp_w"]
+ cp_i=self.const["cp_i"]
+ cp_s=self.const["cp_s"]
+ cp_solution=self.const["cp_solution"]
+
+ # soluton concentration and freezing-point depression
+ solid_fraction=self.const["solid_fraction"]
+ T_m=self.const["T_eq"]+273.15
+ k_f=self.const["k_f"]
+ M_s=self.const["M_s"]
+ depression=self.const["depression"]
+ T_eq_l=T_m-depression
+
+ # nucleation kinetics
+ kb=self.const["kb"]
+ b=self.const["b"]
+
+ # latent heat of fusion
+ Dh=self.const["Dh"]
+
+ # heat transfer coefficient shelf
+ K_shelf=self.k["s0"]
+
+ # time step
+ dt=0.1# [s]
+
+ # initial temperature
+ T_0=self.opcond.cooling["start"]+273.15
+ T_k=T_0# [K] temperature profile at t = T_old
+ T_new=T_k
+
+ # shelf temperature and coooling rate
+ T_shelf_cool=self.opcond.tempProfile(dt)+273.15
+
+ # preallocation of results array
+ T_product_cooling=np.zeros_like(T_shelf_cool)
+ time_cooling=np.zeros_like(T_shelf_cool)
+
+ # initial temperature
+ T_old=T_0
+
+ # random numbers
+ np.random.seed(seed)
+ F_rand=np.random.random()
+
+ # preallocate the expected number of nuclei
+ E_t,Nt_cool_end,Nt_sol_end=0,None,None
+
+ # solving for the cooling stage
+ fori,T_shelfinenumerate(T_shelf_cool):
+ # computing new temperature
+ T_new=T_old+dt*(A*K_shelf*(T_shelf-T_old))/(
+ cp_solution*mass
+ )
+ # saving results
+ T_product_cooling[i]=T_new
+ time_cooling[i]=dt*i
+ # update temperature
+ T_old=T_new
+ # check if nucleation occurs stochastically
+ J=kb*(T_eq_l-T_new)**b*(T_eq_l-T_new>0)
+ # frequency of nucleation events
+ K_v=J*V
+ # expected number of nuclei
+ E_t+=K_v*dt
+ # CDF of probability
+ F_nuc=1-np.exp(-E_t)
+ ifself.opcond.cnTemp==None:
+ # Monte-Carlo approach
+ ifF_nuc>F_rand:
+ Nt_cool_end=i
+ break
+ else:
+ # controlled nucleation
+ ifT_new<=(self.opcond.cnTemp+273.15):
+ Nt_cool_end=i
+ break
+
+ # check if nucleation occured
+ ifNt_cool_endisNone:
+ raiseValueError("Nucleation did not occur, prolong process time.")
+
+ # final solution arrays for cooling stage
+ T_product_cooling=T_product_cooling[0:Nt_cool_end]
+ T_shelf_cooling=T_shelf_cool[0:Nt_cool_end]
+ time_cooling=time_cooling[0:Nt_cool_end]
+
+ # mass fraction of ice over time
+ w_i_cooling=np.zeros(len(T_shelf_cooling))
+
+ # 2. Nucleation phase
+ t_nuc=Nt_cool_end*dt
+ T_nuc=T_new
+
+ # save nucleation results
+ self._stats["T_nuc"]=T_nuc-273.15
+ self._stats["t_nuc"]=t_nuc/60
+ self._stats["t_sol"]=None
+ self._stats["t_fr"]=None
+
+ # solving the quadratic equation
+ B=-T_m-T_nuc-Dh*mass_water/(cp_solution*mass)
+ C=(
+ Dh*mass_water*T_m/(cp_solution*mass)
+ -mass_solute*(k_f/M_s)*Dh/(cp_solution*mass)
+ +T_m*T_nuc
+ )
+
+ # equilibrium temperatures determined from quadratic equation (A = 1)
+ T_eq=0.5*(-B-np.sqrt(B**2-4*C))
+
+ # computed mass of ice from the quadratic equation
+ m_i_eq=mass_water-mass_solute*(k_f/M_s)/(T_m-T_eq)
+
+ # solution of ice nucleation stage
+ w_i_nucl=m_i_eq/mass
+
+ # preallocation
+ T_product_solid=np.zeros(len(T_shelf_cool[Nt_cool_end:]))
+ w_i_solid=np.zeros(len(T_shelf_cool[Nt_cool_end:]))
+ time_solid=np.zeros(len(T_shelf_cool[Nt_cool_end:]))
+
+ # saving results of the ice nucleation stage
+ w_i_solid[0]=w_i_nucl
+ T_product_solid[0]=T_eq
+
+ # shelf temperatures
+ T_old=T_eq
+ w_i_new=w_i_nucl
+
+ # solidification stage
+ fori,T_shelfinenumerate(T_shelf_cool[Nt_cool_end:]):
+ # thermal properties
+ cp=cp_s*w_s+cp_i*w_i_new+cp_w*(1-w_s-w_i_new)# [J/kgK]
+ # temperature computation
+ T_new=T_old+(
+ dt
+ *(A*K_shelf*(T_shelf-T_old))
+ *(
+ cp*rho_l*V
+ +(Dh*k_f*mass_solute/M_s)*(1/(T_m-T_old)**2)
+ )
+ **-1
+ )
+ # update temperatures
+ T_old=T_new
+ # save results
+ T_product_solid[i]=T_new
+ time_solid[i]=t_nuc+dt*i
+ # new ice mass fraction
+ w_i_new=(mass_water-(k_f*mass_solute/M_s)/(T_m-T_new))/mass
+ # solidification sigma
+ sigma_new=w_i_new*mass/(mass-mass_solute)
+ # save solidification time
+ if(self._stats["t_sol"]isNone)&(sigma_new>=0.9):
+ Nt_sol_end=i
+ self._stats["t_sol"]=dt*i/60
+ self._stats["t_fr"]=(t_nuc+dt*i)/60
+
+ # check if solidification is completed
+ ifNt_sol_endisNone:
+ raiseValueError("Solidification is not completed, prolong process time.")
+
+ # ice mass fraction evolution
+ w_i_solid=(
+ mass_water-(k_f*mass_solute/M_s)/(T_m-T_product_solid)
+ )/mass
+
+ # time
+ time_results=np.concatenate((time_cooling,time_solid),axis=0)/3600
+
+ # temperatures
+ temp_results=(
+ np.concatenate((T_product_cooling,T_product_solid),axis=0)-273.15
+ )
+ shelf_results=T_shelf_cool-273.15
+
+ # ice mass fraction
+ ice_results=np.concatenate((w_i_cooling,w_i_solid),axis=0)
+
+ # save results
+ self._domain=0
+ self._prop=(T_eq_l,w_s)
+ self._time=time_results
+ self._shelfTemp=shelf_results
+ self._temp=temp_results
+ self._iceMassFraction=ice_results
+
+ return[
+ self._stats["T_nuc"],
+ self._stats["t_nuc"],
+ self._stats["t_sol"],
+ self._stats["t_fr"],
+ ]
+
+ # -------------------------------------------------------------------------- #
+ # 1D implementation of the freezing model
+ # -------------------------------------------------------------------------- #
+
+ # spatial model in 1 dimension
+ def_run_1D(self,seed=0):
+ # write into local vars for convenience: geometry
+ height=self.const["height"]
+ A=self.const["A"]
+ V=self.const["V"]
+
+ # density
+ rho_l=self.const["rho_l"]
+
+ # mass of individual components
+ mass=self.const["mass"]
+ mass_water=self.const["mass_water"]
+ mass_solute=self.const["mass_solute"]
+ w_s=mass_solute/mass
+
+ # heat conductivities
+ lambda_w=self.const["lambda_w"]
+ lambda_i=self.const["lambda_i"]
+ lambda_s=self.const["lambda_s"]
+
+ # specific heat capacities
+ cp_w=self.const["cp_w"]
+ cp_i=self.const["cp_i"]
+ cp_s=self.const["cp_s"]
+ cp_solution=self.const["cp_solution"]
+
+ # soluton concentration and freezing-point depression
+ solid_fraction=self.const["solid_fraction"]
+ T_m=self.const["T_eq"]+273.15
+ k_f=self.const["k_f"]
+ M_s=self.const["M_s"]
+ depression=self.const["depression"]
+ T_eq_l=T_m-depression
+
+ # nucleation kinetics
+ kb=self.const["kb"]
+ b=self.const["b"]
+
+ # general
+ k_B=self.const["k_B"]
+
+ # latent heat of fusion
+ Dh=self.const["Dh"]
+
+ # heat transfer coefficient shelf
+ K_shelf=self.k["s0"]
+
+ # additional VISF parameters
+ ifself.const["configuration"]=="VISF":
+ # vacuum pressure
+ p_vac=self.const["p_vac"]
+ # evaporation coefficient
+ kappa=self.const["kappa"]
+ # latent heat of vaporization for water [J/kg]
+ dHe=self.const["Dh_evaporation"]
+ # mass of water molecule [kg]
+ m_water=self.const["m_water"]
+ # time for vacuum start [h]
+ t_vac_start=self.const["t_vac_start"]
+ # duration of the vacuum [h]
+ t_vac_duration=self.const["t_vac_duration"]
+
+ # [W/mK] effective heat conductivity
+ lambda_eff=solid_fraction*lambda_s+(1-solid_fraction)*lambda_w
+ # heat diffusivity
+ diffusivity_cooling=lambda_eff/(cp_solution*rho_l)
+
+ Nz=30# [/] number of spatial grid points
+ dz=height/Nz# [m] size of spatial step
+ z=np.linspace(0,height,Nz)# [m] z-position array
+
+ # limiting value of alpha for computing dt from CFL condition
+ alpha_max=lambda_i/(cp_i*rho_l)
+
+ # temporal discretization
+ dt=(
+ 0.4*dz**2/alpha_max
+ )# [s] size of time step (determined by the CFL cond.)
+
+ Nt_exp=(
+ int(np.ceil(self.opcond.t_tot/dt))+1
+ )# the total number of timesteps
+
+ # Preallocation of arrays for saving results from the cooling stage.
+ N_save_cool=10000
+ cooling_temp_results=np.zeros((N_save_cool,Nz))
+ cooling_ice_results=np.zeros((N_save_cool,Nz))
+ cooling_shelf_results=np.zeros(N_save_cool)
+ cooling_time_results=np.zeros(N_save_cool)
+ i_save=0
+
+ # initial temperature
+ T_0=self.opcond.cooling["start"]+273.15
+ T_k=np.zeros(Nz)+T_0# [K] temperature profile at t = T_old
+
+ # shelf temperature and coooling rate
+ T_shelf_cool=self.opcond.tempProfile(dt)+273.15
+
+ # random numbers
+ np.random.seed(seed)
+ F_rand=np.random.random()
+
+ # preallocate the expected number of nuclei
+ E_t,Nt_cool_end,Nt_sol_end=0,None,None
+
+ fori,T_shelfinenumerate(T_shelf_cool):
+ # overall heat flux
+ q_overall=K_shelf*(T_shelf-T_k[0])
+ # boundary condition on the bottom of the vial: ghost grid point
+ T_bottom_BC=T_k[0]+q_overall*dz/lambda_eff
+
+ ifself.const["configuration"]=="VISF":
+ # temperature at the top of the liquid
+ T_l=T_k[-1]
+ # vapour temperature
+ T_v=T_k[-1]
+ # vapour pressure
+ p_vap=Utils.vapour_pressure_liquid(T_l)
+
+ # evaporative heat flux (only in the span of vacuum, hold_2)
+ if(dt*i>t_vac_start*3600)&(
+ dt*i<(t_vac_start+t_vac_duration)*3600
+ ):
+ # vapour flux
+ N_w=Utils.vapour_flux(kappa,m_water,k_B,p_vac,p_vap,T_l,T_v)
+ # evaporative heat flux
+ q_e=-N_w*dHe
+ else:
+ q_e=0
+
+ else:
+ q_e=0
+
+ # boundary condition in the top of the vial
+ T_top_BC=T_k[Nz-1]+q_e*dz/lambda_eff
+
+ # updating the temperature at the bottom
+ T_bottom=T_k[0]+(diffusivity_cooling*dt/dz**2)*(
+ T_k[1]-2*T_k[0]+T_bottom_BC
+ )
+
+ # updating the bulk temperatures
+ T_center=T_k[1:Nz-1]+(diffusivity_cooling*dt/dz**2)*(
+ T_k[2:Nz]-2*T_k[1:Nz-1]+T_k[0:Nz-2]
+ )
+
+ # Top point.
+ # updating the top temperature
+ T_top=T_k[Nz-1]+(diffusivity_cooling*dt/dz**2)*(
+ T_top_BC-2*T_k[Nz-1]+T_k[Nz-2]
+ )
+
+ # update temperatures
+ T_k=np.concatenate((T_bottom,T_center,T_top),axis=None)
+
+ # nucleation rate (global aproach)
+ superCooledMask=T_k<T_eq_l
+ J_z=np.zeros(z.shape)
+ J_z[superCooledMask]=kb*(T_eq_l-T_k[superCooledMask])**b
+
+ # compute the frequency of nucleation
+ K_v=A*simps(J_z,z)
+
+ # check the range of LCS = LocalSupercooling
+ LCS_i=T_k<T_eq_l
+ # the not so supercool part
+ LCS_i_r=~LCS_i
+
+ # sparsely saving results
+ ifnp.mod(i,np.ceil(Nt_exp/N_save_cool))==0:
+ # saving temperature results
+ cooling_temp_results[i_save,:]=T_k-273.15
+ # saving shelf temperature results
+ cooling_shelf_results[i_save]=T_shelf-273.15
+ # ice
+ cooling_ice_results[i_save]=0*T_k
+ # saving sparse time array
+ cooling_time_results[i_save]=dt*i
+ # update saving index
+ i_save=i_save+1
+
+ # expected number of nuclei
+ E_t+=K_v*dt
+ # CDF of probability
+ F_nuc=1-np.exp(-E_t)
+ ifself.opcond.cnTemp==None:
+ # Monte-Carlo approach
+ ifF_nuc>F_rand:
+ Nt_cool_end=i
+ t_nuc=dt*i
+ T_nuc=T_k
+ break
+ else:
+ # controlled nucleation
+ ifT_k.any()<=(self.opcond.cnTemp+273.15):
+ Nt_cool_end=i
+ t_nuc=dt*i
+ T_nuc=T_k
+ break
+
+ # check if nucleation occured
+ ifNt_cool_endisNone:
+ raiseValueError("Nucleation did not occur, prolong process time.")
+
+ # kinetic mean nucleation temperature
+ T_nuc_kin=(A/K_v)*simps(T_nuc*J_z,z)
+
+ # save nucleation temperatures
+ self._stats={
+ "T_nuc_min":T_nuc.min()-273.15,
+ "T_nuc_kin":T_nuc_kin-273.15,
+ "T_nuc_mean":np.mean(T_nuc)-273.15,
+ "T_nuc_max":T_nuc.max()-273.15,
+ }
+
+ # save nucleation time
+ self._stats["t_nuc"]=t_nuc/60
+ self._stats["t_sol"]=None
+ self._stats["t_fr"]=None
+
+ " 2. Ice nucleation stage. "
+
+ # solving the quadratic equation (vectorized (A = 1)
+ B=-T_m-T_nuc-Dh*mass_water/(cp_solution*mass)
+ C=(
+ Dh*mass_water*T_m/(cp_solution*mass)
+ -mass_solute*(k_f/M_s)*Dh/(cp_solution*mass)
+ +T_m*T_nuc
+ )
+
+ # lcoal supercooling
+ LCS_i=T_k<T_eq_l
+ # the not so supercool part
+ LCS_i_r=~LCS_i
+
+ # equilibrium temperatures determined from quadratic equation (A = 1)
+ T_eq_sol_1=0.5*(-B-np.sqrt(B**2-4*C))
+
+ # temperaure field after nucleation
+ T_after_nuc=T_k*LCS_i_r+T_eq_sol_1*LCS_i
+
+ # computed mass of ice from the quadratic equation (T_k used just to get
+ # the same dimensions)
+ m_i_nucl_sol_1=mass_water-mass_solute*(k_f/M_s)/(T_m-T_eq_sol_1)
+ m_i_eq=0*T_k*LCS_i_r+m_i_nucl_sol_1*LCS_i
+
+ # save temperature profile after nucleation
+ cooling_temp_results[i_save,:]=T_after_nuc-273.15
+ # saving shelf temperature results
+ cooling_shelf_results[i_save]=T_shelf-273.15
+ # ice
+ cooling_ice_results[i_save]=m_i_eq/(mass_water+mass_solute)
+ # saving sparse time array
+ cooling_time_results[i_save]=dt*i
+
+ i_end=i
+ i_save_end=i_save
+
+ " 3. Solidification stage. "
+
+ N_save_solid=10000
+
+ # arrays of temperature
+ T_k=T_after_nuc# [K] temperature profile at t = T_old
+
+ # array of ice mass fractions
+ w_i_k=m_i_eq/(mass_water+mass_solute)
+ Nt_solid=Nt_exp-i_end
+
+ # Preallocation of arrays for saving results from the cooling stage.
+ solid_temp_results=np.zeros((N_save_solid,Nz))
+ solid_ice_results=np.zeros((N_save_solid,Nz))
+ solid_shelf_results=np.zeros(N_save_solid)
+ solid_time_results=np.zeros(N_save_solid)
+ i_save=0
+
+ # FDM scheme for solving the heat conduction equation
+ fori,T_shelfinenumerate(T_shelf_cool[i_end:]):
+ # local supercooling
+ LCS_i=T_k<T_eq_l
+ # the not so supercool part
+ LCS_i_r=~LCS_i
+
+ " Thermal properties. "
+ # heat capacity
+ cp_solution=(
+ cp_s*solid_fraction
+ +cp_i*w_i_k
+ +cp_w*(1-solid_fraction-w_i_k)
+ )
+ # heat conductivity
+ lambda_eff=lambda_i*w_i_k+lambda_w*(1-w_i_k)
+ # nonlinear capacitance term
+ beta=Dh*k_f*mass_solute/(M_s*rho_l*V*cp_solution)
+ # total nonlinear capacitance term
+ BETA=np.ones(Nz)*LCS_i_r+(1+beta/(T_k-T_m)**2)*LCS_i
+
+ " Boundary condition at the bottom, z = 0."
+ # overall heat flux
+ q_overall=K_shelf*(T_shelf-T_k[0])
+ # boundary condition on the bottom of the vial: ghost grid point
+ T_bottom_BC=T_k[0]+q_overall*dz/lambda_eff[0]
+
+ ifself.const["configuration"]=="VISF":
+ # temperature at the top of the liquid
+ T_l=T_k[-1]
+ # vapour temperature
+ T_v=T_k[-1]
+ # vapour pressure
+ p_vap=Utils.vapour_pressure_solid(T_l)
+
+ # evaporative heat flux (only in the span of vacuum, hold_2)
+ if(t_nuc+dt*i>t_vac_start*3600)&(
+ t_nuc+dt*i<(t_vac_start+t_vac_duration)*3600
+ ):
+ # vapour flux
+ N_w=Utils.vapour_flux(kappa,m_water,k_B,p_vac,p_vap,T_l,T_v)
+ # evaporative heat flux
+ q_e=-N_w*dHe
+ else:
+ q_e=0
+
+ else:
+ q_e=0
+
+ " Boundary condition at the top, z = height."
+ # boundary condition in the top of the vial
+ T_top_BC=T_k[Nz-1]+q_e*dz/lambda_eff[Nz-1]
+
+ " Bottom point. "
+ # update the point at the bottom
+ T_bottom=T_k[0]+(dt/(cp_solution[0]*rho_l))*(
+ (lambda_eff[1]-lambda_eff[0])*(T_k[1]-T_bottom_BC)/(4*dz**2)
+ +lambda_eff[0]*(T_k[1]-2*T_k[0]+T_bottom_BC)/dz**2
+ )*BETA[0]**(-1)
+
+ " Center points. "
+ # update the bulk points
+ T_center=T_k[1:Nz-1]+(dt/(cp_solution[1:Nz-1]*rho_l))*(
+ (lambda_eff[2:Nz]-lambda_eff[0:Nz-2])
+ *(T_k[2:Nz]-T_k[0:Nz-2])
+ /(4*dz**2)
+ +lambda_eff[1:Nz-1]
+ *(T_k[2:Nz]-2*T_k[1:Nz-1]+T_k[0:Nz-2])
+ /dz**2
+ )*BETA[1:Nz-1]**(-1)
+
+ " Top point. "
+ # update the temperature at the top
+ T_top=T_k[Nz-1]+(dt/(cp_solution[Nz-1]*rho_l))*(
+ (lambda_eff[Nz-1]-lambda_eff[Nz-2])
+ *(T_top_BC-T_k[Nz-2])
+ /(4*dz**2)
+ +lambda_eff[Nz-1]
+ *(T_top_BC-2*T_k[Nz-1]+T_k[Nz-2])
+ /dz**2
+ )*BETA[Nz-1]**(-1)
+
+ " Processing the results I: temperatures. "
+ # update temperatures
+ T_k=np.concatenate((T_bottom,T_center,T_top),axis=None)
+
+ " Check for local supercooling. "
+ # check the range of LCS
+ LCS_i=np.where(T_k<T_eq_l,1,0)
+ # the not so supercool region
+ LCS_i_r=np.where(LCS_i,0,1)
+
+ " Processing the results II: ice mass fractions. "
+ # ice mass distribution
+ m_ice=(
+ np.zeros(Nz)*LCS_i_r
+ +(mass_water-mass_solute*(k_f/M_s)/(T_m-T_k))*LCS_i
+ )
+ # ice mass fraction distribution
+ w_i_k=m_ice/mass
+
+ " Saving results of the cooling stage. "
+ # sparse saving
+
+ ifnp.mod(i,np.ceil(Nt_solid/N_save_solid))==0:
+ # temperature results
+ solid_temp_results[i_save,:]=T_k-273.15
+ # ice mass fraction results
+ solid_ice_results[i_save,:]=w_i_k
+ # shelf temperature results
+ solid_shelf_results[i_save]=T_shelf-273.15
+ # saving sparse time array
+ solid_time_results[i_save]=t_nuc+dt*i
+ # updating saving index
+ i_save=i_save+1
+
+ # solidification sigma
+ sigma_new=(1/z[-1])*simps(m_ice,z)/(mass-mass_solute)
+ # save solidification time
+ if(self._stats["t_sol"]isNone)&(sigma_new>=0.9):
+ Nt_sol_end=i
+ self._stats["t_sol"]=dt*i/60
+ self._stats["t_fr"]=(t_nuc+dt*i)/60
+
+ # check if solidification is completed
+ ifNt_sol_endisNone:
+ raiseValueError("Solidification is not completed, prolong process time.")
+
+ # temperature distributions in the product
+ temp_results=np.concatenate(
+ (
+ cooling_temp_results[:i_save_end+1,:],
+ solid_temp_results[:i_save-1,:],
+ ),
+ axis=0,
+ )
+
+ # ice mass fraction results
+ ice_results=np.concatenate(
+ (
+ cooling_ice_results[:i_save_end+1,:],
+ solid_ice_results[:i_save-1,:],
+ ),
+ axis=0,
+ )
+
+ # total time array
+ time_results=(
+ np.concatenate(
+ (
+ cooling_time_results[:i_save_end+1],
+ solid_time_results[:i_save-1],
+ ),
+ axis=0,
+ )
+ /3600
+ )
+
+ # total shelf temperature array
+ shelf_results=np.concatenate(
+ (
+ cooling_shelf_results[:i_save_end+1],
+ solid_shelf_results[:i_save-1],
+ ),
+ axis=0,
+ )
+
+ # save results
+ self._domain=z
+ self._nuc=(t_nuc,T_nuc)
+ self._prop=(T_eq_l,w_s)
+ self._time=time_results
+ self._shelfTemp=shelf_results
+ self._temp=temp_results
+ self._iceMassFraction=ice_results
+
+ return[
+ self._stats["T_nuc_min"],
+ self._stats["T_nuc_kin"],
+ self._stats["T_nuc_mean"],
+ self._stats["T_nuc_max"],
+ self._stats["t_nuc"],
+ self._stats["t_sol"],
+ self._stats["t_fr"],
+ ]
+
+ # -------------------------------------------------------------------------- #
+ # 2D implementation of the freezing model
+ # -------------------------------------------------------------------------- #
+
+ # spatial model in 2 dimensions
+ def_run_2D(self,seed=0):
+ # write into local vars for convenience: geometry
+ height=self.const["height"]
+ diameter=self.const["diameter"]
+ V=self.const["V"]
+ radius=diameter/2
+
+ # density
+ rho_l=self.const["rho_l"]
+
+ # mass of individual components
+ mass=self.const["mass"]
+ mass_water=self.const["mass_water"]
+ mass_solute=self.const["mass_solute"]
+ w_s=mass_solute/mass
+
+ # heat conductivities
+ lambda_w=self.const["lambda_w"]
+ lambda_i=self.const["lambda_i"]
+ lambda_s=self.const["lambda_s"]
+
+ # specific heat capacities
+ cp_w=self.const["cp_w"]
+ cp_i=self.const["cp_i"]
+ cp_s=self.const["cp_s"]
+ cp_solution=self.const["cp_solution"]
+
+ # soluton concentration and freezing-point depression
+ solid_fraction=self.const["solid_fraction"]
+ T_m=self.const["T_eq"]+273.15
+ k_f=self.const["k_f"]
+ M_s=self.const["M_s"]
+ depression=self.const["depression"]
+ T_eq_l=T_m-depression
+
+ # nucleation kinetics
+ kb=self.const["kb"]
+ b=self.const["b"]
+
+ # general
+ k_B=self.const["k_B"]
+
+ # latent heat of fusion
+ Dh=self.const["Dh"]
+
+ # shelf heat transfer coefficient
+ K_shelf=self.k["s0"]
+
+ # additional VISF parameters
+ ifself.const["configuration"]=="VISF":
+ # vacuum pressure
+ p_vac=self.const["p_vac"]
+ # evaporation coefficient
+ kappa=self.const["kappa"]
+ # latent heat of vaporization for water [J/kg]
+ dHe=self.const["Dh_evaporation"]
+ # mass of water molecule [kg]
+ m_water=self.const["m_water"]
+ # time for vacuum start [h]
+ t_vac_start=self.const["t_vac_start"]
+ # duration of the vacuum [h]
+ t_vac_duration=self.const["t_vac_duration"]
+
+ # additional jacket-ramped freezing parameters
+ ifself.const["configuration"]=="jacket":
+ # vacuum pressure
+ air_gap=self.const["air_gap"]
+ # evaporation coefficient
+ lambda_air=self.const["lambda_air"]
+ # compute the wall heat transfer coefficient
+ K_wall=(1/K_shelf+air_gap/lambda_air)**(-1)# [W/m2K]
+
+ # [W/mK] effective heat conductivity
+ k_eff=solid_fraction*lambda_s+(1-solid_fraction)*lambda_w
+
+ # heat diffusivity
+ alpha=k_eff/(cp_solution*rho_l)
+
+ Nz=30# [/] number of spatial grid points
+ dz=height/Nz# [m] size of spatial step
+ z=np.linspace(0,height,Nz)# [m] z-position array
+
+ # spatal discretization in the radial direction
+ Nr=15# [/] number of spatial grid
+ dr=radius/Nr# [m] size of spatial step
+ r=np.linspace(0,radius,Nr)# [m] z-position array
+
+ # limiting value of alpha for computing dt from CFL condition
+ alpha_max=lambda_i/(cp_i*rho_l)
+
+ # temporal discretization
+ dt=(
+ (0.4/alpha_max)*(dz**2*dr**2)/(dr**2+dz**2)
+ )# [s] size of time step
+
+ Nt_exp=(
+ int(np.ceil(self.opcond.t_tot/dt))+1
+ )# the total number of timesteps
+
+ # Preallocation of arrays for saving results from the cooling stage.
+ N_save_cool=10000
+ cooling_temp_results=np.zeros((N_save_cool,Nz,Nr))
+ cooling_ice_results=np.zeros((N_save_cool,Nz,Nr))
+ cooling_shelf_results=np.zeros(N_save_cool)
+ cooling_time_results=np.zeros(N_save_cool)
+ i_save=0
+
+ # initial temperature
+ T_0=self.opcond.cooling["start"]+273.15
+ T_k=np.zeros((Nz,Nr))+T_0# [K] temperature profile at t = T_old
+ T_new=T_k
+
+ # shelf temperature and coooling rate
+ T_shelf_cool=self.opcond.tempProfile(dt)+273.15
+
+ # random numbers
+ np.random.seed(seed)
+ F_rand=np.random.random()
+
+ # preallocate the expected number of nuclei
+ E_t,Nt_cool_end,Nt_sol_end=0,None,None
+
+ # finite difference method for solving the heat conduction equation
+ fori,T_shelfinenumerate(T_shelf_cool):
+ "Boundary condition at the bottom, z = 0."
+ # overall heat flux
+ q_overall=K_shelf*(T_shelf-T_k[0,:])
+ # boundary condition on the bottom of the vial: ghost grid point
+ T_bottom=T_k[0,:]+q_overall*dz/k_eff
+
+ ifself.const["configuration"]=="VISF":
+ # temperature at the top of the liquid
+ T_l=T_k[-1,:]
+ # vapour temperature
+ T_v=T_k[-1,:]
+ # vapour pressure
+ p_vap=Utils.vapour_pressure_solid(T_l)
+
+ # evaporative heat flux (only in the span of vacuum, hold_2)
+ if(dt*i>t_vac_start*3600)&(
+ dt*i<(t_vac_start+t_vac_duration)*3600
+ ):
+ # vapour flux
+ N_w=Utils.vapour_flux(kappa,m_water,k_B,p_vac,p_vap,T_l,T_v)
+ # evaporative heat flux
+ q_e=-N_w*dHe
+ else:
+ q_e=0
+
+ else:
+ q_e=0
+
+ # boundary condition in the top of the vial
+ T_top=T_k[Nz-1,:]+q_e*dz/k_eff
+
+ " Boundary conditions at the edge, r = R. "
+
+ ifself.const["configuration"]=="jacket":
+ # jacket heat flux
+ q_jacket=K_wall*(T_shelf-T_k[:,Nr-1])
+
+ else:
+ # no cooling from jacket
+ q_jacket=0
+
+ # boundary condition in the top of the vial
+ T_edge=T_k[:,Nr-1]+q_jacket*dz/k_eff
+
+ " Boundary condition in the centre, r = 0. "
+ # axial symmetry
+ T_center=T_k[:,0]
+
+ " Bottom boundary points: z = 0. "
+ # bottom-centre: r = 0
+ T_new[0,0]=T_k[0,0]+alpha*dt*(
+ 2*(T_k[0,1]-2*T_k[0,0]+T_center[0])/dr**2
+ +(T_k[1,0]-2*T_k[0,0]+T_bottom[0])/dz**2
+ )
+
+ # bottom corner (left or right - symmetry): r = R
+ T_new[0,Nr-1]=T_k[0,Nr-1]+alpha*dt*(
+ (1/r[Nr-1])*(T_edge[0]-T_k[0,Nr-2])/(2*dr)
+ +(T_edge[0]-2*T_k[0,Nr-1]+T_k[0,Nr-2])/dr**2
+ +(T_k[1,Nr-1]-2*T_k[0,Nr-1]+T_bottom[Nr-1])/dz**2
+ )
+
+ # the rest of the bottom: 0 < r < R
+ T_new[0,1:Nr-1]=T_k[0,1:Nr-1]+alpha*dt*(
+ (1/r[1:Nr-1])*(T_k[0,2:Nr]-T_k[0,0:Nr-2])/(2*dr)
+ +(T_k[0,2:Nr]-2*T_k[0,1:Nr-1]+T_k[0,0:Nr-2])/dr**2
+ +(T_k[1,1:Nr-1]-2*T_k[0,1:Nr-1]+T_bottom[1:Nr-1])
+ /dz**2
+ )
+
+ " Top boundary points: z = H. "
+ # top centre: r = 0
+ T_new[Nz-1,0]=T_k[Nz-1,0]+alpha*dt*(
+ 2*(T_k[Nz-1,1]-2*T_k[Nz-1,0]+T_center[-1])/dr**2
+ +(T_top[0]-2*T_k[Nz-1,0]+T_k[Nz-2,0])/dz**2
+ )
+
+ # top corner (left and right - symmetry): r = R
+ T_new[Nz-1,Nr-1]=T_k[Nz-1,Nr-1]+alpha*dt*(
+ (1/r[Nr-1])*(T_edge[Nz-1]-T_k[Nz-1,Nr-2])/(2*dr)
+ +(T_edge[Nz-1]-2*T_k[Nz-1,Nr-1]+T_k[Nz-1,Nr-2])
+ /dr**2
+ +(T_top[Nr-1]-2*T_k[Nz-1,Nr-1]+T_k[Nz-2,Nr-1])
+ /dz**2
+ )
+
+ # the rest of the vial's top: 0 < r < R
+ T_new[Nz-1,1:Nr-1]=T_k[Nz-1,1:Nr-1]+alpha*dt*(
+ (1/r[1:Nr-1])
+ *(T_k[Nz-1,2:Nr]-T_k[Nz-1,0:Nr-2])
+ /(2*dr)
+ +(
+ T_k[Nz-1,2:Nr]
+ -2*T_k[Nz-1,1:Nr-1]
+ +T_k[Nz-1,0:Nr-2]
+ )
+ /dr**2
+ +(
+ T_top[1:Nr-1]
+ -2*T_k[Nz-1,1:Nr-1]
+ +T_k[Nz-2,1:Nr-1]
+ )
+ /dz**2
+ )
+
+ " Edge boundary points: r = R. "
+ # curved cylinder surface in contact with air: 0 < z < H
+ T_new[1:Nz-1,Nr-1]=T_k[1:Nz-1,Nr-1]+alpha*dt*(
+ (1/r[Nr-1])
+ *(T_edge[1:Nz-1]-T_k[1:Nz-1,Nr-2])
+ /(2*dr)
+ +(
+ T_edge[1:Nz-1]
+ -2*T_k[1:Nz-1,Nr-1]
+ +T_k[1:Nz-1,Nr-2]
+ )
+ /dr**2
+ +(
+ T_k[2:Nz,Nr-1]
+ -2*T_k[1:Nz-1,Nr-1]
+ +T_k[0:Nz-2,Nr-1]
+ )
+ /dz**2
+ )
+
+ " Center boundary points: r = 0. "
+ # symmetry boundary condition in the centre of the cylinder: 0 < z < H
+ T_new[1:Nz-1,0]=T_k[1:Nz-1,0]+alpha*dt*(
+ 2
+ *(T_k[1:Nz-1,1]-2*T_k[1:Nz-1,0]+T_center[1:Nz-1])
+ /dr**2
+ +(T_k[2:Nz,0]-2*T_k[1:Nz-1,0]+T_k[0:Nz-2,0])/dz**2
+ )
+
+ " Bulk of the vial points: : 0 < z < H & 0 < r < R. "
+ # all the other points in the domain
+ T_new[1:Nz-1,1:Nr-1]=T_k[1:Nz-1,1:Nr-1]+alpha*dt*(
+ (1/r[1:Nr-1])
+ *(T_k[1:Nz-1,2:Nr]-T_k[1:Nz-1,0:Nr-2])
+ /(2*dr)
+ +(
+ T_k[1:Nz-1,2:Nr]
+ -2*T_k[1:Nz-1,1:Nr-1]
+ +T_k[1:Nz-1,0:Nr-2]
+ )
+ /dr**2
+ +(
+ T_k[2:Nz,1:Nr-1]
+ -2*T_k[1:Nz-1,1:Nr-1]
+ +T_k[0:Nz-2,1:Nr-1]
+ )
+ /dz**2
+ )
+
+ " Processing the results of the time iteration. "
+ # update temperatures
+ T_k=T_new
+
+ # Stochastic nucleation for the entire vial.
+ # nucleation rate (global aproach)
+ superCooledMask=T_k<T_eq_l
+ J_z_r=np.zeros((Nz,Nr))
+ J_z_r[superCooledMask]=kb*(T_eq_l-T_k[superCooledMask])**b
+
+ # frequency of nucleation events itegrated over z
+ K_r=2*np.pi*simps(r*J_z_r,r)
+ # frequency of nucleation events itegrated over r
+ K_v=simps(K_r,z)
+
+ # Check for local supercooling.
+ # check the range of LCS = LocalSupercooling
+ LCS_i=T_k<T_eq_l
+ # the not so supercool part
+ LCS_i_r=~LCS_i
+
+ # Saving results of the cooling stage.
+ # sparsely saving results
+ ifnp.mod(i,np.ceil(Nt_exp/N_save_cool))==0:
+ # saving temperature results
+ cooling_temp_results[i_save,:,:]=T_k-273.15
+ # saving shelf temperature results
+ cooling_shelf_results[i_save]=T_shelf-273.15
+ # ice
+ cooling_ice_results[i_save,:,:]=0*T_k
+ # saving sparse time array
+ cooling_time_results[i_save]=dt*i
+ # update saving index
+ i_save=i_save+1
+
+ # expected number of nuclei
+ E_t+=K_v*dt
+ # CDF of probability
+ F_nuc=1-np.exp(-E_t)
+ ifself.opcond.cnTemp==None:
+ # Monte-Carlo approach
+ ifF_nuc>F_rand:
+ Nt_cool_end=i
+ t_nuc=dt*i
+ T_nuc=T_k
+ break
+ else:
+ # controlled nucleation
+ ifT_k.any()<=(self.opcond.cnTemp+273.15):
+ Nt_cool_end=i
+ t_nuc=dt*i
+ T_nuc=T_k
+ break
+
+ # check if nucleation occured
+ ifNt_cool_endisNone:
+ raiseValueError("Nucleation did not occur, prolong process time.")
+
+ # compute the mean kinetic mean nucleation temperature
+ f_z=2*np.pi*simps(r*T_nuc*J_z_r,r)
+ T_nuc_kin=(1/K_v)*simps(f_z,z)
+
+ # save nucleation temperatures
+ self._stats={
+ "T_nuc_min":T_nuc.min()-273.15,
+ "T_nuc_kin":T_nuc_kin-273.15,
+ "T_nuc_mean":np.mean(T_nuc)-273.15,
+ "T_nuc_max":T_nuc.max()-273.15,
+ }
+
+ # save nucleation time
+ self._stats["t_nuc"]=t_nuc/60
+ self._stats["t_sol"]=None
+ self._stats["t_fr"]=None
+
+ # solving the quadratic equation (vectorized (A = 1)
+ B=-T_m-T_nuc-Dh*mass_water/(cp_solution*mass)
+ C=(
+ Dh*mass_water*T_m/(cp_solution*mass)
+ -mass_solute*(k_f/M_s)*Dh/(cp_solution*mass)
+ +T_m*T_nuc
+ )
+
+ # lcoal supercooling
+ LCS_i=T_k<T_eq_l
+ # the not so supercool part
+ LCS_i_r=~LCS_i
+
+ # equilibrium temperatures determined from quadratic equation (A = 1)
+ T_eq_sol_1=0.5*(-B-np.sqrt(B**2-4*C))
+
+ # temperaure field after nucleation
+ T_after_nuc=T_k*LCS_i_r+T_eq_sol_1*LCS_i
+
+ # computed mass of ice from the quadratic equation (T_k used just to get
+ # the same dimensions)
+ m_i_nucl_sol_1=mass_water-mass_solute*(k_f/M_s)/(T_m-T_eq_sol_1)
+ m_i_eq=0*T_k*LCS_i_r+m_i_nucl_sol_1*LCS_i
+
+ cooling_temp_results[i_save,:,:]=T_after_nuc-273.15
+ # saving shelf temperature results
+ cooling_shelf_results[i_save]=T_shelf-273.15
+ # ice
+ cooling_ice_results[i_save]=m_i_eq/(mass_water+mass_solute)
+ # saving sparse time array
+ cooling_time_results[i_save]=dt*i
+ i_end=i
+ i_save_end=i_save
+ N_save_solid=10000
+
+ # arrays of temperature
+ T_k=T_after_nuc# [K] temperature profile at t = T_old
+
+ # array of ice mass fractions
+ w_i_new=m_i_eq/(mass_water+mass_solute)
+
+ Nt_solid=Nt_exp-i_end
+
+ # Preallocation of arrays for saving results from the cooling stage.
+ solid_temp_results=np.zeros((N_save_solid,Nz,Nr))
+ solid_ice_results=np.zeros((N_save_solid,Nz,Nr))
+ solid_shelf_results=np.zeros(N_save_solid)
+ solid_time_results=np.zeros(N_save_solid)
+ i_save=0
+
+ # FDM scheme for solving the heat conduction equation
+ fori,T_shelfinenumerate(T_shelf_cool[i_end:]):
+ # thermal properties
+ # heat capacity
+ cp_eff=(
+ cp_s*solid_fraction
+ +cp_i*w_i_new
+ +cp_w*(1-solid_fraction-w_i_new)
+ )
+ # heat conductivity
+ k_eff=lambda_i*w_i_new+lambda_w*(1-w_i_new)
+ # heat diffusivity
+ alpha=k_eff/(cp_eff*rho_l)
+ # beta parameter beta = beta(t,z,r)
+ beta=Dh*k_f*mass_solute/(M_s*rho_l*V*cp_eff)
+ # the total nonlinear capacitance term
+ BETA=np.ones((Nz,Nr))*LCS_i_r+(1+beta/(T_k-T_m)**2)*LCS_i
+
+ # Boundary condition at the bottom, z = 0.
+ # overall heat flux
+ q_overall=K_shelf*(T_shelf-T_k[0,:])
+ # q_overall = sigma_B*(T_shelf**4 - T_old[0,:]**4)
+ # boundary condition on the bottom of the vial: ghost grid point
+ T_bottom=T_k[0,:]+q_overall*dz/k_eff[0,:]
+
+ ifself.const["configuration"]=="VISF":
+ # temperature at the top of the liquid
+ T_l=T_k[-1,:]
+ # vapour temperature
+ T_v=T_k[-1,:]
+ # vapour pressure
+ p_vap=Utils.vapour_pressure_solid(T_l)
+
+ # evaporative heat flux (only in the span of vacuum, hold_2)
+ if(t_nuc+dt*i>t_vac_start*3600)&(
+ t_nuc+dt*i<(t_vac_start+t_vac_duration)*3600
+ ):
+ # vapour flux
+ N_w=Utils.vapour_flux(kappa,m_water,k_B,p_vac,p_vap,T_l,T_v)
+ # evaporative heat flux
+ q_e=-N_w*dHe
+ else:
+ q_e=0
+
+ else:
+ q_e=0
+
+ # boundary condition in the top of the vial
+ T_top=T_k[Nz-1,:]+q_e*dz/k_eff[Nz-1,:]
+
+ # boundary condition at the side of the vial
+ ifself.const["configuration"]=="jacket":
+ # jacket heat flux
+ q_jacket=K_wall*(T_shelf-T_k[:,Nr-1])
+
+ else:
+ # no cooling from jacket
+ q_jacket=0
+
+ # boundary condition in the top of the vial
+ T_edge=T_k[:,Nr-1]+q_jacket*dz/k_eff[:,Nr-1]
+
+ # Boundary condition in the centre, r = 0.
+ # axial symmetry
+ T_center=T_k[:,0]
+
+ # Bottom boundary points: z = 0.
+ # bottom-centre: r = 0
+ T_new[0,0]=T_k[0,0]+(dt/(cp_eff[0,0]*rho_l))*(
+ 2*k_eff[0,0]*(T_k[0,1]-2*T_k[0,0]+T_center[0])/dr**2
+ +(k_eff[0,1]-k_eff[0,0])*(T_k[0,1]-T_center[0])/(4*dr**2)
+ +(k_eff[1,0]-k_eff[0,0])*(T_k[1,0]-T_bottom[0])/(4*dz**2)
+ +k_eff[0,0]*(T_k[1,0]-2*T_k[0,0]+T_bottom[0])/dz**2
+ )*BETA[0,0]**(-1)
+
+ # bottom corner (left or right - symmetry): r = R
+ T_new[0,Nr-1]=T_k[0,Nr-1]+(dt/(cp_eff[0,Nr-1]*rho_l))*(
+ (k_eff[0,Nr-1]/r[Nr-1])*(T_edge[0]-T_k[0,Nr-2])/(2*dr)
+ +(k_eff[0,Nr-1]-k_eff[0,Nr-2])
+ *(T_edge[0]-T_k[0,Nr-1])
+ /(4*dr**2)
+ +(k_eff[1,Nr-1]-k_eff[0,Nr-1])
+ *(T_k[1,Nr-1]-T_bottom[Nr-1])
+ /(4*dz**2)
+ +k_eff[0,Nr-1]
+ *(T_edge[0]-2*T_k[0,Nr-1]+T_k[0,Nr-2])
+ /dr**2
+ +k_eff[0,Nr-1]
+ *(T_k[1,Nr-1]-2*T_k[0,Nr-1]+T_bottom[Nr-1])
+ /dz**2
+ )*BETA[0,Nr-1]**(-1)
+
+ # the rest of the bottom: 0 < r < R
+ T_new[0,1:Nr-1]=T_k[0,1:Nr-1]+(
+ dt/(cp_eff[0,1:Nr-1]*rho_l)
+ )*(
+ (k_eff[0,1:Nr-1]/r[1:Nr-1])
+ *(T_k[0,2:Nr]-T_k[0,0:Nr-2])
+ /(2*dr)
+ +(k_eff[0,2:Nr]-k_eff[0,0:Nr-2])
+ *(T_k[0,2:Nr]-T_k[0,0:Nr-2])
+ /(4*dr**2)
+ +(k_eff[1,1:Nr-1]-k_eff[0,1:Nr-1])
+ *(T_k[1,1:Nr-1]-T_bottom[1:Nr-1])
+ /(4*dz**2)
+ +k_eff[0,1:Nr-1]
+ *(T_k[0,2:Nr]-2*T_k[0,1:Nr-1]+T_k[0,0:Nr-2])
+ /dr**2
+ +k_eff[0,1:Nr-1]
+ *(T_k[1,1:Nr-1]-2*T_k[0,1:Nr-1]+T_bottom[1:Nr-1])
+ /dz**2
+ )*BETA[
+ 0,1:Nr-1
+ ]**(
+ -1
+ )
+
+ # Top boundary points: z = H.
+ # top centre: r = 0
+ T_new[Nz-1,0]=T_k[Nz-1,0]+(dt/(cp_eff[Nz-1,0]*rho_l))*(
+ 2
+ *k_eff[Nz-1,0]
+ *(T_k[Nz-1,1]-2*T_k[Nz-1,0]+T_center[Nz-1])
+ /dr**2
+ +(k_eff[Nz-1,1]-k_eff[Nz-1,0])
+ *(T_k[Nz-1,1]-T_center[Nz-1])
+ /(4*dr**2)
+ +(k_eff[Nz-1,0]-k_eff[Nz-2,0])
+ *(T_top[0]-T_k[Nz-2,0])
+ /(4*dz**2)
+ +k_eff[Nz-1,0]
+ *(T_top[0]-2*T_k[Nz-1,0]+T_k[Nz-2,0])
+ /dz**2
+ )*BETA[Nz-1,0]**(-1)
+
+ # top edge corner (left or right - symmetry): r = R
+ T_new[Nz-1,Nr-1]=T_k[Nz-1,Nr-1]+(
+ dt/(cp_eff[Nz-1,Nr-1]*rho_l)
+ )*(
+ (k_eff[Nz-1,Nr-1]/r[Nr-1])
+ *(T_edge[Nz-1]-T_k[Nz-1,Nr-2])
+ /(2*dr)
+ +(k_eff[Nz-1,Nr-1]-k_eff[Nz-1,Nr-2])
+ *(T_edge[Nz-1]-T_k[Nz-1,Nr-2])
+ /(4*dr**2)
+ +(k_eff[Nz-1,Nr-1]-k_eff[Nz-2,Nr-1])
+ *(T_top[Nr-1]-T_k[Nz-2,Nr-1])
+ /(4*dz**2)
+ +k_eff[Nz-1,Nr-1]
+ *(T_edge[Nz-1]-2*T_k[Nz-1,Nr-1]+T_k[Nz-1,Nr-2])
+ /dr**2
+ +k_eff[Nz-1,Nr-1]
+ *(T_top[Nr-1]-2*T_k[Nz-1,Nr-1]+T_k[Nz-2,Nr-1])
+ /dz**2
+ )*BETA[
+ Nz-1,Nr-1
+ ]**(
+ -1
+ )
+
+ # the rest of the vial's top: 0 < r < R
+ T_new[Nz-1,1:Nr-1]=T_k[Nz-1,1:Nr-1]+(
+ dt/(cp_eff[Nz-1,1:Nr-1]*rho_l)
+ )*(
+ (k_eff[Nz-1,1:Nr-1]/r[1:Nr-1])
+ *(T_k[Nz-1,2:Nr]-T_k[Nz-1,0:Nr-2])
+ /(2*dr)
+ +(k_eff[Nz-1,2:Nr]-k_eff[Nz-1,0:Nr-2])
+ *(T_k[Nz-1,2:Nr]-T_k[Nz-1,0:Nr-2])
+ /(4*dr**2)
+ +(k_eff[Nz-1,1:Nr-1]-k_eff[Nz-2,1:Nr-1])
+ *(T_top[1:Nr-1]-T_k[Nz-2,1:Nr-1])
+ /(4*dz**2)
+ +k_eff[Nz-1,1:Nr-1]
+ *(
+ T_k[Nz-1,2:Nr]
+ -2*T_k[Nz-1,1:Nr-1]
+ +T_k[Nz-1,0:Nr-2]
+ )
+ /dr**2
+ +k_eff[Nz-1,1:Nr-1]
+ *(
+ T_top[1:Nr-1]
+ -2*T_k[Nz-1,1:Nr-1]
+ +T_k[Nz-2,1:Nr-1]
+ )
+ /dz**2
+ )*BETA[
+ Nz-1,1:Nr-1
+ ]**(
+ -1
+ )
+
+ # Outer boundary points: r = R.
+ # curved cylinder surface in contact with air: 0 < z < H
+ T_new[1:Nz-1,Nr-1]=T_k[1:Nz-1,Nr-1]+(
+ dt/(cp_eff[1:Nz-1,Nr-1]*rho_l)
+ )*(
+ (k_eff[1:Nz-1,Nr-1]/r[Nr-1])
+ *(T_edge[1:Nz-1]-T_k[1:Nz-1,Nr-2])
+ /(2*dr)
+ +(k_eff[1:Nz-1,Nr-1]-k_eff[1:Nz-1,Nr-2])
+ *(T_edge[1:Nz-1]-T_k[1:Nz-1,Nr-2])
+ /(4*dr**2)
+ +(k_eff[2:Nz,Nr-1]-k_eff[0:Nz-2,Nr-1])
+ *(T_k[2:Nz,Nr-1]-T_k[0:Nz-2,Nr-1])
+ /(4*dz**2)
+ +k_eff[1:Nz-1,Nr-1]
+ *(
+ T_edge[1:Nz-1]
+ -2*T_k[1:Nz-1,Nr-1]
+ +T_k[1:Nz-1,Nr-2]
+ )
+ /dr**2
+ +k_eff[1:Nz-1,Nr-1]
+ *(
+ T_k[2:Nz,Nr-1]
+ -2*T_k[1:Nz-1,Nr-1]
+ +T_k[0:Nz-2,Nr-1]
+ )
+ /dz**2
+ )*BETA[
+ 1:Nz-1,Nr-1
+ ]**(
+ -1
+ )
+
+ # Center boundary points: r = R.
+ # symmetry boundary condition in the centre of the cylinder: 0 < z < H
+ T_new[1:Nz-1,0]=T_k[1:Nz-1,0]+(
+ dt/(cp_eff[1:Nz-1,0]*rho_l)
+ )*(
+ 2
+ *k_eff[1:Nz-1,0]
+ *(T_k[1:Nz-1,1]-2*T_k[1:Nz-1,0]+T_center[1:Nz-1])
+ /dr**2
+ +(k_eff[1:Nz-1,1]-k_eff[1:Nz-1,0])
+ *(T_k[1:Nz-1,1]-T_center[1:Nz-1])
+ /(4*dr**2)
+ +(k_eff[2:Nz,0]-k_eff[0:Nz-2,0])
+ *(T_k[2:Nz,0]-T_k[0:Nz-2,0])
+ /(4*dz**2)
+ +k_eff[1:Nz-1,0]
+ *(T_k[2:Nz,0]-2*T_k[1:Nz-1,0]+T_k[0:Nz-2,0])
+ /dz**2
+ )*BETA[
+ 1:Nz-1,0
+ ]**(
+ -1
+ )
+
+ # Bulk of the vial points: : 0 < z < H & 0 < r < R.
+ # all the other points in the domain
+ T_new[1:Nz-1,1:Nr-1]=T_k[1:Nz-1,1:Nr-1]+(
+ dt/(cp_eff[1:Nz-1,1:Nr-1]*rho_l)
+ )*(
+ (k_eff[1:Nz-1,1:Nr-1]/r[1:Nr-1])
+ *(T_k[1:Nz-1,2:Nr]-T_k[1:Nz-1,0:Nr-2])
+ /(2*dr)
+ +(k_eff[1:Nz-1,2:Nr]-k_eff[1:Nz-1,0:Nr-2])
+ *(T_k[1:Nz-1,2:Nr]-T_k[1:Nz-1,0:Nr-2])
+ /(4*dr**2)
+ +(k_eff[2:Nz,1:Nr-1]-k_eff[0:Nz-2,1:Nr-1])
+ *(T_k[2:Nz,1:Nr-1]-T_k[0:Nz-2,1:Nr-1])
+ /(4*dz**2)
+ +k_eff[1:Nz-1,1:Nr-1]
+ *(
+ T_k[1:Nz-1,2:Nr]
+ -2*T_k[1:Nz-1,1:Nr-1]
+ +T_k[1:Nz-1,0:Nr-2]
+ )
+ /dr**2
+ +k_eff[1:Nz-1,1:Nr-1]
+ *(
+ T_k[2:Nz,1:Nr-1]
+ -2*T_k[1:Nz-1,1:Nr-1]
+ +T_k[0:Nz-2,1:Nr-1]
+ )
+ /dz**2
+ )*BETA[
+ 1:Nz-1,1:Nr-1
+ ]**(
+ -1
+ )
+
+ # update temperatures
+ T_k=T_new
+
+ # Check for local supercooling.
+ LCS_i=T_k<T_eq_l
+ # the not so supercool part
+ LCS_i_r=~LCS_i
+
+ # ice mass distribution
+ m_ice=(
+ np.zeros((Nz,Nr))*LCS_i_r
+ +(mass_water-mass_solute*(k_f/M_s)/(T_m-T_new))*LCS_i
+ )
+ # ice mass fraction distribution
+ w_i_new=m_ice/(mass_water+mass_solute)
+
+ # sparse saving
+ ifnp.mod(i,np.ceil(Nt_solid/N_save_solid))==0:
+ # saving temperature profiles
+ solid_temp_results[i_save,:,:]=T_new-273.15
+ # saving ice mass fractions
+ solid_ice_results[i_save,:,:]=w_i_new
+ # shelf temperature
+ solid_shelf_results[i_save]=T_shelf-273.15
+ # saving sparse time array
+ solid_time_results[i_save]=t_nuc+dt*i
+ # update saving index
+ i_save=i_save+1
+
+ # solidification sigma
+ w_i_r=2*np.pi*simps(r*w_i_new,r)
+ w_i=simps(w_i_r,z)
+ sigma_new=(
+ (1/(np.pi*r[-1]**2*z[-1]))*w_i*mass/(mass-mass_solute)
+ )
+
+ # save solidification time
+ if(self._stats["t_sol"]isNone)&(sigma_new>=0.9):
+ Nt_sol_end=i
+ self._stats["t_sol"]=dt*i/60
+ self._stats["t_fr"]=(t_nuc+dt*i)/60
+
+ # check if solidification is completed
+ ifNt_sol_endisNone:
+ raiseValueError("Solidification is not completed, prolong process time.")
+
+ # temperature distributions in the product
+ temp_results=np.concatenate(
+ (
+ cooling_temp_results[:i_save_end+1,:],
+ solid_temp_results[:i_save-1,:],
+ ),
+ axis=0,
+ )
+
+ # ice mass fraction results
+ ice_results=np.concatenate(
+ (
+ cooling_ice_results[:i_save_end+1,:],
+ solid_ice_results[:i_save-1,:],
+ ),
+ axis=0,
+ )
+
+ # total time array
+ time_results=(
+ np.concatenate(
+ (
+ cooling_time_results[:i_save_end+1],
+ solid_time_results[:i_save-1],
+ ),
+ axis=0,
+ )
+ /3600
+ )
+
+ # total shelf temperature array
+ shelf_results=np.concatenate(
+ (
+ cooling_shelf_results[:i_save_end+1],
+ solid_shelf_results[:i_save-1],
+ ),
+ axis=0,
+ )
+
+ # save reuslts
+ self._domain=(z,r)
+ self._prop=(T_eq_l,w_s)
+ self._nuc=(t_nuc,T_nuc)
+ self._time=time_results
+ self._shelfTemp=shelf_results
+ self._temp=temp_results
+ self._iceMassFraction=ice_results
+
+ return[
+ self._stats["T_nuc_min"],
+ self._stats["T_nuc_kin"],
+ self._stats["T_nuc_mean"],
+ self._stats["T_nuc_max"],
+ self._stats["t_nuc"],
+ self._stats["t_sol"],
+ self._stats["t_fr"],
+ ]
+
+ # -------------------------------------------------------------------------- #
+ # Plotting functions
+ # -------------------------------------------------------------------------- #
+
+ # plot distribution of the desired variable
+
[docs]defplot_cdf(self,what:str="t_nuc",comp:bool=False):
+ """Create CDF plots for Snowing object.
+
+ Args:
+ what (str, optional): What to plot, i.e., keys of the stats dict.
+ Valid options are t_nuc, T_nuc, t_sol, t_fr.
+ Defaults to "t_nuc". Also accepting t_nucleation, T_nucleation
+ and t_solidification.
+ comp (bool, optional): If complementary CDF is plotted. Defaults to False.
+
+ Args:
+ what (str, optional): What to plot, i.e., keys of the stats dict.
+ Valid options are t_nuc, T_nuc, t_sol, t_fr.
+ Defaults to "t_nuc". Also accepting t_nucleation, T_nucleation
+ and t_solidification.
+
+ Raises:
+ NotImplementedError: This plot not available for single vial simulation.
+ """
+
+ # plot settings
+ plt.rcParams.update(
+ {
+ "axes.axisbelow":True,
+ "xtick.direction":"in",
+ "ytick.direction":"in",
+ "xtick.bottom":True,
+ "xtick.top":True,
+ "ytick.left":True,
+ "ytick.right":True,
+ }
+ )
+
+ # enable this plot only for Nrep > 1
+ ifself.Nrep<=1:
+ raiseNotImplementedError(
+ "This plot only available for multiple simulations. "
+ +"Run with Nrep > 1."
+ )
+
+ # plot settings
+ plt.rcParams.update(
+ {
+ "lines.linewidth":3,
+ "axes.axisbelow":True,
+ "xtick.direction":"in",
+ "ytick.direction":"in",
+ "xtick.bottom":True,
+ "xtick.top":True,
+ "ytick.left":True,
+ "ytick.right":True,
+ }
+ )
+
+ # align the names with Snowfall and Snowflake
+ ifwhat=="t_nucleation":
+ what="t_nuc"
+ elifwhat=="T_nucleation":
+ what="T_nuc"
+ elifwhat=="t_solidification":
+ what="t_sol"
+
+ # if nucleation temperatures requested change complementary to True
+ ifwhat=="T_nuc":
+ comp=True
+
+ # get data
+ self._statsMultiple_df=self.results
+
+ # plot attributes
+ plt.figure()
+ ifwhat=="T_nuc":
+ plt.xlabel("nucleation temperature, $T^{nuc}$ [°C]")
+ plt.ylabel("nucleation temperature CDF, $F_{nT}$ [K$^{-1}$]")
+ elifwhat=="t_nuc":
+ plt.xlabel("nucleation time, $t^{nuc}$ [min]")
+ plt.ylabel("nucleation time CDF, $F_{nt}$ [min$^{-1}$]")
+ elifwhat=="t_sol":
+ plt.xlabel("solidification time, $t^{sol}$ [min]")
+ plt.ylabel("solidification time CDF, $F_{sol}$ [min$^{-1}$]")
+ elifwhat=="t_fr":
+ plt.xlabel("complete freezing time, $t^{fr}$ [min]")
+ plt.ylabel("frrozen product time CDF, $F_{fr}$ [min$^{-1}$]")
+ else:
+ raiseValueError(
+ 'Property to plot incorrectly specified. Use "T_nuc", "t_nuc", "t_sol" or "t_fr" instead.'
+ )
+
+ # plot different nucleation temperatures for spatial models
+ if(self.const["dimensionality"]!="homogeneous")&(what=="T_nuc"):
+ self._data_df_Tnuc=self._statsMultiple_df[
+ ["T_nuc_min","T_nuc_kin","T_nuc_mean","T_nuc_max"]
+ ]
+ sns.ecdfplot(data=self._data_df_Tnuc,complementary=comp)
+ # in other cases, just plot the desired quantity
+ else:
+ sns.ecdfplot(data=self._statsMultiple_df,x=what,complementary=comp)
+ plt.grid(axis="both",color="lightgray",linestyle="solid")
+ plt.ylim(-0.05,1.05)
+ plt.show()
+
+ # plot the stats (categorical plots) of whichever desired variable
+
[docs]defplot(self,what:str="t_nuc",kind:str="box"):
+ """Create categorical plots for Snowing object.
+
+ Args:
+ what (str, optional): What to plot, i.e., keys of the stats dict.
+ Valid options are t_nuc, T_nuc, t_sol, t_fr.
+ Defaults to "t_nuc". Also accepting t_nucleation, T_nucleation
+ and t_solidification.
+ kind (str, optional): Any sns.catplot 'kind' input is allowed.
+ Defaults to "box".
+ """
+
+ # enable this plot only for Nrep > 1
+ ifself.Nrep<=1:
+ raiseNotImplementedError(
+ "This plot is only available for multiple simulations. "
+ +"Run with Nrep > 1."
+ )
+
+ # plot settings
+ plt.rcParams.update(
+ {
+ "axes.axisbelow":True,
+ "xtick.direction":"in",
+ "ytick.direction":"in",
+ "xtick.bottom":True,
+ "xtick.top":True,
+ "ytick.left":True,
+ "ytick.right":True,
+ }
+ )
+
+ # align the names with Snowfall and Snowflake
+ ifwhat=="t_nucleation":
+ what="t_nuc"
+ elifwhat=="T_nucleation":
+ what="T_nuc"
+ elifwhat=="t_solidification":
+ what="t_sol"
+
+ # get data
+ self._statsMultiple_df=self.results
+
+ # plot different nucleation temperatures for spatial models
+ if(self.const["dimensionality"]!="homogeneous")&(what=="T_nuc"):
+ self._data_df_Tnuc=self._statsMultiple_df[
+ ["T_nuc_min","T_nuc_kin","T_nuc_mean","T_nuc_max"]
+ ]
+ sns.catplot(data=self._data_df_Tnuc,kind=kind)
+ # in other cases, just plot the desired quantity
+ else:
+ sns.catplot(data=self._statsMultiple_df,x=what,kind=kind)
+
+ # plot the temperature evolution
+ def_plot_temperature_evolution(self):
+ # temperature evolution plot
+ plt.figure()
+ # plot nucleation time and equilibrum temperature
+ plt.plot(
+ [self._stats["t_nuc"]/60,self._stats["t_nuc"]/60],
+ [self._shelfTemp.min()-5,self._shelfTemp.max()+5],
+ "--",
+ color="red",
+ linewidth=1,
+ label="$t = t^{nuc}$",
+ )
+ # plot time of frozen product
+ ifself._stats["t_fr"]isnotNone:
+ plt.plot(
+ [self._stats["t_fr"]/60,self._stats["t_fr"]/60],
+ [self._shelfTemp.min()-5,self._shelfTemp.max()+5],
+ "--",
+ color="indigo",
+ linewidth=1,
+ label="$t = t^{fr}$",
+ )
+ plt.plot(
+ [-0.05*self._time.max(),1.05*self._time.max()],
+ [self._prop[0]-273.15,self._prop[0]-273.15],
+ "--",
+ color="magenta",
+ linewidth=1,
+ label="$T = T^{eq}_{l}$",
+ )
+ # plot shelf temperature evolution
+ plt.plot(self._time,self._shelfTemp,"b-",linewidth=2.5,label="shelf T")
+
+ ifself.const["dimensionality"]=="homogeneous":
+ # plot temperature evolution for homogeneous model
+ plt.plot(self._time,self._temp,"k-",linewidth=2.5,label="product T")
+
+ # access legend objects automatically created from data
+ handles,labels=plt.gca().get_legend_handles_labels()
+
+ # plot product temperature evolution for 1D case
+ ifself.const["dimensionality"]=="spatial_1D":
+ # colors
+ colors,my_cmap,norm=Utils.colormap(self._domain)
+ # plot temperature evolution for different vertical positions
+ foriinrange(0,len(self._domain),round(len(self._domain)/11)):
+ plt.plot(self._time,self._temp[:,i],color=colors[i],linewidth=2.5)
+ # colorbar
+ sm=plt.cm.ScalarMappable(cmap=my_cmap,norm=norm)
+ cbar=plt.colorbar(sm,ax=plt.gca(),aspect=35)
+ cbar.set_label("dimensionless vertical position, z/H [-]",rotation=90)
+ # add colormap to the legend entries
+ cmaps_gradients=my_cmap(np.linspace(0,1,30))
+ # combine patches into a list and add to the handles and labels
+ cmap_patches=[
+ patches.Patch(facecolor=c,edgecolor=c)forcincmaps_gradients
+ ]
+ handles.append(cmap_patches)
+ labels.append("product T")
+
+ # plot product temperature evolution for 2D case
+ elifself.const["dimensionality"]=="spatial_2D":
+ # colors
+ colors,my_cmap,norm=Utils.colormap(self._domain[0])
+ # plot temperature evolution for different vertical positions
+ foriinrange(0,len(self._domain[0]),round(len(self._domain[0])/11)):
+ plt.plot(
+ self._time,self._temp[:,i,-1],color=colors[i],linewidth=2.5
+ )
+ # colorbar
+ sm=plt.cm.ScalarMappable(cmap=my_cmap,norm=norm)
+ cbar=plt.colorbar(sm,ax=plt.gca(),aspect=35)
+ cbar.set_label(
+ "dimensionless vertical position (at r = 0), z/H [-]",rotation=90
+ )
+ # add colormap to the legend entries
+ cmaps_gradients=my_cmap(np.linspace(0,1,30))
+ # combine patches into a list and add to the handles and labels
+ cmap_patches=[
+ patches.Patch(facecolor=c,edgecolor=c)forcincmaps_gradients
+ ]
+ handles.append(cmap_patches)
+ labels.append("product T")
+
+ # plot labels and attributes
+ plt.xlabel("process time, $t$ [h]")
+ plt.ylabel("temperature, $T$ [°C]",labelpad=8)
+ plt.grid(axis="both",color="lightgray",linestyle="solid")
+ plt.xlim(-0.05*self._time.max(),1.05*self._time.max())
+ plt.ylim(self._shelfTemp.min()-5,self._shelfTemp.max()+5)
+ plt.title("Product temperature evolution")
+ plt.legend(
+ handles=handles,
+ labels=labels,
+ handler_map={list:HandlerTuple(ndivide=None,pad=0)},
+ loc="best",
+ )
+ plt.show()
+
+ # plot ice mass fraction evolution
+ def_plot_ice_mass_fraction_evolution(self):
+ # ice mass fraction plot
+ plt.figure()
+ # plot lines for nucleation time and maximum ice mass fraction
+ plt.plot(
+ [self._stats["t_nuc"]/60,self._stats["t_nuc"]/60],
+ [-0.05,1.05],
+ "--",
+ color="red",
+ linewidth=1,
+ label="$t = t^{nuc}$",
+ )
+ # plot time of frozen product
+ ifself._stats["t_fr"]isnotNone:
+ plt.plot(
+ [self._stats["t_fr"]/60,self._stats["t_fr"]/60],
+ [-0.05,1.05],
+ "--",
+ color="indigo",
+ linewidth=1,
+ label="$t = t^{fr}$",
+ )
+ plt.plot(
+ [-0.05*self._time.max(),1.05*self._time.max()],
+ [1-self._prop[1],1-self._prop[1]],
+ "--",
+ color="blue",
+ linewidth=1,
+ label="$w_{i} = 1 - w_{s}$",
+ )
+
+ ifself.const["dimensionality"]=="homogeneous":
+ # plot temperature evolution for homogeneous model
+ plt.plot(
+ self._time,
+ self._iceMassFraction,
+ "k-",
+ linewidth=2.5,
+ label="product $w_{i}$",
+ )
+
+ # access legend objects automatically created from data
+ handles,labels=plt.gca().get_legend_handles_labels()
+
+ # plot product temperature evolution for 1D case
+ ifself.const["dimensionality"]=="spatial_1D":
+ # colors
+ colors,my_cmap,norm=Utils.colormap(self._domain)
+ # colorbar
+ sm=plt.cm.ScalarMappable(cmap=my_cmap,norm=norm)
+ cbar=plt.colorbar(sm,ax=plt.gca(),aspect=35)
+ cbar.set_label("dimensionless vertical position, z/H [-]",rotation=90)
+ # plot temperature evolution for different vertical positions
+ foriinrange(0,len(self._domain),round(len(self._domain)/11)):
+ plt.plot(
+ self._time,
+ self._iceMassFraction[:,i],
+ color=colors[i],
+ linewidth=2.5,
+ )
+ # add colormap to the legend entries
+ cmaps_gradients=my_cmap(np.linspace(0,1,30))
+ # combine patches into a list and add to the handles and labels
+ cmap_patches=[
+ patches.Patch(facecolor=c,edgecolor=c)forcincmaps_gradients
+ ]
+ handles.append(cmap_patches)
+ labels.append("product $w_{i}$")
+
+ # plot product temperature evolution for 2D case
+ elifself.const["dimensionality"]=="spatial_2D":
+ # colors
+ colors,my_cmap,norm=Utils.colormap(self._domain[0])
+ # colorbar
+ sm=plt.cm.ScalarMappable(cmap=my_cmap,norm=norm)
+ cbar=plt.colorbar(sm,ax=plt.gca(),aspect=35)
+ cbar.set_label(
+ "dimensionless vertical position (at r = 0), z/H [-]",rotation=90
+ )
+ # plot temperature evolution for different vertical positions
+ foriinrange(0,len(self._domain[0]),round(len(self._domain[0])/11)):
+ plt.plot(
+ self._time,
+ self._iceMassFraction[:,i,-1],
+ color=colors[i],
+ linewidth=2.5,
+ )
+ # add colormap to the legend entries
+ cmaps_gradients=my_cmap(np.linspace(0,1,30))
+ # combine patches into a list and add to the handles and labels
+ cmap_patches=[
+ patches.Patch(facecolor=c,edgecolor=c)forcincmaps_gradients
+ ]
+ handles.append(cmap_patches)
+ labels.append("product $w_{i}$")
+
+ # plot labels and attributes
+ plt.xlabel("process time, $t$ [h]")
+ plt.ylabel("ice mass fraction, $w_i$ [-]")
+ plt.grid(axis="both",color="lightgray",linestyle="solid")
+ plt.xlim(-0.05*self._time.max(),1.05*self._time.max())
+ plt.ylim(-0.05,1.05)
+ plt.title("Ice mass fraction evolution")
+ plt.legend(
+ handles=handles,
+ labels=labels,
+ handler_map={list:HandlerTuple(ndivide=None,pad=0)},
+ loc="best",
+ )
+ plt.show()
+
+ # plot evolution function
+
[docs]defplot_evolution(self,what:str="temperature"):
+ """Function to plot the spatial evolution of temperature or ice mass fractions.
+
+ Args:
+ what (str, optional): Quantity to be plotted. Defaults to "temperature".
+
+ Raises:
+ NotImplementedError: This plotting functionality is only implemented for single simulations (when Nrep = 1), otherwise, an error is raised.
+ ValueError: Raises error if quantity for plotting is incorrectly specified.
+ """
+
+ # plot settings
+ plt.rcParams.update(
+ {
+ "axes.axisbelow":True,
+ "xtick.direction":"in",
+ "ytick.direction":"in",
+ "xtick.bottom":True,
+ "xtick.top":True,
+ "ytick.left":True,
+ "ytick.right":True,
+ }
+ )
+
+ ifself.Nrep>1:
+ raiseNotImplementedError(
+ "This plot only available for single simulations. "
+ +"Run with Nrep = 1."
+ )
+
+ ifwhat=="temperature":
+ self._plot_temperature_evolution()
+ elifwhat=="ice_mass_fraction":
+ self._plot_ice_mass_fraction_evolution()
+ else:
+ raiseValueError(
+ 'Property to plot incorrectly specified. Use "temperature" or "ice_mass_fraction" instead.'
+ )
+
+ # repr function
+ def__repr__(self)->str:
+ """Return string representation of the Snowing class.
+
+ Returns:
+ str: The Snowing class string representation giving some basic info.
+ """
+ return(
+ f"Snowing([{self.Nrep} Snowing{'s'ifself.Nrep>1else''}, "
+ +f"dimensionality: {self.const['dimensionality']}, "
+ +f"configuration: {self.const['configuration']}])"
+ )
+importnumpyasnp
+frommatplotlib.colorsimportLinearSegmentedColormap,BoundaryNorm
+importmatplotlib.colorsasmc
+importmatplotlib.pyplotasplt
+importcolorsys
+
+# -------------------------------------------------------------------------- #
+# Vacuum-induced surface freezing helper functions
+# -------------------------------------------------------------------------- #
+
+# More information on the following consititutive equations used when
+# simulating VISF can be found in the following literature:
+#
+# Marek, R.; Straub, J. Analysis of the evaporation coefficient and
+# the condensation coefficient of water. International Journal of Heat
+# and Mass Transfer 2001, 44, 39–53.
+# https://doi.org/10.1016/S0017-9310(00)00086-7
+#
+# Murphy, D. M.; Koop, T. Review of the vapour pressures of ice and
+# supercooled water for atmospheric applications. Quarterly Journal of
+# the Royal Meteorological Society 2005, 131, 1539–1565.
+# https://doi.org/10.1256/qj.04.94
+
+
+# vapour pressure estimation (temperature in K, pressure in Pa)
+
[docs]defvapour_pressure_liquid(T_liq):
+ """A function to compute the water vapour pressure above
+ the liquid product in a vial based on the product surface
+ temperature.
+
+ Args:
+ T_liq (np.ndarray): Liquid product temperature at the surface.
+
+ Returns:
+ np.ndarray: Water vapour pressure above the product surface.
+ """
+ # pressure in Pa
+ p_liq=np.exp(
+ 54.842763
+ -6763.22/T_liq
+ -4.210*np.log(T_liq)
+ +0.000367*T_liq
+ +np.tanh(0.0415*(T_liq-218.8))
+ *(53.878-1331.22/T_liq-9.44523*np.log(T_liq)+0.014025*T_liq)
+ )
+ # return vapour pressure
+ returnp_liq
+
+
+# vapour pressure estimation (temperature in K, pressure in Pa)
+
[docs]defvapour_pressure_solid(T_sol):
+ """A function to compute the water vapour pressure above
+ the frozen product in a vial based on the frozen product
+ surface temperature.
+
+ Args:
+ T_sol (np.ndarray): Frozen product temperature at the surface.
+
+ Returns:
+ np.ndarray: Water vapour pressure above the product surface.
+ """
+ # pressure in Pa
+ p_sol=np.exp(
+ 9.550426-5723.265/T_sol+3.53068*np.log(T_sol)-0.00728332*T_sol
+ )
+ # return vapour pressure
+ returnp_sol
+
+
+# vapour pressure estimation (temperature in K, pressure in Pa)
+
[docs]defvapour_flux(kappa,m_water,k_B,p_vac,p_vap,T_l,T_v):
+ """A function to compute the water vapour flux at the top surface,
+ which is exposed to vacuum during VISF.
+
+ Args:
+ kappa (float): Evaporation efficiency.
+ m_water (float): Mass of a water molecule.
+ k_B (float): Boltzmann constant.
+ p_vac (np.ndarray): Chamber vacuum pressure.
+ p_vap (np.ndarray): Vapour pressure.
+ T_l (np.ndarray): Product temperature.
+ T_v (np.ndarray): Vapour temperature.
+
+ Returns:
+ np.ndarray: Flux of water vapour at the top surface of the frozen
+ or liquid product
+ """
+ # calculate flux
+ N_w=(
+ (2/(2-kappa))
+ *np.sqrt(m_water*kappa**2/(2*np.pi*k_B))
+ *(p_vap/np.sqrt(T_l)-p_vac/np.sqrt(T_v))
+ )
+ # return vapour flux
+ returnN_w
[docs]defcolormap(z):
+ """A function to create a colormap used to plot the time evolutions of
+ temperature and ice mass fractions with respect to vertical positions
+ in the vial.
+
+ Args:
+ z (np.ndarray): Discretized spatial domain in the vertical direction.
+ """
+
+ deflighten_color(color,amount=0.5):
+ try:
+ c=mc.cnames[color]
+ except:
+ c=color
+ c=colorsys.rgb_to_hls(*mc.to_rgb(c))
+ returncolorsys.hls_to_rgb(c[0],1-amount*(1-c[1]),c[2])
+
+ my_cmap=LinearSegmentedColormap.from_list(
+ "my_BBG",["black",lighten_color("lightslategray",0.75)]
+ )
+ Nz=len(z)
+ colors=my_cmap(np.linspace(0,1,Nz))
+ # number of lines to plot the evolutions
+ N_lines=11
+ colors=my_cmap(np.linspace(0,1,Nz))
+ cmaplist=[my_cmap(i)foriinrange(my_cmap.N)]
+ my_cmap=LinearSegmentedColormap.from_list("Custom cmap",cmaplist,my_cmap.N)
+ bounds=np.linspace(0,1,N_lines)
+ norm=BoundaryNorm(bounds,my_cmap.N)
+
+ returncolors,my_cmap,norm
+
diff --git a/docs/_sources/api/ethz_snow.rst.txt b/docs/_sources/api/ethz_snow.rst.txt
index 82333fe..8516aa3 100644
--- a/docs/_sources/api/ethz_snow.rst.txt
+++ b/docs/_sources/api/ethz_snow.rst.txt
@@ -36,6 +36,22 @@ ethz\_snow.snowflake module
:undoc-members:
:show-inheritance:
+ethz\_snow.snowing module
+-------------------------
+
+.. automodule:: ethz_snow.snowing
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
+ethz\_snow.utils module
+-----------------------
+
+.. automodule:: ethz_snow.utils
+ :members:
+ :undoc-members:
+ :show-inheritance:
+
Module contents
---------------
diff --git a/docs/_sources/authors.rst.txt b/docs/_sources/authors.rst.txt
index 5763ccc..31ca2ee 100644
--- a/docs/_sources/authors.rst.txt
+++ b/docs/_sources/authors.rst.txt
@@ -3,4 +3,5 @@ Contributors
============
* Leif-Thore Deck
-* Dr. David R. Ochsenbein
+* Andraž Košir
+* Dr. David R. Ochsenbein
diff --git a/docs/_sources/changelog.rst.txt b/docs/_sources/changelog.rst.txt
index f57a306..a096571 100644
--- a/docs/_sources/changelog.rst.txt
+++ b/docs/_sources/changelog.rst.txt
@@ -2,6 +2,19 @@
Changelog
=========
+Version 2.0
+===========
+
+- New implementation of freezing model for a single container with spatial resolution (termed snowing). The model considers the gradients in temperature and in ice mass that form during freezing.
+- Currently supported features of the new spatial model:
+ - Simulation of three different freezing process configuration:
+ - Shelf-ramped freezing (commonly used in freeze-drying)
+ - Vacuum-induced surface freezing (process variation where vacuum is applied to promote freezing through surface evaporation)
+ - Jacket-ramped freezing (process variation where vial is surrounded by a temperature-controlled jacket
+ - Simulation of different heat transfer boundary conditions (natural convection, conduction, thermal radiation, surface evaporation)
+ - Simulation with different dimensionalities (0D, 1D, 2D)
+
+
Version 1.1
===========
@@ -24,4 +37,4 @@ Version 1.0
- Arbitrary cooling protocols (i.e., user may choose cooling rate, integrate holding steps and controlled nucleation)
- Tracking of nucleation times, nucleation temperatures and solidification times for all vials
- Stochastic nucleation in the form of a Monte Carlo approach as well as controlled nucleation in the form of forced initiation of nucleation at a certain point
- - Cubic geometry of vial and rectangular arrangement on the shelf
\ No newline at end of file
+ - Cubic geometry of vial and rectangular arrangement on the shelf
diff --git a/docs/_sources/tutorial.rst.txt b/docs/_sources/tutorial.rst.txt
index 220e2b4..bc27e75 100644
--- a/docs/_sources/tutorial.rst.txt
+++ b/docs/_sources/tutorial.rst.txt
@@ -16,11 +16,11 @@ Or install the latest release version via ``pip install ethz_snow``.
First steps: General information
================================
-The package contains several files, two of which enable the simulation of freezing processes. **snowflake.py** is the basis version of the model, capable of simulating the freezing process of a shelf comprising an arbitrary number of vials exactly once. **snowfall.py** is setup to run such simulation repetitively, namely for **Nrep** times. The idea behind the naming, of course, is that Snowfall in reality comprises a large number of Snowflakes, all of which are unique. In the same way, every freezing process will be unique as a consequence of the stochasticity of primary ice nucleation.
+The package contains several files, three of which enable the simulation of freezing processes. **snowflake.py** is the basis version of the model, capable of simulating the freezing process of a shelf comprising an arbitrary number of vials exactly once. **snowfall.py** is setup to run such simulation repetitively, namely for **Nrep** times. The idea behind the naming, of course, is that Snowfall in reality comprises a large number of Snowflakes, all of which are unique. In the same way, every freezing process will be unique as a consequence of the stochasticity of primary ice nucleation. **snowing.py** (snow-internal gradients) provides stochastic simulations of freezing with spatial resolution within a single vial; it has been integrated into the package in version 2.0.
-Both models require a number of operatingConditions and constants as input parameters, which may be adjusted by the user. We define operating conditions as those input parameters, that are related to the numerics of the model (e.g. time step, number of repetitions) or to technical aspects (e.g. heat transfer parameters, cooling protocol). Constants, on the other hand, refer to formulation specific parameters, such as the volume of a vial, its composition and the required individual physicochemical parameters of the individual components (including ice nucleation kinetics).
+All models require a number of operatingConditions and constants as input parameters, which may be adjusted by the user. We define operating conditions as those input parameters, that are related to the numerics of the model (e.g. time step, number of repetitions) or to technical aspects (e.g. heat transfer parameters, cooling protocol). Constants, on the other hand, refer to formulation specific parameters, such as the volume of a vial, its composition and the required individual physicochemical parameters of the individual components (including ice nucleation kinetics).
-Operating conditions are to be provided directly when calling the **snowflake** / **Snowfall** functions, while we specify the constants separately in a .yaml file. If no values are specified, the default configurations are used; these are stored in the **operatingConditions.py** and the **snowConfig_default.yaml** (see the next subsection for more details on the latter).
+Operating conditions are to be provided directly when calling the **snowflake** / **snowfall** / **snowing** functions, while we specify the constants separately in a .yaml file. If no values are specified, the default configurations are used; these are stored in the **operatingConditions.py** and the **snowConfig_default.yaml** (see the next subsection for more details on the latter).
Constants
---------
@@ -116,7 +116,7 @@ Finally, we may define the Snowfall class. We set the pool_size parameter to the
S = Snowfall(pool_size=8, k=d, Nrep=50, N_vials=(7,7,1), opcond=op)
-We then start the simulation via **S.run()** and may check whether it completed via **S.simulationStatus()**. In case we are only interested in a single repetition, the **Snowflake** class may be used instead. Compared to **Snowfall**, **Snowflake** does not require Nrep or pool_size as input. However, it is able to store information on the thermal evolution of all vials, which is a feature that was removed for **Snowfall** to increase computational performance.
+We then start the simulation via **S.run()** and may check whether it completed via **S.simulationStatus**. In case we are only interested in a single repetition, the **Snowflake** class may be used instead. Compared to **Snowfall**, **Snowflake** does not require Nrep or pool_size as input. However, it is able to store information on the thermal evolution of all vials, which is a feature that was removed for **Snowfall** to increase computational performance.
Simulation output
=================
@@ -151,4 +151,151 @@ Note that due to the geometry applied, the heat transfer settings for a pallet c
initial = {"temp": 20}
op = OperatingConditions(t_tot=6e6,cooling= c, holding =dict(duration=6e6,temp=-8) )
-In order to simulate a constant storage temperature, an arbitrarily small cooling rate may be defined in addition with a holding step. In this way, the temperature is set for the entire process duration to the storage temperature, which is -8°C in this case. Note that due to the large system size, typically longer process durations have to be simulated for pallets compared to shelf freezing.
\ No newline at end of file
+In order to simulate a constant storage temperature, an arbitrarily small cooling rate may be defined in addition with a holding step. In this way, the temperature is set for the entire process duration to the storage temperature, which is -8°C in this case. Note that due to the large system size, typically longer process durations have to be simulated for pallets compared to shelf freezing.
+
+Version 2.0. Freezing simulations with internal gradients
+=========================================================
+
+In version 2.0. new functionalities related to spatial phenomena during freezing are integrated into the package: when simulating the freezing process in a single container of arbitrary size, the model considers gradients of temperature and of ice mass fraction within the container. Such spatial simulation of freezing is currently only available for a single container, i.e., it does not consider thermal interactions with potential neighboring containers. Hence the spatial freezing model in version 2.0. provides complimentary information to the process-scale freezing models introduced earlier. An upcoming scientific publication will discuss the model development, validation, and relevant use cases in detail.
+
+The spatial model (termed Snowing) accounts for different **dimensionalities (0D, 1D and 2D)** of the vial and for different **freezing configurations (shelf-ramped freezing, vacuum-induced surface freezing (VISF) and jacket-ramped freezing)**. These three configurations represent common freezing conditions employed in both commercial manufacturing and in academia. Both model dimensionality and freezing configuration are specified in the yaml file used to specify the constants. To this end, we add additional blocks of parameters used for the spatial simulation to the yaml. An example of the updated yaml file is provided below:
+
+.. code-block:: yaml
+
+ # model dimensionality/temperature resolution within vial (homogeneous, spatial_1D, spatial_2D)
+ dimensionality: spatial_1D
+
+ # freezing configuration to be simulated (shelf, VISF, jacket)
+ configuration: shelf
+
+ # all vial-related parameters
+ vial:
+ # define the geometry
+ geometry:
+ # base shape of vial (snow and snowfall only accepts cube, snowing rewrites it to cylinder)
+ shape: cube
+ # diameter of the vial for cylindrical geometry (only spatial model) [m]
+ diameter: 0.01
+ # height of the filled product for cubic (snow/snowfall) and cylindrical geometry (only spatial model) [m]
+ height: 0.01
+ # length of the vial when shape is cube [m]
+ length: 0.01
+ # width of the vial when shape is cube [m]
+ width: 0.01
+
+ # parameters used in VISF simulation
+ VISF:
+ # vacuum pressure [Pa]
+ p_vac: 100
+ # evaporation coefficient [-]
+ kappa: 0.01
+ # latent heat of vaporization for water [J/kg]
+ Dh_evaporation: 2500.9e3
+ # mass of water molecule [kg]
+ m_water: 2.99e-26
+ # time for vacuum start [h]
+ t_vac_start: 0.75
+ # duration of the vacuum [h]
+ t_vac_duration: 0.1
+
+ # parameters for jacket-ramped freezing
+ jacket:
+ # air gap between the wall and the vial [m]
+ air_gap: 0.001
+ # heat conductivity of air [W/mK]
+ lambda_air: 0.025
+
+ # all parameters related to water
+ water:
+ # specific heat capacity of liquid water [J/kgK]
+ cp_w: 4187
+ # specific heat capacity of ice [J/kgK]
+ cp_i: 2108
+ # heat conductivity of liquid water [W/mK]
+ lambda_w: 0.598
+ # heat conductivity of ice [W/mK]
+ lambda_i: 2.25
+ # latent heat of fusion of water [J/kg]
+ Dh: 333550
+
+ # all solution parameters
+ solution:
+ # solute mass fraction [-]
+ solid_fraction: 0.05
+ # melting temperature of pure water [°C]
+ T_eq: 0
+ # density of liquid water [kg/m3]
+ rho_l: 1000
+ # specific heat capacity of solute [J/kgK]
+ cp_s: 1240
+ # heat conductivity of solute [W/mK]
+ lambda_s: 0.126
+ # cryoscopic constant of water [K/kgmol]
+ k_f: 1.853
+ # molar mass of solute [kg/mol]
+ M_s: 0.3423
+
+ # nucleation kinetics
+ kinetics:
+ # pre-exponential nucleation parameter [1/m3sKb]
+ kb: 1e-9
+ # exponential nucleation parameter [-]
+ b: 12
+
+ general:
+ # Stefan-Boltzmann constant [W/m2]
+ sigma_B: 5.67e-8
+ # Boltzmann constant [JK]
+ k_B: 1.38e-23
+
+Vial geometry, other constants and operating conditions are specified as established in previous versions. We start by importing the new **Snowing** module:
+
+.. code-block:: python
+
+ from ethz_snow.snowing import Snowing
+
+Additionally, we also need to import the operatingConditions, define the heat transfer parameters in a dictionary and constants in a Yaml file linked to Snowing via configPath (same as in previous versions, see above). A sample configuration of the spatial model may in this case be initiated by creating an instance of the **Snowing** class:
+
+.. code-block:: python
+
+ S = Snowing(k=d, opcond=op)
+
+The simulation is then carried out using: ``S.run()`` In this case, simulation is carried out in 1D (considering heat transfer only in the vertical direction) with the freezing configuration being set to the conventional shelf-ramped freezing. The two parameters ``k=d`` and ``opcond=op`` (considering heat transfer and operating conditions are identical to the ones used for Snowfall and Snowflake). Different dimensionalities of the model can be called by varying the **dimensionality** parameter specified in the yaml file, possible keys are: **homogeneous**, **spatial_1D** or **spatial_2D**, while different freezing configurations can be simulated by modifying the **configuration** parameter, options are: **shelf**, **VISF**, **jacket**.
+
+In order to evaluate the variability in nucleation times, temperatures and solidification times due to the stochasticity of nucleation a larger number of the single vial simulations may be carried out. This can be achieved by adding an integer parameter **Nrep** denoting the number of repeated independent simulations:
+
+.. code-block:: python
+
+ S = Snowing(k=d, opcond=op, Nrep = 100)
+
+When **Nrep > 1**, the user can plot the statistics in a form of a cumulative probability function of a desired variable using:
+
+.. code-block:: python
+
+ S.plot_cdf(what = "T_nuc")
+ S.plot_cdf(what = "t_nuc")
+ S.plot_cdf(what = "t_sol")
+ S.plot_cdf(what = "t_fr")
+
+Besides nucleation times (t_nuc), the user can plot the cumulative distribution functions of nucleation temperatures (T_nuc), solidification times (t_sol) and times of complete freezing (t_fr). In case of 1D or 2D model complexity, temperature at the time of nucleation is a field, hence the choice of nucleation temperature is not straightforward. To this end, ``S.plot_cdf(what = "T_nuc")`` plots distributions of 4 different temperatures: minimum, kinetic mean, mean and maximum temperature at nucleation. For more information see the relevant publication regarding the spatial model. Finally, the following command allows the user to get the statistics on all the relevant variables (output is a DataFrame):
+
+.. code-block:: python
+
+ S.results
+
+In case of a single simulation (**Nrep = 1**), the following commands also provide detailed simulation results (time array, shelf temperature profile, temperature and ice mass fraction field evolution):
+
+.. code-block:: python
+
+ time = S.time
+ shelf = S.shelfTemp
+ temp = S.temp
+ ice = S.iceMassFraction
+
+Additionally, the time evolution of product temperature and ice mass fraction may be plotted using:
+
+.. code-block:: python
+
+ S.plot_evolution(what = "temperature")
+ S.plot_evolution(what = "ice_mass_fraction")
+
diff --git a/docs/_static/_sphinx_javascript_frameworks_compat.js b/docs/_static/_sphinx_javascript_frameworks_compat.js
new file mode 100644
index 0000000..8549469
--- /dev/null
+++ b/docs/_static/_sphinx_javascript_frameworks_compat.js
@@ -0,0 +1,134 @@
+/*
+ * _sphinx_javascript_frameworks_compat.js
+ * ~~~~~~~~~~
+ *
+ * Compatability shim for jQuery and underscores.js.
+ *
+ * WILL BE REMOVED IN Sphinx 6.0
+ * xref RemovedInSphinx60Warning
+ *
+ */
+
+/**
+ * select a different prefix for underscore
+ */
+$u = _.noConflict();
+
+
+/**
+ * small helper function to urldecode strings
+ *
+ * See https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/decodeURIComponent#Decoding_query_parameters_from_a_URL
+ */
+jQuery.urldecode = function(x) {
+ if (!x) {
+ return x
+ }
+ return decodeURIComponent(x.replace(/\+/g, ' '));
+};
+
+/**
+ * small helper function to urlencode strings
+ */
+jQuery.urlencode = encodeURIComponent;
+
+/**
+ * This function returns the parsed url parameters of the
+ * current request. Multiple values per key are supported,
+ * it will always return arrays of strings for the value parts.
+ */
+jQuery.getQueryParameters = function(s) {
+ if (typeof s === 'undefined')
+ s = document.location.search;
+ var parts = s.substr(s.indexOf('?') + 1).split('&');
+ var result = {};
+ for (var i = 0; i < parts.length; i++) {
+ var tmp = parts[i].split('=', 2);
+ var key = jQuery.urldecode(tmp[0]);
+ var value = jQuery.urldecode(tmp[1]);
+ if (key in result)
+ result[key].push(value);
+ else
+ result[key] = [value];
+ }
+ return result;
+};
+
+/**
+ * highlight a given string on a jquery object by wrapping it in
+ * span elements with the given class name.
+ */
+jQuery.fn.highlightText = function(text, className) {
+ function highlight(node, addItems) {
+ if (node.nodeType === 3) {
+ var val = node.nodeValue;
+ var pos = val.toLowerCase().indexOf(text);
+ if (pos >= 0 &&
+ !jQuery(node.parentNode).hasClass(className) &&
+ !jQuery(node.parentNode).hasClass("nohighlight")) {
+ var span;
+ var isInSVG = jQuery(node).closest("body, svg, foreignObject").is("svg");
+ if (isInSVG) {
+ span = document.createElementNS("http://www.w3.org/2000/svg", "tspan");
+ } else {
+ span = document.createElement("span");
+ span.className = className;
+ }
+ span.appendChild(document.createTextNode(val.substr(pos, text.length)));
+ node.parentNode.insertBefore(span, node.parentNode.insertBefore(
+ document.createTextNode(val.substr(pos + text.length)),
+ node.nextSibling));
+ node.nodeValue = val.substr(0, pos);
+ if (isInSVG) {
+ var rect = document.createElementNS("http://www.w3.org/2000/svg", "rect");
+ var bbox = node.parentElement.getBBox();
+ rect.x.baseVal.value = bbox.x;
+ rect.y.baseVal.value = bbox.y;
+ rect.width.baseVal.value = bbox.width;
+ rect.height.baseVal.value = bbox.height;
+ rect.setAttribute('class', className);
+ addItems.push({
+ "parent": node.parentNode,
+ "target": rect});
+ }
+ }
+ }
+ else if (!jQuery(node).is("button, select, textarea")) {
+ jQuery.each(node.childNodes, function() {
+ highlight(this, addItems);
+ });
+ }
+ }
+ var addItems = [];
+ var result = this.each(function() {
+ highlight(this, addItems);
+ });
+ for (var i = 0; i < addItems.length; ++i) {
+ jQuery(addItems[i].parent).before(addItems[i].target);
+ }
+ return result;
+};
+
+/*
+ * backward compatibility for jQuery.browser
+ * This will be supported until firefox bug is fixed.
+ */
+if (!jQuery.browser) {
+ jQuery.uaMatch = function(ua) {
+ ua = ua.toLowerCase();
+
+ var match = /(chrome)[ \/]([\w.]+)/.exec(ua) ||
+ /(webkit)[ \/]([\w.]+)/.exec(ua) ||
+ /(opera)(?:.*version|)[ \/]([\w.]+)/.exec(ua) ||
+ /(msie) ([\w.]+)/.exec(ua) ||
+ ua.indexOf("compatible") < 0 && /(mozilla)(?:.*? rv:([\w.]+)|)/.exec(ua) ||
+ [];
+
+ return {
+ browser: match[ 1 ] || "",
+ version: match[ 2 ] || "0"
+ };
+ };
+ jQuery.browser = {};
+ jQuery.browser[jQuery.uaMatch(navigator.userAgent).browser] = true;
+}
diff --git a/docs/_static/basic.css b/docs/_static/basic.css
index 912859b..7d5974c 100644
--- a/docs/_static/basic.css
+++ b/docs/_static/basic.css
@@ -4,7 +4,7 @@
*
* Sphinx stylesheet -- basic theme.
*
- * :copyright: Copyright 2007-2021 by the Sphinx team, see AUTHORS.
+ * :copyright: Copyright 2007-2022 by the Sphinx team, see AUTHORS.
* :license: BSD, see LICENSE for details.
*
*/
@@ -222,7 +222,7 @@ table.modindextable td {
/* -- general body styles --------------------------------------------------- */
div.body {
- min-width: 450px;
+ min-width: 360px;
max-width: 800px;
}
@@ -236,7 +236,6 @@ div.body p, div.body dd, div.body li, div.body blockquote {
a.headerlink {
visibility: hidden;
}
-
a.brackets:before,
span.brackets > a:before{
content: "[";
@@ -247,6 +246,7 @@ span.brackets > a:after {
content: "]";
}
+
h1:hover > a.headerlink,
h2:hover > a.headerlink,
h3:hover > a.headerlink,
@@ -334,13 +334,11 @@ aside.sidebar {
p.sidebar-title {
font-weight: bold;
}
-
div.admonition, div.topic, blockquote {
clear: left;
}
/* -- topics ---------------------------------------------------------------- */
-
div.topic {
border: 1px solid #ccc;
padding: 7px;
@@ -428,10 +426,6 @@ table.docutils td, table.docutils th {
border-bottom: 1px solid #aaa;
}
-table.footnote td, table.footnote th {
- border: 0 !important;
-}
-
th {
text-align: left;
padding-right: 5px;
@@ -615,6 +609,7 @@ ul.simple p {
margin-bottom: 0;
}
+/* Docutils 0.17 and older (footnotes & citations) */
dl.footnote > dt,
dl.citation > dt {
float: left;
@@ -632,6 +627,33 @@ dl.citation > dd:after {
clear: both;
}
+/* Docutils 0.18+ (footnotes & citations) */
+aside.footnote > span,
+div.citation > span {
+ float: left;
+}
+aside.footnote > span:last-of-type,
+div.citation > span:last-of-type {
+ padding-right: 0.5em;
+}
+aside.footnote > p {
+ margin-left: 2em;
+}
+div.citation > p {
+ margin-left: 4em;
+}
+aside.footnote > p:last-of-type,
+div.citation > p:last-of-type {
+ margin-bottom: 0em;
+}
+aside.footnote > p:last-of-type:after,
+div.citation > p:last-of-type:after {
+ content: "";
+ clear: both;
+}
+
+/* Footnotes & citations ends */
+
dl.field-list {
display: grid;
grid-template-columns: fit-content(30%) auto;
@@ -731,8 +753,9 @@ dl.glossary dt {
.classifier:before {
font-style: normal;
- margin: 0.5em;
+ margin: 0 0.5em;
content: ":";
+ display: inline-block;
}
abbr, acronym {
@@ -756,6 +779,7 @@ span.pre {
-ms-hyphens: none;
-webkit-hyphens: none;
hyphens: none;
+ white-space: nowrap;
}
div[class*="highlight-"] {
diff --git a/docs/_static/doctools.js b/docs/_static/doctools.js
index 8cbf1b1..c3db08d 100644
--- a/docs/_static/doctools.js
+++ b/docs/_static/doctools.js
@@ -2,322 +2,263 @@
* doctools.js
* ~~~~~~~~~~~
*
- * Sphinx JavaScript utilities for all documentation.
+ * Base JavaScript utilities for all Sphinx HTML documentation.
*
- * :copyright: Copyright 2007-2021 by the Sphinx team, see AUTHORS.
+ * :copyright: Copyright 2007-2022 by the Sphinx team, see AUTHORS.
* :license: BSD, see LICENSE for details.
*
*/
+"use strict";
-/**
- * select a different prefix for underscore
- */
-$u = _.noConflict();
-
-/**
- * make the code below compatible with browsers without
- * an installed firebug like debugger
-if (!window.console || !console.firebug) {
- var names = ["log", "debug", "info", "warn", "error", "assert", "dir",
- "dirxml", "group", "groupEnd", "time", "timeEnd", "count", "trace",
- "profile", "profileEnd"];
- window.console = {};
- for (var i = 0; i < names.length; ++i)
- window.console[names[i]] = function() {};
-}
- */
-
-/**
- * small helper function to urldecode strings
- *
- * See https://developer.mozilla.org/en-US/docs/Web/JavaScript/Reference/Global_Objects/decodeURIComponent#Decoding_query_parameters_from_a_URL
- */
-jQuery.urldecode = function(x) {
- if (!x) {
- return x
+const _ready = (callback) => {
+ if (document.readyState !== "loading") {
+ callback();
+ } else {
+ document.addEventListener("DOMContentLoaded", callback);
}
- return decodeURIComponent(x.replace(/\+/g, ' '));
};
/**
- * small helper function to urlencode strings
+ * highlight a given string on a node by wrapping it in
+ * span elements with the given class name.
*/
-jQuery.urlencode = encodeURIComponent;
+const _highlight = (node, addItems, text, className) => {
+ if (node.nodeType === Node.TEXT_NODE) {
+ const val = node.nodeValue;
+ const parent = node.parentNode;
+ const pos = val.toLowerCase().indexOf(text);
+ if (
+ pos >= 0 &&
+ !parent.classList.contains(className) &&
+ !parent.classList.contains("nohighlight")
+ ) {
+ let span;
-/**
- * This function returns the parsed url parameters of the
- * current request. Multiple values per key are supported,
- * it will always return arrays of strings for the value parts.
- */
-jQuery.getQueryParameters = function(s) {
- if (typeof s === 'undefined')
- s = document.location.search;
- var parts = s.substr(s.indexOf('?') + 1).split('&');
- var result = {};
- for (var i = 0; i < parts.length; i++) {
- var tmp = parts[i].split('=', 2);
- var key = jQuery.urldecode(tmp[0]);
- var value = jQuery.urldecode(tmp[1]);
- if (key in result)
- result[key].push(value);
- else
- result[key] = [value];
- }
- return result;
-};
+ const closestNode = parent.closest("body, svg, foreignObject");
+ const isInSVG = closestNode && closestNode.matches("svg");
+ if (isInSVG) {
+ span = document.createElementNS("http://www.w3.org/2000/svg", "tspan");
+ } else {
+ span = document.createElement("span");
+ span.classList.add(className);
+ }
-/**
- * highlight a given string on a jquery object by wrapping it in
- * span elements with the given class name.
- */
-jQuery.fn.highlightText = function(text, className) {
- function highlight(node, addItems) {
- if (node.nodeType === 3) {
- var val = node.nodeValue;
- var pos = val.toLowerCase().indexOf(text);
- if (pos >= 0 &&
- !jQuery(node.parentNode).hasClass(className) &&
- !jQuery(node.parentNode).hasClass("nohighlight")) {
- var span;
- var isInSVG = jQuery(node).closest("body, svg, foreignObject").is("svg");
- if (isInSVG) {
- span = document.createElementNS("http://www.w3.org/2000/svg", "tspan");
- } else {
- span = document.createElement("span");
- span.className = className;
- }
- span.appendChild(document.createTextNode(val.substr(pos, text.length)));
- node.parentNode.insertBefore(span, node.parentNode.insertBefore(
+ span.appendChild(document.createTextNode(val.substr(pos, text.length)));
+ parent.insertBefore(
+ span,
+ parent.insertBefore(
document.createTextNode(val.substr(pos + text.length)),
- node.nextSibling));
- node.nodeValue = val.substr(0, pos);
- if (isInSVG) {
- var rect = document.createElementNS("http://www.w3.org/2000/svg", "rect");
- var bbox = node.parentElement.getBBox();
- rect.x.baseVal.value = bbox.x;
- rect.y.baseVal.value = bbox.y;
- rect.width.baseVal.value = bbox.width;
- rect.height.baseVal.value = bbox.height;
- rect.setAttribute('class', className);
- addItems.push({
- "parent": node.parentNode,
- "target": rect});
- }
+ node.nextSibling
+ )
+ );
+ node.nodeValue = val.substr(0, pos);
+
+ if (isInSVG) {
+ const rect = document.createElementNS(
+ "http://www.w3.org/2000/svg",
+ "rect"
+ );
+ const bbox = parent.getBBox();
+ rect.x.baseVal.value = bbox.x;
+ rect.y.baseVal.value = bbox.y;
+ rect.width.baseVal.value = bbox.width;
+ rect.height.baseVal.value = bbox.height;
+ rect.setAttribute("class", className);
+ addItems.push({ parent: parent, target: rect });
}
}
- else if (!jQuery(node).is("button, select, textarea")) {
- jQuery.each(node.childNodes, function() {
- highlight(this, addItems);
- });
- }
+ } else if (node.matches && !node.matches("button, select, textarea")) {
+ node.childNodes.forEach((el) => _highlight(el, addItems, text, className));
}
- var addItems = [];
- var result = this.each(function() {
- highlight(this, addItems);
- });
- for (var i = 0; i < addItems.length; ++i) {
- jQuery(addItems[i].parent).before(addItems[i].target);
- }
- return result;
};
-
-/*
- * backward compatibility for jQuery.browser
- * This will be supported until firefox bug is fixed.
- */
-if (!jQuery.browser) {
- jQuery.uaMatch = function(ua) {
- ua = ua.toLowerCase();
-
- var match = /(chrome)[ \/]([\w.]+)/.exec(ua) ||
- /(webkit)[ \/]([\w.]+)/.exec(ua) ||
- /(opera)(?:.*version|)[ \/]([\w.]+)/.exec(ua) ||
- /(msie) ([\w.]+)/.exec(ua) ||
- ua.indexOf("compatible") < 0 && /(mozilla)(?:.*? rv:([\w.]+)|)/.exec(ua) ||
- [];
-
- return {
- browser: match[ 1 ] || "",
- version: match[ 2 ] || "0"
- };
- };
- jQuery.browser = {};
- jQuery.browser[jQuery.uaMatch(navigator.userAgent).browser] = true;
-}
+const _highlightText = (thisNode, text, className) => {
+ let addItems = [];
+ _highlight(thisNode, addItems, text, className);
+ addItems.forEach((obj) =>
+ obj.parent.insertAdjacentElement("beforebegin", obj.target)
+ );
+};
/**
* Small JavaScript module for the documentation.
*/
-var Documentation = {
-
- init : function() {
- this.fixFirefoxAnchorBug();
- this.highlightSearchWords();
- this.initIndexTable();
- if (DOCUMENTATION_OPTIONS.NAVIGATION_WITH_KEYS) {
- this.initOnKeyListeners();
- }
+const Documentation = {
+ init: () => {
+ Documentation.highlightSearchWords();
+ Documentation.initDomainIndexTable();
+ Documentation.initOnKeyListeners();
},
/**
* i18n support
*/
- TRANSLATIONS : {},
- PLURAL_EXPR : function(n) { return n === 1 ? 0 : 1; },
- LOCALE : 'unknown',
+ TRANSLATIONS: {},
+ PLURAL_EXPR: (n) => (n === 1 ? 0 : 1),
+ LOCALE: "unknown",
// gettext and ngettext don't access this so that the functions
// can safely bound to a different name (_ = Documentation.gettext)
- gettext : function(string) {
- var translated = Documentation.TRANSLATIONS[string];
- if (typeof translated === 'undefined')
- return string;
- return (typeof translated === 'string') ? translated : translated[0];
- },
-
- ngettext : function(singular, plural, n) {
- var translated = Documentation.TRANSLATIONS[singular];
- if (typeof translated === 'undefined')
- return (n == 1) ? singular : plural;
- return translated[Documentation.PLURALEXPR(n)];
- },
-
- addTranslations : function(catalog) {
- for (var key in catalog.messages)
- this.TRANSLATIONS[key] = catalog.messages[key];
- this.PLURAL_EXPR = new Function('n', 'return +(' + catalog.plural_expr + ')');
- this.LOCALE = catalog.locale;
+ gettext: (string) => {
+ const translated = Documentation.TRANSLATIONS[string];
+ switch (typeof translated) {
+ case "undefined":
+ return string; // no translation
+ case "string":
+ return translated; // translation exists
+ default:
+ return translated[0]; // (singular, plural) translation tuple exists
+ }
},
- /**
- * add context elements like header anchor links
- */
- addContextElements : function() {
- $('div[id] > :header:first').each(function() {
- $('\u00B6').
- attr('href', '#' + this.id).
- attr('title', _('Permalink to this headline')).
- appendTo(this);
- });
- $('dt[id]').each(function() {
- $('\u00B6').
- attr('href', '#' + this.id).
- attr('title', _('Permalink to this definition')).
- appendTo(this);
- });
+ ngettext: (singular, plural, n) => {
+ const translated = Documentation.TRANSLATIONS[singular];
+ if (typeof translated !== "undefined")
+ return translated[Documentation.PLURAL_EXPR(n)];
+ return n === 1 ? singular : plural;
},
- /**
- * workaround a firefox stupidity
- * see: https://bugzilla.mozilla.org/show_bug.cgi?id=645075
- */
- fixFirefoxAnchorBug : function() {
- if (document.location.hash && $.browser.mozilla)
- window.setTimeout(function() {
- document.location.href += '';
- }, 10);
+ addTranslations: (catalog) => {
+ Object.assign(Documentation.TRANSLATIONS, catalog.messages);
+ Documentation.PLURAL_EXPR = new Function(
+ "n",
+ `return (${catalog.plural_expr})`
+ );
+ Documentation.LOCALE = catalog.locale;
},
/**
* highlight the search words provided in the url in the text
*/
- highlightSearchWords : function() {
- var params = $.getQueryParameters();
- var terms = (params.highlight) ? params.highlight[0].split(/\s+/) : [];
- if (terms.length) {
- var body = $('div.body');
- if (!body.length) {
- body = $('body');
- }
- window.setTimeout(function() {
- $.each(terms, function() {
- body.highlightText(this.toLowerCase(), 'highlighted');
- });
- }, 10);
- $('