ARTICLE DETAIL

资讯详情

深耕网站视觉设计与运营推广的一线实战洞察。

SymPy 梁弯曲分析完全指南:用奇异函数求解 Beam 与 Beam3D 的支座反力、剪力弯矩与挠度

SymPy 梁弯曲分析完全指南:用奇异函数求解 Beam 与 Beam3D 的支座反力、剪力弯矩与挠度 SymPy 梁弯曲分析完全指南用奇异函数求解 Beam 与 Beam3D 的支座反力、剪力弯矩与挠度【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy本文以 SymPy 的连续介质力学模块sympy.physics.continuum_mechanics中的 Beam 与 Beam3D 文档 为骨架系统讲解如何用奇异函数Singularity Function对二维、三维梁进行符号化弯曲求解。读完本文你将掌握梁模型的定义、各类载荷与支座的施加、反力求解、剪力/弯矩/转角/挠度曲线的导出、组合梁与铰处理、影响线ILD计算以及结果可视化并了解每步背后的源码实现与测试验证。模块定位与数学基础beam.py模块的用途一句话即可概括见 模块 docstring本模块利用奇异函数求解力学中的二维梁弯曲问题。梁Beam是一种主要依靠抵抗弯曲来承受载荷的结构构件其力学行为由三个要素刻画横截面轮廓截面二次矩/惯性矩 I、长度与材料属性弹性模量 E。SymPy 的Beam类把这三者抽象为构造参数并将整个求解过程符号化所有载荷、反力、弯矩、转角、挠度都以关于位置变量x的奇异函数表达式呈现而不是数值表。奇异函数SingularityFunction(x, a, n)的核心价值在于它把集中力、集中力矩、均布载荷、线性/抛物线分布载荷统一成同一套表示从而让「载荷 → 剪力 → 弯矩 → 转角 → 挠度」的四次积分/微分链可以在符号层面一次性完成无需分段列方程。核心类 Beam构造与属性构造函数与参数Beam类的构造函数签名如下见 beam.pyBeam(length, elastic_modulus, second_moment, areaSymbol(A), variableSymbol(x), base_charC, ild_variableSymbol(a))各参数说明参数含义说明length梁长任意 Sympifyable可以是符号elastic_modulus弹性模量 E衡量材料刚度也可以是关于位置的连续函数second_moment截面二次矩 I反映面积相对中性轴的分布可以是关于位置的连续函数也可以传入几何模块的Polygon等形状对象此时假设形状对象的 x 轴与梁弯曲轴对齐模块内部自动计算截面二次矩area横截面积默认Symbol(A)用于剪应力计算variable沿梁的位置变量默认Symbol(x)可修改b.variable y会抛出TypeError因为 setter 强制要求Symbol类型base_char积分常量基字符当边界条件不足以确定积分常数时用它生成顺序符号如C3、C4默认Cild_variable影响线移动载荷位置变量默认Symbol(a)符号约定Sign Convention求解时必须采用一致的符号约定结果会自动遵循所选约定但约定必须满足规则——在梁轴正侧产生正剪力的载荷应给出负弯矩。模块 docstring 对此附有allowed-sign-conventions.png示意图说明正弯矩与正转角的方向。另外Beam3D类假设任何分布载荷/分布力矩都施加于梁的全跨度上见 Beam3D docstring。关键属性b.length、b.elastic_modulus、b.second_moment、b.area构造参数的可读写属性setter 内部一律sympify化。b.variable位置变量默认x。b.boundary_conditions返回形如{bending_moment: [], deflection: [...], shear_force: [], slope: [...]}的字典每类边界条件是(位置, 值)元组列表。可通过bc_deflection、bc_slope、bc_bending_moment、bc_shear_force四个属性分别读写见 boundary_conditions。b.load载荷分布曲线奇异函数表达式只读。b.applied_loads所有已施加载荷的(value, start, order, end)元组列表。b.reaction_loads求解后的反力字典。b.rotation_jumps/b.deflection_jumps转动铰/滑动铰处的转角跳变与挠度跳变字典。b.ild_reactions/b.ild_shear/b.ild_moment影响线计算结果。施加载荷apply_load 与阶次体系apply_load(value, start, order, endNone)是载荷建模的唯一入口见 apply_load。核心思想载荷的大小单位是[力/(距离^(n1))]其中 n 是载荷阶次order。载荷类型单位order集中力矩kN·m-2集中力点载荷kN-1均布载荷kN/m0线性分布ramp载荷kN/m²1抛物线分布载荷kN/m³2更高阶分布载荷……start载荷起点对集中力矩与集中力而言就是作用位置。end可选分布载荷在梁内的终止点。源码明确要求提供end时order不能为负否则抛出ValueError(If end is provided the order of the load cannot be negative...)。实现上apply_load先在起点累加value*SingularityFunction(x, start, order)再从end点减去f value*(x-start)**order的前order1阶 Taylor 补偿项coeff*SingularityFunction(x, end, i)/factorial(i)从而把「到端为止」的分布载荷精确截断。内部维护_applied_loads与_original_load两份记录前者供remove_load/绘图使用后者是未代入反力时的原始载荷方程用于影响线计算。remove_load(value, start, order, endNone)精确移除载荷若载荷不存在则抛出ValueError(No such load distribution exists on the beam object.)见 remove_load。完整示例节选自 Beam docstring一根 4 米梁从半长处到末端受 6 N/m 均布载荷两端简支、端部挠度受约束约定向下为正 from sympy.physics.continuum_mechanics.beam import Beam from sympy import symbols, Piecewise E, I symbols(E, I) R1, R2 symbols(R1, R2) b Beam(4, E, I) b.apply_load(R1, 0, -1) # 支座反力未知量先当载荷施加 b.apply_load(6, 2, 0) # 从 x2 开始的均布载荷 b.apply_load(R2, 4, -1) b.bc_deflection [(0, 0), (4, 0)] b.boundary_conditions {bending_moment: [], deflection: [(0, 0), (4, 0)], shear_force: [], slope: []} b.load R1*SingularityFunction(x, 0, -1) R2*SingularityFunction(x, 4, -1) 6*SingularityFunction(x, 2, 0) b.solve_for_reaction_loads(R1, R2) b.load -3*SingularityFunction(x, 0, -1) 6*SingularityFunction(x, 2, 0) - 9*SingularityFunction(x, 4, -1) b.shear_force() 3*SingularityFunction(x, 0, 0) - 6*SingularityFunction(x, 2, 1) 9*SingularityFunction(x, 4, 0) b.bending_moment() 3*SingularityFunction(x, 0, 1) - 3*SingularityFunction(x, 2, 2) 9*SingularityFunction(x, 4, 1) b.slope() (-3*SingularityFunction(x, 0, 2)/2 SingularityFunction(x, 2, 3) - 9*SingularityFunction(x, 4, 2)/2 7)/(E*I) b.deflection() (7*x - SingularityFunction(x, 0, 3)/2 SingularityFunction(x, 2, 4)/4 - 3*SingularityFunction(x, 4, 3)/2)/(E*I) b.deflection().rewrite(Piecewise) (7*x - Piecewise((x**3, x 0), (0, True))/2 - 3*Piecewise(((x - 4)**3, x 4), (0, True))/2 Piecewise(((x - 2)**4, x 2), (0, True))/4)/(E*I)注意rewrite(Piecewise)把奇异函数翻译成分段函数便于理解、代入求值或转数值计算。支座与反力求解apply_support声明式支座apply_support(loc, typefixed)会在指定位置施加支座并返回未知反力符号见 apply_supporttypefixed0 个自由度同时施加反力与反力矩返回(R_loc, M_loc)两个符号typepin1 个自由度返回反力符号typeroller2 个自由度返回反力符号非法类型抛出ValueError(Invalid support type. Choose from pin, roller, or fixed.)。实现上pin/roller会自动apply_load(reaction_load, loc, -1)并把(loc, 0)追加到bc_deflectionfixed还会施加阶次 -2 的反力矩并把(loc, 0)追加到bc_slope。因此使用apply_support后反力符号被内部跟踪后续solve_for_reaction_loads()可以不传任何参数直接求解。solve_for_reaction_loads反力/跳变的统一求解solve_for_reaction_loads(*reactions)同时求解全部反力、转动铰转角跳变与滑动铰挠度跳变见 solve_for_reaction_loads。其底层流程清晰可见用limit(self.shear_force(), x, l)与limit(self.bending_moment(), x, l)建立梁端剪力/弯矩平衡方程遍历shear_force、bending_moment、slope、deflection四类边界条件分别生成方程剪力/弯矩处还会过滤含无穷大的奇异项引入积分常量C3、C4把所有方程梁端平衡 边界条件 积分路径交给linsolve一次性解出(C3, C4) 反力 转角跳变 挠度跳变将解写回self._reaction_loads、self._rotation_jumps、self._deflection_jumps并回代到self._load中。三种等效用法见 solve_for_reaction_loads 示例 from sympy.physics.continuum_mechanics.beam import Beam from sympy import symbols E, I symbols(E, I) R1, R2 symbols(R1, R2) # 方式一用 apply_load 手动施加未知反力求解时必须显式传入符号 b Beam(30, E, I) b.apply_load(-8, 0, -1); b.apply_load(R1, 10, -1); b.apply_load(R2, 30, -1); b.apply_load(120, 30, -2) b.bc_deflection [(10, 0), (30, 0)] b.solve_for_reaction_loads(R1, R2) b.reaction_loads {R1: 6, R2: 2} # 方式二用 apply_support 施加支座并传入返回的符号 b Beam(30, E, I) b.apply_load(-8, 0, -1) R1 b.apply_support(10, pin); R2 b.apply_support(30, pin) b.apply_load(120, 30, -2) b.solve_for_reaction_loads(R1, R2) b.reaction_loads {R_10: 6, R_30: 2} # 方式三用 apply_support 施加支座后直接空参求解符号被内部跟踪 b.solve_for_reaction_loads() b.reaction_loads {R_10: 6, R_30: 2}另外重复传入相同符号会抛出ValueError(Duplicate Symbols passed into solve_for_reaction_loads())对应测试 test_for_duplicate_symbols_to_solve_for_reaction_loads。全符号求解示例任意长度 L 的梁docstring 还给出了一个完全符号化的算例见 beam.py梁长 L受向下集中力 P1L/4 处、向上集中力 P2L/8 处、逆时针力矩 M1L/2 处、顺时针力矩 M23L/4 处、向下的均布 q1L/2→3L/4、向上的均布 q23L/4→L。符号载荷不需要任何假设即可求解但把 L 声明为positiveTrue有助于算法化简 E, I symbols(E, I) L symbols(L, positiveTrue) P1, P2, M1, M2, q1, q2 symbols(P1, P2, M1, M2, q1, q2) R1, R2 symbols(R1, R2) b Beam(L, E, I) b.apply_load(R1, 0, -1); b.apply_load(R2, L, -1) b.apply_load(P1, L/4, -1); b.apply_load(-P2, L/8, -1) b.apply_load(M1, L/2, -2); b.apply_load(-M2, 3*L/4, -2) b.apply_load(q1, L/2, 0, 3*L/4); b.apply_load(-q2, 3*L/4, 0, L) b.bc_deflection [(0, 0), (L, 0)] b.solve_for_reaction_loads(R1, R2) print(b.reaction_loads[R1]) (-3*L**2*q1 L**2*q2 - 24*L*P1 28*L*P2 - 32*M1 32*M2)/(32*L) print(b.reaction_loads[R2]) (-5*L**2*q1 7*L**2*q2 - 8*L*P1 4*L*P2 32*M1 - 32*M2)/(32*L)这是「设计公式符号化推导」的典型场景反力表达式直接以 L、P、M、q 的符号形式给出可直接用于后续设计验算。结果曲线剪力、弯矩、转角与挠度四条核心曲线的计算关系源码实现见 shear_force、bending_moment、slope、deflection剪力 V(x) -∫ q(x) dx 对载荷积分带负号 弯矩 M(x) ∫ V(x) dx 转角 θ(x) -∫ M(x)/(E·I) dx C3 挠度 y(x) ∫ θ(x) dx C4当未设置任何转角/挠度边界条件时slope()/deflection()会保留积分常量用base_char生成的C3、C4或C3:5系列符号设置边界条件后则通过linsolve消去常量。对组合梁I为Piecewise且_joined_beamTrueslope()/deflection()走分段积分路径逐段以Piecewise的各I值积分并利用SingularityFunction(x, prev_end, 0)拼接段间连续性确保不同惯性矩交界处转角/挠度连续见 beam.py 与 beam.py。shear_stress()即shear_force()/area给出平均剪应力曲线。特征量提取极值、反弯点与最大挠度max_shear_force()返回(位置, 最大剪力绝对值)。算法先收集奇异函数位置点逐区间求剪力零点与区间端点、边界极限的绝对值最大值对剪力斜率为零/常数的区间走NotImplementedError兜底分支见 max_shear_force。符号化载荷同样支持对应测试 test_max_shear_force_symbolic_load。max_bmoment()同上逻辑作用于弯矩曲线返回(位置, 最大弯矩绝对值)见 max_bmoment。max_deflection()在(0, length)内求解slope()0的点挠度极值点取这些点处deflection()绝对值最大者无解返回None见 max_deflection。point_cflexure()返回弯矩变号正→负或负→正的反弯点集合。实现先过滤掉阶次 0 的奇异项得到「无奇异」弯矩表达式按奇异点切分区间后用solveset求根再以根两侧中点弯矩之积是否取负判断变号见 point_cflexure。组合梁join 与变截面join(beam, viafixed)把另一根梁刚性连接到当前梁右端返回新的组合梁对象用于处理弹性模量或截面二次矩不连续的情况见 joinviafixed刚性连接默认viahinge通过转动铰连接新梁会在连接处自动施加转动铰其他取值抛出ValueError(Invalid joining method. Choose from fixed or hinge.)两根梁弹性模量不同会抛NotImplementedError当前未实现异 E 拼接截面二次矩不同时新梁的second_moment自动构造为Piecewise((I1, xL1), (I2, xL1L2))。典型用例docstring 示例见 beam.py悬臂梁总长 4 米前 2 米惯性矩1.5*I、后 2 米为I自由端受 20 N 集中力 from sympy.physics.continuum_mechanics.beam import Beam from sympy import symbols E, I symbols(E, I) R1, R2 symbols(R1, R2) b1 Beam(2, E, 1.5*I) b2 Beam(2, E, I) b b1.join(b2, fixed) b.apply_load(20, 4, -1) b.apply_load(R1, 0, -1) b.apply_load(R2, 0, -2) b.bc_slope [(0, 0)] b.bc_deflection [(0, 0)] b.solve_for_reaction_loads(R1, R2) b.load 80*SingularityFunction(x, 0, -2) - 20*SingularityFunction(x, 0, -1) 20*SingularityFunction(x, 4, -1)组合梁相关回归由 test_composite_beam 覆盖。铰与滑动铰处理不连续跳变apply_rotation_hinge(loc)在loc处施加转动铰使转角发生跳变。返回转角跳变符号P_loc物理量纲为转角加载时以E*I*rotation_jump作为阶次 -3 的奇异函数施加并在该处强制弯矩边界条件为 0见 apply_rotation_hinge。求解后可用b.rotation_jumps读取跳变值例如示例中的{P_12: -1875/(16*E*I), P_5: 9625/(24*E*I)}。apply_sliding_hinge(loc)施加滑动铰使挠度发生跳变。返回挠度跳变符号W_loc以E*I*deflection_jump作为阶次 -4 的奇异函数施加并在该处强制剪力边界条件为 0见 apply_sliding_hinge。示例结果为{W_8: 85/24}。两个方法的 docstring 均给出含多个支座 铰的超静定算例见 beam.py 与 beam.py对应测试 test_apply_rotation_hinge 与 test_apply_sliding_hinge。影响线ILD移动载荷下的响应函数影响线Influence Line Diagram描述单位移动载荷沿梁移动时某支座反力、某点剪力或弯矩随移动载荷位置a的变化规律。Beam内置完整支持位置变量默认Symbol(a)即构造参数的ild_variable。solve_for_ild_reactions(value, *reactions)求解反力影响线方程结果存入b.ild_reactions。内部利用_original_load未代入反力的原始载荷叠加value*SingularityFunction(x, a, -1)后积分出剪力/弯矩再结合平衡与边界条件联立求解见 solve_for_ild_reactions。solve_for_ild_shear(distance, value, *reactions)/solve_for_ild_moment(distance, value, *reactions)分别求指定截面剪力、弯矩的影响线方程存入b.ild_shear/b.ild_moment见 solve_for_ild_shear 与 solve_for_ild_moment。绘图方法plot_ild_reactions()、plot_ild_shear()、plot_ild_moment()必须在对应 solve 之后调用否则抛出ValueError提示先生成方程见 plot_ild_reactions。⚠️重要警告docstring 明确标注影响线方程在代入a 0或a ll 为梁长时可能给出错误结果使用与绘图时应避开这两个端点。12 米简支梁在距起点 8 米处设第二支座的剪力影响线示例见 solve_for_ild_shear 示例 from sympy import symbols from sympy.physics.continuum_mechanics.beam import Beam E, I symbols(E, I) R_0, R_8 symbols(R_0, R_8) b Beam(12, E, I) p0 b.apply_support(0, roller) p8 b.apply_support(8, roller) b.solve_for_ild_reactions(1, R_0, R_8) b.solve_for_ild_shear(4, 1, R_0, R_8) b.ild_shear -(-SingularityFunction(a, 0, 0) SingularityFunction(a, 12, 0) 2)*SingularityFunction(a, 4, 0) - SingularityFunction(-a, 0, 0) - SingularityFunction(a, 0, 0) SingularityFunction(a, 0, 1)/8 SingularityFunction(a, 12, 0)/2 - SingularityFunction(a, 12, 1)/8 1对应测试见 test_solve_for_ild_reactions、test_solve_for_ild_shear、test_solve_for_ild_moment。可视化plot 系列与 draw绘图方法统一签名plot_xxx(subsNone)其中subs是以{符号: 数值}形式传入的字典表达式中任何除位置变量外的符号若未在subs中给出会抛出ValueError(Value of %s was not passed.)。这些绘图依赖 matplotlib模块顶部的__doctest_requires__声明了这些 doctest 需要 matplotlib见 beam.py。plot_shear_force()绿色V曲线plot_bending_moment()蓝色M曲线plot_slope()品红色θ曲线plot_deflection()红色δ曲线plot_shear_stress()红色τ曲线需要构造时传入areaplot_loading_results()返回PlotGrid(4, 1, ...)将剪力、弯矩、转角、挠度四张图纵向排布。docstring 中的综合示例8 米梁、E200 GPa、I400e-6 m⁴、集中力 5 kN2m、均布 10 kN/m4→8m、两端简支见 beam.py from sympy.physics.continuum_mechanics.beam import Beam from sympy import symbols R1, R2 symbols(R1, R2) b Beam(8, 200*(10**9), 400*(10**-6)) b.apply_load(5000, 2, -1) b.apply_load(R1, 0, -1) b.apply_load(R2, 8, -1) b.apply_load(10000, 4, 0, end8) b.bc_deflection [(0, 0), (8, 0)] b.solve_for_reaction_loads(R1, R2) b.plot_shear_force() Plot object containing: [0]: cartesian line: 13750*SingularityFunction(x, 0, 0) - 5000*SingularityFunction(x, 2, 0) - 10000*SingularityFunction(x, 4, 1) 31250*SingularityFunction(x, 8, 0) 10000*SingularityFunction(x, 8, 1) for x over (0.0, 8.0)draw(pictorialTrue)则返回梁的示意图依赖 numpy棕色梁身、黑色竖直箭头表示集中力与反力、圆弧箭头表示力矩、阴影区域表示分布载荷、支座符号pin/roller/fixed以及铰白点/白竖线都会绘制正载荷画在梁上方、负荷载画在下方同方向叠加的分布载荷会自动合并见 draw。注意pictorialTrue是缩放的示意图pictorialFalse要求所有分布载荷为数值符号化分布载荷会抛ValueError若梁长为符号表达式绘图时所有符号默认取 10 代值无法判断方向的符号化载荷会以warnings.warn提示示意图可能与计算采用的符号约定不一致。三维扩展Beam3DBeam3D(Beam)处理任意空间方向加载、且各轴截面二次矩不同的梁见 Beam3D。构造签名Beam3D(length, elastic_modulus, shear_modulus, second_moment, area, variableSymbol(x))second_moment可传[I_y, I_z]列表分别绕 y、z 轴单值则两轴相同新增shear_modulus剪切模量 G内部polar_moment()返回极惯性矩非列表时为2*I列表时为I_yI_z载荷施加带方向参数apply_load(value, start, order, diry)、apply_moment_load(value, start, order, diry)dir取x/y/zshear_force()、bending_moment()、slope()、deflection()均返回三个元素的列表分别对应 x/y/z 方向另有axial_force()、axial_stress()、torsional_moment()、angular_deflection()solve_slope_deflection()采用 Timoshenko 梁理论联立求解x 方向用轴向变形方程Derivative(E*A*Derivative(defl(x),x),x)load_x0y/z 方向用含剪切变形的耦合方程组经dsolve与linsolve消元得到三段解析解见 solve_slope_deflectionsolve_for_torsion()处理沿 x 轴出/入梁截面的扭转力矩得到angular_deflection()角位移随 x 变化绘图与极值方法与 2D 版一致但带方向维度plot_shear_force(dirall)、max_shear_force()、max_bending_moment()别名max_bmoment、max_deflection()均返回按方向排列的列表。docstring 完整示例l 米长梁两端固支受 y 向均布 q、z 向均布力矩 m见 beam.py from sympy.physics.continuum_mechanics.beam import Beam3D from sympy import symbols, simplify, collect, factor l, E, G, I, A symbols(l, E, G, I, A) b Beam3D(l, E, G, I, A) x, q, m symbols(x, q, m) b.apply_load(q, 0, 0, diry) b.apply_moment_load(m, 0, -1, dirz) b.shear_force() [0, -q*x, 0] b.bending_moment() [0, 0, -m*x q*x**2/2] b.bc_slope [(0, [0, 0, 0]), (l, [0, 0, 0])] b.bc_deflection [(0, [0, 0, 0]), (l, [0, 0, 0])] b.solve_slope_deflection() factor(b.slope()) [0, 0, x*(-l x)*(-A*G*l**3*q 2*A*G*l**2*q*x - 12*E*I*l*q - 72*E*I*m 24*E*I*q*x)/(12*E*I*(A*G*l**2 12*E*I))]注意Beam3D的边界条件值与 2D 不同bc_slope/bc_deflection的值是沿三轴的三元素列表/元组如[(0, [0, 0, 0]), (l, [0, 0, 0])]。测试与使用前提模块的回归保障由 test_beam.py 提供共 33 个测试函数覆盖基础求解test_Beam、边界条件不足test_insufficient_bconditions、超静定test_statically_indeterminate、带单位符号test_beam_units、变位置变量test_variable_moment、组合梁、反弯点单根/多根/慢用例、载荷移除、apply_support、铰与滑动铰、三类极值、ILD反力/剪力/弯矩/带铰、奇异的 ramp 起点大于终点、符号化 start/end、Beam3D 全家极惯性矩、抛物线载荷、扭矩、三维极值等。使用提示与限制均来自 docstring 与源码可直接作为事实依据必须保持一致的符号约定Beam3D假设分布载荷/分布力矩施加于全跨度apply_load的end参数仅适用于order 0的分布载荷负阶次传end会抛ValueError影响线方法在a 0、a l处结果可能错误draw需要 numpy绘图方法需要 matplotlib组合梁要求两段弹性模量相同异 E 尚未实现plot_xxx系列必须为所有非位置变量提供subs数值否则抛ValueError。按「建梁 → 施载 → 设边界 → 求反力 → 取曲线 → 提取征量/绘图」的标准流程上述 API 可以完整覆盖工程梁弯曲问题的符号求解与验证闭环。【免费下载链接】sympyA computer algebra system written in pure Python项目地址: https://gitcode.com/GitHub_Trending/sy/sympy创作声明:本文部分内容由AI辅助生成(AIGC),仅供参考
返回列表