guanlan-src/scripts/sim_hub/pain_pitch_imbalance.py

# -*- coding: utf-8 -*-
"""变桨/叶轮角度不一致 -> 气动不平衡 -> 1P 激励 -> 载荷抬升。
① 机理: 三叶片推力随方位角, 可正演(相对量)。
② 实测: 一个真实机群的塔筒 1P 筛查结果 (机群相对, 去共模 z 分数)。台号重编, 不含场名与厂商。
③ 后果: 电量与寿命 —— 运行数据无载荷/应变通道, 本模型不自算, 全部引文献。"""
import os, json, 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"

H,R,ALPHA,K = 90.,60.,0.20,0.060
psi=np.linspace(0,4*np.pi,900)
def blade_moment(dt):
    My=np.zeros_like(psi)
    for i in range(3):
        a=psi+i*2*np.pi/3; z=H+R*np.cos(a)
        My+=(np.maximum(z,1.)/H)**(2*ALPHA)*(1.0-K*(dt if i==0 else 0.))*np.cos(a)
    return My,(np.abs(np.fft.rfft(My[:600]))[1])/300.

M=json.load(open("/Users/yuanying/wind-analytics/outputs/blade_imbalance/mingyang_tower1p_screen.json",
                 encoding="utf-8"))
sc=np.array([t["score_dc"] for t in M]); tp=[t["type"] for t in M]
o=np.argsort(-sc)

fig=plt.figure(figsize=(14.6,5.9), facecolor=BG)
gs=fig.add_gridspec(1,3,left=.055,right=.985,top=.665,bottom=.165,wspace=.30)

ax=fig.add_subplot(gs[0]); ax.set_facecolor("#fff")
for d,lab,col in [(0.,"三叶片一致",OK),(2.,"一片偏 2°",BAD)]:
    My,_=blade_moment(d); ax.plot(psi/(2*np.pi),My,color=col,lw=2.6,label=lab)
ax.set_xlim(0,2); ax.set_xlabel("风轮转过的圈数",color=MUT,fontsize=10.5)
ax.set_ylabel("轮毂弯矩 (相对)",color=MUT,fontsize=10.5)
ax.set_title("① 机理：一片偏了，每转一圈晃一次",color=INK,fontsize=13,loc="left",pad=10)
ax.legend(loc="lower center",ncol=2,fontsize=10.5,frameon=False,labelcolor=INK2)

ax2=fig.add_subplot(gs[1]); ax2.set_facecolor("#fff")
CMAP={"Y轴lean":BAD,"Z轴lean":WARN,"混合":"#9fc4d4"}
ax2.bar(range(len(sc)),sc[o],color=[CMAP.get(tp[i],"#9fc4d4") for i in o],width=.84)
ax2.axhline(3,color=BAD,ls="--",lw=1.6); ax2.axhline(2,color=WARN,ls=":",lw=1.5)
ax2.text(len(sc)*.98,3.15,"z = 3",color=BAD,fontsize=10.5,ha="right")
ax2.text(len(sc)*.98,2.12,"z = 2",color=WARN,fontsize=10.5,ha="right")
ax2.set_xlabel(f"{len(sc)} 台机组（台号已重编）",color=MUT,fontsize=10)
ax2.set_ylabel("塔筒 1P 去共模 z 分数",color=MUT,fontsize=10.5)
ax2.set_title("② 实测：42 台里只挑出 2 台",color=INK,fontsize=13,loc="left",pad=10)
hs=[plt.Rectangle((0,0),1,1,color=v) for v in CMAP.values()]
ax2.legend(hs,[f"{k} {tp.count(k)} 台" for k in CMAP],fontsize=9.4,frameon=False,
           labelcolor=INK2,loc="upper right")
ax2.text(.30,.60,f"z≥3： {int((sc>=3).sum())} 台\nz≥2： {int((sc>=2).sum())} 台\n机群中位 {np.median(sc):+.2f}",
         transform=ax2.transAxes,fontsize=11,color=INK2,va="top")

ax3=fig.add_subplot(gs[2]); ax3.set_facecolor("#fff"); ax3.axis("off")
ax3.set_title("③ 后果：电量与寿命双损",color=INK,fontsize=13,loc="left",pad=10)
rows=[("年发电量损失","0.6 %","Δθ = 2° 时；最不利 1.2 %",BAD),
      ("剩余寿命损失","放大 2.9 – 3.2 倍","非线性放大；非旋转部件受影响最大",BAD),
      ("设计容差","± 0.3° / 极差 0.6°","GL 1989–2010 · DIBt 1993 · DS 472",OK),
      ("行业超标比例","35.3 %","1100+ 台激光实测筛 195 台，欧美 2013–2021",WARN)]
for i,(k,v,d,col) in enumerate(rows):
    y=.90-i*.235
    ax3.add_patch(plt.Rectangle((0.01,y-.185),0.97,.215,transform=ax3.transAxes,
                  facecolor="#f2fafc",edgecolor=LINE,lw=1))
    ax3.text(.05,y-.035,k,transform=ax3.transAxes,fontsize=11,color=INK2,va="top")
    ax3.text(.05,y-.078,v,transform=ax3.transAxes,fontsize=14.5,weight="bold",color=col,va="top")
    ax3.text(.05,y-.150,d,transform=ax3.transAxes,fontsize=9.0,color=MUT,va="top")
ax3.text(.01,-.055,"③ 引自 Saathoff 等 (2021)，Wind Energy Science 6:1079-1087。",
         transform=ax3.transAxes,fontsize=9.6,color=MUT)

for a in (ax,ax2):
    a.tick_params(colors=MUT,labelsize=9.4)
    for s in a.spines.values(): s.set_color(LINE)
    a.grid(color=LINE,lw=.7,alpha=.55)

fig.text(.055,.930,"变桨角度不一致：肉眼看不出，机器每转一圈晃一次",color=INK,fontsize=20,weight="bold")
fig.text(.055,.882,"三个叶片桨距不一致，气动推力就不再三重对称。合力矩多出一个「每转一圈来一次」的分量，"
                   "直接打在轮毂、主轴和塔架上。",color=INK2,fontsize=12.5)
fig.text(.055,.843,"偏差两三度，站在塔下看不出来 —— 但它同时吃电量和寿命，寿命那头还是非线性放大的。",
         color=MUT,fontsize=11)
fig.text(.055,.030,"① 为气动不平衡正演（切变 + 推力随桨距变化，相对量）；② 为真实机群塔筒 1P 筛查实测"
                   "（台号已重编，不含场名与厂商）；③ 运行数据无载荷与应变通道，本模型不出绝对载荷与寿命，数字全部引文献。",
         color=MUT,fontsize=9.3)
out=os.path.join(os.path.dirname(__file__),"assets","pain_pitch.png")
fig.savefig(out,dpi=150,facecolor=BG); print("SAVED",out)
print(f"真实: {len(sc)} 台 | z>=3: {int((sc>=3).sum())} | z>=2: {int((sc>=2).sum())} | 中位 {np.median(sc):+.2f} | 最大 {sc.max():.2f}")