直接压片单元模型开发教程（11）：成形性模型——Ryshkewitch-Duckworth 与 Kuentz-Leuenberger 方程

直接压片单元模型开发教程（11）：成形性模型——Ryshkewitch-Duckworth 与 Kuentz-Leuenberger 方程

版本声明块

• 工具/软件：本系列自研TabletPress单元类的成形性拟合模块（本系列二次开发示例类，非官方 API）；压缩表征术语依 USP⟨1062⟩，抗张强度与破碎力表述依 USP⟨1217⟩。
• 语言/环境：Python 3.9 + numpy（<2）+ scipy（<1.14）+ matplotlib
• 本文目标：把"抗张强度随固体分数近似指数增大"这句定性话，变成两个可用scipy.optimize.curve_fit标定、可外推、可报误差的成形性模型。

一句话结论：Ryshkewitch–Duckworth（R–D）方程把"抗张强度–孔隙率"写成指数关系，其三种常见等价写法σ = σ0·e^(−k·ε)、log(TS) = k·SF + A、ln(TS/TS0) = −k_b·ε只差一个对数底，标定时必须统一底（混用 ln 与 log10 会使 k 相差 ln10≈2.3026 倍）；k_b<10 易压、k_b>10 难压是文献经验阈值、以原文为准；Kuentz–Leuenberger（K–L）相对密度形式σ_t,p = σ_max·[1 − (1−ρ_r)/(1−ρ_c,σ)·e^(ρ_r−ρ_c,σ)]用σ_max（全密度抗张强度）与ρ_c,σ（强度临界密度）两个参数给出"低于临界密度则无强度"的物理图像，但其精确代数形式仅得二次引用，须以原始论文为准。

〇、本篇要解决的认知问题

1. 同样是描述"抗张强度–孔隙率"，为什么会出现 R–D 的三种写法？它们到底是不是同一个方程？
2. σ0、k、k_b、A分别代表什么物理量？键合能力k_b的"<10 易压 / >10 难压"能当常数用吗？
3. Kuentz–Leuenberger 方程里的σ_max与强度临界密度ρ_c,σ是什么？它比 R–D 多说了什么？
4. compressibility（可压性）、compactibility（成形性）、tabletability（压片性）三个词在建模里各对应哪条曲线，为什么不能混用？
5. 给定一组"固体分数–抗张强度"实测点，怎么用 numpy/scipy 稳健地拟合出参数并给出可信度？

一、机制解析

1.1 compactibility 是"密度→强度"的映射，不是"压力→强度"

要理解成形性模型，先把三个极易混用的术语钉死（底账第 6.6 节，术语依 USP⟨1062⟩）：

关键点在于自变量不同：可压性问"给多大压力能压多密"，成形性问"压到某个密度能有多硬"。这两条曲线的自变量分别是压力与密度。R–D 与 K–L 都属于成形性模型——它们的横轴是固体分数 SF（或孔隙率 ε），纵轴是抗张强度 TS。USP⟨1062⟩ 明确：抗张强度随固体分数增大而近似指数增加，粉体成形性常用 R–D 方程良好描述。

这就解释了一个反直觉现象：两块片剂如果固体分数相同，成形性模型会预测它们抗张强度相同，而不管你是"一次压到"还是"两段压到"。这既是 S–D/K–L 的便利，也是它们的局限——真正的压片性还依赖压力历史与应变速率（后续第 13、14、18 篇接力）。

压实路径（可压性 Heckel） 成形性（R–D / K–L） 压力 P ─────────────► 固体分数 SF ─────────────► 抗张强度 TS (固体分数–压力曲线) (抗张强度–固体分数曲线) ▲ ▲ └────── 压片性 tabletability ──┘ (抗张强度–压力，两段复合)

这张 ASCII 图是本篇与前后篇的分工地图：第 8 篇出 SF§，第 11 篇出 TS(SF)，两者串起来才是 TS§。

1.2 R–D 的三种写法：同一方程的不同坐标

底账第 6.1 节登记了 R–D 的三种等价写法：

它们其实是同一件事在"不同对数底 + 不同自变量（孔隙率 ε 还是固体分数 SF）"下的表达。把指数形式两边取自然对数：

ln σ = ln σ0 − k·ε ⟹ ln(TS/TS0) = −k·ε

与第三种写法逐字一致，故k_b ≡ k（在都用自然对数时）。再代入SF = 1 − ε：

log(TS) = k·(1 − ε) + A = (k + A) − k·ε

与第一种写法也自洽，只是截距被平移了。唯一的陷阱是对数底：若第一种用e（自然对数隐含）、第二种在文献里用log10，那两条直线的斜率会差一个ln 10 ≈ 2.3026的因子。这是标定时最常见的翻车点，第 2 节代码里会显式验证。

1.3 键合能力 k_b：经验阈值，不是物理常数

R–D 的核心输出之一是把材料分为"好压"与"难压"两档。底账明确记录：键合能力判别k_b<10 易压、k_b>10 难压属文献经验阈值，须以原文为准。这句话有三层含义，建模时必须听懂：

1. 它是一个经验分档，来自特定文献在特定固体分数窗口内的拟合结果，不是从第一性原理推出的普适常数。
2. k_b的数值依赖拟合窗口——你在 SF=0.7–0.95 上拟合与在 SF=0.5–1.0 上拟合，得到的k_b会不同。这正是K/Py类参数"材料特异、须实测拟合"的普遍规律（底账第 9 项"不可断言项"）。
3. 因此教程里可以打印"该批数据 k_b=…，按文献经验分档属易压/难压"，但不得把 10 这个阈值写成绝对判据，也不能据此对未测材料下结论。

1.4 Kuentz–Leuenberger：从相对密度直接给强度

K–L 相对密度形式为：

σ_t,p = σ_max · [ 1 − (1 − ρ_r)/(1 − ρ_c,σ) · e^(ρ_r − ρ_c,σ) ]

其中σ_max是全密度时的抗张强度，ρ_r是片剂相对密度，ρ_c,σ是强度临界密度。这个方程的物理图像很干净，取两个端点看：

• 当ρ_r = ρ_c,σ：括号内1 − 1·e^0 = 1 − 1 = 0，即σ=0。低于临界密度时片剂根本没有强度（现实里表现为"压不成型/一碰就散"）。
• 当ρ_r = 1：(1−ρ_r)=0，括号内1 − 0 = 1，即σ=σ_max，回到全密度强度上限。

也就是说，K–L 用一个"临界密度"把 R–D 的纯指数关系截断成"有阈值、再单调上升"的 S 形起步曲线，比 R–D 多描述了一件事：材料存在一个成形下限密度。对高孔隙率段（压得很松）的预测，K–L 通常比 R–D 更贴近"压不成型"的真实行为。

但必须诚实标注（底账"不可断言项"第 3 项）：K–L 方程的精确代数形式仅得二次引用，指数项精确写法各文献转录略有差异（部分文献指数含 γ 系数），本篇采用 Su 等（2018）可检索的二次引用形式，建议核对原始论文。同理，Leuenberger 的压力型 compactibility 方程σ = σ_max·[1 − e^(−γ·…·P)]（γ 为压缩敏感度）指数具体形式各文献表述不一，以原始论文为准。教程代码里我会把这些二次引用形式全部写进注释并用"以原始论文为准"收尾。

1.5 三模型的分工与本篇在系列中的位置

到这里可以把成形性模型的"零件"列清楚：

• R–D：无阈值、纯指数，刻画"密度越高越强"的整体趋势，参数少（σ0、k）好标定。
• K–L：带强度临界密度 ρ_c,σ，描述"低于某密度不成型"，对松片段更真实，但形式来源需核对。
• Tye 等（2005）：用 R–D 型指数关系描述"表压速度—抗张强度—固体分数"，并未提出与 R–D 不同的独立公式（底账明确）——所以不要把它当成第四个方程。

本篇是整个系列"TS(SF) 曲线"的生产者；它吃进的是第 12 篇由破碎力换算来的抗张强度、第 10 篇算出的固体分数，吐出的是可供第 19 篇 QbD 响应面直接调用的 compactibility 参数。

1.6 为什么"密度→强度"会长成指数：物理直觉与配套关系

R–D 由 Ryshkewitch（1953）与 Duckworth（1953）提出，被广泛用于拟合"抗张强度–孔隙率"以获得完整的 compactibility 曲线（底账第 6.1 节）。指数形式背后有一条可讲的物理直觉：片剂的强度来自颗粒之间有效键合点的数量与承力固相截面积。孔隙率每下降一点，单位断面上的固相承力面积与有效键合点都近似按比例增加，于是强度随密度下降呈近似指数上升——这就是σ = σ0·e^(−k·ε)的由来，也是σ0被解释为"零孔隙率（全密度）抗张强度"的原因。

这条直觉同时给出模型的适用边界：在中等—高固体分数段，承力网络已经形成，R–D 拟合通常很好；但在极低密度段（片剂接近散沙、承力结构尚未建立），R–D 会高估强度，因为纯指数外推假设"键合随密度连续增长"，而现实里低密度根本没有稳定的键合骨架。这正是 K–L 引入强度临界密度ρ_c,σ想补的洞——K–L 让曲线在ρ_c,σ处归零，而不是像 R–D 那样永远给正值。

配套地，USP⟨1062⟩ 还指出：许多药物材料在常见固体分数范围内，可用log(compression pressure) = a − b×(solid fraction)描述（a、b 为经验常数）。请务必分清：这条是"压力–固体分数"（可压性侧）的对数线性化，与 R–D 的"强度–固体分数"是两条不同的曲线，自变量与因变量都不是一回事。但在建模流水线里它们是紧邻的两步：先用前者从压力得到密度，再用后者从密度得到强度。

R–D 与 K–L 的选型可归纳成一张简表：

这张表的最后一行的取舍很关键：如果一篇报告要接受审计、且不允许出现"二次引用形式"的模糊，就优先用 R–D；K–L 更适合作为"解释低密度阈值"的补充模型，并明确标注来源等级。这是本系列"真实性优先于形式丰富度"的一贯取向。

1.7 参数标定的实操顺序（一份可复用的清单）

把"怎么把参数正确地拟合出来"固化成清单，能避免绝大多数返工：

1. 数据准备：固体分数统一为 0–1 小数，孔隙率ε = 1 − SF；抗张强度统一为 MPa（由第 12 篇的破碎力换算而来）。
2. 结构选择：先做线性化看趋势（ln(TS)对ε是否近似直线），判断单指数是否成立；若低密度段明显塌陷、出现"压不成型"阈值，再考虑 K–L 的带阈值形式。
3. 统一对数底：全文只用自然对数或只用 log10，并在报告里注明；任何跨底换算写显式因子ln10。
4. 初值与边界：用线性化结果反推p0；K–L 给ρ_c∈(0,1)的bounds；对指数形式设maxfev上限防止不收敛。
5. 误差与窗口：报告参数1σ（np.sqrt(np.diag(pcov))），并同时给出拟合所用的 SF 窗口——参数只在窗口内可信，窗口外是外推。
6. 分档只作标注：k_b的"易压/难压"只在报告里以"文献经验分档（阈值以原文为准）"出现，绝不当绝对判据。

这份清单会在第 17 篇（参数估计与标定）被完整工程化：那时我们会把curve_fit、加权拟合、AIC/R² 选模与误差传播串成流水线，并把标定结果回灌到TabletPress单元。本篇先埋下接口——你只要记住"统一底、给初值、报窗口、标来源"这十二个字，就能避开本篇绝大多数坑。

二、完整代码与逐行剖析

下面是一段可直接运行的脚本：先用三路拟合验证 R–D 三种写法的等价性（含对数底换算），再拟合 K–L 的两个参数，最后把两条 compactibility 曲线画在一张图上。

importnumpyasnpfromscipy.optimizeimportcurve_fitimportmatplotlib.pyplotasplt plt.rcParams["font.sans-serif"]=["SimHei","Microsoft YaHei","DejaVu Sans"]plt.rcParams["axes.unicode_minus"]=False# ============================================================# 0. 演示数据# 形态依据：USP⟨1062⟩"抗张强度随固体分数增大而近似指数增加"。# 具体数值以自己实测（第 10 篇 SF + 第 12 篇 TS）为准，此处仅作可运行演示。# ============================================================SF=np.array([0.70,0.75,0.80,0.85,0.90,0.95])# 固体分数 0–1TS=np.array([0.45,0.85,1.60,2.90,5.10,8.60])# 抗张强度 MPaeps=1.0-SF# 孔隙率 = 1 − 固体分数# ============================================================# 1. R–D 三种写法（函数形式）# ============================================================defrd_exp(eps,sigma0,k):"""σ = σ0·e^(−k·ε) 自然对数为底"""returnsigma0*np.exp(-k*eps)defrd_ln(eps,TS0,k_b):"""ln(TS/TS0) = −k_b·ε ⟹ TS = TS0·e^(−k_b·ε)"""returnTS0*np.exp(-k_b*eps)defrd_log10(SF,k,A):"""log10(TS) = k·SF + A ⟹ TS = 10^(k·SF + A)"""return10.0**(k*SF+A)# ---- 1a. 走"自然对数"路线：用 ln(TS) 对 ε 线性回归 ----slope_ln,intercept_ln=np.polyfit(eps,np.log(TS),1)# 斜率 = −k_bk_b=-slope_ln TS0=np.exp(intercept_ln)print(f"[自然对数] ln(TS)=−k_b·ε+lnTS0 → k_b={k_b:.3f}TS0={TS0:.3f}MPa")# ---- 1b. 走"log10"路线：log10(TS) 对 SF 线性回归 ----slope_l10,intercept_l10=np.polyfit(SF,np.log10(TS),1)k_log10,A_log10=slope_l10,intercept_l10print(f"[以 10 为底] log10(TS)={k_log10:.3f}·SF+{A_log10:.3f}")# ---- 1c. 等价性验证：混用对数底会差 ln10 倍 ----print(f"验证：k_log10 × ln10 ={k_log10*np.log(10):.3f}vs k_b ={k_b:.3f}")assertabs(k_log10*np.log(10)-k_b)<0.05,"对数底不一致！检查 log10 与 ln 混用"# ---- 1d. curve_fit 直接拟合指数形式（带初值与协方差） ----p0=(TS[0]/np.exp(-1.0*eps[0]),10.0)# 初值：σ0 取首点外推、k 取 10 量级popt_rd,pcov_rd=curve_fit(rd_exp,eps,TS,p0=p0,maxfev=20000)sigma0_fit,k_fit=popt_rd err_sigma0,err_k=np.sqrt(np.diag(pcov_rd))print(f"[curve_fit] σ0={sigma0_fit:.3f}±{err_sigma0:.3f}MPa, k={k_fit:.3f}±{err_k:.3f}")# 用 k_b 做"易压/难压"文献经验分档（阈值以原文为准）grade="易压（k_b<10）"ifk_fit<10else"难压（k_b>10）"print(f"文献经验分档（阈值以原文为准）：{grade}")# ============================================================# 2. Kuentz–Leuenberger 相对密度形式# 注：精确代数形式仅得二次引用（Su 等 2018 可检索形式），# 部分文献指数项含 γ 系数 —— 以原始论文为准。# ============================================================defkuentz_leuenberger(rho_r,sigma_max,rho_c):"""σ = σ_max·[1 − (1−ρ_r)/(1−ρ_c)·e^(ρ_r−ρ_c)]"""returnsigma_max*(1.0-(1.0-rho_r)/(1.0-rho_c)*np.exp(rho_r-rho_c))rho_r=SF.copy()# 相对密度 ≈ 固体分数（演示口径）p0_kl=(max(TS)*1.5,0.6)# 初值：σ_max 取上限外推、ρ_c 取偏低popt_kl,pcov_kl=curve_fit(kuentz_leuenberger,rho_r,TS,p0=p0_kl,bounds=([0.0,0.0],[np.inf,0.999]),maxfev=20000)sigma_max,rho_c=popt_kl err_smax,err_rhoc=np.sqrt(np.diag(pcov_kl))print(f"[K–L] σ_max={sigma_max:.3f}±{err_smax:.3f}MPa, ρ_c,σ={rho_c:.3f}±{err_rhoc:.3f}")print(f" ρ_r=ρ_c,σ 时 σ→0，即强度临界密度 ρ_c,σ={rho_c:.3f}")# ============================================================# 3. 绘制 compactibility 曲线# ============================================================sf_grid=np.linspace(SF.min(),1.0,200)eps_grid=1.0-sf_grid plt.figure(figsize=(7.2,5.0))plt.scatter(SF,TS,s=55,color="#c0392b",zorder=3,label="实测点（演示）")plt.plot(sf_grid,rd_exp(eps_grid,sigma0_fit,k_fit),lw=2,label=f"R–D: σ0={sigma0_fit:.2f}MPa, k={k_fit:.1f}")plt.plot(sf_grid,kuentz_leuenberger(sf_grid,sigma_max,rho_c),lw=2,ls="--",label=f"K–L: σ_max={sigma_max:.1f}MPa, ρ_c,σ={rho_c:.2f}")plt.axvline(rho_c,color="#7f8c8d",ls=":",lw=1.2)plt.text(rho_c,max(TS)*0.35," 强度临界密度 ρ_c,σ",color="#7f8c8d",fontsize=9)plt.xlabel("固体分数 SF");plt.ylabel("抗张强度 TS [MPa]")plt.title("compactibility 曲线：R–D 与 K–L（虚线为二次引用形式，以原始论文为准）")plt.legend(frameon=False);plt.tight_layout()plt.savefig("compactibility_rd_kl.png",dpi=150)plt.show()print("\n完成。R–D 三种写法等价（差对数底），K–L 多了强度临界密度。")

逐段剖析：

• 数据段：eps = 1 − SF是全文反复出现的关系（底账第 6.3 节：固体分数 = 表观密度/真密度 = 1 − 孔隙率）。SF用小数（0–1）而不是百分数，是本系列约定；若你手上数据是 85%，务必先除 100。
• 三路拟合：1a用np.polyfit(eps, ln(TS), 1)得到斜率−k_b与截距ln TS0；1b用log10(TS)对SF回归得k_log10。1c的assert abs(k_log10·ln10 − k_b) < 0.05是本篇最重要的"反直觉自检"——它把"三种写法等价"从口头结论变成一条会报警的断言。
• curve_fit 而非手工线性化：1d之所以再用curve_fit拟rd_exp，是因为线性化（先取对数再回归）会改变残差权重——大 TS 点的对数误差被压缩、小 TS 点被放大。直接对原始 TS 做非线性最小二乘能得到更符合直观的拟合，代价是需要p0初值。两者结论应接近，若差很多说明数据在该窗口不适合单指数。
• 误差输出：np.sqrt(np.diag(pcov))给出参数的 1σ 误差，这是第 17 篇参数标定的标准收尾；k_fit ± err_k直接决定"易压/难压"分档是否稳。
• K–L 的边界处理：curve_fit里用bounds=([0,0],[inf,0.999])把ρ_c,σ限制在 (0,1)，避免出现ρ_c≥1导致1−ρ_c≤0的除零/符号错乱（K–L 只在ρ_r>ρ_c,σ区间有物理意义）。注释里再次声明"二次引用形式，以原始论文为准"。
• 绘图：sf_grid一直画到 1.0，是为了让 R–D 与 K–L 在"全密度端"的差异显现——R–D 会单调爬升，K–L 则在 1.0 处正好收敛到σ_max。

反直觉默认值与坑（表格）

三、常见报错与排查

现象一：assert报"对数底不一致"或 k 与文献数值差 2.3 倍。
根因：一处用np.log（自然对数）、一处用np.log10（常用对数），两条直线的斜率天然差ln10≈2.3026。解法：在脚本开头统一约定"所有线性化一律用自然对数"，并在换算处显式写k_log10 * np.log(10)；本系列后续所有拟合都沿用此约定。

现象二：curve_fit抛RuntimeError: Optimal parameters not found或返回荒谬参数。
根因：指数形式对σ0、k高度非线性，p0给得太远（例如k=1而真值 20）会发散或陷入局部极小；数据跨度太窄（SF 只覆盖 0.90–0.95）时参数不可辨识。解法：p0用线性化结果反推（本脚本就是这么给的），必要时加bounds与maxfev，并检查 SF 跨度是否足够。

现象三：K–L 拟合给出ρ_c,σ > 1或负强度。
根因：1−ρ_c出现在分母，当ρ_c→1计算爆炸；且原式仅在ρ_r>ρ_c,σ时有物理意义，ρ_r<ρ_c,σ会算出负值。解法：bounds把ρ_c卡在 (0,1)，绘图/预测时用np.clip把自变量限制在拟合窗口内，不做远距离外推。

现象四：k_b=10被当成"合格线"直接判废。
根因：把文献经验阈值误当官方标准。真相：底账明确该阈值"属文献经验阈值，须以原文为准"，且k_b依赖拟合窗口。解法：只在报告里写"按某文献经验分档属易压/难压（阈值以原文为准）"，并同时列出拟合窗口与 1σ 误差；不得据此对未测物料下结论。

现象五：K–L 的σ_max与 R–D 的σ0当成同一个量。
根因：两者都可叫"零孔隙率抗张强度"，但 K–L 的σ_max是ρ_r=1端点的渐近值，R–D 的σ0是ε=0处的截距，二者在有限数据上外推结果不同。解法：报告里对同一批数据同时给出两个参数并注明定义域，不要相互替换。

四、动手练习

1. 对数底等价性自证：用本篇脚本，把1b的np.log10改成np.log，重跑。判定：改完后k_log10会变成原来的 2.3026 倍左右，assert因乘积不再等于k_b而失败——你要能解释"为什么换个底，斜率就变了"。
2. 一键分档：把演示数据整体乘 0.3（模拟"难压"物料），重跑并读取grade。判定：k_fit变化后能正确落入k_b>10或<10档，且脚本在报告里同时打印拟合窗口SF.min()–SF.max()与 1σ 误差，符合"经验阈值须标注"的纪律。
3. K–L 与 R–D 的分歧区：在 SF∈[0.55,0.75] 生成一组"低密度段明显压不成型"的数据（TS≈0），重跑两模型。判定：K–L 能在ρ_r→ρ_c,σ处把 TS 拉向 0，而 R–D 仍给正值；你能指出这正是"K–L 多一个强度临界密度"的收益，同时能说清该式来源属二次引用。

五、小结与下一篇预告

本篇把"成形性"这条曲线做成了可标定的模型：R–D 三种写法本质是同一指数关系的不同对数底与坐标，拟合时必须统一底，k_b只是文献经验阈值；K–L 相对密度形式用σ_max与强度临界密度ρ_c,σ补上了"低于某密度无强度"的物理图像，但其精确代数形式以原始论文为准。你还钉死了 compressibility / compactibility / tabletability 三条曲线的自变量差异——这是后面把"压力→密度→强度"串成链条的接口。

这段代码的输入是"固体分数 + 抗张强度"。固体分数来自第 10 篇的几何量计算，而抗张强度则要从破碎力换算而来——这正是下一篇的主题。

下一篇：第 12 篇《破碎力与抗张强度——USP⟨1217⟩ 单位换算、加载速率与径向破碎公式》——讲透"硬度"实为破碎力的误称、kp/N/SCU 的换算与 SCU 的不可靠性、加载速率 ≤20 N/s 或 ≤3.5 mm/s 的测试口径，以及径向破碎σ=2F/(πDt)与凸面片扩展式的可运行实现。

本篇认知问题回显（FAQ）

Q1：Ryshkewitch–Duckworth 方程的三种写法是不是同一个方程？为什么会有三种？
A1：是同一个方程的三种等价表达。σ=σ0·e^(−k·ε)（自然对数为底、自变量为孔隙率）、ln(TS/TS0)=−k_b·ε（同上，参数改名）、log(TS)=k·SF+A（常用对数、自变量换成固体分数 SF=1−ε）——它们只差对数底与坐标平移，故k_b≡k（同为自然对数时）。标定时必须统一对数底，否则 k 会相差 ln10≈2.3026 倍。

Q2：σ0、k、k_b、A 各代表什么？"k_b<10 易压 / >10 难压"能当常数用吗？
A2：σ0、TS0 是零孔隙率（全密度）抗张强度的外推值；k、k_b 是表征"键合能力"的指数衰减常数（二者等价）；A 是 log 形式的截距（与 k、σ0 相关）。"<10 易压、>10 难压"是文献经验阈值、须以原文为准，且 k_b 依赖拟合的固体分数窗口，不是普适常数，不能据此对未测物料下结论。

Q3：Kuentz–Leuenberger 方程里的 σ_max 与强度临界密度 ρ_c,σ 是什么？它比 R–D 多了什么？
A3：σ_t,p = σ_max·[1 − (1−ρ_r)/(1−ρ_c,σ)·e^(ρ_r−ρ_c,σ)]中，σ_max 是 ρ_r=1（全密度）时的抗张强度上限，ρ_c,σ 是强度临界密度。当 ρ_r=ρ_c,σ 时 σ=0（压不成型），ρ_r=1 时 σ=σ_max。相比纯指数的 R–D，K–L 多描述了"材料存在成形下限密度"这一特征。其精确代数形式仅得二次引用（部分文献指数含 γ 系数），以原始论文为准。

Q4：compressibility、compactibility、tabletability 在建模里各对应哪条曲线？
A4：可压性 compressibility = 固体分数–压缩压力（Heckel/Kawakita）；成形性 compactibility = 抗张强度–固体分数（孔隙率，即 R–D/K–L）；压片性 tabletability = 抗张强度–压缩压力。前两者自变量分别是压力与密度，后者是两者复合，三者不可混用（术语依 USP⟨1062⟩）。

Q5：给定一组"固体分数–抗张强度"数据，怎么稳健拟合并给出可信度？
A5：先按ε=1−SF转换自变量、统一用自然对数做线性化初值，再scipy.optimize.curve_fit对原式做非线性拟合；用np.sqrt(np.diag(pcov))取参数 1σ 误差，报告时同时给出拟合窗口。K–L 需加bounds（0<ρ_c<1）并只在 ρ_r>ρ_c,σ 区间解释。曲线来源若有二次引用，须注明"以原始论文为准"。