(積分)例題3:等分布荷重を受ける単純梁＜【② 構造力学の「構造」例題】様を参考にsympyで実行した。

Last updated at 2025-07-28Posted at 2025-01-16

前　 (積分)例題2:等分布荷重を受ける片持ち梁へ
次　 (積分)例題4:等分布荷重を受ける固定-ローラ梁へ
　　私のsympyの目次へ

オリジナル

Youtube 仮想仕事の原理の教育に関する研究会様　
例題3:等分布荷重を受ける単純梁 (9:00〜10:50)
https://youtu.be/NMP0Q78Wl9Q?t=540
＜目次(18:50)
https://youtu.be/NMP0Q78Wl9Q?t=1130
＜フル【② 構造力学の「構造」例題】(0:00〜19:21)
https://youtu.be/NMP0Q78Wl9Q

sympyで

・ver0.2
　nth_integral_with_conditions_multi_points_debug で

# ver0.2
from sympy import *
import sympy as sp
def nth_integral_with_conditions_multi_points_debug(expr, var, n, conditions):
    C = sp.symbols(f'C1:{n+1}')  # 積分定数
    F = expr
    for i in range(n):
        F = sp.integrate(F, var) + C[i]
    eqs = []
    for (point, deriv_order, val) in conditions:
        dF = F.diff(var, deriv_order)
        eq = dF.subs(var, point) - val
        # print(f"条件: F^{deriv_order}({point}) = {val} → 式: {eq}")
        eqs.append(eq)
    sol = sp.solve(eqs, C, dict=True)
    if not sol:
        raise ValueError(f"\n条件が不十分または矛盾しています。\n"
                         f"定数 {C} に対する方程式系:\n" +
                         "\n".join([str(eq) for eq in eqs]))    
    F_final = F.subs(sol[0])
    return sp.simplify(F_final)
var('x,w,L,EI')
conds = [
    (0, 0, 0),   # 幾何学的境界条件  
    (0, 2, 0),   # 力学的境界条件
    (L, 0, 0),   # 幾何学的境界条件
    (L, 2, 0)    # 力学的境界条件
]
v = nth_integral_with_conditions_multi_points_debug(w/EI, x, 4, conds).expand()
print("#",     v)
print("#",     v   .subs({x:L/2}))
print("#",diff(v,x).subs({x:0}))
# L**3*w*x/(24*EI) - L*w*x**3/(12*EI) + w*x**4/(24*EI)
# 5*L**4*w/(384*EI)
# L**3*w/(24*EI)

・ver0.1

# ver0.1
# (積分)例題3:等分布荷重を受ける単純梁 ＜【② 構造力学の「構造」例題】を参考にsympyで実行した。
from sympy import *
var('x,C1,C2,C3,C4,w,L,EI')
def myMujigennka(v,com):
    return ( str(com)+"*("
            +str((v/com).expand())                \
                 .replace("x**2/L**2","(x/L)**2") \
                 .replace("x**3/L**3","(x/L)**3") \
                 .replace("x**4/L**4","(x/L)**4") 
            +")").replace(" ","")
def myKakeruPw(v,Pw):
    return Pw+"*"+str(v).replace("*"+Pw,"")
EIv4=w                                     #;print(EIv4)
EIv3=integrate(EIv4,x)+w*L   *C1           #;print(EIv3)
EIv2=integrate(EIv3,x)+w*L**2*C2           #;print(EIv2)
EIv1=integrate(EIv2,x)+w*L**3*C3           #;print(EIv1)
EIv =integrate(EIv1,x)+w*L**4*C4           #;print(EIv )
sol =solve([Eq(1/EI*EIv .subs({x:0}),0),   # 幾何的境界条件
            Eq(1/EI*EIv .subs({x:L}),0),
            Eq(    -EIv2.subs({x:0}),0),   # 力学的境界条件
            Eq(    -EIv2.subs({x:L}),0)           
           ],[C1,C2,C3,C4])                # ;print(sol)
v   =(1/EI*(EIv.subs({C1:sol[C1],C2:sol[C2],C3:sol[C3],C4:sol[C4]}))).factor()    #;print(w)
print("#",myKakeruPw(myMujigennka(v,w*L**4/(24*EI)),"w"))                   
print("#",     v   .subs({x:L/2}))
print("#",diff(v,x).subs({x:0}))
# w*L**4/(24*EI)*(x/L-2*(x/L)**3+(x/L)**4)
# 5*L**4*w/(384*EI)
# L**3*w/(24*EI)
#
rep={w:1,L:1,EI:1}
plot (diff(v,x,3).subs(rep),(x,0,1))                   # (+)
plot (diff(v,x,2).subs(rep),(x,0,1),ylim=(-0.2,0.1))   # (+)下に凸

plotで作図

・値は、適当に1.0です。yscaleは、ylim=(,)で対応です。

SFD (-) 例題3:等分布荷重を受ける単純梁

BMD (+)下に凸例題3:等分布荷重を受ける単純梁

sympyで(Beamで）

・ver0.3

# ver0.3
from sympy.physics.continuum_mechanics.beam import Beam
from sympy import *
var('q,E,I,x')             
var('RA,RB')         
var('l',positive=True)    # positive が,あればいいです。
var('RA_2,RB_2')         
def my_Simple_Uniform_Def():
    b=Beam(l,E,I)
    b.apply_load(RA,0,-1)
    b.apply_load(RB,l,-1)
    b.apply_load(q ,0,0,l)
    b.bc_deflection=[(0,0),(l,0)]
    b.solve_for_reaction_loads(RA,RB)
    print("#",b.shear_force())
    print("#",b.bending_moment())
    print("#",b.slope().subs({x:0}))
    return 
# -----------------------------------------------
def my_Simple_Uniform_Draw():
    q_2,l_2,E_2,I_2=1,1,1,1   # ←←←適当に?書いてます。
    b_2=Beam(l_2,E_2,I_2)
    b_2.apply_load(RA_2,0  ,-1)
    b_2.apply_load(RB_2,l_2,-1)
    b_2.apply_load(q_2 ,0,0,l_2)
    b_2.bc_deflection=[(0,0),(l_2,0)]
    b_2.solve_for_reaction_loads(RA_2,RB_2)
    b_2.plot_shear_force()
    b_2.plot_bending_moment()
    # 
    b_2.remove_load  (RA_2,0  ,-1)
    b_2.remove_load  (RB_2,l_2,-1)
    b_2.apply_support(0   ,"pin")
    b_2.apply_support(l_2 ,"roller")
    b_2.draw().show()
    return 
my_Simple_Uniform_Def()
my_Simple_Uniform_Draw()
# l*q*SingularityFunction(x, 0, 0)/2 + l*q*SingularityFunction(x, l, 0)/2 - q*SingularityFunction(x, 0, 1) + q*SingularityFunction(x, l, 1)
# l*q*SingularityFunction(x, 0, 1)/2 + l*q*SingularityFunction(x, l, 1)/2 - q*SingularityFunction(x, 0, 2)/2 + q*SingularityFunction(x, l, 2)/2
# l**3*q/(24*E*I)