Add atmospheric settling velocity validation

src/utopia/preprocessing/dry_deposition_MS.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 # ---- Constants dependants of Temperature, Pressure, altitud. Have to be adapted to specific regions ------------
-muAir = 1.789e-5  # [kg−1.s−1], dynamic viscosity of the air (Heigth, Temp dependent)
+muAir = 1.789e-5  # [kg.m−1.s−1], dynamic viscosity of the air (Heigth, Temp dependent)
 rhoAir = 1.1438  # [kg.m-3], density of air (Heigth, Temp dependent)
 g0 = 9.81  # [m.s-2] gravitational acceleration on earth (Heigth dependent)
 

tests/settlingV_air/settlingVtest_XZ.py

@@ -0,0 +1,358 @@
+'''
+Author: Xiaoyu Zhang (xiaoyu.zhang@aces.su.se)
+Date: 2026-02-11 15:26:29
+LastEditTime: 2026-02-13 13:07:51
+Description: this script is for testing the atmospheric settling velocity,
+             and comparing the settling velocity & Reynolds number with measured data.
+             (test data extracted from the paper: https://pubs.acs.org/doi/10.1021/acsestair.5c00250)
+'''
+
+import pandas as pd
+import numpy as np
+import matplotlib.pyplot as plt
+import re
+import warnings
+warnings.filterwarnings('ignore')
+
+
+# ==================== 1. Original Function from dry_deposition_MS.py ====================
+
+# ---- Constants dependants of Temperature, Pressure, altitud. Have to be adapted to specific regions ------------
+muAir = 1.82e-5 #1.789e-5  # [kg.m−1.s−1], dynamic viscosity of the air (Heigth, Temp dependent)
+rhoAir = 1.205 #1.1438  # [kg.m-3], density of air (Heigth, Temp dependent)
+#!!! 260211 XZ testing: revised muAir and rhoAir to the values at 20°C and ori. vaules are commented behind
+
+g0 = 9.81  # [m.s-2] gravitational acceleration on earth (Heigth dependent)
+
+
+# Constant for Cunningham correction, here from Jennings, S, 1988
+A1 = 1.257
+A2 = 0.400
+A3 = 1.100
+mfpAir = 6.635e-8  # [m] Mean free path on dry air, from Jennings, S, 1988
+
+
+# Calculate the Reynolds Number (based on settling velocity of Stokes,
+# may be optimised with self consistancy)
+def ReynoldsNumberFromStokes(d, rho):
+    vg = np.multiply(d, d)
+    vg = vg * g0 / (18.0 * muAir)
+    vg = np.multiply(vg, (rho - rhoAir))
+
+    Rep = np.multiply(vg, rhoAir)
+    Rep = np.multiply(Rep, d)
+    Rep = np.multiply(Rep, (1.0 / muAir))
+
+    return Rep
+
+
+# Placeholder function to calculate Reynolds number
+def ReynoldsNumberFromVg(d, rho, vg):
+    # Placeholder implementation
+    # This function should calculate the Reynolds number based on the settling velocity
+    Rep = np.multiply(vg, rhoAir)
+    Rep = np.multiply(Rep, d)
+    Rep = np.multiply(Rep, (1.0 / muAir))
+
+    return Rep
+
+
+# Calculate the Drag coefficient number
+def dragCoefficient(d, rho, Rep):
+    # Rep = ReynoldsNumberFromStokes (d, rho)
+    # Set Rep to a specific value for test (e.g., 13.4) !!!!!!!!!!! to remoove after test
+
+    # Rep = 13.35
+    Cd = np.zeros_like(Rep)
+
+    Rep_conditions = [
+        Rep <= 1,
+        (1 < Rep) & (Rep <= 1000.0),
+        (1000.0 < Rep) & (Rep <= 2.0e5),
+        Rep > 2.0e5,
+    ]
+
+    Cd_functions = [
+        lambda x: 24.0 / x,
+        lambda x: 24.0 / x + 5.0 / np.power(x, 0.6) + 0.44,
+        0.44,
+        0.10,
+    ]
+
+    Cd = np.piecewise(Rep, Rep_conditions, Cd_functions)
+
+    # Cd = 24.0 * (1.0 + (0.15 * Rep**0.687)) / Rep
+    # Cd = Cd + 0.42 / (1.0 + (42500.0/ (Rep**1.16)))
+    return Cd
+
+
+# Calculate rate constants of dry settling based on Newton regime
+# (for big particles, generating a turbulent flow)
+def kineticCstdrySettlingNewtonSphere(d, rho, Rep):
+    d = np.array(d)
+    rho = np.array(rho)
+
+    # Calculate the drag coefficient Cd
+    Cd = dragCoefficient(d, rho, Rep)
+    # print("Cd sphere=", Cd)
+    # Calculate the settling velocity
+    v = np.sqrt(4.0 * d * g0 * (rho - rhoAir) / (3.0 * Cd * rhoAir))
+    # v = np.sqrt(4.0 * d * g0 * (rho) / (3.0 * Cd * rhoAir))
+
+    # elt = v / shape  # [m.s-1] / [m] The Boundary layer heigth
+    elt = v
+    return elt
+
+
+# Placeholder function to calculate settling velocity
+def calculate_settling(reynolds_number):
+    # Placeholder implementation
+    # This function should calculate a new estimate for the settling velocity based on the Reynolds number
+    return new_settling_velocity
+
+
+# Subroutine to calculate settling velocity iteratively
+def get_settling(initial_Settling, d, rho, initial_Rep):
+
+    settling_old = initial_Settling
+
+    # Set convergence threshold
+    tolerance = 0.001
+
+    # Maximum number of iterations
+    max_iterations = 20
+
+    # Initialize iteration counter
+    iteration = 0
+
+    # Iterate until convergence or maximum iterations reached
+    while iteration < max_iterations:
+        # Calculate Reynolds number using current settling velocity
+        reynolds = ReynoldsNumberFromVg(d, rho, settling_old)
+
+        # Use Reynolds number to calculate new estimate for settling velocity
+        settling_new = kineticCstdrySettlingNewtonSphere(d, rho, reynolds)
+
+        # Check for convergence
+        cvg = abs((settling_new - settling_old) / settling_new)
+        # print("convergence <", tolerance, "?=", cvg)
+        if np.all(cvg < tolerance):
+            break
+
+        # Update settling velocity for next iteration
+        settling_old = settling_new
+
+        # Increment iteration counter
+        iteration += 1
+
+    # Use final settling velocity for further calculations
+    final_settling_velocity = settling_new
+    reynolds = ReynoldsNumberFromVg(d, rho, final_settling_velocity)
+
+    return final_settling_velocity
+
+
+
+
+# ==================== 2. Read data and process ====================
+def load_and_prepare(csv_path):
+    df = pd.read_csv(csv_path, encoding='utf-8-sig')
+
+    # Fuzzy match key column names (adapt to your latest headers)
+    def find_col(patterns):
+        for col in df.columns:
+            col_lower = str(col).lower()
+            if any(p.lower() in col_lower for p in patterns):
+                return col
+        return None
+
+    col_map = {
+        'diameter': find_col(['diameter', '(mm)']),
+        'density': find_col(['density']),
+        'vs_meas': find_col(['settling velocity', 'v (m s-1)']),
+        're_meas': find_col(['rep']),
+        'shape': find_col(['shape']),
+        'length': find_col(['l (mm)']),
+        'aspect': find_col(['aspect ratio']),
+        'vs_sd': find_col(['settling velocity sd'])
+    }
+
+    print("Automatically detected columns:")
+    for k, v in col_map.items():
+        print(f"  {k:10} -> {v}")
+
+    # Extract main diameter value (handle parentheses)
+    def extract_number(val):
+        if pd.isna(val):
+            return np.nan
+        s = str(val).strip()
+        match = re.search(r'([-+]?\d*\.?\d+)', s)
+        return float(match.group(1)) if match else np.nan
+
+    df['diameter_mm'] = df[col_map['diameter']].apply(extract_number)
+    df['diameter_m'] = df['diameter_mm'] / 1000.0
+    df['density_kgm3'] = pd.to_numeric(df[col_map['density']], errors='coerce') * 1000.0
+    df['vs_measured'] = pd.to_numeric(df[col_map['vs_meas']], errors='coerce')
+    df['vs_sd'] = pd.to_numeric(df[col_map['vs_sd']], errors='coerce')
+    df['Re_measured'] = pd.to_numeric(df[col_map['re_meas']], errors='coerce')
+
+    # delete rows with missing critical data
+    initial_len = len(df)
+    df.dropna(subset=['diameter_m', 'density_kgm3', 'vs_measured', 'Re_measured'], inplace=True)
+    df.reset_index(drop=True, inplace=True)
+    print(f"Valid data: {len(df)} rows (originally {initial_len} rows)")
+
+    return df, col_map
+
+# ==================== 3. Calculate settling velocity and Reynolds number ====================
+def compute_settling(df):
+    vs_calc_list = []
+    Re_calc_list = []
+
+    for idx, row in df.iterrows():
+        d = np.array([row['diameter_m']])
+        rho_p = np.array([row['density_kgm3']])
+
+        vs_stokes = (d**2 * g0 * (rho_p - rhoAir)) / (18.0 * muAir)
+        Re_stokes = ReynoldsNumberFromStokes(d, rho_p)
+
+        vs_final = get_settling(vs_stokes, d, rho_p, Re_stokes)
+        vs_final = vs_final[0] if isinstance(vs_final, np.ndarray) else vs_final
+
+        Re_final = ReynoldsNumberFromVg(d, rho_p, np.array([vs_final]))[0]
+
+        vs_calc_list.append(vs_final)
+        Re_calc_list.append(Re_final)
+
+    df['vs_calc'] = vs_calc_list
+    df['Re_calc'] = Re_calc_list
+
+    # Relative error
+    df['vs_error_%'] = (df['vs_calc'] - df['vs_measured']) / df['vs_measured'] * 100.0
+    df['Re_error_%'] = (df['Re_calc'] - df['Re_measured']) / df['Re_measured'] * 100.0
+
+    return df
+
+# ==================== 4. Save results to CSV ====================
+def save_results(df, col_map):
+    output_cols = [
+        'Number', 'Microplastics', 'Material',
+        col_map['density'], col_map['shape'], 'diameter_mm',
+        col_map['length'], col_map['aspect'],
+        col_map['vs_meas'], col_map['vs_sd'],
+        'vs_calc', 'vs_error_%',
+        col_map['re_meas'], 'Re_calc', 'Re_error_%'
+    ]
+    output_cols = [c for c in output_cols if c in df.columns]
+    df_out = df[output_cols].copy()
+    df_out.to_csv('settling_velocity_results.csv', index=False, encoding='utf-8-sig')
+    print("Result saved to settling_velocity_results.csv")
+
+# ==================== 5. Plotting (Settling Velocity with Error Bars) ====================
+def plot_results(df, col_map):
+    shape_colors = {
+        'sphere': '#1f78b4', 'disk': '#5a9755', 'fragment': '#e31a1c',
+        'fiber': '#ff7f00', 'film': '#6a3d9a'
+    }
+    df['shape_lower'] = df[col_map['shape']].str.lower()
+    shapes = df['shape_lower'].unique()
+
+    # ----- Figure 1: Measured vs Calculated (1x2) -----
+    plt.figure(figsize=(14, 6))
+
+    # Subplot 1: Settling Velocity (with error bars)
+    plt.subplot(1, 2, 1)
+    for shape in shapes:
+        subset = df[df['shape_lower'] == shape].sort_values('diameter_mm')
+        # Measured values: hollow circles + error bars
+        plt.errorbar(subset['diameter_mm'], subset['vs_measured'],
+                     yerr=subset['vs_sd'],
+                     fmt='o',
+                     markerfacecolor='none',
+                     markeredgecolor=shape_colors.get(shape, 'gray'),
+                     markeredgewidth=1.5,
+                     ecolor=shape_colors.get(shape, 'gray'),
+                     elinewidth=1,
+                     capsize=3,
+                     capthick=1,
+                     label=f'{shape.title()} (Measured)')
+        # Calculated values: crosses
+        plt.scatter(subset['diameter_mm'], subset['vs_calc'],
+                    color=shape_colors.get(shape, 'gray'),
+                    marker='x',
+                    label=f'{shape.title()} (Calc)')
+    plt.xlabel('Diameter (mm)')
+    plt.ylabel('Settling Velocity (m/s)')
+    plt.title('Settling Velocity: Measured vs Calculated')
+    plt.grid(True, linestyle='--', alpha=0.7)
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=8)
+
+    # Subplot 2: Reynolds Number (no error bars)
+    plt.subplot(1, 2, 2)
+    for shape in shapes:
+        subset = df[df['shape_lower'] == shape]
+        plt.scatter(subset['diameter_mm'], subset['Re_measured'],
+                    edgecolors=shape_colors.get(shape, 'gray'), facecolors='none',
+                    label=f'{shape.title()} (Measured)', marker='o')
+        plt.scatter(subset['diameter_mm'], subset['Re_calc'],
+                    color=shape_colors.get(shape, 'gray'),
+                    label=f'{shape.title()} (Calc)', marker='x')
+    plt.xlabel('Diameter (mm)')
+    plt.ylabel('Reynolds Number')
+    plt.title('Reynolds Number: Measured vs Calculated')
+    plt.grid(True, linestyle='--', alpha=0.7)
+    plt.xscale('log')
+    plt.yscale('log')
+    plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', fontsize=8)
+
+    plt.tight_layout()
+    plt.savefig('settling_velocity_comparison.png', dpi=300)
+    plt.show()
+
+    # ----- Figure 2: Error scatter plots (1x2, no error bars) -----
+    plt.figure(figsize=(14, 6))
+
+    # Settling velocity error
+    plt.subplot(1, 2, 1)
+    for shape in shapes:
+        subset = df[df['shape_lower'] == shape].dropna(subset=['vs_error_%'])
+        plt.scatter(subset['diameter_mm'], subset['vs_error_%'],
+                    color=shape_colors.get(shape, 'gray'),
+                    label=shape.title(), s=80)
+    plt.axhline(y=0, color='k', linestyle='--', alpha=0.5)
+    plt.xlabel('Diameter (mm)')
+    plt.ylabel('Error (%)')
+    plt.title('Settling Velocity Error')
+    plt.grid(True, linestyle='--', alpha=0.7)
+    plt.xscale('log')
+    plt.ylim(-100, 200)
+    plt.legend()
+
+    # Reynolds number error
+    plt.subplot(1, 2, 2)
+    for shape in shapes:
+        subset = df[df['shape_lower'] == shape].dropna(subset=['Re_error_%'])
+        plt.scatter(subset['diameter_mm'], subset['Re_error_%'],
+                    color=shape_colors.get(shape, 'gray'),
+                    label=shape.title(), s=80)
+    plt.axhline(y=0, color='k', linestyle='--', alpha=0.5)
+    plt.xlabel('Diameter (mm)')
+    plt.ylabel('Error (%)')
+    plt.title('Reynolds Number Error')
+    plt.grid(True, linestyle='--', alpha=0.7)
+    plt.xscale('log')
+    plt.ylim(-100, 200)
+    plt.legend()
+
+    plt.tight_layout()
+    plt.savefig('settling_velocity_errors.png', dpi=300)
+    plt.show()
+
+# ==================== 6. Main Program ====================
+if __name__ == '__main__':
+    df, col_map = load_and_prepare('settlingdata_re.csv')
+    df = compute_settling(df)
+    save_results(df, col_map)
+    plot_results(df, col_map)

tests/settlingV_air/settling_velocity_results.csv

@@ -0,0 +1,24 @@
+Number,Microplastics,Material,Density (g cm-3),Shape,diameter_mm,L (mm),Aspect Ratio,Settling Velocity v (m s-1),Settling Velocity SD(m s-1),vs_calc,vs_error_%,Rep,Re_calc,Re_error_%
+1,Bead0.083,PE,0.95,sphere,0.083,0.083,1.0,0.34,0.11,0.18281945999531735,-46.22957058961255,1.84,1.0046532028259159,-45.39928245511327
+2,Bead0.116,PE,0.95,sphere,0.116,0.116,1.0,0.42,0.09,0.2882539946223951,-31.368096518477355,3.48,2.2138540312262855,-36.38350484981938
+3,Bead0.165,PE,0.95,sphere,0.165,0.165,1.0,0.66,0.12,0.5102837806816257,-22.68427565429914,7.12,5.574569928243091,-21.70547853591165
+4,Bead0.196,PE,0.95,sphere,0.196,0.196,1.0,0.95,0.15,0.6617748932730574,-30.339484918625537,12.17,8.587801884243447,-29.434659948698048
+5,Bead0.231,PE,0.95,sphere,0.231,0.231,1.0,1.19,0.24,0.8370790596566882,-29.65722187758923,17.97,12.80248031048008,-28.756370002893263
+6,Disk0.188,PET,1.38,disk,0.19,0.30(0.04),1.6,0.91,0.22,0.8584227403906463,-5.667830726302612,11.18,10.79867507760651,-3.410777481158219
+7,Disk0.247,PET,1.38,disk,0.25,0.33(0.05),1.35,0.94,0.26,1.243648894661229,32.30307390013075,15.18,20.585122500917326,35.606867594975796
+8,Disk0.418,PET,1.38,disk,0.42,1.15(0.06),2.76,1.55,0.31,2.2506745608137724,45.20481037508208,42.35,62.58606567185991,47.78291776118043
+9,Disk0.577,PET,1.38,disk,0.58,1.22(0.10),2.12,1.99,0.4,3.0528737059485325,53.41073899238856,75.05,117.23370511469396,56.207468507253786
+10,Disk0.674,PET,1.38,disk,0.67,2.16(0.09),3.2,1.42,0.55,3.4538691860093227,143.23022436685372,62.55,153.2132575453092,144.94525586780048
+11,Disk0.787,PET,1.38,disk,0.79,2.39(0.12),3.03,2.08,0.59,3.942940677836301,89.56445566520676,106.99,206.2352955091356,92.761281904043
+12,Fragment0.510,PA,1.15,fragment,0.51,1.10(0.38),2.1,1.64,0.33,2.414739980330067,47.240242703052864,54.67,81.53727774240895,49.14446267131689
+13,Fragment0.605,PA,1.15,fragment,0.61,1.83(0.51),2.98,1.94,0.45,2.846465863132925,46.725044491387884,82.42,114.9612490492229,39.482224034485434
+14,Flock0.194,PA,1.15,fiber,0.194,2,10.31,0.4,0.12,0.7624865311333263,90.62163278333158,5.07,9.793762438628445,93.17085677768134
+15,Flock0.222,PA,1.15,fiber,0.222,3,13.51,0.38,0.13,0.9222467465646966,142.69651225386752,5.51,13.555506987556155,146.0164607541952
+16,Flock0.298,PA,1.15,fiber,0.298,4,13.4,0.49,0.13,1.3504516923179228,175.6023861873312,9.54,26.64470869200236,179.29464037738325
+17,Fiber0.640,PES,1.38,fiber,0.64,1.47(1.34),2.0,0.46,0.2,3.323332890240265,622.4636717913619,19.24,140.82166620622488,631.9213420281959
+18,Fiber0.909,PES,1.38,fiber,0.91,2.81(1.41),3.02,0.53,0.21,4.390203706041307,728.3403218945862,31.49,264.5097732889888,739.9802263861188
+19,Fiber1.102,PES,1.38,fiber,1.1,4.20(2.26),4.07,0.68,0.25,5.03202190933213,640.0032219606073,48.98,366.4804967483373,648.2247789880305
+20,Fiber1.345,PES,1.38,fiber,1.35,8.50(3.59),6.19,0.69,0.33,5.7829810957148595,738.113202277516,60.66,516.8936537062721,752.11614524608
+21,Film0.650,PE,0.95,film,0.65,1.39(0.53),2.09,0.67,0.34,2.6700646116523847,298.5171062167738,28.46,114.9081377514687,303.7531192953924
+22,Film0823,PE,0.95,film,0.82,1.55(0.97),1.86,1.0,0.33,3.248002541321009,224.8002541321009,53.79,176.3379841252357,227.82670408112233
+23,Film1.234,PE,0.95,film,1.23,2.83(1.51),2.23,1.14,0.41,4.40265575404312,286.1978731616772,91.95,358.53825416785776,289.92741073176484