直接压片单元模型开发教程4第一个压缩循环仿真——线性/方波/修正正弦三类位移-时间曲线版本声明块工具/软件Python numpy matplotlib纯 numpy 实现位移与速度曲线生成。语言/环境Python 3.9numpy2、scipy1.14、matplotlib。本文目标写出本系列的第一个真仿真函数make_displacement_profile把冲头位移从一句描述变成可计算、可对比、可画速度剖面的数组。一句话结论压缩模拟中冲头位移-时间曲线分三类——线性saw-tooth、方波square、旋转机特有的修正正弦modified sinusoidalUSP⟨1062⟩ 转述用make_displacement_profile(kind, rpm, points)生成后可见方波在平台段速度趋零但边沿加速度发散、线性曲线全程无零速点、修正正弦速度平滑过零三者本质上是在驻留时间与加速度代价之间做不同取舍旋转机位移曲线起始斜率大、最大位移附近斜率趋 0使最大压力附近单位时间体积变化缓慢利于排气与成型。〇、本篇要解决的认知问题前面两篇把机器讲清了现在到了第一次真正算数的时刻。四个问题为什么整个仿真要从位移-时间曲线开始没有它后面哪一步都做不了吗线性saw-tooth、方波square、修正正弦modified sinusoidal三类曲线差别只是形状好看吗它们各自隐含什么物理假设为什么说修正正弦是旋转压片机的天然形状它的起始斜率大、峰值斜率→0是怎么和排气挂上钩的压缩窗口到底该取多长为什么说它是毫秒级而不是一个转台转一圈三类曲线的取舍标准是什么——要不要驻留要驻留就得付出什么代价一、机制解析1.1 位移曲线为什么是第一块砖一句话位移是一切的下游源头。回顾第 1 篇的骨架压力 → 固体分数 → 强度——你可能会问压力不是直接给定的吗不是。压片机调的是力而压辊把力施加到粉床上的时机与速率是由冲头位移曲线决定的位移曲线 y(t) │ ├─ 一阶导 d(Δy)/dt → 冲头瞬时速度 v(t) → 应变速率 → 应变速率敏感性 SRS第 18 篇 │ ├─ 峰值附近平的时长 → 驻留时间 dwell第 7 篇 │ └─ 位移 → 体积变化 → 压缩压力 P(t) → 固体分数 SF → 抗张强度 → 破碎力/CQAUSP⟨1245⟩ 把压缩模拟列为药典级的表征手段正是因为位移-时间是可定义、可复现、可跨设备比较的归一化语言第 7 篇会看到 dwell 如何成为放大归一化参数。所以第一步必须先能生成位移曲线。1.2 三类曲线三种物理假设USP⟨1062⟩本系列转述指出压缩模拟中的冲头位移-时间曲线可分三类线性saw-tooth、方波square、旋转机特有的修正正弦modified sinusoidal并可按单侧或双侧压缩细分。三者不是画法风格的差异而是三种不同的理想化曲线形状隐含的物理假设主要用途线性 saw-tooth等速升降的三角或陡降的锯齿冲头速度恒定升到最大位移后立即反向最简单基准用来看完全没有驻留会怎样方波 square迅速到位 → 保持 → 迅速退回理想凸轮瞬时到位、平台保持上界情形看驻留拉满的极限修正正弦 modified sinusoidal平滑的类正弦过渡压辊滚压产生的连续、有限加速度运动旋转压片机的真实近似先记住一句警告“修正正弦近似是工程上常见做法但不得当作唯一标准。真正的旋转机压辊位移一般基于 Rippie 和 Danielson 提出的方程计算USP⟨1245⟩而该方程的显式表达式本次调研未获原文本系列不给出其形式——凡需要显式的场合一律以 Rippie Danielson 原文或 USP⟨1245⟩ 所引公式为准。本篇的三类曲线是可计算的参照系”不是标准答案。1.3 为什么修正正弦是旋转机的天然形状第 3 篇已经给出关键观察伍善根、沈震曲线斜率 d(Δy)/dt 即冲头瞬时速度旋转压片机位移曲线起始斜率大冲头刚接触压辊时速度大在最大位移 Δy 附近斜率趋近 0而单冲机在最大压力附近斜率仍较大气体来不及排出。现在把这条观察翻译成该选哪条曲线方波确实在平台段速度为零——这最符合峰值附近速度趋于零但它的边沿速度无限大、加速度发散。真实机械承受不了无限加速度所以方波只能是上界参照。线性曲线全程速度恒定、没有零速点对应最大压力附近单位时间体积变化快、气体来不及排出那一类——更适合描述单冲式行为。修正正弦半正弦过渡是两者的折中起点处斜率高对应起始速度快峰值附近速度平滑穿过 0对应最大位移附近斜率→0且加速度有限真实机械可实现。这就是它被称为旋转机天然形状的原因——它同时满足了第 3 篇里那条起始快、峰值慢的双条件又没有方波的加速度发散。这三者的取舍可以写成一张权衡表曲线峰值处速度边沿加速度驻留物理可实现性线性 saw-tooth不为 0反向跳变有限但有跳变无中加速度不连续方波 square严格为 0平台发散理想化最大平台整体低理想化上界修正正弦平滑过 0仅峰值瞬时有限且连续峰值处瞬时过零、占比低高旋转机近似这张表就是本篇的核心结论。它解释了为什么后续所有动力学量驻留时间、应变速率、SRS都必须挂在同一条位移曲线上算而不能各算各的。需要额外点破一层有驻留和驻留长是两件事。方波是驻留长的极端线性是完全没有驻留修正正弦是有一个平滑的、瞬时性的驻留。真实压片机上人们说的驻留时间是毫秒级的一个真实时长第 7 篇给出 7–33 ms 这类实测值它来自压辊与冲头头部平面接触这一几何事实而不是来自曲线形状本身。本篇用一个|v| ≤ 5%·峰值的粗阈值去量化它只是为了让你在三条曲线之间建立谁更接近驻留的直觉——它不能替代第 7 篇的精确公式。1.4 单侧与双侧压缩USP⟨1062⟩ 转述指出压缩模拟可分单侧/双侧。这个区分直接影响位移曲线的振幅分配单侧压缩只有一侧冲头运动另一侧固定。总压缩位移 Δy 全部由运动侧承担。双侧压缩两侧相向运动。这正是第 3 篇里旋转压片机的常态——上辊驱动上冲、下辊驱动下冲压缩阶段冲头运动对称USP⟨1245⟩。注意对称只描述运动学两侧位移对称不等于两侧力对称力传递比 R 最大下冲力/最大上冲力 常小于 1第 14 篇详述。因此在生成位移曲线时“双侧意味着同一条曲线被复制到两侧而力的分配要另开一套模型。本篇的make_displacement_profile生成的是一条单侧等效位移”双侧只是把它镜像到另一侧。1.5 压缩窗口为什么是毫秒一个必须算清的账新手最容易犯的错是把一个压缩循环当成转台转一圈。这两者的时长差了三个数量级。转台转一圈的周期是T_rev 60 / rpm。以 50 min⁻¹ 为例T_rev 1.2 s。而冲头真正处于压缩状态被压辊压下、粉床被压缩的时间是压辊作用的那一小段圆弧所对应的时间——用第 7 篇的口径说这是驻留时间在毫秒量级的来源Natoli 白皮书、PMC 论文给出的实测 dwell 在 7–33 ms 这个量级。做一个简单的量级校核取压缩事件窗口t_comp 10 ms、总压缩位移Δy 3 mm则修正正弦的峰值速度约为π·Δy / t_comp π×3 mm / 10 ms ≈ 942 mm/s ≈ 0.94 m/s。而 USP⟨1245⟩ 给出的实机数据是大型高速旋转压片机冲头垂直速度可达约 1 m/s。两者在同一量级。这说明只有把压缩窗口取在毫秒尺度约 1 m/s这个实机数字才会自然出现如果你把压缩窗口取成 1.2 s峰值速度只有 7.85 mm/s与实机差了两个数量级。所以本篇的函数里t_comp_ms是一个必须显式给出的参数而rpm的作用是让你看到压缩窗口只占一转的极小比例这个事实。这个比例本身就是一个很有冲击力的数字10 ms / 1.2 s ≈0.83%。也就是说压片机每转一圈真正在压的只有不到百分之一的时间。剩下 99% 的时间在干什么在充填、在输送、在脱模、在等下一个循环。这句话值得反复念。关于压缩窗口的精确取值从哪里来它由压辊几何辊径、位置、转台 PCD 与转速共同决定其位移形式回到 Rippie Danielson 方程。具体取值与显式方程以官方文档与原始文献为准本篇用毫秒量级的窗口 修正正弦形状作为工程近似并在代码中把该假设写在注释与打印里。二、完整代码与逐行剖析下面这份脚本实现make_displacement_profile(kind, rpm, points)并画出三类曲线的位移剖面与速度剖面。# -*- coding: utf-8 -*- 第 4 篇配套脚本生成三类冲头位移-时间曲线并对比速度剖面。 三类曲线依据 USP1062 转述 - linear : 线性 saw-tooth本实现对采用等速升降的三角近似 - square : 方波迅速到位 - 保持 - 迅速退回边沿用有限斜坡近似 - sinusoidal : 修正正弦本实现采用半正弦近似常见工程做法非唯一标准 重要声明 1) 旋转机真实压辊位移一般基于 Rippie 和 Danielson 提出的方程USP1245 其显式表达式未获原文本脚本不作复现一律以原文与 USP1245 所引公式为准。 2) 修正正弦仅为常见近似不得当作唯一标准。 importnumpyasnpimportmatplotlib.pyplotaspltdefmake_displacement_profile(kind,rpm,points2000,delta_y_mm3.0,t_comp_ms10.0): 生成一条压缩事件的冲头位移-时间曲线及其速度剖面。 参数 ---- kind : str linear | square | sinusoidal。 rpm : float 转台转速min^-1。用于计算“压缩窗口占一转的比例”。 points : int 采样点数。 delta_y_mm : float 最大压缩位移 Δymm示意值真实值由压辊几何决定。 t_comp_ms : float 一个压缩事件的时间窗口ms毫秒量级与 dwell 的 ms 尺度一致。 返回 ---- dict: {t_ms, y_mm, v_mm_s, t_rev_s, duty, kind} ifkindnotin(linear,square,sinusoidal):raiseValueError(kind 必须是 linear / square / sinusoidal 之一)unp.linspace(0.0,1.0,points)# 归一化相位 0..1dfloat(delta_y_mm)ifkindlinear:# 等速升至峰值u0.5再等速返回全程速度恒定峰值处速度反向跳变无零速点ynp.where(u0.5,2.0*d*u,2.0*d*(1.0-u))elifkindsquare:# 有限斜坡近似方波平台段速度为零对应最大驻留边沿加速度较大edge0.04xpnp.array([0.0,0.25-edge/2,0.25edge/2,0.75-edge/2,0.75edge/2,1.0])fpnp.array([0.0,0.0,d,d,0.0,0.0])ynp.interp(u,xp,fp)else:# sinusoidal修正正弦的半正弦近似# 起点斜率最大起始速度快峰值处速度平滑穿过 0最大位移附近斜率-0yd*np.sin(np.pi*u)t_comp_st_comp_ms/1000.0# 压缩窗口st_su*t_comp_s v_mm_snp.gradient(y,t_s)# 速度剖面mm/st_rev_s60.0/rpm# 一转周期sdutyt_comp_s/t_rev_s# 压缩窗口占一转比例return{t_ms:t_s*1000.0,y_mm:y,v_mm_s:v_mm_s,t_rev_s:t_rev_s,duty:duty,kind:kind}defdwell_fraction(prof,thresh_ratio0.05):计算“准驻留时间占比”|v| 小于 (thresh_ratio * 峰值|v|) 的时间占比。vnp.abs(prof[v_mm_s])v_peakv.max()ifv_peak0:return1.0maskvthresh_ratio*v_peakreturnfloat(np.mean(mask))defplot_profiles(profiles,out_pngdisplacement_profiles.png):绘制三类曲线的位移剖面与速度剖面对比图。plt.rcParams[font.sans-serif][SimHei,Microsoft YaHei,DejaVu Sans]plt.rcParams[axes.unicode_minus]Falsefig,(ax1,ax2)plt.subplots(1,2,figsize(13,5.0))colors{linear:#1f4e79,square:#c00000,sinusoidal:#2e8b57}forprofinprofiles:kprof[kind]ax1.plot(prof[t_ms],prof[y_mm],colorcolors[k],lw2.0,labelk)ax2.plot(prof[t_ms],prof[v_mm_s],colorcolors[k],lw2.0,labelk)ax1.set_xlabel(时间 t (ms))ax1.set_ylabel(冲头压入位移 Δy (mm))ax1.set_title(三类位移-时间曲线Δy3 mm压缩窗口 10 ms示意)ax1.legend(fontsize9)ax1.grid(alpha0.3)ax2.set_xlabel(时间 t (ms))ax2.set_ylabel(冲头速度 v d(Δy)/dt (mm/s))ax2.set_title(速度剖面斜率即瞬时速度)ax2.legend(fontsize9)ax2.grid(alpha0.3)ax2.annotate(线性的峰值处速度反向跳变\n无零速点,xy(5.0,0.0),xytext(2.2,320),arrowpropsdict(arrowstyle-,color#1f4e79,lw1.1),fontsize8,color#1f4e79)ax2.annotate(修正正弦平滑过 0\n峰值斜率→0,xy(5.0,0.0),xytext(5.2,-520),arrowpropsdict(arrowstyle-,color#2e8b57,lw1.1),fontsize8,color#2e8b57)fig.tight_layout()fig.savefig(out_png,dpi140)plt.show()returnout_pngdefmain():rpm50.0kinds[linear,square,sinusoidal]profiles[make_displacement_profile(k,rpmrpm,points4000)forkinkinds]print(f转速 rpm {rpm:.0f}min^-1 - 一转周期 T_rev {profiles[0][t_rev_s]:.3f}s)print(f压缩窗口 10 ms仅占一转的{profiles[0][duty]*100:.2f}%\n)print(f{曲线:12}{峰值|v| (mm/s):18}{峰值|v| (m/s):15}{准驻留占比:12})print(-*57)forprofinprofiles:vpnp.abs(prof[v_mm_s]).max()print(f{prof[kind]:12}{vp:18.1f}{vp/1000.0:15.3f}{dwell_fraction(prof):12.3f})print(\n注修正正弦峰值速度约 0.94 m/s与 USP1245 所述实机约 1 m/s 同量级)print( 该一致性来自“压缩窗口取毫秒量级”这一设定而非曲线本身。)outplot_profiles(profiles)print(f图已保存:{out})if__name____main__:main()逐段剖析u np.linspace(0, 1, points)而不是直接构造时间数组。先把三类曲线都定义在归一化相位u上最后一次性映射到t u·t_comp。这样做的好处是修改压缩窗口t_comp_ms时曲线形状不变只有时间轴伸缩——这正是形状与时长解耦的工程写法。linear用np.where构造三角。2d·u在u0.5时从 0 升到 d2d·(1−u)在u≥0.5时从 d 降到 0。速度恒为±2d/t_comp 0.6 m/s在峰值处从 600 跳到 −600这是线性曲线无零速点的数学体现。文献里也叫 saw-tooth锯齿另一种写法是缓慢升、瞬时降的锯齿本实现取等速升降的三角因为它更能说明没有零速点这个要害。square用np.interp配有限斜坡。如果写成np.where(u0.25, d, 0)边沿会得到d/dt量级的速度——在 10 ms、4000 点的采样下可达上万 m/s是个纯粹的数值伪影。所以这里用edge0.04的斜坡把边沿软化成有限速度。这个细节本身就是一条教学结论理想方波在物理上不可实现它的加速度发散。sinusoidal用半正弦d·sin(πu)。关键在于它的导数π·d·cos(πu)/t_comp。在u0起点处导数最大 π·d/t_comp起始速度大在u0.5最大位移处导数为 0峰值斜率→0。一条公式同时满足了第 3 篇那条起始快、峰值慢的双条件——这就是它被称为旋转机天然形状的根本原因。再次提醒这是常见近似不是唯一标准。np.gradient(y, t_s)求速度。注意第二个参数传的是真实时间步秒不是索引。若误传索引得到的速度单位就变成了 mm/步数值会因采样点数不同而变——这是数值微分的经典陷阱第 6 篇会专门处理平滑与差分。dwell_fraction是准驻留的粗定义。它统计|v| ≤ 5%·峰值的时间占比。之所以叫准是因为真正的驻留时间定义第 7 篇是压辊与冲头头部平面接触的时间与峰值附近速度小有关但不等价。把它标成准驻留是为了不越界。main()里那句该一致性来自设定而非曲线是刻意的。0.94 m/s 与 USP⟨1245⟩ 的约 1 m/s 同量级是压缩窗口取毫秒量级这一设定的产物如果你把窗口改成 1.2 s修正正弦的峰值速度会掉到 7.85 mm/s与实机差两个数量级。不要把这个巧合当成模型验证通过。反直觉的默认值提醒np.where(u 0.5, ...)在u恰为 0.5 时走的是else分支结果2d×(1−0.5)d与if分支一致所以不产生跳变。但如果你把阈值改成非对称值例如u 0.4峰值位置与平台位置就会错位速度剖面上会出现一个哑平台。做曲线时分段函数的接缝处一定要打印被定义域余数捕获的临界点。三、常见报错与排查1现象方波的速度剖面出现上万 m/s 的尖峰图像被压扁成一条线。根因理想方波边沿的速度是Δy/Δt当Δt由采样步长决定时尖峰大小随采样点数变化是数值伪影而非物理速度。解法用有限斜坡本篇edge0.04软化边沿并在图上标注方波边沿为理想化、此处用有限斜坡近似或者干脆不把方波放进同一张速度图改用准驻留占比这类标量指标对比。2现象np.gradient算出的速度数量级离谱例如比预期大 1000 倍。根因把np.gradient(y)直接调用第二参数缺省为 1于是速度实际上是每步的位移。解法必须写成np.gradient(y, t_s)且t_s用秒。单位核对mm / s mm/s。3现象三类曲线的峰值速度算出来都远低于 1 m/s。根因压缩窗口取成了一转周期而不是毫秒级压缩事件窗口。解法显式传入t_comp_ms并用duty t_comp / T_rev自查——在 rpm50 时该比例约 0.83%。若你算出 duty 接近 100%说明窗口取错了。4现象直接套用某篇文献给出的修正正弦公式结果与自己的数据对不上怀疑公式错。根因“修正正弦是常见近似不同来源的具体形式与参数并不统一旋转机的权威来源是 Rippie Danielson 方程而其显式表达式本次调研未获原文。解法不要声称自己的公式是标准式”。凡出现具体公式标注以 Rippie Danielson 原文或 USP⟨1245⟩ 所引公式为准把所用近似明确写成本系列采用的工程近似。5现象把双侧压缩写成两条独立的位移曲线。根因混淆了运动学与力学。双侧指两侧相向、位移对称USP⟨1245⟩上辊驱动上冲、下辊驱动下冲但两侧力不一定相等。解法双侧用同一条位移曲线镜像到两侧表达力的分配另开模型力传递比 R第 14 篇。四、动手练习练习 1基础跑通三类曲线并核对量级。保存脚本为displacement_profiles.py运行。成功标准量化终端打印的duty为0.83%rpm50修正正弦峰值 |v| 落在930–955 mm/s即 0.93–0.955 m/s方波准驻留占比最高、线性最低。练习 2进阶做一次转速扫描并画 dwell-转速曲线。对 rpm ∈ {20, 30, 50, 70, 100} 计算duty并确认duty与 rpm 成正比。成功标准量化输出的duty%随 rpm 单调递增且duty × rpm近似为常数各点相对偏差 1e-6用一句话解释为什么转速越高压缩窗口占一转的比例越大提示窗口时长固定一转时长随转速缩短。练习 3挑战给三类曲线各加一个形状参数并做权衡分析。给square加斜坡宽度参数edge∈{0.01,0.04,0.10}给sinusoidal加形变指数p用y d·(sin(πu))**pp∈{1,2,3}。成功标准量化输出一张 3×3 的表格列为峰值 |v|“与准驻留占比”要求edge从 0.01 增到 0.10 时方波峰值 |v| 单调下降p从 1 增到 3 时修正正弦准驻留占比单调上升。并在结论里明确指出驻留变好总是伴随峰值速度或加速度代价。五、小结与下一篇预告本篇完成了本系列的第一次真仿真make_displacement_profile(kind, rpm, points)生成三类位移曲线np.gradient(y, t_s)给出速度剖面。三类曲线的本质差异不在形状而在驻留与加速度的取舍——线性曲线全程无零速点、方波平台驻留最大但边沿加速度发散、修正正弦平滑过零且加速度有限因而最贴近旋转机的真实行为。请务必记住两条红线修正正弦是常见近似、不得当作唯一标准Rippie Danielson 位移方程的显式形式以原文与 USP⟨1245⟩ 所引公式为准本系列不作复现。还有一条最容易忽略的量级结论压缩窗口是毫秒级只占一转的不到 1%rpm50、10 ms 时约 0.83%正是这个设定让修正正弦的峰值速度约 0.94 m/s与 USP⟨1245⟩ 所述实机的约 1 m/s 落在同一量级。这不是验证是量级自洽。下一篇第 5 篇把今天生成的位移曲线装进一个可复用的容器定义自研的TabletPress单元类本系列二次开发示例类非官方 API把冲头直径、冲数、PCD、转速、预压/主压设定封装为构造入参用step()输出结果对象并演示如何把它串进一条最小 flowsheet 拓扑。第 6 篇将把本篇的数值微分做扎实速度/加速度与平滑第 7 篇将给出驻留时间的精确定义与四套公式第 13 篇详述双辊预压/主压第 14 篇详述力传递比。本篇认知问题回显FAQQ1为什么整个仿真要从位移-时间曲线开始A因为位移是下游一切的源头其一阶导给出冲头瞬时速度斜率 d(Δy)/dt 即瞬时速度峰值附近的平给出驻留时间位移经体积变化给出压缩压力进而决定固体分数与强度USP⟨1062⟩ 也将位移-时间曲线作为压缩模拟的基础表征。Q2线性、方波、修正正弦三类曲线的差别只是形状吗A不是。线性saw-tooth全程速度恒定、峰值处反向跳变、无零速点方波平台段速度为零、驻留最大但边沿加速度发散理想化修正正弦起始斜率大、峰值平滑过零、加速度有限是最贴近旋转机的近似。Q3为什么修正正弦是旋转机的天然形状它怎么和排气挂钩A旋转机位移曲线起始斜率大、在最大位移附近斜率趋近 0伍善根、沈震——半正弦恰好同时满足这两个条件斜率趋零意味着最大压力附近单位时间体积变化缓慢为排气与成型留出时间窗单冲机峰值附近斜率仍大气体来不及排出。Q4压缩窗口该取多长为什么是毫秒级A压缩窗口对应压辊作用的那一小段圆弧与 dwell 同属毫秒量级实测 dwell 如 Natoli NP-500 为 17 ms、Fette 1200i 为 7–33 ms取 10 ms、Δy3 mm 时修正正弦峰值速度约 0.94 m/s与 USP⟨1245⟩ 所述实机约 1 m/s 同量级取成一转周期rpm50 时 1.2 s则峰值速度仅 7.85 mm/s差两个数量级。Q5三类曲线的取舍标准是什么A核心是驻留 vs 加速度代价的权衡——方波驻留最大但加速度发散线性无驻留修正正弦居中且加速度有限此外修正正弦只是常见工程近似不得当作唯一标准权威形式应回到 Rippie Danielson 原文与 USP⟨1245⟩ 所引公式。