| 12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182 |
- # -*- coding: utf-8 -*-
- """涡激振荡 (VIV) 原理演示 —— 为什么停机比运行更危险。
- 工程经验模型: 涡脱频率 f_s = St*U/D; 锁频带; Scruton 数 Sc = 2*m*zeta/(rho*D^2)。
- **不是 CFD, 不是气弹耦合。** 参数为通用示意值, 用于说明机理, 不对应任何具体塔筒。"""
- import os, numpy as np, matplotlib
- matplotlib.use("Agg"); import matplotlib.pyplot as plt
- from matplotlib import font_manager as fm
- for c in ["PingFang SC","Hiragino Sans GB","Songti SC","STHeiti"]:
- try:
- fm.findfont(fm.FontProperties(family=c), fallback_to_default=False)
- matplotlib.rcParams["font.sans-serif"]=[c,"DejaVu Sans"]; break
- except Exception: pass
- matplotlib.rcParams["axes.unicode_minus"]=False
- BG="#edf7f9"; INK="#133047"; INK2="#28536b"; MUT="#61788b"; LINE="#cbe2eb"
- SPEC="#0e8898"; MEAS="#2077a8"; WARN="#d48818"; BAD="#d94b43"; OK="#23875b"
- St=0.20; D=4.0; f1=0.30; rho=1.225; m=3000.
- U=np.linspace(0,25,600); fs=St*U/D
- U_cr=f1*D/St
- LOCK=(0.85,1.30) # 锁频带 (经验区间)
- def sc(z): return 2*m*z/(rho*D**2)
- def amp(U,z):
- """锁频带内振幅抬升, 峰值 ~ C/Sc (经验形); 带外按 1/失谐 迅速衰减。"""
- r=U/U_cr; s=sc(z); pk=min(0.9, 1.0/max(s,0.35))
- band=np.exp(-((r-1.03)/0.16)**2)
- return pk*band
- Z_STOP, Z_RUN = 0.005, 0.05
- fig=plt.figure(figsize=(14.6,5.5), facecolor=BG)
- gs=fig.add_gridspec(1,3, left=.055,right=.985,top=.70,bottom=.14,wspace=.28)
- ax=fig.add_subplot(gs[0]); ax.set_facecolor("#fff")
- ax.plot(U,fs,color=MEAS,lw=2.4,label="涡脱频率 $f_s=St\\,U/D$")
- ax.axhline(f1,color=BAD,lw=2.0,ls="--",label=f"塔筒一阶频率 {f1:.2f} Hz")
- ax.axvspan(U_cr*LOCK[0],U_cr*LOCK[1],color=WARN,alpha=.16)
- ax.plot([U_cr],[f1],marker="o",ms=9,color=BAD,zorder=5)
- ax.annotate(f"临界风速 {U_cr:.1f} m/s",(U_cr,f1),xytext=(U_cr+3.2,f1*1.42),
- arrowprops=dict(arrowstyle="->",color=BAD,lw=1.4),color=BAD,fontsize=12,weight="bold")
- ax.set_xlim(0,25); ax.set_ylim(0,0.95)
- ax.set_xlabel("风速 (m/s)",color=MUT,fontsize=10.5); ax.set_ylabel("频率 (Hz)",color=MUT,fontsize=10.5)
- ax.set_title("① 涡脱频率撞上塔筒固有频率",color=INK,fontsize=13.5,loc="left",pad=10)
- ax.legend(loc="upper left",fontsize=10,frameon=False,labelcolor=INK2)
- ax2=fig.add_subplot(gs[1]); ax2.set_facecolor("#fff")
- a_stop=np.array([amp(u,Z_STOP) for u in U]); a_run=np.array([amp(u,Z_RUN) for u in U])
- ax2.axvspan(U_cr*LOCK[0],U_cr*LOCK[1],color=WARN,alpha=.16)
- ax2.fill_between(U,0,a_stop,color=BAD,alpha=.16)
- ax2.plot(U,a_stop,color=BAD,lw=2.6,label=f"停机 / 吊装期 阻尼 {Z_STOP*100:.1f}%")
- ax2.plot(U,a_run,color=OK,lw=2.4,label=f"正常运行 阻尼 {Z_RUN*100:.0f}%")
- ax2.set_xlim(0,25); ax2.set_ylim(0,max(a_stop)*1.25)
- ax2.set_xlabel("风速 (m/s)",color=MUT,fontsize=10.5)
- ax2.set_ylabel("横风向振幅 / 塔筒直径",color=MUT,fontsize=10.5)
- ax2.set_title("② 同样的风，差别全在阻尼",color=INK,fontsize=13.5,loc="left",pad=10)
- ax2.legend(loc="upper right",fontsize=10,frameon=False,labelcolor=INK2)
- ax2.annotate(f"相差约 {max(a_stop)/max(a_run):.0f} 倍",(U_cr,max(a_stop)*.55),
- color=BAD,fontsize=12.5,weight="bold",ha="center")
- ax3=fig.add_subplot(gs[2]); ax3.set_facecolor("#fff")
- zz=np.linspace(0.002,0.06,300); pk=[amp(U_cr*1.03,z) for z in zz]
- ax3.plot(sc(zz),pk,color=SPEC,lw=2.6)
- for z,lab,col in [(Z_STOP,"停机 / 吊装",BAD),(Z_RUN,"正常运行",OK)]:
- ax3.plot([sc(z)],[amp(U_cr*1.03,z)],marker="o",ms=10,color=col,zorder=5)
- ax3.annotate(lab,(sc(z),amp(U_cr*1.03,z)),xytext=(12,10),textcoords="offset points",
- color=col,fontsize=11.5,weight="bold")
- ax3.set_xlabel("Scruton 数 $Sc=2m\\zeta/(\\rho D^2)$ ← 越小越危险",color=MUT,fontsize=10.5)
- ax3.set_ylabel("锁频峰值振幅 / 直径",color=MUT,fontsize=10.5)
- ax3.set_title("③ 叶轮不转，气动阻尼就没了",color=INK,fontsize=13.5,loc="left",pad=10)
- for a in (ax,ax2,ax3):
- a.tick_params(colors=MUT,labelsize=9.5)
- for s in a.spines.values(): s.set_color(LINE)
- a.grid(color=LINE,lw=.7,alpha=.6)
- fig.text(.055,.935,"涡激振荡：为什么停着的机组比转着的更危险",color=INK,fontsize=20,weight="bold")
- fig.text(.055,.885,"风吹过塔筒会周期性脱落旋涡。脱落频率一旦撞上塔筒自身频率就进入「锁频」——"
- "振动不再随风速变化，而是自己把自己放大。",color=INK2,fontsize=12.5)
- fig.text(.055,.845,"关键在两点：临界风速通常落在很常见的风速段，不是极端天气；而最危险的时刻是吊装期和不带叶片的停机 —— "
- "叶轮转动提供的气动阻尼没了，只剩结构那点阻尼。",color=MUT,fontsize=11)
- fig.text(.055,.035,"工程经验模型（涡脱频率 + 锁频带 + Scruton 数），非 CFD、非气弹耦合；参数为通用示意值，"
- "不对应任何具体塔筒。用于说明机理与量级关系，不作设计依据。",color=MUT,fontsize=9.3)
- out=os.path.join(os.path.dirname(__file__),"assets","pain_viv.png")
- fig.savefig(out,dpi=150,facecolor=BG); print("SAVED",out)
- print(f"U_cr={U_cr:.2f} m/s | Sc停机={sc(Z_STOP):.2f} Sc运行={sc(Z_RUN):.2f} | 振幅比={max(a_stop)/max(a_run):.1f}x")
|
