Линч - Витены бутлалтын загвар (Lynch - Whiten)

Энэхүү python программ нь 1977 онд Линчийн боловсруулсан матрицын аргаар бутлалтын циклыг загварчлалыг гүйцэтгэнэ.

# ---------- Fresh feed ----------
Q_fresh = 1812.0      # tph
F80_feed_mm = 45.0    # for information

# Example feed PSD (replace with real data!)
feed_psd = np.array([0.02, 0.04, 0.09, 0.20, 0.30, 0.35], dtype=float)
feed_psd /= feed_psd.sum()

print("Check feed P80 (mm):", calc_P80(feed_psd))

def build_breakage_matrix(size_mid):
    """
    Simple lower-triangular Whiten-style breakage matrix.
    You can calibrate this later from testwork.
    """
    n = len(size_mid)
    B = np.zeros((n, n))
    for i in range(n):
        B[i, i] = 0.4
        if i >= 1:
            B[i, i-1] += 0.3
        if i >= 2:
            B[i, i-2] += 0.2
        if i >= 3:
            B[i, i-3] += 0.1
        B[i] /= B[i].sum()
    return B
def build_classification_matrix(size_mid, css_mm=18.0):
    """
    Lynch selection function (probability of breakage).
    """
    n = len(size_mid)
    C = np.zeros((n, n))
    for i, d in enumerate(size_mid):
        if d > 2.5 * css_mm:
            pb = 0.95
        elif d > css_mm:
            pb = 0.70
        elif d > 0.5 * css_mm:
            pb = 0.40
        else:
            pb = 0.05
        C[i, i] = pb
    return C

B = build_breakage_matrix(size_mid)
C = build_classification_matrix(size_mid, css_mm=18.0)  # guess CSS

def cone_crusher(feed_psd, B, C):
    """
    Whiten cone crusher:
        P = F (I - C) + F C B
    """
    F = feed_psd.reshape(1, -1)
    I = np.eye(len(feed_psd))
    P = F @ (I - C) + F @ C @ B
    P = P.ravel()
    P = np.maximum(P, 0)
    P /= P.sum()
    return P

def build_partition_curve(size_mid, d50=14.0, sharpness=2.0):
    """
    Partition curve: probability to undersize.
      p_u(d) = 1 / (1 + (d/d50)^sharpness)
    """
    d = size_mid
    return 1.0 / (1.0 + (d / d50) ** sharpness)

p_u = build_partition_curve(size_mid, d50=14.0, sharpness=2.0)

def screen_step(feed_mass_by_class, p_u):
    """
    feed_mass_by_class: tph per size class
    """
    mass_U = feed_mass_by_class * p_u
    mass_O = feed_mass_by_class - mass_U

    Q_U = mass_U.sum()
    Q_O = mass_O.sum()

    psd_U = mass_U / Q_U if Q_U > 0 else np.zeros_like(mass_U)
    psd_O = mass_O / Q_O if Q_O > 0 else np.zeros_like(mass_O)

    return mass_U, mass_O, Q_U, Q_O, psd_U, psd_O
    def simulate_circuit(Q_fresh, feed_psd, B, C, p_u,
                     max_iter=1000, tol=1e-4):
    n = len(feed_psd)

    Q_recycle = 0.0
    psd_recycle = np.zeros(n)

    for it in range(max_iter):
        Q_cr_feed = Q_fresh + Q_recycle

        if Q_cr_feed > 0:
            psd_cr_feed = (Q_fresh * feed_psd + Q_recycle * psd_recycle) / Q_cr_feed
        else:
            psd_cr_feed = feed_psd.copy()

        psd_cr_prod = cone_crusher(psd_cr_feed, B, C)

        mass_screen_feed = Q_cr_feed * psd_cr_prod
        (mass_U, mass_O,
         Q_U, Q_O,
         psd_U, psd_O) = screen_step(mass_screen_feed, p_u)

        # screen feed PSD for plotting
        psd_screen_feed = mass_screen_feed / mass_screen_feed.sum()

        Q_recycle_new = Q_O
        psd_recycle_new = psd_O.copy()

        rel_change = abs(Q_recycle_new - Q_recycle) / max(Q_cr_feed, 1e-6)
        Q_recycle, psd_recycle = Q_recycle_new, psd_recycle_new

        if rel_change < tol:
            print(f"Converged in {it+1} iterations.")
            break

    results = {
        # rates
        "Q_fresh": Q_fresh,
        "Q_crusher_feed": Q_cr_feed,
        "Q_recycle": Q_recycle,
        "Q_product": Q_U,
        # PSDs
        "psd_fresh": feed_psd,
        "psd_crusher_feed": psd_cr_feed,
        "psd_crusher_product": psd_cr_prod,
        "psd_screen_feed": psd_screen_feed,
        "psd_screen_oversize": psd_O,
        "psd_screen_undersize": psd_U,
    }
    return results
    
    res = simulate_circuit(Q_fresh, feed_psd, B, C, p_u)

Q_prod = res["Q_product"]
P80_prod = calc_P80(res["psd_screen_undersize"])

print("\n=== Steady-state summary ===")
print(f"Fresh feed rate          : {res['Q_fresh']:.1f} tph")
print(f"Crusher feed rate        : {res['Q_crusher_feed']:.1f} tph")
print(f"Recycle (screen oversize): {res['Q_recycle']:.1f} tph")
print(f"Final product rate       : {Q_prod:.1f} tph")
print(f"Final product P80 (mm)   : {P80_prod:.2f} mm")