オリジナル
基本例題8.33 P289 (pdf)＜「材料力学演習(20221021)一括(ver.3.3)」を勉強したい。#sympy
？？？タブレット等で、pdfを開く事ができないかも。？？？
sympyで。単位荷重法で。
・(2024-12-16追加)
# p289 例題8.33
from sympy import *
EI,P1,P2,a,b,x =symbols('EI,P1,P2,a,b,x')
def myCollect(stｒ,P1,P2):
return -collect((sympify(str).expand().simplify()),[P1,P2])
MAB_str="-a*P1-(a+b)*P2+(P1+P2)*x" ;print("#",MAB_str)
MBC_str=MAB_str +"-P1*(x-a)" ;print("#",MBC_str)
MAB =myCollect(MAB_str,P1,P2) #;print("#",MAB_str)
MBC =myCollect(MBC_str,P1,P2) #;print("#",MBC_str)
δ1 = 1/EI*Integral(-MAB*(-MAB.subs({P1:1,P2:0})),(x,0,a)) ;print("#",δ1)
δ2 =str(1/EI*Integral(-MAB*(-MAB.subs({P1:0,P2:1})),(x,0,a))) + "+" + \
str(1/EI*Integral(-MBC*(-MBC.subs({P1:0,P2:1})),(x,a,a+b))) ;print("#",δ2)
# -a*P1-(a+b)*P2+(P1+P2)*x
# -a*P1-(a+b)*P2+(P1+P2)*x-P1*(x-a)
# Integral((-a + x)*(P1*(-a + x) + P2*(-a - b + x)), (x, 0, a))/EI
# Integral((P1*(-a + x) + P2*(-a - b + x))*(-a - b + x), (x, 0, a))/EI+Integral(P2*(-a - b + x)**2, (x, a, a + b))/EI
sympyで。Castigliano's second theoremで。
・出力の順番。注意です。
# ver0.1
from sympy import *
EI,P1,P2,a,b,x =symbols('EI,P1,P2,a,b,x')
MAB_str="-a*P1-(a+b)*P2+(P1+P2)*x" #;print("#",MAB_str)
MBC_str=MAB_str +"-P1*(x-a)" #;print("#",MBC_str)
U = "1/(2*EI)*integrate((" + MAB_str+")**2,(x,0,a ))+" \
+"1/(2*EI)*integrate((" + MBC_str+")**2,(x,a,a+b))" #;print("#",U)
print("#",collect((diff(sympify(U),P1).expand()),[P1,P2]))
print("#",collect((diff(sympify(U),P2).expand()),[P1,P2]))
# P2*(a**3/(3*EI) + a**2*b/(2*EI)) + P1*a**3/(3*EI)
# P1*(a**3/(3*EI) + a**2*b/(2*EI)) + P2*(a**3/(3*EI) + a**2*b/EI + a*b**2/EI + b**3/(3*EI))
sympyで。モーメントのつりあいより。
微分方程式の直接解法で
The Myosotis Methodで
# ver0.3
from sympy import *
var('EI,P1,P2,a,b')
δ1 =(P2*b)*a**2/(2*EI)+(P1+P2)*a**3/(3*EI) ;print("#",δ1)
θ1 =(P2*b)*a / EI +(P1+P2)*a**2/(2*EI) ;print("#",θ1)
δ2d=P2*b**3/(3*EI) #;print("#",δ2d)
δ2 =factor(δ1+θ1*b+δ2d) ;print("#",δ2)
# P2*a**2*b/(2*EI) + a**3*(P1 + P2)/(3*EI)
# P2*a*b/EI + a**2*(P1 + P2)/(2*EI)
# (2*P1*a**3 + 3*P1*a**2*b + 2*P2*a**3 + 6*P2*a**2*b + 6*P2*a*b**2 + 2*P2*b**3)/(6*EI)
Castigliano's second theoremで。
# ver0.4
from sympy import *
EI,P1,P2,a,b,x =symbols('EI,P1,P2,a,b,x')
M01_str="-(P1*(a-x)+P2*(a+b-x))" ;print("#",M01_str)
M12_str= "-P2*(a+b-x)" ;print("#",M12_str)
U = "1/(2*EI)*integrate((" + M01_str+")**2,(x,0,a ))+" \
+"1/(2*EI)*integrate((" + M12_str+")**2,(x,a,a+b))" ;print("#",U)
δ1=collect((diff(sympify(U),P1).expand()),[P1,P2]) ;print("#",δ1)
δ2=collect((diff(sympify(U),P2).expand()),[P1,P2]) ;print("#",δ2)
# P2*(a**3/(3*EI) + a**2*b/(2*EI)) + P1*a**3/(3*EI)
# P1*(a**3/(3*EI) + a**2*b/(2*EI)) + P2*(a**3/(3*EI) + a**2*b/EI + a*b**2/EI + b**3/(3*EI))
割り込みでテスト2例
・以下は単体で、実行できません。ver0.4 の続き
# ver0.4 の続き
var('P,L')
rep={P1:P,a:L,P2:0,b:0}
print("#",δ1.subs(rep),diff(δ1.subs(rep),L))
rep={P1:4*P,a:L,P2:P,b:3*L}
print("#",δ1.subs(rep),diff(δ1.subs(rep),L))
print("#",δ2.subs(rep),diff(δ2.subs(rep),L))
# L**3*P/(3*EI) L**2*P/EI
# 19*L**3*P/(6*EI) 19*L**2*P/(2*EI)
# 86*L**3*P/(3*EI) 86*L**2*P/EI
申し訳ありません。作図は適当にしました。sympyのBeamで。
from sympy import *
from sympy.physics.continuum_mechanics.beam import Beam
from sympy import symbols
E,I,RO,RA,PB=symbols('E,I,RO,RA,PB')
a,b =symbols('a,b',positive=True)
var("x,P1,P2")
# def myBeam_Katamoti_RP_Tawami(a,b,E,I,PB):
# b=Beam(a+b,E,I)
# b.apply_load( RO,0, -1)
# b.apply_load( RA,a ,-1)
# # b.apply_load(-PB,a+b,-1)
# b.apply_load(-PB,2,-1)
# b.bc_deflection.append((0,0))
# b.bc_deflection.append((a,0))
# b.bc_slope .append((0,0))
# b.solve_for_reaction_loads(RO,RA)
# return b.deflection()
def myBeam_Katamoti_P1P2_Sakuzu(a,b,E,I,P1,P2):
b=Beam(a+b,E,I)
# b.apply_load(PB,a+b,-1) # error
b.apply_load (P1,0.75,-1)
b.apply_load (P2,1.00,-1)
b.apply_support(0,'fixed')
p = b.draw()
p.show()
return
# δ =myBeam_Katamoti_RP_Tawami(a,b,E,I,PB)
myBeam_Katamoti_P1P2_Sakuzu (0.75,0.25,1,1,1,1)
いつもの? sympyの実行環境と 参考のおすすめです。
いつもと違うおすすめです。
「弾性曲線方程式の様々な表示」 の表＜ 各諸量とたわみの関係 ＜ wikipedia
・上記の英語版は、見つける事ができませんでした。
参考文献
＞一部に等分布荷重を受ける単純支持はり
＞JSME p89
＞単純支持ばり 7)
＞構造力学公式集 p136