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

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

前　 (積分)例題1:集中荷重を受ける片持ち梁へ
次　 (積分)例題3:等分布荷重を受ける単純梁へ
　　私のsympyの目次へ

オリジナル

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

sympyで

・ver0.4
マトリックス変位法
　片持ちはり。等分布荷重なので、等価節点力を追加しました。
　未知数を定義は、リスト同士の足し算と引き算です。

# ver0.4
from sympy import *
var('EI,L')
var('Pyi,Mi,Pyj,Mj')
var('vi ,θi,vj ,θj')
var('q')
def beam_element_stiffness(EI, L):
    """Euler-Bernoulli 梁要素剛性マトリクス（曲げのみ, 4×4）"""
    return EI / L**3 * Matrix([
        [12,     6*L,   -12,    6*L],
        [6*L,  4*L**2, -6*L,  2*L**2],
        [-12,   -6*L,  12  ,   -6*L],
        [6*L,  2*L**2, -6*L,  4*L**2],
    ])
f   =Matrix([[Pyi],[Mi],[Pyj],[Mj]]) 
u   =Matrix([[ vi],[θi],[ vj],[θj]])         
rep ={vi:0,θi:0,Pyj:0,Mj:0}
vars=list( set( [x[0] for x in f.tolist()]   \
               +[x[0] for x in u.tolist()])  
          -set( list(rep.keys())          ))          
print("#",solve( Eq(f+q*L/12*Matrix([[6],[L],[6],[-L]]),                \
                    beam_element_stiffness(EI, L)*u).subs(rep),vars))
# {Mi: -L**2*q/2, Pyi: -L*q, vj: L**4*q/(8*EI), θj: L**3*q/(6*EI)}

・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, 1, 0),   # 幾何学的境界条件
    (L, 2, 0),   # 力学的境界条件
    (L, 3, 0)    # 力学的境界条件
]
v = nth_integral_with_conditions_multi_points_debug(w/EI, x, 4, conds).expand()
print("#",     v)
print("#",     v   .subs({x:L}))
print("#",diff(v,x).subs({x:L}))
# L**2*w*x**2/(4*EI) - L*w*x**3/(6*EI) + w*x**4/(24*EI)
# L**4*w/(8*EI)
# L**3*w/(6*EI)

・ver0.1

# ver0.1
# 例題2:等分布荷重を受ける片持ち梁 ＜【② 構造力学の「構造」例題】を参考に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*EIv1.subs({x:0}),0),
            Eq(    -EIv3.subs({x:L}),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("#",myKakeruPw(myMujigennka(v,w*L**4/(24*EI)),"w"))                   
print("#",myKakeruPw(     v   .subs({x:L})         ,"w"))
print("#",myKakeruPw(diff(v,x).subs({x:L})         ,"w"))
# w*L**4/(24*EI)*(6*(x/L)**2-4*(x/L)**3+(x/L)**4)
# w*L**4/(8*EI)
# w*L**3/(6*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))  # (-)上に凸

plotで作図

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

SFD (+) 例題2:等分布荷重を受ける片持ち梁

BMD (-)上に凸例題2:等分布荷重を受ける片持ち梁

sympyで(Beamで）

・ver0.3

from sympy.physics.continuum_mechanics.beam import Beam
from sympy import *
var('q,E,I,x')             
var('RA,MA')         
var('l',positive=True)    # positive が,あればいいです。
var('RA_2,MA_2')         
def my_Cantilever_Def():
    b=Beam(l,E,I)
    b.apply_load(RA,0,-1)
    b.apply_load(MA,0,-2)
    b.apply_load(q ,0,0,l)
    b.bc_deflection=[(0,0)]
    b.bc_slope     =[(0,0)]
    b.solve_for_reaction_loads(RA,MA)
    print("#",b.deflection().subs({x:l}))
    print("#",b.slope     ().subs({x:l}))
    return 
# -----------------------------------------------
def my_Cantilever_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(MA_2,0,-2)
    b_2.apply_load(q_2 ,0,0,l_2)
    b_2.bc_deflection=[(0,0)]
    b_2.bc_slope     =[(0,0)]
    b_2.solve_for_reaction_loads(RA_2,MA_2)
    b_2.plot_shear_force()
    b_2.plot_bending_moment()
    # 
    b_2.remove_load(RA_2,0,-1)
    b_2.remove_load(MA_2,0,-2)
    b_2.apply_support(0  ,"fixed")
    b_2.draw().show()
    return 
my_Cantilever_Def()
my_Cantilever_Draw()
# l**4*q/(8*E*I)
# l**3*q/(6*E*I)