# 斜拉桥三个教学模型：公式、接口与验证 本目录的 `mechanics.js` 是无外部依赖的 ES module。模型用于解释索角、调索耦合和多跨传力；不是某座真实桥的设计验算，也不能据此决定现场张拉量。 统一单位：力 kN，长度和位移 m，弯矩 kN·m，刚度 kN/m，抗弯刚度 EI 为 kN·m²，角度计算用 rad。图形若放大变形，必须注明放大倍数或“变形示意”。 ## 1. Harp、fan、semi-fan：固定比较条件 `cableLayout(options)` 展示同一跨径的右半跨。塔在 x=0，梁面 z=0，总跨径 `span=120`，塔高 `towerHeight=40`，默认每半跨十根索，梁锚点 x=6、12、…、60 m。`count` 可改为 6、10、14 等正整数。每根索默认承担相同竖向分力 V=1000 kN，轴向刚度参数 EA=2×10⁶ kN。左右两侧可分别调用并独立改变 V，再把两侧结果传给 `towerResponse`。 塔端锚点高程： - harp：zᵢ=h·xᵢ/(L/2)，因此各索平行。 - fan：zᵢ=h，所有索锚于理想化塔顶一点。 - semi-fan：zᵢ=h[0.55+0.45·xᵢ/(L/2)]，锚点分布在塔上部。 计算 θ=atan2(z,x)，索长 ℓ=√(x²+z²)，并逐索输出： ```text T = V / sinθ H = V cotθ k_v,elastic = (EA / ℓ) sin²θ ``` 这里的 `verticalStiffness` 只包含索轴向伸长投影形成的弹性刚度，塔端固定、梁端水平位移固定，忽略初始拉力的几何刚度、索垂度、梁塔变形和锚固区细节。不能把它称作全桥竖向刚度。高塔锚点既改变 θ 也改变 ℓ，因此“fan 的每根索弹性竖向刚度都比 harp 大”并不是这些公式保证的结论。 等 V、等梁锚点且不降低塔锚点时，θ 增大使 T 与 H 减小。`totals.T` 只是各索力大小的算术和，不能解释为一个合力；实际合力应按方向相加。右半跨的水平分力大小 H 指向塔。 返回值：`{type,span,towerHeight,count,verticalForce,EA,cables,totals}`；每条 `cables` 含 `{index,x,z,length,angle,angleDeg,V,H,T,verticalStiffness}`。 ## 2. 索力优化：塔位支承与无塔支承对照 基础梁参数为 L=120 m，EI=8×10⁸ kN·m²，q=80 kN/m，塔位 x=60 m。各索竖向分力 Vᵢ 是可调的向上集中力。为保留旧对照，`DEFAULT_MODEL` 仍是四索、`DENSE_MODEL` 仍是十二索，二者 API 默认 `towerSupport:false`。**当前交互采用 `towerSupport:true`：塔位竖向支承主梁，主梁在该处连续，转角不被锁定，也没有插入梁内铰。** 可直接使用 `SUPPORTED_DENSE_MODEL`，或 `{...DEFAULT_MODEL,towerSupport:true}`。 塔的侧移、斜索伸长和初始无应力长度、轴弯耦合、施工阶段、徐变收缩、二阶效应均不进入梁的响应求解。塔的竖向支承视为不沉降的理想约束；另算的塔侧移不反馈给索力。 ### 当前塔位支承：用力法保证 w(a)=0 取 a=`towerX`，先在没有塔位竖向约束的简支基系上求当前索力作用下的 `w_free(x)`。f(x,a) 为该简支基系的单位向下力柔度。塔对梁的支反力向上为正： ```text R塔 = w_free(a) / f(a,a) w(x) = w_free(x) − R塔 f(x,a) M(x) = M_free(x) − R塔 m(x,a) R左 = R左,free − R塔(L−a)/L R右 = R右,free − R塔 a/L ``` m(x,a)=x(L−a)/L−max(0,x−a) 为单位向下力弯矩影响。塔位弯矩通常不为零；这里没有把连续梁截断成两根独立简支梁。 优化使用消去 R塔 后的响应，不能先按无塔支承优化再只在显示时把塔位挠度强制设零。单位索托力的影响函数变为 ```text f_supported(x,b) = f(x,b) − f(a,b) f(x,a) / f(a,a) m_supported(x,b) = m(x,b) − f(a,b) m(x,a) / f(a,a) ``` 均布无索基线也按同一支承条件求解：R塔,0=w_q(a)/f(a,a)，w₀=w_q−R塔,0 f(x,a)，M₀=M_q−R塔,0 m(x,a)。目标归一化使用这组 w₀、M₀；积分区间额外在塔位分段。导出的 `pointLoadFlexibility` 和 `pointLoadMoment` 仍代表简支基系，支承修正由 `analyzeBeam` 内部施加。 分担宽度初值把塔点纳入相邻支点：四索为 `[1920,1440,1440,1920]` kN；十二索为 `[640,640,640,640,640,800,800,640,640,640,640,640]` kN。这些只是优化起点，R塔 仍由位移协调计算。 当前120 m塔位支承模型复算表（最大值为121点采样）： | 索数与状态 | 左端反力 / kN | 塔反力 / kN | 右端反力 / kN | 挠度 RMS / mm | 弯曲能 / kN·m | 最大 \|M\| / kN·m | 最大 \|w\| / mm | |---|---:|---:|---:|---:|---:|---:|---:| | 四索初值 | 889.920000 | 1100.160000 | 889.920000 | 0.663332689 | 0.636263424 | 4949.120000 | 0.926625567 | | 四索挠度目标 | 721.566557 | 698.489517 | 721.566557 | 0.067931300 | 0.328846577 | 5722.402636 | 0.136819237 | | 四索弯曲能目标 | 747.000805 | 552.018842 | 747.000805 | 0.095499682 | 0.320371200 | 5111.980684 | 0.206445958 | | 十二索初值 | 304.533333 | 990.933333 | 304.533333 | 0.016479238 | 0.008166400 | 928.000000 | 0.030634667 | | 十二索挠度目标 | 255.149071 | 997.849741 | 255.149071 | 0.001919842 | 0.005980979 | 1013.517031 | 0.005752514 | | 十二索弯曲能目标 | 252.469161 | 1018.079951 | 252.469161 | 0.002140739 | 0.005939364 | 1076.170347 | 0.006707034 | 四索挠度目标为外侧各2312.762970、内侧各1416.425714 kN；弯曲能目标为外侧各2237.999382、内侧各1538.990392 kN。十二索挠度目标左半侧约为 `[715.649,630.033,639.494,650.988,571.969,837.793]`，弯曲能目标约为 `[725.328,618.307,641.788,654.494,580.005,818.569]`，右半侧镜像。两类优化均已收敛；目标下降不等于弯矩峰值下降，本表的峰值 M 即可作反例。 无索时，两等跨连续梁的独立闭式检验为 R塔=5qL/8=6000 kN、R左=R右=3qL/16=1800 kN、M塔=−qL²/32=−36000 kN·m。总反力9600 kN，塔位位移为零。若将这个 R塔 作为桥面处向下力交给塔，塔底轴压为6000 kN。 ### 无塔位支承对照的响应公式 以梁挠度向下、正弯矩为正： ```text R_A = qL/2 − Σ Vᵢ(L−aᵢ)/L R_B = qL/2 − Σ Vᵢaᵢ/L M(x) = R_A x − qx²/2 + Σ Vᵢ max(0,x−aᵢ) w₀(x) = qx(L³−2Lx²+x³)/(24EI) w(x) = w₀(x) − Σ f(x,aᵢ)Vᵢ ``` 单位向下集中力的挠度影响函数 f，以 b=L−a 记： ```text x≤a: f(x,a)= b x [L²−b²−x²]/(6LEI) x>a: f(x,a)= a(L−x)[L²−a²−(L−x)²]/(6LEI) ``` 调一根索会改变整根梁的 M(x)、w(x)，也会改变其他索作用点的梁位移。此处的 V 是外力控制变量，所以“其他索作用点位移改变”不等于已模拟真实被动索的张力重分配；真实张力重分配还须联立索、塔、梁变形协调。 ### 两个目标与边界 无塔位支承的四索初值 `equalTributaryForces()` 按相邻支点中线间的分担宽度取 Vᵢ=q·24=1920 kN；启用塔支承时按前节修正分担宽度。初值不是最优解。 `objective:'deflection'` 为默认目标： ```text J_w = [∫w² dx / ∫w₀² dx] + (λ/n)Σ[(Vᵢ−Vᵢ,base)/Vᵢ,base]² ``` `objective:'energy'` 为弯曲能目标： ```text U = ∫ M²/(2EI) dx J_U = U/U₀ + (λ/n)Σ[(Vᵢ−Vᵢ,base)/Vᵢ,base]² ``` w₀、U₀ 是同一根梁、同一恒载下没有斜索支承力时的响应。λ=10⁻⁵；小正则项抑制无必要的力偏离，不代表真实桥造价。两个目标都是四个 Vᵢ 的凸二次函数，包含交叉项，不能逐索独立求解。 四索默认约束 0≤Vᵢ≤3200 kN，十二索为0≤Vᵢ≤2000 kN。边界可为标量或与索数相同的数组。**这是竖向分力边界，不是轴向索力 T 边界，也没有支座抗拔或应力约束。** 手动输入较大 V 可能使端支座或塔位支座反力为负，`supportUplift` 检查所有三个反力；这种状态对于不承受拉力的支座不可行。默认两种优化结果的三个反力均为正。 `optimizationStep` 执行一轮精确坐标下降：依次固定其他变量，求当前变量二次函数的最小点，再投影到它的允许区间。一轮中间结果可能不对称，这是扫描次序造成的；收敛后的对称工况应恢复对称。`optimizeForces` 重复扫描到投影梯度小于10⁻¹⁵；默认上限为无塔支承200000轮、有塔支承800000轮。返回 `converged`，不应忽略失败标志。十二索有塔支承的挠度目标从初值需要约381381轮，测试机器约0.6秒；近相关的影响列使其明显慢于能量目标。 每个相邻作用点之间采用五点 Gauss 积分。M² 为四次多项式、w² 为八次多项式，积分对本模型在浮点误差范围内精确。曲线默认以 1 m 间隔输出；`maxAbsMoment` 和 `maxAbsDeflection` 是这 121 个采样点上的最大绝对值，而不是经过连续极值搜索的严格峰值。 ### JavaScript 接口 ```js import { SUPPORTED_DENSE_MODEL, equalTributaryForces, analyzeBeam, optimizationStep, optimizeForces, cableLayout, solveMultiSpan, towerResponse, } from './mechanics.js'; let forces = equalTributaryForces(SUPPORTED_DENSE_MODEL); let current = analyzeBeam(forces, SUPPORTED_DENSE_MODEL, { objective: 'deflection' }); forces = optimizationStep(forces, SUPPORTED_DENSE_MODEL, { objective: 'deflection' }).forces; const optimized = optimizeForces(forces, SUPPORTED_DENSE_MODEL, { objective: 'energy' }); const response = analyzeBeam(optimized.forces, SUPPORTED_DENSE_MODEL, { objective: 'energy' }); ``` `analyzeBeam` 返回 `x,moment,deflection,forces,angles,tensions,horizontalComponents,reactions,energy,rmsDeflection,maxAbsMoment,maxAbsDeflection,cableDeflections,towerDeflection,towerMoment,objective,objectiveComponents,objectiveType,bounds,withinBounds,supportUplift,equilibriumResidual,model`。`reactions` 为 `{left,right,tower}`，未启用塔位支承时 tower=0。`equilibriumResidual` 检查三个反力与全部索托力的整体竖向平衡。 优化函数返回 `forces,before,after,improvement,iterations,converged,projectedGradient`；其中 before、after 是数值目标值，不是分析结果对象。完整优化额外返回 `history`；单步额外返回 `maxChange`。如果输入超出边界，先投影，before 是投影后起点的目标值。 保留的**无塔位支承四索对照**结果（不能作为当前交互默认结果）： | 状态 | 四根索的 V / kN | 最大采样挠度绝对值 / mm | 挠度 RMS / mm | 弯曲能 / kN·m | |---|---|---:|---:|---:| | 分担宽度初值 | 1920，1920，1920，1920 | 8.7264 | 6.3524 | 1.327104 | | 挠度目标最优 | 2195.503，1856.135，1856.135，2195.503 | 0.177994 | 0.081513 | 0.350284 | | 弯曲能目标最优 | 2172.469，1869.573，1869.573，2172.469 | 0.217864 | 0.100801 | 0.349238 | 这些很小的残余挠度来自理想线弹性恒载调平问题；不能解读为真实桥施工精度。 ### 无塔位支承十二索对照 `DENSE_MODEL` 保持 L、EI、q 和塔位置不变，把梁上作用点改为 `[8,16,24,32,40,48,72,80,88,96,104,112]` m。左右各六根索，每根 V 独立，边界0～2000 kN。这里的 `DENSE_MODEL` 默认没有塔位支承；当前有支承交互使用 `SUPPORTED_DENSE_MODEL`。无塔支承的分担宽度基线为 `[640,640,640,640,640,1280,1280,640,640,640,640,640]` kN；靠塔的两根分担较宽的梁区段。 从基线开始的实测收敛结果： | 目标 | 左半侧六个 V / kN（右半侧镜像） | 主梁最大采样 \|M\| / kN·m | 最大采样 \|w\| / mm | 坐标扫描轮数 | |---|---|---:|---:|---:| | 挠度平方积分 | 637.475，620.736，594.355，582.609，610.955，1429.296 | 3015.994 | 0.149755 | 98372 | | 弯曲能 | 720.008，626.105，622.507，754.253，105.970，1716.907 | 3264.369 | 0.138575 | 9287 | 两个返回值均为 `converged:true`；不对称起点也已测试。挠度平方积分最小不等于最大挠度最小，弯曲能最小不等于最大弯矩最小，因此不能根据目标名称把峰值大小关系说反。最优索力接近镜像是对称荷载与对称模型的结果，并没有把左右索力强制绑定。 ### 主塔响应：预先给定索力后的受力复算 `towerResponse({leftCables,rightCables,verticalLoads:[],towerHeight,EI?,samples:101})` 把主塔视作底端固定的线弹性悬臂杆。每根输入索需含 `{z,H,V}`：z 为索在塔上的锚固点距**塔固定底端**的高度，H、V 均为非负大小；左侧索对塔施加 −H，右侧索对塔施加 +H，两侧均向下施加 −V。EI 省略时只返回内力，位移返回 `null`；给定 EI 时计算右偏为正的位移。 `verticalLoads` 接受 `{z,V}`，V向下为正，允许为负。把塔位支反力传给塔时，使用 `verticalLoads:[{z:15,V:beam.reactions.tower}]`：梁所受向上反力对应塔所受向下力。负R塔对应塔所受向上力，需在界面提示支座抗拔。该附加竖向力沿塔轴作用，只改变轴力、底端压缩合量和竖向反力，不改变本一阶模型的水平剪力、弯矩或侧移。未设置该参数时旧API行为保持不变。 以每根索的有符号水平作用 Fᵢ=±Hᵢ 记： ```text Q(z) = Σ_(aᵢ≥z) Fᵢ M(z) = Σ Fᵢ max(aᵢ−z,0) N_c(z) = Σ_(aᵢ≥z) Vᵢ （压缩为正） M(0) = Σ右 Hᵢaᵢ − Σ左 Hᵢaᵢ u(z) = Σ Fᵢ f_c(z,aᵢ) z≤a: f_c = z²(3a−z)/(6EI) z>a: f_c = a²(3z−a)/(6EI) ``` 轴力求和还包括 `verticalLoads` 中位于截面及其上方的有符号V；因此有塔位支承时，桥面以下塔段的轴力不仅是两侧索竖向分力之和，还要加R塔。剪力和轴力在集中力位置取紧贴锚点下方截面的值。此处Q、M按锚点以上外力合量定义符号，底端支承反力、反力矩与对应外力合量反向。返回：`{z,shear,moment,axial,deflection,bottom:{shear,moment,compression},topDeflection,maxAbsMoment,maxAbsDeflection,loads,EI,towerHeight,reactions}`。曲线坐标额外纳入所有锚点；因此可见弯矩斜率变化。塔自重、局部锚固、索竖向分力偏心、二阶轴压效应均未计入。 交互图把塔底放在桥面以下 15 m，桥面以上最高锚点为 40 m。`cableLayout` 中 z 相对于桥面，调用前应平移： ```js const tower = towerResponse({ leftCables: left.map(c => ({ ...c, z: 15 + c.z })), rightCables: right.map(c => ({ ...c, z: 15 + c.z })), verticalLoads: [{ z: 15, V: response.reactions.tower }], towerHeight: 55, EI: 2e8, // 示例值；不代表具体桥塔截面 }); ``` 若从梁优化模型取数据，`analyzeBeam` 的 `horizontalComponents` 给出 H 的大小；按 `positions[i]