化：0.71L處彎矩降91.6%的Matlab實(shí)現(xiàn))
做懸挑結(jié)構(gòu)設(shè)計(jì)的人應(yīng)該都體會(huì)過(guò)那種尷尬梁伸出去越長(zhǎng)根部彎矩越大截面只能越做越厚材料大半都在“伺候”那一小段危險(xiǎn)截面。如果現(xiàn)場(chǎng)條件允許在梁下面加一道中間支撐情況會(huì)立刻不一樣——可這道支撐到底放在哪里最劃算很多人會(huì)憑直覺(jué)往自由端附近放我一開(kāi)始也是這么干的。直到我把“最大彎矩隨支座位置變化”的曲線畫(huà)出來(lái)才發(fā)現(xiàn)最優(yōu)位置并不在自由端而是在離固定端約0.71倍梁長(zhǎng)處最大彎矩能降到原來(lái)的十分之一不到。這篇內(nèi)容源自一個(gè)非常經(jīng)典的工程優(yōu)化命題均布載荷懸臂梁的支座位置優(yōu)化。目標(biāo)是在梁上加一個(gè)可移動(dòng)的簡(jiǎn)支支座通過(guò)調(diào)整它的位置來(lái)降低整根梁的最大彎矩最終減小截面尺寸或提高安全冗余。文章會(huì)從力學(xué)建模講起推導(dǎo)支座反力、固定端彎矩和梁上彎矩分布再給出一套完整的Matlab腳本包含符號(hào)推導(dǎo)驗(yàn)證、數(shù)值掃描尋優(yōu)和彎矩圖繪制。適合學(xué)結(jié)構(gòu)力學(xué)、材料力學(xué)的學(xué)生做支架、吊車梁、懸挑平臺(tái)、機(jī)械臂支撐底座的工程師以及所有想用Matlab做參數(shù)化優(yōu)化的人。說(shuō)明一下代碼只用了Matlab基礎(chǔ)函數(shù)不依賴任何工具箱從R2016b到R2026b都能直接跑。下面我按自己實(shí)際做這個(gè)問(wèn)題的思路完整走一遍。1. 懸臂梁的彎矩困境與支座優(yōu)化的基本思路先看原始問(wèn)題。一根長(zhǎng)度為 $L$、抗彎剛度為 $EI$ 的懸臂梁左端固定右端自由全梁受均布載荷 $q$ 作用。這是材料力學(xué)教科書(shū)最經(jīng)典的工況根部固定端彎矩為$$ M_{\max}\frac{qL^2}{2} $$這個(gè)值大家都很熟悉。問(wèn)題在于懸臂梁越靠自由端彎矩越小大部分材料并沒(méi)有充分發(fā)揮作用根部卻早已接近承載極限?，F(xiàn)實(shí)工程里一般有兩種應(yīng)對(duì)思路一是增加截面高度或換材料配合減小懸伸長(zhǎng)度二是加支撐改變結(jié)構(gòu)的傳力路徑。加支撐就是本題的核心思路。在懸臂梁下方加一個(gè)豎向簡(jiǎn)支支座原來(lái)的純懸臂梁就變成“固定端 中間簡(jiǎn)支 自由端懸挑”的混合結(jié)構(gòu)。這個(gè)結(jié)構(gòu)有兩個(gè)關(guān)鍵好處中間支座提供向上的反力相當(dāng)于幫固定端分擔(dān)了一部分外荷載力矩同時(shí)整個(gè)結(jié)構(gòu)的彎矩分布會(huì)被重新分配原本集中在根部的巨大負(fù)彎矩會(huì)被削掉一大塊。但問(wèn)題也隨之而來(lái)這個(gè)中間支座放在哪里如果放在離固定端太近的位置支座反力雖然能幫固定端分擔(dān)力矩但力臂太短效果有限如果放在靠近自由端的位置懸挑段太長(zhǎng)自由端附近的負(fù)彎矩又會(huì)變得很大。直覺(jué)上很多人會(huì)認(rèn)為“越靠近自由端越好”因?yàn)橹坞x荷載作用重心越近越能直接托住載荷。但實(shí)際計(jì)算會(huì)發(fā)現(xiàn)最優(yōu)位置既不靠固定端也不靠自由端而是大致在 $a0.71L$ 的地方。這個(gè)結(jié)論乍看有點(diǎn)反直覺(jué)但可以從彎矩分配的角度解釋。加支座后梁的最大彎矩通常出現(xiàn)在三個(gè)候選位置固定端負(fù)彎矩、支座右側(cè)懸挑段的負(fù)彎矩、AB段內(nèi)部可能出現(xiàn)的正彎矩峰值。這三個(gè)值都隨支座位置 $a$ 變化最優(yōu)位置就是讓其中最大的那個(gè)盡量小。結(jié)構(gòu)優(yōu)化里最樸素的原則就是“削峰填谷”讓幾個(gè)局部峰值的絕對(duì)值相等整體最大值就最小。這個(gè)原則也是后面Matlab數(shù)值優(yōu)化目標(biāo)函數(shù)的設(shè)計(jì)基礎(chǔ)。我不建議只看理論公式想當(dāng)然最好把整個(gè)梁的彎矩分布算出來(lái)再對(duì)每個(gè)可能的支座位置取一遍全梁絕對(duì)值最大值這樣最穩(wěn)妥。2. 力學(xué)建模柔度法求支座反力與三個(gè)彎矩峰值2.1 結(jié)構(gòu)模型與超靜定次數(shù)模型設(shè)置如下梁長(zhǎng) $L$抗彎剛度 $EI$均布載荷 $q$。左端 $A$ 固定右端 $C$ 自由。在距離左端 $a$ 的位置加一個(gè)豎向簡(jiǎn)支支座 $B$也就是限制該點(diǎn)豎向位移但允許轉(zhuǎn)動(dòng)。在支座位置變化的過(guò)程中梁始終保持連續(xù)。這個(gè)結(jié)構(gòu)實(shí)際上是一次超靜定固定端有豎向反力和約束彎矩兩個(gè)未知量中間支座有一個(gè)豎向反力總共三個(gè)未知力而平面一般力系只有兩個(gè)獨(dú)立平衡方程。所以需要補(bǔ)充一個(gè)變形協(xié)調(diào)條件。2.2 去掉B支座回到基本體系用結(jié)構(gòu)力學(xué)里的力法柔度法求解。先把B支座去掉得到一個(gè)左端固定、右端自由的懸臂梁這就是基本體系。現(xiàn)在做兩件事第一算均布載荷 $q$ 單獨(dú)作用下$xa$ 處的豎向撓度。對(duì)懸臂梁距固定端 $x$ 處的彎矩為 $M(\xi)q(L-\xi)^2/2$代入撓曲線微分方程 $EIwM(\xi)$兩次積分并利用固定端轉(zhuǎn)角、撓度均為零的邊界條件可以得到$$ w_q(a)\frac{qa^2}{24EI}\left(6L^2-4Laa^2\right) $$第二在 $xa$ 處加一個(gè)豎直向上的單位力計(jì)算該點(diǎn)的撓度。這里有個(gè)容易繞進(jìn)去的點(diǎn)單位力作用在距離固定端 $a$ 的位置它對(duì)該點(diǎn)產(chǎn)生的撓度其實(shí)等同于一根長(zhǎng)度為 $a$ 的懸臂梁端部受集中力 $P1$ 時(shí)的端部撓度也就是$$ w_1(a)\frac{a^3}{3EI} $$為什么會(huì)等同于“長(zhǎng)度為a的懸臂梁”因?yàn)?$xa$ 右側(cè)那段梁并沒(méi)有參與受力傳遞單位力只讓 $xa$ 截面到固定端這一段發(fā)生彎曲。這一點(diǎn)很多初學(xué)的人會(huì)算錯(cuò)我后面還會(huì)提到。2.3 變形協(xié)調(diào)求支座反力B點(diǎn)是實(shí)際支座真實(shí)位移必須為零。所以基本體系下均布載荷產(chǎn)生的向下?lián)隙缺仨毐恢ё戳?$R_B$ 產(chǎn)生的向上撓度完全抵消$$ w_q(a)-R_B,w_1(a)0 $$解得$$ R_B\frac{q(6L^2-4Laa^2)}{8a} $$這個(gè)公式是后面所有計(jì)算的核心。可以做一個(gè)快速校核當(dāng) $aL$也就是把支座放到自由端時(shí)結(jié)構(gòu)變成“一端固定、一端簡(jiǎn)支”的固簡(jiǎn)梁此時(shí)$$ R_B\frac{q(6L^2-4L^2L^2)}{8L}\frac{3qL}{8} $$這正是固簡(jiǎn)梁受均布載荷時(shí)簡(jiǎn)支端反力的經(jīng)典結(jié)果。公式可信。2.4 固定端反力與彎矩有了 $R_B$整梁豎向平衡可以求固定端豎向反力$$ R_AqL-R_B $$再對(duì)A點(diǎn)取矩得到固定端約束彎矩以大小表示$$ M_A\frac{qL^2}{2}-R_B a\frac{q(4La-2L^2-a^2)}{8} $$注意這個(gè) $M_A$ 是“固定端負(fù)彎矩的大小”實(shí)際截面彎矩寫(xiě)作 $-M_A$。當(dāng) $a$ 大于 $(2-\sqrt{2})L\approx0.586L$ 時(shí)$M_A0$固定端依舊是上側(cè)受拉如果 $a$ 小于這個(gè)值說(shuō)明支座離固定端太近固定端彎矩甚至可能反向。優(yōu)化時(shí)一般只關(guān)心合理區(qū)間。2.5 全梁彎矩分布與三個(gè)候選峰值現(xiàn)在寫(xiě)整根梁的彎矩表達(dá)式。符號(hào)約定采用材料力學(xué)常用約定正彎矩使梁下側(cè)受拉負(fù)彎矩使梁上側(cè)受拉。BC段$a\le x\le L$是懸挑段取截面右側(cè)隔離體只有均布載荷因此$$ M(x)-\frac{q(L-x)^2}{2} $$這個(gè)式子說(shuō)明支座B右側(cè)截面的彎矩并不是0而是$$ M_B-\frac{q(L-a)^2}{2} $$這個(gè)細(xì)節(jié)非常重要。很多人會(huì)想當(dāng)然認(rèn)為“簡(jiǎn)支支座處彎矩為0”那是指梁端鉸支座或者是梁在此處斷開(kāi)的鉸。這里B支座是連續(xù)梁的中間支座梁在B截面處是連續(xù)的左右彎矩必須相等而這個(gè)值由懸挑段決定并不為零。AB段$0\le x\le a$的彎矩取截面左側(cè)隔離體$$ M(x)-M_AR_A x-\frac{qx^2}{2} $$把 $M_A$ 和 $R_A$ 代入后$xa$ 處剛好得到 $-\frac{q(L-a)^2}{2}$與BC段連續(xù)說(shuō)明推導(dǎo)自洽。從這兩段表達(dá)式能提煉出三個(gè)可能的幅值峰值固定端負(fù)彎矩大小 $M_A\frac{q(4La-2L^2-a^2)}{8}$支座右側(cè)懸挑段最大負(fù)彎矩大小 $M_{BC}\frac{q(L-a)^2}{2}$AB段內(nèi)部可能出現(xiàn)的正彎矩峰值先求剪力零點(diǎn)位置 $x_0R_A/q$如果 $x_0\in(0,a)$則峰值約為$$ M_{ABp}\max\left(0,\ -M_A\frac{R_A^2}{2q}\right) $$這三個(gè)峰值隨 $a$ 的變化趨勢(shì)不同。固定端負(fù)彎矩在 $a$ 從0.586L增加到L的過(guò)程中從0增大到 $qL^2/8$懸挑段負(fù)彎矩則從 $0.0857qL^2$ 單調(diào)降到0AB段內(nèi)部正彎矩峰值通常比兩端都小但在某些 $a$ 值下不能忽略。優(yōu)化目標(biāo)函數(shù)其實(shí)就是$$ M_{\max}(a)\max\left(M_A,\ M_{BC},\ M_{ABp}\right) $$我們要找的就是讓這個(gè)函數(shù)最小的 $a$。2.6 解析平衡式與最優(yōu)位置如果暫時(shí)忽略AB段內(nèi)部正彎矩峰值最后驗(yàn)證它不控制最核心的優(yōu)化原則是“固定端負(fù)彎矩與懸挑段負(fù)彎矩峰值相等”$$ \frac{q(4La-2L^2-a^2)}{8}\frac{q(L-a)^2}{2} $$化簡(jiǎn)得$$ 5a^2-12La6L^20 $$解這個(gè)一元二次方程舍去大于 $L$ 的根得到$$ a_{opt}\frac{6-\sqrt{6}}{5}L\approx0.7101L $$此時(shí)固定端負(fù)彎矩$$ M_A\frac{qL^2}{8}\left(4\cdot0.7101-2-0.7101^2\right)\approx0.042qL^2 $$懸挑段負(fù)彎矩$$ M_{BC}\frac{qL^2}{2}(1-0.7101)^2\approx0.042qL^2 $$二者剛好相等。再檢查AB段內(nèi)部正彎矩峰值此時(shí) $R_A0.3549qL$$x_00.3549L$峰值約為 $0.021qL^2$確實(shí)小于兩端不控制。所以最優(yōu)位置理論解成立。對(duì)比原始懸臂梁的 $0.5qL^2$最大彎矩下降到約 $0.042qL^2$降幅達(dá)到91.6%左右。這個(gè)數(shù)據(jù)相當(dāng)可觀。3. Matlab實(shí)現(xiàn)從符號(hào)推導(dǎo)、掃描尋優(yōu)到成圖3.1 代碼設(shè)計(jì)思路我實(shí)現(xiàn)這個(gè)優(yōu)化時(shí)沒(méi)有只寫(xiě)一段“悶頭優(yōu)化”的代碼而是分成四步用符號(hào)計(jì)算驗(yàn)證解析公式確保理論推導(dǎo)沒(méi)有筆誤寫(xiě)一個(gè)目標(biāo)函數(shù)calc_Mmax(a)輸入支座位置返回整根梁絕對(duì)值最大的彎矩在物理合理區(qū)間內(nèi)掃描 $a$畫(huà)出“支座位置-最大彎矩”曲線看全貌用fminbnd在掃描到的最優(yōu)鄰域做精確尋優(yōu)再畫(huà)最優(yōu)彎矩圖。目標(biāo)函數(shù)里不直接寫(xiě)“平衡兩個(gè)峰值”的簡(jiǎn)化公式而是老老實(shí)實(shí)離散整根梁、逐點(diǎn)算彎矩再取絕對(duì)值最大值。這樣最笨但最穩(wěn)能把AB段內(nèi)部正彎矩峰值也自動(dòng)包含進(jìn)去不會(huì)出現(xiàn)“只盯著兩端、漏掉中間峰值”的問(wèn)題。3.2 核心函數(shù)與變量說(shuō)明變量含義如下表變量含義L梁長(zhǎng)無(wú)量綱取1q均布載荷集度無(wú)量綱取1EI抗彎剛度無(wú)量綱取1a支座B到固定端A的距離R_B支座B的豎向反力向上為正R_A固定端豎向反力向上為正M_A固定端負(fù)彎矩大小N全梁離散點(diǎn)數(shù)取1001由于問(wèn)題線性最優(yōu)位置和 $q、EI$ 無(wú)關(guān)所以都把參數(shù)歸一化成1只關(guān)注 $a/L$ 和 $M/(qL^2)$ 這種無(wú)量綱結(jié)果。3.3 完整Matlab代碼下面是我實(shí)際跑通過(guò)的完整腳本注釋比較詳細(xì)可以直接復(fù)制到.m文件運(yùn)行%% 均布載荷懸臂梁的最優(yōu)支座位置優(yōu)化 % 模型A端固定B點(diǎn)簡(jiǎn)支可移動(dòng)C端自由 % 全梁受均布載荷 q求使最大彎矩最小的支座位置 a clear; clc; close all; %% 1. 基本參數(shù)無(wú)量綱化 L 1; % 梁長(zhǎng) q 1; % 均布載荷集度 EI 1; % 抗彎剛度 %% 2. 符號(hào)推導(dǎo)驗(yàn)證可選 syms a_sym real assume(0 a_sym L); % 基本體系去掉B支座的懸臂梁 % 均布載荷在 xa 處產(chǎn)生的向下?lián)隙?delta_q q*a_sym^2*(6*L^2 - 4*L*a_sym a_sym^2)/(24*EI); % 單位力在 xa 處產(chǎn)生的向上撓度 delta_1 a_sym^3/(3*EI); % 變形協(xié)調(diào)w_q - R_B*w_1 0 R_B_sym simplify(delta_q / delta_1); fprintf(支座反力 R_B ); pretty(R_B_sym); % 固定端負(fù)彎矩大小 M_A_sym simplify(q*L^2/2 - R_B_sym*a_sym); fprintf(固定端負(fù)彎矩 M_A ); pretty(M_A_sym); %% 3. 定義目標(biāo)函數(shù)給定 a計(jì)算全梁最大彎矩絕對(duì)值 function Mmax calc_Mmax(a, L, q) % 由柔度法求出的支座反力 R_B q*(6*L^2 - 4*L*a a^2)/(8*a); R_A q*L - R_B; % 固定端負(fù)彎矩大小上側(cè)受拉 M_A q*L^2/2 - R_B*a; % 離散整根梁 N 1001; x linspace(0, L, N); M zeros(size(x)); for i 1:N if x(i) a % AB段固定端負(fù)彎矩 豎向反力 均布載荷 M(i) -M_A R_A*x(i) - q*x(i)^2/2; else % BC段懸挑段均布載荷 M(i) -q*(L - x(i))^2/2; end end Mmax max(abs(M)); end %% 4. 掃描畫(huà)曲線看全局趨勢(shì) a_scan linspace(0.30*L, 0.95*L, 500); M_scan zeros(size(a_scan)); for k 1:length(a_scan) M_scan(k) calc_Mmax(a_scan(k), L, q); end figure(Color, white, Position, [100 100 900 700]); subplot(2,1,1); plot(a_scan/L, M_scan/(q*L^2), b-, LineWidth, 1.6); grid on; hold on; xlabel(支座位置 a/L); ylabel(最大彎矩 / (qL^2)); title(最大彎矩隨支座位置的變化); %% 5. 用 fminbnd 精確優(yōu)化 % 從掃描曲線看最優(yōu)在0.7附近且目標(biāo)函數(shù)單峰給區(qū)間 [0.6, 0.8] [a_opt, M_opt] fminbnd((a) calc_Mmax(a, L, q), 0.6*L, 0.8*L); fprintf(\n 優(yōu)化結(jié)果 \n); fprintf(最優(yōu)支座位置 a/L %.4f\n, a_opt/L); fprintf(對(duì)應(yīng)最大彎矩 Mmax/(qL^2) %.4f\n, M_opt/(q*L^2)); % 解析解對(duì)照 a_theory (6 - sqrt(6))/5 * L; fprintf(理論解析解 a/L %.4f\n, a_theory/L); % 在掃描圖上標(biāo)記最優(yōu)點(diǎn) plot(a_opt/L, M_opt/(q*L^2), ro, MarkerSize, 8, MarkerFaceColor, r); legend(Mmax(a), 最優(yōu)點(diǎn), Location, North); %% 6. 繪制最優(yōu)支座位置下的彎矩圖與無(wú)支座懸臂梁對(duì)比 subplot(2,1,2); x linspace(0, L, 401); M_opt_curve zeros(size(x)); M_no_support zeros(size(x)); R_B_opt q*(6*L^2 - 4*L*a_opt a_opt^2)/(8*a_opt); R_A_opt q*L - R_B_opt; M_A_opt q*L^2/2 - R_B_opt*a_opt; for i 1:length(x) if x(i) a_opt M_opt_curve(i) -M_A_opt R_A_opt*x(i) - q*x(i)^2/2; else M_opt_curve(i) -q*(L - x(i))^2/2; end % 無(wú)支座的純懸臂梁彎矩 M_no_support(i) -q*(L - x(i))^2/2; end plot(x/L, M_opt_curve/(q*L^2), b-, LineWidth, 1.8); hold on; plot(x/L, M_no_support/(q*L^2), k--, LineWidth, 1.3); grid on; xlabel(x/L); ylabel(彎矩 M / (qL^2)); title(最優(yōu)支座位置下的彎矩分布對(duì)比); legend(加支座且 aa_{opt}, 無(wú)支座純懸臂梁, Location, SouthEast);3.4 運(yùn)行結(jié)果與關(guān)鍵輸出腳本運(yùn)行后會(huì)輸出類似下面的命令行信息支座反力 R_B q*(6*L^2 - 4*L*a a^2)/(8*a) 固定端負(fù)彎矩 M_A q*(4*L*a - 2*L^2 - a^2)/8 優(yōu)化結(jié)果 最優(yōu)支座位置 a/L 0.7101 對(duì)應(yīng)最大彎矩 Mmax/(qL^2) 0.0420 理論解析解 a/L 0.7101第一張圖是最大彎矩隨支座位置變化的曲線整體呈“U形”最小值出現(xiàn)在0.71附近。第二張圖是加支座后的最優(yōu)彎矩分布與無(wú)支座懸臂梁的對(duì)比能清楚看到根部負(fù)彎矩被大幅削減彎矩圖在兩個(gè)區(qū)域峰值接近相等。3.5 為什么用掃描加fminbnd而不是純符號(hào)求解理論上解析解已經(jīng)拿到了是不是不需要數(shù)值優(yōu)化實(shí)際做研究或工程時(shí)我非常推薦還是把掃描和數(shù)值優(yōu)化寫(xiě)上。原因有三第一解析解只針對(duì)“均布載荷固定端單個(gè)中間支座”這一種理想模型。如果載荷變成集中力、附加端部力或者支座變成彈性支撐解析式立刻失效但目標(biāo)函數(shù)幾乎不用改。第二解析求解過(guò)程中容易出偽根。比如我看到有人直接solve(M_A M_BC, a)會(huì)得到一個(gè) $a1.6899L$ 的根明顯超出梁長(zhǎng)。如果沒(méi)有物理區(qū)間過(guò)濾很容易被這個(gè)偽根帶偏。第三數(shù)值掃描能讓你看到整個(gè)可行域的趨勢(shì)。曲線長(zhǎng)什么樣、哪里是谷底、哪里可能出現(xiàn)多峰這些信息對(duì)工程判斷非常有用遠(yuǎn)超一個(gè)孤立的最優(yōu)數(shù)值。4. 結(jié)果對(duì)比與工程化擴(kuò)展4.1 不同支座位置的效果對(duì)比下面這張表列出了幾種典型支座位置下的關(guān)鍵結(jié)果無(wú)量綱化載荷 $q1$梁長(zhǎng) $L1$支座位置 $a/L$支座反力 $R_B/(qL)$固定端負(fù)彎矩 $M_A/(qL^2)$懸挑段峰值 $M_{BC}/(qL^2)$AB段內(nèi)正峰值 $M_{ABp}/(qL^2)$全梁最大彎矩 $M_{max}/(qL^2)$無(wú)支座00.5000000.50000.501.06250.03130.125000.12500.71最優(yōu)0.64510.04200.04200.02100.04200.900.44580.09880.00500.05460.0988幾個(gè)值得注意的結(jié)論支座放到0.5L附近時(shí)懸挑段負(fù)彎矩達(dá)到0.125是控制值雖然固定端彎矩已經(jīng)很小但整體最優(yōu)性不如0.71L。支座放到0.9L時(shí)懸挑段幾乎消失了但固定端負(fù)彎矩回升到0.0988AB段內(nèi)部正彎矩也有0.0546最大彎矩仍然比最優(yōu)值大兩倍多。最優(yōu)0.71L時(shí)固定端和懸挑段兩個(gè)峰值幾乎完全相等這就是“削峰填谷”的典型效果。從反力角度看最優(yōu)狀態(tài)下支座反力約為 $0.645qL$接近總荷載的三分之二。這意味著削減梁彎矩的代價(jià)是支撐結(jié)構(gòu)要承受更大的力。設(shè)計(jì)支架、立柱、地基時(shí)必須按這個(gè)反力來(lái)校核不能只看梁本身。4.2 優(yōu)化前后降幅量化無(wú)支座時(shí)全梁最大彎矩$$ M_{\max,0}\frac{qL^2}{2}0.5qL^2 $$加最優(yōu)支座后$$ M_{\max,opt}\approx0.042qL^2 $$降幅約91.6%。換句話說(shuō)在同樣的截面和許用應(yīng)力下懸臂長(zhǎng)度理論上的承載能力可以大幅提升或者說(shuō)同樣的外載荷下梁截面可以明顯減小。這就是這類支座優(yōu)化的工程價(jià)值。4.3 向更復(fù)雜工況擴(kuò)展這套建模思路和Matlab腳本并不是只能算這一道題。我實(shí)際用過(guò)它做了下面幾個(gè)方向的擴(kuò)展都只需要小改目標(biāo)函數(shù)多點(diǎn)支座優(yōu)化。如果現(xiàn)場(chǎng)可以加兩道甚至三道支撐結(jié)構(gòu)變成更高次超靜定。這時(shí)不要讓每個(gè)支座反力都走力法手推太痛苦。直接把支座位置寫(xiě)成向量 $[a_1, a_2, \dots]$每個(gè)位置對(duì)應(yīng)一個(gè)豎向約束用桿件有限元或矩陣位移法組裝剛度矩陣然后目標(biāo)函數(shù)依然取“全梁最大彎矩絕對(duì)值”用fminsearch或多變量粒子群優(yōu)化即可。最優(yōu)結(jié)果會(huì)比單支座更低但要注意支座反力分配更復(fù)雜甚至可能出現(xiàn)個(gè)別支座反力為負(fù)必須加約束。彈性支座。現(xiàn)實(shí)中支撐結(jié)構(gòu)不可能是絕對(duì)剛性的。如果B支座的支撐剛度為 $k$變形協(xié)調(diào)方程從“位移等于零”改成“位移等于反力除以剛度”$$ w_q(a)-R_B w_1(a)\frac{R_B}{k} $$解出來(lái)的反力會(huì)變小最優(yōu)位置也會(huì)向自由端方向移動(dòng)。這個(gè)修正只需要把R_B的表達(dá)式換掉目標(biāo)函數(shù)結(jié)構(gòu)不用動(dòng)。集中力或復(fù)雜載荷。如果梁上除了均布載荷還有集中力、集中力偶基本體系的撓度公式改用疊加法或者直接在代碼里用數(shù)值積分求撓度。目標(biāo)函數(shù)里只要把每一段彎矩表達(dá)式寫(xiě)對(duì)優(yōu)化流程完全一致。無(wú)支座懸臂梁的根部彎矩其實(shí)也可以通過(guò)改變固定端位置來(lái)優(yōu)化。比如一根托架要伸到指定位置固定端可以沿墻上下移動(dòng)本質(zhì)上也是“位置優(yōu)化”問(wèn)題。把固定端坐標(biāo)設(shè)成變量重新寫(xiě)幾何關(guān)系同樣用這個(gè)方法掃一遍即可。4.4 哪些量會(huì)隨EI和q變化哪些不會(huì)有一個(gè)很實(shí)用的判斷在線性彈性范圍內(nèi)最優(yōu)的 $a/L$ 只與梁長(zhǎng)和支座約束形式有關(guān)與 $q$ 和 $EI$ 無(wú)關(guān)。因?yàn)橹ё戳蛷澗乇磉_(dá)式里$q$ 和 $EI$ 都以線性因子的形式出現(xiàn)優(yōu)化目標(biāo)兩側(cè)可以約掉。但如果把最大撓度也作為約束條件情況就不同了。撓度與 $EI$ 直接相關(guān)而支座位置會(huì)同時(shí)影響撓曲線形狀。如果設(shè)計(jì)規(guī)范限定了最大撓度可能需要把撓度寫(xiě)入目標(biāo)函數(shù)變成多目標(biāo)優(yōu)化最優(yōu)位置往往就不再是0.71L了。這也是我建議你把腳本保留下來(lái)、隨時(shí)可以加約束的原因。5. 實(shí)操中的坑與我的幾點(diǎn)補(bǔ)充5.1 最容易踩的坑把中間支座當(dāng)成“梁端鉸”這是我必須單獨(dú)拿出來(lái)說(shuō)的一個(gè)坑。很多人做這道題時(shí)會(huì)想當(dāng)然地認(rèn)為B點(diǎn)是簡(jiǎn)支支座所以B點(diǎn)彎矩一定等于零然后畫(huà)AB段的彎矩圖就畫(huà)錯(cuò)了甚至連解析推導(dǎo)都會(huì)出矛盾。實(shí)際上B點(diǎn)是連續(xù)梁的中間支座不是梁的端部。梁在B截面處并沒(méi)有斷開(kāi)所以彎矩是連續(xù)的大小由懸挑段決定等于$$ M_B-\frac{q(L-a)^2}{2} $$只有當(dāng)你把梁理解成“AE段和BC段是兩根梁在B點(diǎn)用鉸連接”B點(diǎn)彎矩才是零。那本質(zhì)上是一個(gè)不同的結(jié)構(gòu)模型不是本題的連續(xù)梁模型。我建議你做這類問(wèn)題前先明確“支座處的梁截面是否連續(xù)”再?zèng)Q定用哪個(gè)公式。5.2 解析求解的偽根問(wèn)題有人喜歡直接用符號(hào)計(jì)算求解“兩峰相等”的方程。我用solve跑過(guò)一次結(jié)果除了0.7101L還得到一個(gè) $a1.6899L$ 的根。這個(gè)根顯然是偽根因?yàn)槌^(guò)梁長(zhǎng)了。問(wèn)題在于符號(hào)求解不會(huì)自動(dòng)幫你過(guò)濾物理約束。因此解析解拿到后一定要做兩件事一是限定 $0aL$二是代入AB段內(nèi)部正彎矩峰值公式驗(yàn)證它確實(shí)不是控制值。如果省略第二步在某些參數(shù)組合下你可能會(huì)得到一個(gè)看似平衡、實(shí)際被第三個(gè)峰值控制的結(jié)果。5.3 數(shù)值優(yōu)化中的單位與離散點(diǎn)我在代碼里把 $L、q、EI$ 全部設(shè)為1。因?yàn)檫@個(gè)問(wèn)題是無(wú)量綱的最優(yōu)位置不受影響。但如果是實(shí)際工程結(jié)構(gòu)不要把單位搞混尤其是 $q$ 如果是“每米多少千?！?L$ 用“米”最后彎矩單位就是“千牛·米”反力單位是“千?！?。我見(jiàn)過(guò)有人把 $q$ 當(dāng)成總荷載結(jié)果反力大了一倍所有結(jié)果全錯(cuò)。離散點(diǎn)數(shù)量一般取1000到2000就足夠。太少會(huì)讓 $M_{\max}$ 曲線出現(xiàn)毛刺太多則沒(méi)有必要地拖慢掃描。如果你在目標(biāo)函數(shù)里用向量化寫(xiě)法替代for循環(huán)速度還能更快但可讀性會(huì)差一些我保留了循環(huán)寫(xiě)法因?yàn)楦菀卓炊?.4 如果算出來(lái)支座反力為負(fù)說(shuō)明方案不成立支座反力為負(fù)意味著該支座實(shí)際上需要承受拉力也就是支座“脫空”了。比如你把中間支座放到0.5L甚至更靠前的位置$R_A$ 會(huì)變成負(fù)值固定端不僅要提供向上的反力還要往下“拉住”梁這在普通支座設(shè)計(jì)里是不允許的。所以腳本里輸出的R_A和R_B都要檢查。如果出現(xiàn)負(fù)數(shù)說(shuō)明這個(gè)布置方案不能直接使用要么調(diào)整位置要么改用能承受拉力的支座要么改為兩端均固定。這個(gè)約束條件可以在優(yōu)化時(shí)顯式加入比如if R_A 0 || R_B 0 Mmax 1e6; % 罰函數(shù) return; end5.5 關(guān)于這個(gè)項(xiàng)目我最想保留的一點(diǎn)經(jīng)驗(yàn)我最初做這個(gè)優(yōu)化時(shí)憑直覺(jué)認(rèn)為支撐越靠近自由端越好結(jié)果算出來(lái)0.71L時(shí)自己都有點(diǎn)意外。后來(lái)想明白了懸挑段負(fù)彎矩和固定端負(fù)彎矩是兩個(gè)此消彼長(zhǎng)的量只有讓它們相等整體才算“最優(yōu)”。這種“峰值平衡”的思路在做支座位置、加勁肋位置、加強(qiáng)板位置等各類結(jié)構(gòu)優(yōu)化時(shí)都會(huì)反復(fù)用到。一個(gè)建議是把這段腳本里的calc_Mmax函數(shù)當(dāng)成一個(gè)通用模板。以后遇到類似的“移動(dòng)支撐找最優(yōu)位置”問(wèn)題改一改彎矩表達(dá)式就能用。再配合掃描畫(huà)圖先看全局趨勢(shì)、再優(yōu)化整個(gè)過(guò)程非?？煽?。對(duì)我來(lái)說(shuō)這種問(wèn)題的魅力不在于一個(gè)0.71L的答案而在于它逼你把超靜定結(jié)構(gòu)的傳力路徑真正想清楚。算完之后再回頭畫(huà)彎矩圖你會(huì)明顯覺(jué)得自己對(duì)“結(jié)構(gòu)如何分配內(nèi)力”的理解上了一個(gè)臺(tái)階。 懸臂梁支座優(yōu)化：0.71L處彎矩降91.6%的Matlab實(shí)現(xiàn)
做懸挑結(jié)構(gòu)設(shè)計(jì)的人應(yīng)該都體會(huì)過(guò)那種尷尬梁伸出去越長(zhǎng)根部彎矩越大截面只能越做越厚材料大半都在“伺候”那一小段危險(xiǎn)截面。如果現(xiàn)場(chǎng)條件允許在梁下面加一道中間支撐情況會(huì)立刻不一樣——可這道支撐到底放在哪里最劃算很多人會(huì)憑直覺(jué)往自由端附近放我一開(kāi)始也是這么干的。直到我把“最大彎矩隨支座位置變化”的曲線畫(huà)出來(lái)才發(fā)現(xiàn)最優(yōu)位置并不在自由端而是在離固定端約0.71倍梁長(zhǎng)處最大彎矩能降到原來(lái)的十分之一不到。這篇內(nèi)容源自一個(gè)非常經(jīng)典的工程優(yōu)化命題均布載荷懸臂梁的支座位置優(yōu)化。目標(biāo)是在梁上加一個(gè)可移動(dòng)的簡(jiǎn)支支座通過(guò)調(diào)整它的位置來(lái)降低整根梁的最大彎矩最終減小截面尺寸或提高安全冗余。文章會(huì)從力學(xué)建模講起推導(dǎo)支座反力、固定端彎矩和梁上彎矩分布再給出一套完整的Matlab腳本包含符號(hào)推導(dǎo)驗(yàn)證、數(shù)值掃描尋優(yōu)和彎矩圖繪制。適合學(xué)結(jié)構(gòu)力學(xué)、材料力學(xué)的學(xué)生做支架、吊車梁、懸挑平臺(tái)、機(jī)械臂支撐底座的工程師以及所有想用Matlab做參數(shù)化優(yōu)化的人。說(shuō)明一下代碼只用了Matlab基礎(chǔ)函數(shù)不依賴任何工具箱從R2016b到R2026b都能直接跑。下面我按自己實(shí)際做這個(gè)問(wèn)題的思路完整走一遍。1. 懸臂梁的彎矩困境與支座優(yōu)化的基本思路先看原始問(wèn)題。一根長(zhǎng)度為 $L$、抗彎剛度為 $EI$ 的懸臂梁左端固定右端自由全梁受均布載荷 $q$ 作用。這是材料力學(xué)教科書(shū)最經(jīng)典的工況根部固定端彎矩為$$ M_{\max}\frac{qL^2}{2} $$這個(gè)值大家都很熟悉。問(wèn)題在于懸臂梁越靠自由端彎矩越小大部分材料并沒(méi)有充分發(fā)揮作用根部卻早已接近承載極限?，F(xiàn)實(shí)工程里一般有兩種應(yīng)對(duì)思路一是增加截面高度或換材料配合減小懸伸長(zhǎng)度二是加支撐改變結(jié)構(gòu)的傳力路徑。加支撐就是本題的核心思路。在懸臂梁下方加一個(gè)豎向簡(jiǎn)支支座原來(lái)的純懸臂梁就變成“固定端 中間簡(jiǎn)支 自由端懸挑”的混合結(jié)構(gòu)。這個(gè)結(jié)構(gòu)有兩個(gè)關(guān)鍵好處中間支座提供向上的反力相當(dāng)于幫固定端分擔(dān)了一部分外荷載力矩同時(shí)整個(gè)結(jié)構(gòu)的彎矩分布會(huì)被重新分配原本集中在根部的巨大負(fù)彎矩會(huì)被削掉一大塊。但問(wèn)題也隨之而來(lái)這個(gè)中間支座放在哪里如果放在離固定端太近的位置支座反力雖然能幫固定端分擔(dān)力矩但力臂太短效果有限如果放在靠近自由端的位置懸挑段太長(zhǎng)自由端附近的負(fù)彎矩又會(huì)變得很大。直覺(jué)上很多人會(huì)認(rèn)為“越靠近自由端越好”因?yàn)橹坞x荷載作用重心越近越能直接托住載荷。但實(shí)際計(jì)算會(huì)發(fā)現(xiàn)最優(yōu)位置既不靠固定端也不靠自由端而是大致在 $a0.71L$ 的地方。這個(gè)結(jié)論乍看有點(diǎn)反直覺(jué)但可以從彎矩分配的角度解釋。加支座后梁的最大彎矩通常出現(xiàn)在三個(gè)候選位置固定端負(fù)彎矩、支座右側(cè)懸挑段的負(fù)彎矩、AB段內(nèi)部可能出現(xiàn)的正彎矩峰值。這三個(gè)值都隨支座位置 $a$ 變化最優(yōu)位置就是讓其中最大的那個(gè)盡量小。結(jié)構(gòu)優(yōu)化里最樸素的原則就是“削峰填谷”讓幾個(gè)局部峰值的絕對(duì)值相等整體最大值就最小。這個(gè)原則也是后面Matlab數(shù)值優(yōu)化目標(biāo)函數(shù)的設(shè)計(jì)基礎(chǔ)。我不建議只看理論公式想當(dāng)然最好把整個(gè)梁的彎矩分布算出來(lái)再對(duì)每個(gè)可能的支座位置取一遍全梁絕對(duì)值最大值這樣最穩(wěn)妥。2. 力學(xué)建模柔度法求支座反力與三個(gè)彎矩峰值2.1 結(jié)構(gòu)模型與超靜定次數(shù)模型設(shè)置如下梁長(zhǎng) $L$抗彎剛度 $EI$均布載荷 $q$。左端 $A$ 固定右端 $C$ 自由。在距離左端 $a$ 的位置加一個(gè)豎向簡(jiǎn)支支座 $B$也就是限制該點(diǎn)豎向位移但允許轉(zhuǎn)動(dòng)。在支座位置變化的過(guò)程中梁始終保持連續(xù)。這個(gè)結(jié)構(gòu)實(shí)際上是一次超靜定固定端有豎向反力和約束彎矩兩個(gè)未知量中間支座有一個(gè)豎向反力總共三個(gè)未知力而平面一般力系只有兩個(gè)獨(dú)立平衡方程。所以需要補(bǔ)充一個(gè)變形協(xié)調(diào)條件。2.2 去掉B支座回到基本體系用結(jié)構(gòu)力學(xué)里的力法柔度法求解。先把B支座去掉得到一個(gè)左端固定、右端自由的懸臂梁這就是基本體系。現(xiàn)在做兩件事第一算均布載荷 $q$ 單獨(dú)作用下$xa$ 處的豎向撓度。對(duì)懸臂梁距固定端 $x$ 處的彎矩為 $M(\xi)q(L-\xi)^2/2$代入撓曲線微分方程 $EIwM(\xi)$兩次積分并利用固定端轉(zhuǎn)角、撓度均為零的邊界條件可以得到$$ w_q(a)\frac{qa^2}{24EI}\left(6L^2-4Laa^2\right) $$第二在 $xa$ 處加一個(gè)豎直向上的單位力計(jì)算該點(diǎn)的撓度。這里有個(gè)容易繞進(jìn)去的點(diǎn)單位力作用在距離固定端 $a$ 的位置它對(duì)該點(diǎn)產(chǎn)生的撓度其實(shí)等同于一根長(zhǎng)度為 $a$ 的懸臂梁端部受集中力 $P1$ 時(shí)的端部撓度也就是$$ w_1(a)\frac{a^3}{3EI} $$為什么會(huì)等同于“長(zhǎng)度為a的懸臂梁”因?yàn)?$xa$ 右側(cè)那段梁并沒(méi)有參與受力傳遞單位力只讓 $xa$ 截面到固定端這一段發(fā)生彎曲。這一點(diǎn)很多初學(xué)的人會(huì)算錯(cuò)我后面還會(huì)提到。2.3 變形協(xié)調(diào)求支座反力B點(diǎn)是實(shí)際支座真實(shí)位移必須為零。所以基本體系下均布載荷產(chǎn)生的向下?lián)隙缺仨毐恢ё戳?$R_B$ 產(chǎn)生的向上撓度完全抵消$$ w_q(a)-R_B,w_1(a)0 $$解得$$ R_B\frac{q(6L^2-4Laa^2)}{8a} $$這個(gè)公式是后面所有計(jì)算的核心。可以做一個(gè)快速校核當(dāng) $aL$也就是把支座放到自由端時(shí)結(jié)構(gòu)變成“一端固定、一端簡(jiǎn)支”的固簡(jiǎn)梁此時(shí)$$ R_B\frac{q(6L^2-4L^2L^2)}{8L}\frac{3qL}{8} $$這正是固簡(jiǎn)梁受均布載荷時(shí)簡(jiǎn)支端反力的經(jīng)典結(jié)果。公式可信。2.4 固定端反力與彎矩有了 $R_B$整梁豎向平衡可以求固定端豎向反力$$ R_AqL-R_B $$再對(duì)A點(diǎn)取矩得到固定端約束彎矩以大小表示$$ M_A\frac{qL^2}{2}-R_B a\frac{q(4La-2L^2-a^2)}{8} $$注意這個(gè) $M_A$ 是“固定端負(fù)彎矩的大小”實(shí)際截面彎矩寫(xiě)作 $-M_A$。當(dāng) $a$ 大于 $(2-\sqrt{2})L\approx0.586L$ 時(shí)$M_A0$固定端依舊是上側(cè)受拉如果 $a$ 小于這個(gè)值說(shuō)明支座離固定端太近固定端彎矩甚至可能反向。優(yōu)化時(shí)一般只關(guān)心合理區(qū)間。2.5 全梁彎矩分布與三個(gè)候選峰值現(xiàn)在寫(xiě)整根梁的彎矩表達(dá)式。符號(hào)約定采用材料力學(xué)常用約定正彎矩使梁下側(cè)受拉負(fù)彎矩使梁上側(cè)受拉。BC段$a\le x\le L$是懸挑段取截面右側(cè)隔離體只有均布載荷因此$$ M(x)-\frac{q(L-x)^2}{2} $$這個(gè)式子說(shuō)明支座B右側(cè)截面的彎矩并不是0而是$$ M_B-\frac{q(L-a)^2}{2} $$這個(gè)細(xì)節(jié)非常重要。很多人會(huì)想當(dāng)然認(rèn)為“簡(jiǎn)支支座處彎矩為0”那是指梁端鉸支座或者是梁在此處斷開(kāi)的鉸。這里B支座是連續(xù)梁的中間支座梁在B截面處是連續(xù)的左右彎矩必須相等而這個(gè)值由懸挑段決定并不為零。AB段$0\le x\le a$的彎矩取截面左側(cè)隔離體$$ M(x)-M_AR_A x-\frac{qx^2}{2} $$把 $M_A$ 和 $R_A$ 代入后$xa$ 處剛好得到 $-\frac{q(L-a)^2}{2}$與BC段連續(xù)說(shuō)明推導(dǎo)自洽。從這兩段表達(dá)式能提煉出三個(gè)可能的幅值峰值固定端負(fù)彎矩大小 $M_A\frac{q(4La-2L^2-a^2)}{8}$支座右側(cè)懸挑段最大負(fù)彎矩大小 $M_{BC}\frac{q(L-a)^2}{2}$AB段內(nèi)部可能出現(xiàn)的正彎矩峰值先求剪力零點(diǎn)位置 $x_0R_A/q$如果 $x_0\in(0,a)$則峰值約為$$ M_{ABp}\max\left(0,\ -M_A\frac{R_A^2}{2q}\right) $$這三個(gè)峰值隨 $a$ 的變化趨勢(shì)不同。固定端負(fù)彎矩在 $a$ 從0.586L增加到L的過(guò)程中從0增大到 $qL^2/8$懸挑段負(fù)彎矩則從 $0.0857qL^2$ 單調(diào)降到0AB段內(nèi)部正彎矩峰值通常比兩端都小但在某些 $a$ 值下不能忽略。優(yōu)化目標(biāo)函數(shù)其實(shí)就是$$ M_{\max}(a)\max\left(M_A,\ M_{BC},\ M_{ABp}\right) $$我們要找的就是讓這個(gè)函數(shù)最小的 $a$。2.6 解析平衡式與最優(yōu)位置如果暫時(shí)忽略AB段內(nèi)部正彎矩峰值最后驗(yàn)證它不控制最核心的優(yōu)化原則是“固定端負(fù)彎矩與懸挑段負(fù)彎矩峰值相等”$$ \frac{q(4La-2L^2-a^2)}{8}\frac{q(L-a)^2}{2} $$化簡(jiǎn)得$$ 5a^2-12La6L^20 $$解這個(gè)一元二次方程舍去大于 $L$ 的根得到$$ a_{opt}\frac{6-\sqrt{6}}{5}L\approx0.7101L $$此時(shí)固定端負(fù)彎矩$$ M_A\frac{qL^2}{8}\left(4\cdot0.7101-2-0.7101^2\right)\approx0.042qL^2 $$懸挑段負(fù)彎矩$$ M_{BC}\frac{qL^2}{2}(1-0.7101)^2\approx0.042qL^2 $$二者剛好相等。再檢查AB段內(nèi)部正彎矩峰值此時(shí) $R_A0.3549qL$$x_00.3549L$峰值約為 $0.021qL^2$確實(shí)小于兩端不控制。所以最優(yōu)位置理論解成立。對(duì)比原始懸臂梁的 $0.5qL^2$最大彎矩下降到約 $0.042qL^2$降幅達(dá)到91.6%左右。這個(gè)數(shù)據(jù)相當(dāng)可觀。3. Matlab實(shí)現(xiàn)從符號(hào)推導(dǎo)、掃描尋優(yōu)到成圖3.1 代碼設(shè)計(jì)思路我實(shí)現(xiàn)這個(gè)優(yōu)化時(shí)沒(méi)有只寫(xiě)一段“悶頭優(yōu)化”的代碼而是分成四步用符號(hào)計(jì)算驗(yàn)證解析公式確保理論推導(dǎo)沒(méi)有筆誤寫(xiě)一個(gè)目標(biāo)函數(shù)calc_Mmax(a)輸入支座位置返回整根梁絕對(duì)值最大的彎矩在物理合理區(qū)間內(nèi)掃描 $a$畫(huà)出“支座位置-最大彎矩”曲線看全貌用fminbnd在掃描到的最優(yōu)鄰域做精確尋優(yōu)再畫(huà)最優(yōu)彎矩圖。目標(biāo)函數(shù)里不直接寫(xiě)“平衡兩個(gè)峰值”的簡(jiǎn)化公式而是老老實(shí)實(shí)離散整根梁、逐點(diǎn)算彎矩再取絕對(duì)值最大值。這樣最笨但最穩(wěn)能把AB段內(nèi)部正彎矩峰值也自動(dòng)包含進(jìn)去不會(huì)出現(xiàn)“只盯著兩端、漏掉中間峰值”的問(wèn)題。3.2 核心函數(shù)與變量說(shuō)明變量含義如下表變量含義L梁長(zhǎng)無(wú)量綱取1q均布載荷集度無(wú)量綱取1EI抗彎剛度無(wú)量綱取1a支座B到固定端A的距離R_B支座B的豎向反力向上為正R_A固定端豎向反力向上為正M_A固定端負(fù)彎矩大小N全梁離散點(diǎn)數(shù)取1001由于問(wèn)題線性最優(yōu)位置和 $q、EI$ 無(wú)關(guān)所以都把參數(shù)歸一化成1只關(guān)注 $a/L$ 和 $M/(qL^2)$ 這種無(wú)量綱結(jié)果。3.3 完整Matlab代碼下面是我實(shí)際跑通過(guò)的完整腳本注釋比較詳細(xì)可以直接復(fù)制到.m文件運(yùn)行%% 均布載荷懸臂梁的最優(yōu)支座位置優(yōu)化 % 模型A端固定B點(diǎn)簡(jiǎn)支可移動(dòng)C端自由 % 全梁受均布載荷 q求使最大彎矩最小的支座位置 a clear; clc; close all; %% 1. 基本參數(shù)無(wú)量綱化 L 1; % 梁長(zhǎng) q 1; % 均布載荷集度 EI 1; % 抗彎剛度 %% 2. 符號(hào)推導(dǎo)驗(yàn)證可選 syms a_sym real assume(0 a_sym L); % 基本體系去掉B支座的懸臂梁 % 均布載荷在 xa 處產(chǎn)生的向下?lián)隙?delta_q q*a_sym^2*(6*L^2 - 4*L*a_sym a_sym^2)/(24*EI); % 單位力在 xa 處產(chǎn)生的向上撓度 delta_1 a_sym^3/(3*EI); % 變形協(xié)調(diào)w_q - R_B*w_1 0 R_B_sym simplify(delta_q / delta_1); fprintf(支座反力 R_B ); pretty(R_B_sym); % 固定端負(fù)彎矩大小 M_A_sym simplify(q*L^2/2 - R_B_sym*a_sym); fprintf(固定端負(fù)彎矩 M_A ); pretty(M_A_sym); %% 3. 定義目標(biāo)函數(shù)給定 a計(jì)算全梁最大彎矩絕對(duì)值 function Mmax calc_Mmax(a, L, q) % 由柔度法求出的支座反力 R_B q*(6*L^2 - 4*L*a a^2)/(8*a); R_A q*L - R_B; % 固定端負(fù)彎矩大小上側(cè)受拉 M_A q*L^2/2 - R_B*a; % 離散整根梁 N 1001; x linspace(0, L, N); M zeros(size(x)); for i 1:N if x(i) a % AB段固定端負(fù)彎矩 豎向反力 均布載荷 M(i) -M_A R_A*x(i) - q*x(i)^2/2; else % BC段懸挑段均布載荷 M(i) -q*(L - x(i))^2/2; end end Mmax max(abs(M)); end %% 4. 掃描畫(huà)曲線看全局趨勢(shì) a_scan linspace(0.30*L, 0.95*L, 500); M_scan zeros(size(a_scan)); for k 1:length(a_scan) M_scan(k) calc_Mmax(a_scan(k), L, q); end figure(Color, white, Position, [100 100 900 700]); subplot(2,1,1); plot(a_scan/L, M_scan/(q*L^2), b-, LineWidth, 1.6); grid on; hold on; xlabel(支座位置 a/L); ylabel(最大彎矩 / (qL^2)); title(最大彎矩隨支座位置的變化); %% 5. 用 fminbnd 精確優(yōu)化 % 從掃描曲線看最優(yōu)在0.7附近且目標(biāo)函數(shù)單峰給區(qū)間 [0.6, 0.8] [a_opt, M_opt] fminbnd((a) calc_Mmax(a, L, q), 0.6*L, 0.8*L); fprintf(\n 優(yōu)化結(jié)果 \n); fprintf(最優(yōu)支座位置 a/L %.4f\n, a_opt/L); fprintf(對(duì)應(yīng)最大彎矩 Mmax/(qL^2) %.4f\n, M_opt/(q*L^2)); % 解析解對(duì)照 a_theory (6 - sqrt(6))/5 * L; fprintf(理論解析解 a/L %.4f\n, a_theory/L); % 在掃描圖上標(biāo)記最優(yōu)點(diǎn) plot(a_opt/L, M_opt/(q*L^2), ro, MarkerSize, 8, MarkerFaceColor, r); legend(Mmax(a), 最優(yōu)點(diǎn), Location, North); %% 6. 繪制最優(yōu)支座位置下的彎矩圖與無(wú)支座懸臂梁對(duì)比 subplot(2,1,2); x linspace(0, L, 401); M_opt_curve zeros(size(x)); M_no_support zeros(size(x)); R_B_opt q*(6*L^2 - 4*L*a_opt a_opt^2)/(8*a_opt); R_A_opt q*L - R_B_opt; M_A_opt q*L^2/2 - R_B_opt*a_opt; for i 1:length(x) if x(i) a_opt M_opt_curve(i) -M_A_opt R_A_opt*x(i) - q*x(i)^2/2; else M_opt_curve(i) -q*(L - x(i))^2/2; end % 無(wú)支座的純懸臂梁彎矩 M_no_support(i) -q*(L - x(i))^2/2; end plot(x/L, M_opt_curve/(q*L^2), b-, LineWidth, 1.8); hold on; plot(x/L, M_no_support/(q*L^2), k--, LineWidth, 1.3); grid on; xlabel(x/L); ylabel(彎矩 M / (qL^2)); title(最優(yōu)支座位置下的彎矩分布對(duì)比); legend(加支座且 aa_{opt}, 無(wú)支座純懸臂梁, Location, SouthEast);3.4 運(yùn)行結(jié)果與關(guān)鍵輸出腳本運(yùn)行后會(huì)輸出類似下面的命令行信息支座反力 R_B q*(6*L^2 - 4*L*a a^2)/(8*a) 固定端負(fù)彎矩 M_A q*(4*L*a - 2*L^2 - a^2)/8 優(yōu)化結(jié)果 最優(yōu)支座位置 a/L 0.7101 對(duì)應(yīng)最大彎矩 Mmax/(qL^2) 0.0420 理論解析解 a/L 0.7101第一張圖是最大彎矩隨支座位置變化的曲線整體呈“U形”最小值出現(xiàn)在0.71附近。第二張圖是加支座后的最優(yōu)彎矩分布與無(wú)支座懸臂梁的對(duì)比能清楚看到根部負(fù)彎矩被大幅削減彎矩圖在兩個(gè)區(qū)域峰值接近相等。3.5 為什么用掃描加fminbnd而不是純符號(hào)求解理論上解析解已經(jīng)拿到了是不是不需要數(shù)值優(yōu)化實(shí)際做研究或工程時(shí)我非常推薦還是把掃描和數(shù)值優(yōu)化寫(xiě)上。原因有三第一解析解只針對(duì)“均布載荷固定端單個(gè)中間支座”這一種理想模型。如果載荷變成集中力、附加端部力或者支座變成彈性支撐解析式立刻失效但目標(biāo)函數(shù)幾乎不用改。第二解析求解過(guò)程中容易出偽根。比如我看到有人直接solve(M_A M_BC, a)會(huì)得到一個(gè) $a1.6899L$ 的根明顯超出梁長(zhǎng)。如果沒(méi)有物理區(qū)間過(guò)濾很容易被這個(gè)偽根帶偏。第三數(shù)值掃描能讓你看到整個(gè)可行域的趨勢(shì)。曲線長(zhǎng)什么樣、哪里是谷底、哪里可能出現(xiàn)多峰這些信息對(duì)工程判斷非常有用遠(yuǎn)超一個(gè)孤立的最優(yōu)數(shù)值。4. 結(jié)果對(duì)比與工程化擴(kuò)展4.1 不同支座位置的效果對(duì)比下面這張表列出了幾種典型支座位置下的關(guān)鍵結(jié)果無(wú)量綱化載荷 $q1$梁長(zhǎng) $L1$支座位置 $a/L$支座反力 $R_B/(qL)$固定端負(fù)彎矩 $M_A/(qL^2)$懸挑段峰值 $M_{BC}/(qL^2)$AB段內(nèi)正峰值 $M_{ABp}/(qL^2)$全梁最大彎矩 $M_{max}/(qL^2)$無(wú)支座00.5000000.50000.501.06250.03130.125000.12500.71最優(yōu)0.64510.04200.04200.02100.04200.900.44580.09880.00500.05460.0988幾個(gè)值得注意的結(jié)論支座放到0.5L附近時(shí)懸挑段負(fù)彎矩達(dá)到0.125是控制值雖然固定端彎矩已經(jīng)很小但整體最優(yōu)性不如0.71L。支座放到0.9L時(shí)懸挑段幾乎消失了但固定端負(fù)彎矩回升到0.0988AB段內(nèi)部正彎矩也有0.0546最大彎矩仍然比最優(yōu)值大兩倍多。最優(yōu)0.71L時(shí)固定端和懸挑段兩個(gè)峰值幾乎完全相等這就是“削峰填谷”的典型效果。從反力角度看最優(yōu)狀態(tài)下支座反力約為 $0.645qL$接近總荷載的三分之二。這意味著削減梁彎矩的代價(jià)是支撐結(jié)構(gòu)要承受更大的力。設(shè)計(jì)支架、立柱、地基時(shí)必須按這個(gè)反力來(lái)校核不能只看梁本身。4.2 優(yōu)化前后降幅量化無(wú)支座時(shí)全梁最大彎矩$$ M_{\max,0}\frac{qL^2}{2}0.5qL^2 $$加最優(yōu)支座后$$ M_{\max,opt}\approx0.042qL^2 $$降幅約91.6%。換句話說(shuō)在同樣的截面和許用應(yīng)力下懸臂長(zhǎng)度理論上的承載能力可以大幅提升或者說(shuō)同樣的外載荷下梁截面可以明顯減小。這就是這類支座優(yōu)化的工程價(jià)值。4.3 向更復(fù)雜工況擴(kuò)展這套建模思路和Matlab腳本并不是只能算這一道題。我實(shí)際用過(guò)它做了下面幾個(gè)方向的擴(kuò)展都只需要小改目標(biāo)函數(shù)多點(diǎn)支座優(yōu)化。如果現(xiàn)場(chǎng)可以加兩道甚至三道支撐結(jié)構(gòu)變成更高次超靜定。這時(shí)不要讓每個(gè)支座反力都走力法手推太痛苦。直接把支座位置寫(xiě)成向量 $[a_1, a_2, \dots]$每個(gè)位置對(duì)應(yīng)一個(gè)豎向約束用桿件有限元或矩陣位移法組裝剛度矩陣然后目標(biāo)函數(shù)依然取“全梁最大彎矩絕對(duì)值”用fminsearch或多變量粒子群優(yōu)化即可。最優(yōu)結(jié)果會(huì)比單支座更低但要注意支座反力分配更復(fù)雜甚至可能出現(xiàn)個(gè)別支座反力為負(fù)必須加約束。彈性支座。現(xiàn)實(shí)中支撐結(jié)構(gòu)不可能是絕對(duì)剛性的。如果B支座的支撐剛度為 $k$變形協(xié)調(diào)方程從“位移等于零”改成“位移等于反力除以剛度”$$ w_q(a)-R_B w_1(a)\frac{R_B}{k} $$解出來(lái)的反力會(huì)變小最優(yōu)位置也會(huì)向自由端方向移動(dòng)。這個(gè)修正只需要把R_B的表達(dá)式換掉目標(biāo)函數(shù)結(jié)構(gòu)不用動(dòng)。集中力或復(fù)雜載荷。如果梁上除了均布載荷還有集中力、集中力偶基本體系的撓度公式改用疊加法或者直接在代碼里用數(shù)值積分求撓度。目標(biāo)函數(shù)里只要把每一段彎矩表達(dá)式寫(xiě)對(duì)優(yōu)化流程完全一致。無(wú)支座懸臂梁的根部彎矩其實(shí)也可以通過(guò)改變固定端位置來(lái)優(yōu)化。比如一根托架要伸到指定位置固定端可以沿墻上下移動(dòng)本質(zhì)上也是“位置優(yōu)化”問(wèn)題。把固定端坐標(biāo)設(shè)成變量重新寫(xiě)幾何關(guān)系同樣用這個(gè)方法掃一遍即可。4.4 哪些量會(huì)隨EI和q變化哪些不會(huì)有一個(gè)很實(shí)用的判斷在線性彈性范圍內(nèi)最優(yōu)的 $a/L$ 只與梁長(zhǎng)和支座約束形式有關(guān)與 $q$ 和 $EI$ 無(wú)關(guān)。因?yàn)橹ё戳蛷澗乇磉_(dá)式里$q$ 和 $EI$ 都以線性因子的形式出現(xiàn)優(yōu)化目標(biāo)兩側(cè)可以約掉。但如果把最大撓度也作為約束條件情況就不同了。撓度與 $EI$ 直接相關(guān)而支座位置會(huì)同時(shí)影響撓曲線形狀。如果設(shè)計(jì)規(guī)范限定了最大撓度可能需要把撓度寫(xiě)入目標(biāo)函數(shù)變成多目標(biāo)優(yōu)化最優(yōu)位置往往就不再是0.71L了。這也是我建議你把腳本保留下來(lái)、隨時(shí)可以加約束的原因。5. 實(shí)操中的坑與我的幾點(diǎn)補(bǔ)充5.1 最容易踩的坑把中間支座當(dāng)成“梁端鉸”這是我必須單獨(dú)拿出來(lái)說(shuō)的一個(gè)坑。很多人做這道題時(shí)會(huì)想當(dāng)然地認(rèn)為B點(diǎn)是簡(jiǎn)支支座所以B點(diǎn)彎矩一定等于零然后畫(huà)AB段的彎矩圖就畫(huà)錯(cuò)了甚至連解析推導(dǎo)都會(huì)出矛盾。實(shí)際上B點(diǎn)是連續(xù)梁的中間支座不是梁的端部。梁在B截面處并沒(méi)有斷開(kāi)所以彎矩是連續(xù)的大小由懸挑段決定等于$$ M_B-\frac{q(L-a)^2}{2} $$只有當(dāng)你把梁理解成“AE段和BC段是兩根梁在B點(diǎn)用鉸連接”B點(diǎn)彎矩才是零。那本質(zhì)上是一個(gè)不同的結(jié)構(gòu)模型不是本題的連續(xù)梁模型。我建議你做這類問(wèn)題前先明確“支座處的梁截面是否連續(xù)”再?zèng)Q定用哪個(gè)公式。5.2 解析求解的偽根問(wèn)題有人喜歡直接用符號(hào)計(jì)算求解“兩峰相等”的方程。我用solve跑過(guò)一次結(jié)果除了0.7101L還得到一個(gè) $a1.6899L$ 的根。這個(gè)根顯然是偽根因?yàn)槌^(guò)梁長(zhǎng)了。問(wèn)題在于符號(hào)求解不會(huì)自動(dòng)幫你過(guò)濾物理約束。因此解析解拿到后一定要做兩件事一是限定 $0aL$二是代入AB段內(nèi)部正彎矩峰值公式驗(yàn)證它確實(shí)不是控制值。如果省略第二步在某些參數(shù)組合下你可能會(huì)得到一個(gè)看似平衡、實(shí)際被第三個(gè)峰值控制的結(jié)果。5.3 數(shù)值優(yōu)化中的單位與離散點(diǎn)我在代碼里把 $L、q、EI$ 全部設(shè)為1。因?yàn)檫@個(gè)問(wèn)題是無(wú)量綱的最優(yōu)位置不受影響。但如果是實(shí)際工程結(jié)構(gòu)不要把單位搞混尤其是 $q$ 如果是“每米多少千?！?L$ 用“米”最后彎矩單位就是“千牛·米”反力單位是“千?！?。我見(jiàn)過(guò)有人把 $q$ 當(dāng)成總荷載結(jié)果反力大了一倍所有結(jié)果全錯(cuò)。離散點(diǎn)數(shù)量一般取1000到2000就足夠。太少會(huì)讓 $M_{\max}$ 曲線出現(xiàn)毛刺太多則沒(méi)有必要地拖慢掃描。如果你在目標(biāo)函數(shù)里用向量化寫(xiě)法替代for循環(huán)速度還能更快但可讀性會(huì)差一些我保留了循環(huán)寫(xiě)法因?yàn)楦菀卓炊?.4 如果算出來(lái)支座反力為負(fù)說(shuō)明方案不成立支座反力為負(fù)意味著該支座實(shí)際上需要承受拉力也就是支座“脫空”了。比如你把中間支座放到0.5L甚至更靠前的位置$R_A$ 會(huì)變成負(fù)值固定端不僅要提供向上的反力還要往下“拉住”梁這在普通支座設(shè)計(jì)里是不允許的。所以腳本里輸出的R_A和R_B都要檢查。如果出現(xiàn)負(fù)數(shù)說(shuō)明這個(gè)布置方案不能直接使用要么調(diào)整位置要么改用能承受拉力的支座要么改為兩端均固定。這個(gè)約束條件可以在優(yōu)化時(shí)顯式加入比如if R_A 0 || R_B 0 Mmax 1e6; % 罰函數(shù) return; end5.5 關(guān)于這個(gè)項(xiàng)目我最想保留的一點(diǎn)經(jīng)驗(yàn)我最初做這個(gè)優(yōu)化時(shí)憑直覺(jué)認(rèn)為支撐越靠近自由端越好結(jié)果算出來(lái)0.71L時(shí)自己都有點(diǎn)意外。后來(lái)想明白了懸挑段負(fù)彎矩和固定端負(fù)彎矩是兩個(gè)此消彼長(zhǎng)的量只有讓它們相等整體才算“最優(yōu)”。這種“峰值平衡”的思路在做支座位置、加勁肋位置、加強(qiáng)板位置等各類結(jié)構(gòu)優(yōu)化時(shí)都會(huì)反復(fù)用到。一個(gè)建議是把這段腳本里的calc_Mmax函數(shù)當(dāng)成一個(gè)通用模板。以后遇到類似的“移動(dòng)支撐找最優(yōu)位置”問(wèn)題改一改彎矩表達(dá)式就能用。再配合掃描畫(huà)圖先看全局趨勢(shì)、再優(yōu)化整個(gè)過(guò)程非?？煽?。對(duì)我來(lái)說(shuō)這種問(wèn)題的魅力不在于一個(gè)0.71L的答案而在于它逼你把超靜定結(jié)構(gòu)的傳力路徑真正想清楚。算完之后再回頭畫(huà)彎矩圖你會(huì)明顯覺(jué)得自己對(duì)“結(jié)構(gòu)如何分配內(nèi)力”的理解上了一個(gè)臺(tái)階。