From 34b1f89d69aebd73ae3359b39a84e93addbe2053 Mon Sep 17 00:00:00 2001 From: Luis Aleixo Date: Tue, 14 Sep 2021 14:25:23 +0200 Subject: [PATCH] deposition fraction for a breathing model --- cara/plot_output.py | 2 +- cara/results_paper.py | 98 ++++++++++++++++++++++++++----------------- 2 files changed, 60 insertions(+), 40 deletions(-) diff --git a/cara/plot_output.py b/cara/plot_output.py index a231d25f..8797554c 100644 --- a/cara/plot_output.py +++ b/cara/plot_output.py @@ -39,7 +39,7 @@ print('\n') #present_vl_er_histograms(activity='Heavy exercise', mask='No mask') ############ CDFs for comparing the QR-Values in different scenarios ############ -generate_cdf_curves() +#generate_cdf_curves() ############ Deposition Fraction Graph ############ print('\n<<<<<<<<<<< Deposition Fraction for Breathing, seated >>>>>>>>>>>') diff --git a/cara/results_paper.py b/cara/results_paper.py index 237aa727..32b12257 100644 --- a/cara/results_paper.py +++ b/cara/results_paper.py @@ -491,64 +491,84 @@ def calculate_cunningham_slip_factor(d:int): return 1 + (((2*_lambda)/d) * (A1+(A2*(math.e**((-A3*d)/_lambda))))) -# Deposition factor for a breathing model +# Deposition factors for a breathing model def calculate_deposition_factor(): - exposure_mc = breathing_exposure(activity='Heavy exercise', mask='No mask') - exposure_model = exposure_mc.build_model(size=SAMPLE_SIZE) + fig = plt.figure() + ax = fig.add_subplot(1, 1, 1) + + ############ Breathing models ############ + br_seated = breathing_seated_exposure() + br_seated_model = br_seated.build_model(size=SAMPLE_SIZE) + + br_light_activity = breathing_light_activity_exposure() + br_light_activity_model = br_light_activity.build_model(size=SAMPLE_SIZE) + + br_heavy_exercise = breathing_heavy_exercise_exposure() + br_heavy_exercise_model = br_heavy_exercise.build_model(size=SAMPLE_SIZE) rho_p = 1000 mu_air = 1.8*10**-5 FRC = 0.003 Vt = 0.0004 g = 9.8 - BRk = exposure_model.exposed.activity.inhalation_rate k = 1.38*10**-23 T = 300 - - diameters = np.linspace(0.01, 100, 200) #particle diameter (multiply later by 10**(-6)) - fractions = [] - for d in diameters: - d_μm = d*10**(-6) - cunningham_slip_factor = calculate_cunningham_slip_factor(d_μm) + diameters = np.linspace(0.001, 100, 200) #particle diameter (multiply later by 10**(-6)) + activity_fractions = [] + for scenario in (br_seated_model, br_light_activity_model, br_heavy_exercise_model): + BRk = scenario.exposed.activity.inhalation_rate + fractions = [] + for d_μm in diameters: + d = d_μm*10**(-6) + cunningham_slip_factor = calculate_cunningham_slip_factor(d) - f_dep_ine = 0.08 + 0.92 / ( - 1 + (4.09*10**-6 * ( - (((cunningham_slip_factor*rho_p*d_μm**2*(BRk/3600))/mu_air*FRC)**0.8) + ( - 0.01*( - ((cunningham_slip_factor*g*rho_p*d_μm**2*FRC**(2/3))/(mu_air*(BRk/3600))**0.4) * ( - (Vt/FRC)**0.8 - ) + f_dep_ine = 0.08 + 0.92 / ( + 1 + (4.09*10**-6 * ( + (((cunningham_slip_factor*rho_p*d**2*(BRk/3600))/mu_air*FRC)**0.8) + ( + 0.01*( + ((cunningham_slip_factor*g*rho_p*d**2*FRC**(2/3))/(mu_air*(BRk/3600))**0.4) * ( + (Vt/FRC)**0.8 + ) ) ) )**(-2.06) )) - f_dep_diff = 1 - (1/ - ((7380*( - ((k * T * cunningham_slip_factor)/3*math.pi*mu_air*d_μm) * - ((Vt**(1/3)) / (BRk/3600)))**0.539) * ((Vt/FRC)**0.884) + 1)) + f_dep_diff = 1 - 1 / ( + 7380*(((k * T * cunningham_slip_factor)/(3 * math.pi * mu_air * d)*(Vt**(1/3))/(BRk/3600))**0.539 * ((Vt/FRC)**0.884)) + 1) - f_dep_sed = (0.431 * f_dep_ine) + (0.541 * f_dep_diff) + (1.060 * f_dep_ine**2) + (0.685 * f_dep_diff**2) - (1.521 * f_dep_ine * f_dep_diff) - - # Eq. S.1 - if d < 0.3: - f_dep = f_dep_diff - elif d >= 0.3 and d <= 1.: - f_dep = f_dep_sed - elif d > 1: - f_dep = f_dep_ine - - fractions.append(f_dep) - - fig = plt.figure() - ax = fig.add_subplot(1, 1, 1) - ax.plot(diameters, fractions) - + f_dep_sed = (0.431 * f_dep_ine) + (0.541 * f_dep_diff) + (1.060 * f_dep_ine**2) + (0.685 * f_dep_diff**2) - (1.521 * f_dep_ine * f_dep_diff) + + # Eq. S.1 + if d_μm < 0.1: + f_dep = f_dep_diff + elif d_μm > 1: + f_dep = f_dep_ine + else: + f_dep = f_dep_sed + + fractions.append(f_dep) + activity_fractions.append(fractions) + + # Colors can be changed here + colors = ['skyblue', 'royalblue', 'navy'] + + breathing_rates = ['Seated', 'Light activity', 'Heavy activity'] + + lines = [mlines.Line2D([], [], color=color, markersize=15, label=label) + for color, label in zip(colors, breathing_rates)] + + for i in range(3): + ax.plot(diameters, activity_fractions[i], c=colors[i]) + + ax.grid(linestyle='--') + ax.legend(handles=lines, loc='upper left') + ax.set_xscale('log') - plt.ylabel( - 'Deposition fraction (f$_{dep}$)', fontsize=14) plt.xlabel('Particle diameter (μm)', fontsize=14) + plt.ylabel('Deposition fraction (f$_{dep}$)', fontsize=14) + plt.show() ######### Auxiliar functions #########