deposition_fraction_plot with correct values
This commit is contained in:
Luis Aleixo
parent
cb7d8abeb8
commit
9810ee3f1f
1 changed files with 38 additions and 85 deletions
|
|
@ -470,107 +470,60 @@ def generate_cdf_curves():
|
|||
plt.show()
|
||||
|
||||
############ Deposition Fraction Graph ############
|
||||
# :D particle diameter - varies between 0.3*10**(-6) and 30*10**(-6)
|
||||
# :Pp mass density - 1000 kg/m-3
|
||||
# :Uair dynamic viscosity of air - 1.8*10**(-5)
|
||||
# :FRC functional residual capacity of lungs - 0.003m3
|
||||
# :Vt tidal volume - 0.0004m3
|
||||
# :g gravity of earth - 9.8m/s**(-2)
|
||||
# :BRk mask - inhalation rate
|
||||
# Cs(D) Cunningham slip factor
|
||||
# lambda = 0.065
|
||||
# A1 = 1.246
|
||||
# A2 = 0.42
|
||||
# A3 = 0.87
|
||||
|
||||
def calculate_cunningham_slip_factor(d:int):
|
||||
_lambda=0.065*10**-6
|
||||
A1 = 1.246
|
||||
A2 = 0.42
|
||||
A3 = 0.87
|
||||
|
||||
return 1 + (((2*_lambda)/d) * (A1+(A2*(math.e**((-A3*d)/_lambda)))))
|
||||
|
||||
# Deposition factors for a breathing model
|
||||
#
|
||||
def calculate_deposition_factor():
|
||||
|
||||
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.81
|
||||
k = 1.38*10**-23
|
||||
T = 300
|
||||
|
||||
diameters = np.linspace(0.001, 0.01, 1000) #particle diameter (multiply later by 10**(-6))
|
||||
diameters = np.append(diameters, np.linspace(0.01, 100, 400))
|
||||
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)
|
||||
diameters = np.linspace(0.001, 0.01, 1000) #particle diameter (μm)
|
||||
diameters = np.append(diameters, np.linspace(0.01, 0.1, 1000))
|
||||
diameters = np.append(diameters, np.linspace(0.1, 1., 1000))
|
||||
diameters = np.append(diameters, np.linspace(1., 10., 1000))
|
||||
diameters = np.append(diameters, np.linspace(10, 100, 1000))
|
||||
|
||||
stk_ = (cunningham_slip_factor*rho_p*(d**2)*(BRk/3600)) / (mu_air * FRC)
|
||||
t_ = (cunningham_slip_factor*g*rho_p*(d**2)*(FRC**(2/3))) / (mu_air * (BRk / 3600))
|
||||
fractions_et = []
|
||||
fractions_tb = []
|
||||
fractions_al = []
|
||||
fractions_df = []
|
||||
for d in diameters:
|
||||
|
||||
IF = 1 - 0.5 * (1 - (1 / (1 + (0.00076*(d**2.8)))))
|
||||
DF_et = IF * (
|
||||
(1 / ( 1 + np.exp(6.84 + 1.183 * np.log(d)))) + (1 / (1 + np.exp(0.924 - 1.885 * np.log(d))))
|
||||
)
|
||||
fractions_et.append(DF_et)
|
||||
|
||||
f_dep_ine = 0.08 + (0.92 / (
|
||||
1 + ((4.09*10**(-6) * (
|
||||
(stk_**0.8) + 0.01 * (t_**0.4) * ((Vt/FRC)**0.8)
|
||||
) ** -2.06))
|
||||
))
|
||||
DF_tb = (0.00352/d) * (np.exp(-0.234*((np.log(d) + 3.40)**2)) + (63.9 * np.exp(-0.819*((np.log(d) - 1.61)**2))))
|
||||
fractions_tb.append(DF_tb)
|
||||
|
||||
X_ = ((k * T * cunningham_slip_factor) / ((3 * math.pi * mu_air * d)) * (((Vt**(1/3)) / (BRk / 3600))))
|
||||
y_ = Vt / FRC
|
||||
|
||||
f_dep_diff = 1 - (
|
||||
1 / ((7380 * (X_ **0.539) * y_**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_μm < 0.1:
|
||||
f_dep = f_dep_diff
|
||||
elif d_μm >= 0.1 and d_μm < 1.:
|
||||
f_dep = f_dep_sed
|
||||
elif d_μm > 1:
|
||||
f_dep = f_dep_ine
|
||||
DF_al = (0.0155/d) * (np.exp(-0.416*((np.log(d) + 2.84)**2)) + (19.11*np.exp(-0.482 * ((np.log(d) - 1.362)**2))))
|
||||
fractions_al.append(DF_al)
|
||||
|
||||
fractions.append(f_dep)
|
||||
activity_fractions.append(fractions)
|
||||
DF = IF * (0.0587 + (0.911/(1 + np.exp(4.77 + 1.485 * np.log(d)))) + (0.943/(1 + np.exp(0.508 - 2.58 * np.log(d)))))
|
||||
fractions_df.append(DF)
|
||||
|
||||
# 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.plot(diameters, fractions_df, label = 'Total')
|
||||
ax.plot(diameters, fractions_et, label = 'Extrathoracic (ET)', ls='--', lw=0.9)
|
||||
ax.plot(diameters, fractions_tb, label = 'Tracheobronchial (TB)', ls='--', lw=0.9)
|
||||
ax.plot(diameters, fractions_al, label = 'Alveolar-interstitial (AI)', ls='--', lw=0.9)
|
||||
|
||||
ax.grid(linestyle='--')
|
||||
ax.legend(handles=lines, loc='upper left')
|
||||
ax.set_xscale('log')
|
||||
ax.margins(x=0, y=0)
|
||||
plt.legend(bbox_to_anchor=(1.04,1), loc="upper left")
|
||||
|
||||
y_ticks = [0., 0.2, 0.4, 0.6, 0.8, 1]
|
||||
x_ticks = [0.001, 0.01, 0.1, 1, 10, 100]
|
||||
plt.yticks(ticks=y_ticks, labels=[
|
||||
str(i) for i in y_ticks])
|
||||
plt.xticks(ticks=x_ticks, labels=[
|
||||
str(i) for i in x_ticks])
|
||||
plt.xlabel('Particle diameter (μm)', fontsize=14)
|
||||
plt.ylabel('Deposition fraction (f$_{dep}$)', fontsize=14)
|
||||
|
||||
|
||||
fig.set_figwidth(10)
|
||||
plt.tight_layout()
|
||||
plt.show()
|
||||
|
||||
######### Auxiliar functions #########
|
||||
|
|
|
|||
Loading…
Reference in a new issue