# -*- 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")