Python
sympy
物理

sympyで不静定はりを解く

概要

sympy1.1でsympy.physics.continuum_mechanics.beamなる梁問題のソルバーが追加されていたらしいです。
Release Notes for 1.1 · sympy/sympy Wiki · GitHub

  • Implementation of Singularity Functions to solve Beam Bending problems (Sampad Kumar Saha's GSoC project).
  • ...
  • Support for Singularity Functions and 2D Beam Problem Solver.

Continuum Mechanics — SymPy 1.1.1 documentation

荷重条件(apply_load)や境界条件の設定方法(bc_slope, bc_deflection)は気に入ったのですが、肝心のsolve_for_reaction_loadsでは力の釣り合いとモーメントの釣り合いだけで解こうとしていて静定問題しか解けないものです。
そこでたわみ角とたわみに境界条件をいい感じにキメて不静定な問題も解けるソルバーが欲しかったので、作りました。

実装

適当な例題として追加例題(6.16(連続はり))を解いたものです。

statically_indeterminate_beam.py
from sympy.core import Symbol, Eq
from sympy.integrals import integrate
from sympy.series import limit
from sympy.solvers import linsolve

from sympy.physics.continuum_mechanics.beam import Beam


def solve_for_reaction_loads(self: Beam, *reactions: Symbol) -> None:
    x = self.variable
    l = self.length
    shear_curve = limit(self.shear_force(), x, l)
    moment_curve = limit(self.bending_moment(), x, l)

    C1 = Symbol('C1')
    slope_curve = integrate(self.bending_moment(), x) + C1

    C2 = Symbol('C2')
    deflection_curve = integrate(slope_curve, x) + C2

    bc_eqs = []
    for position, value in self.boundary_conditions['slope']:
        eqs = Eq(slope_curve.subs(x, position), value)
        bc_eqs.append(eqs)
    for position, value in self.boundary_conditions['deflection']:
        eqs = Eq(deflection_curve.subs(x, position), value)
        bc_eqs.append(eqs)

    equations = [Eq(shear_curve), Eq(moment_curve), *bc_eqs]
    solutions = next(iter(linsolve(equations, [*reactions, C1, C2])))
    self._reaction_loads = dict(zip(reactions, solutions))

    # solution = solve([Eq(shear_curve), Eq(moment_curve), *bc_eqs], [*reactions, C1, C2])
    # self._reaction_loads = {symbol: solution[symbol] for symbol in reactions}
    self._load = self._load.subs(self._reaction_loads)


Beam.solve_for_reaction_loads = solve_for_reaction_loads


def _main():
    from sympy.core import symbols
    from sympy.plotting import plot

    l = symbols('l', positive=True, real=True)
    RA, RB, RC, RE, RF = symbols('R_A, R_B, R_C, R_E, R_F')
    MA, MF = symbols('M_A, M_F')
    P = symbols('P')
    E, I = symbols('E, I')
    loads = [
        (-RA, 0, -1),
        (-RB, l, -1),
        (-RC, (1 + 2) * l, -1),
        (-RE, (1 + 2 + 2 + 2) * l, -1),
        (-RF, (1 + 2 + 2 + 2 + 3) * l, -1),
        (P, (1 + 2 + 2) * l, -1),
        (MA, 0, -2),
        (MF, 10 * l, -2),
    ]

    b = Beam(10 * l, E, I)
    for load in loads:
        b.apply_load(*load)
    b.bc_deflection = [
        (0, 0),
        (l, 0),
        ((1 + 2) * l, 0),
        ((1 + 2 + 2 + 2) * l, 0),
        ((1 + 2 + 2 + 2 + 3) * l, 0),
    ]
    b.bc_slope = [
        (0, 0),
        ((1 + 2 + 2 + 2 + 3) * l, 0),
    ]

    b.solve_for_reaction_loads(RA, RB, RC, RE, RF, MA, MF)
    print(b.reaction_loads)

    # こことんでもなく雑な数字
    constants = {l: 1, P: 1, E: 1, I: 1}
    plot(
        # b.shear_force().subs(constants),
        b.bending_moment().subs(constants),
        # b.slope().subs(constants),
        # b.deflection().subs(constants),
        (b.variable, 0, b.length.subs(constants))
    )

if __name__ == '__main__':
    _main()

あとがき

いろいろ問題を入れて元のsolve_for_reaction_loadsの実装との互換性や、不静定梁の問題の整合性はある程度確認はしていますが、万が一間違った結果が出て単位を落としたり設計したものが破損しても免責とさせて頂きます。

参考文献

Singularity function - Wikipedia
「材料力学」-はりのコツ- 石田良平 平成20年7月1日
特異関数を用いてたわみ曲線を求める