Go to Qiita Advent Calendar 2024 Top
LoginSignup
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

@mrrclb48z

simple beam,a downward point load はりの応力とSFD，BMD 基本例題4.16 P95 「材料力学演習(20221021)」をsympyでやってみたい。

Last updated at Posted at 2023-06-01

パイソニスタの方へ　教えて下さい。
・symbol,varの重複宣言は間違っていますか？
(2024-11-16　ver0.3追加)

・ver0.2 勉強の成果という事で。
・ver0.1
　sympyのPiecewise関数とmaximum関数でチャレンジ。
　??? sympyのPiecewise関数で最大値は計算できませんでした。
　sympyのPiecewise関数で作図plotできました。

オリジナル

基本例題4.16 p95＜「材料力学演習(20221021)一括(ver.3.3)」を勉強したい。#sympy
　　　　　　　　　　？？？タブレット等で、pdfを開く事ができないかも。？？？　

sympyで

・ver0.2
・(未)最大値の計算は、ユーザー定義関数にする予定です。
・【Python3】Decimal.quantizeを用いて小数点第○位以下を四捨五入

# ver0.2
# 基本例題4.16 P95
from sympy import *
from decimal import Decimal, ROUND_HALF_UP

var('FAC,x,a,MAC,FCB,P,MCB,l,RB',real=True)
# var('RA'      ,real=True)
RA = Symbol('RA')
def myRound1(x, d=1):
    p = Decimal(str(x)).quantize(Decimal(str(1/10**d)), rounding=ROUND_HALF_UP)
    p = float(p)
    return p
def myDic (myStrKey,myValue):
    return {sympify(myStrKey):myValue}
FAC=solve(-RA  +FAC        ,FAC)[0]   ;print("# FAC =",FAC)
MAC=solve(-RA*x+MAC        ,MAC)[0]   ;print("# MAC =",MAC)
FCB=solve(-RA  +P+FCB      ,FCB)[0]   ;print("# FAC =",FCB)
MCB=solve(-RA*x+P*(x-a)+MCB,MCB)[0]   ;print("# MCB =",poly(MCB,P).args[0])
RA =solve(MCB.subs({x:l})  ,RA )[0]   ;print("# RA  =",RA)           
RB =solve( RA+RB-P         ,RB )[0]   ;print("# RB  =",RB)                  ;print()
# 
rep = {l:1.0,a:0.4,P:1000}
rep |= myDic("RA",myRound1(RA.subs(rep)))                               
print("# ",rep )
# 
print("# RA  =",myRound1(RA .subs(rep)))
print("# RB  =",myRound1(RB .subs(rep)))
print("# FAC =",myRound1(FAC.subs(rep)))
print("# FCB =",myRound1(FCB.subs(rep)))
print("# MAC =",MAC.subs(rep))
print("# MCB =",MCB.subs(rep))
print("# Fmax =",myRound1(max    (abs(FAC.subs(rep)),abs(FAC.subs(rep)))))
print("# Mmax =",myRound1(maximum(MAC.subs(rep),x,Interval(0.0,0.4)))    )
# FAC = RA
# MAC = RA*x
# FAC = -P + RA
# MCB = P*(a - x) + RA*x
# RA  = P*(-a + l)/l
# RB  = P*a/l

#  {l: 1.0, a: 0.4, P: 1000, RA: 600.0}
# RA  = 600.0
# RB  = 400.0
# FAC = 600.0
# FCB = -400.0
# MAC = 600.0*x
# MCB = 400.0 - 400.0*x
# Fmax = 600.0
# Mmax = 240.0

・ver0.1
　??? エラーです。

# ver0.1
# 基本例題4.16 P95

from sympy import *
l,a,x  =symbols('l a x'  ,real=True)
P,RA,RB=symbols('P RA RB',real=True)
FAC,MAC=symbols('FAC MAC',real=True)
FCB,MCB=symbols('FCB MCB',real=True)
x      =symbols('x'      ,real=True)
#
myF=Eq(-RA+FAC  ,0)
myM=Eq(-RA*x+MAC,0)          
ans=solve([myF,myM],[FAC,MAC])
#
FAC_ans,MAC_ans=ans[FAC],ans[MAC]
myF=Eq(-RA+P+FCB        ,0)
myM=Eq(-RA*x+P*(x-a)+MCB,0)          
ans=solve([myF,myM],[FCB,MCB])
FCB_ans,MCB_ans=ans[FCB],ans[MCB]
RA_ans=solve(MCB_ans.subs({x:l}),RA)[0]
RB_ans=P-RA_ans
#
FAC_ans=FAC_ans.subs({RA:RA_ans})
FCB_ans=FCB_ans.subs({RA:RA_ans})
MAC_ans=MAC_ans.subs({RA:RA_ans})
MCB_ans=MCB_ans.subs({RA:RA_ans})
# print("# RA_ans ",RA_ans )
# print("# RB_ans ",RB_ans )
# print("# FAC_ans",FAC_ans)
# print("# FCB_ans",FCB_ans)
# print("# MAC_ans",MAC_ans)
# print("# MCB_ans",MCB_ans)
# print("################################")
P_v,l_v,a_v=1000,1.0,0.4
mySubs ={P:P_v,l:l_v,a:a_v}
#
print("# RA_ans ",round(RA_ans .subs(mySubs)))
print("# RB_ans ",round(RB_ans .subs(mySubs)))
print("# FAC_ans",round(FAC_ans.subs(mySubs)))
print("# FCB_ans",round(FCB_ans.subs(mySubs)))
print("# MAC_ans",      MAC_ans.subs(mySubs) )
print("# MCB_ans",      MCB_ans.subs(mySubs) )
BMD=Piecewise(( MAC_ans.subs(mySubs) ,x< a),(MCB_ans.subs(mySubs), True)).subs(mySubs)
myMaxBMD=maximum(BMD,x,Interval(0,l_v))  
# plot (BMD,(x,0,l_v))
print("#    BMD ",     BMD)
print("# ????? MaxBMD ",myMaxBMD)
# RA_ans  600
# RB_ans  400
# FAC_ans 600
# FCB_ans -400
# MAC_ans 600.0*x
# MCB_ans 400.0 - 400.0*x
#    BMD  Piecewise((600.0*x, x < 0.4), (400.0 - 400.0*x, True))
#  ????? MaxBMD  0

BMDを作図(sympyのplot moduleで)

0png.png

(勉強中) Mathematicaで

(勉強中) MATLAB 追加のオプションSymbolic Math Toolbox で

sympyのBeamで

・本命案　勉強中

# ver0.3
# 基本例題4.16 P95
# 作図には、具体的な数値が必要です。
from sympy import *
from sympy.physics.continuum_mechanics.beam import Beam
from sympy import symbols
from decimal import Decimal, ROUND_HALF_UP
E,I,P= symbols('E,I,P')
var("myE,myI,myP,myl,myla")
var("R,x,M,A1,A2,M1,M2,R1,R2")
def remove_last_term(expr):
    # 式を項ごとに分解し、最後の項を除去
    terms = expr.as_ordered_terms()[:-1]  # 最後の項を取り除く
    # 残りの項を使って式を再構成
    return Add(*terms)
def MYshear_stress(expr):
    var('A')
    return remove_last_term(sympify(expr.subs({A:1})))
def myRound1(x, d=1):
    p = Decimal(str(x)).quantize(Decimal(str(1/10**d)), rounding=ROUND_HALF_UP)
    p = float(p)
    return p
def myBeam_Tanjun_reaction_loads(l,E,I,P,la):
    E ,I =symbols('E ,I ')
    b1=Beam(l,E,I)
    E ,I,R1,R2 =symbols('E ,I,R1,R2 ')
    b1.apply_load(R1, 0,-1)
    b1.apply_load(R2, l,-1)
    b1.apply_load(P ,la,-1) 
    # b1.bc_deflection.append([(0,0),(l,0)]) # あってはいけない
    b1.solve_for_reaction_loads(R1,R2)
    R1R2=b1.reaction_loads
    return R1R2[R1],R1R2[R2]
def myBeam_Tanjun_Sakuzu_draw   (l,E,I,P,la):
    E ,I =symbols('E ,I ')
    b2=Beam(l,E,I)
    E ,I,R1,R2 =symbols('E ,I,R1,R2 ')
    b2.apply_load(R1, 0,-1)
    b2.apply_load(R2, l,-1)
    b2.apply_load(P ,la,-1) 
    b2.apply_support(0 ,'pin')
    b2.apply_support(l ,'roller')
    p = b2.draw()  
    p.show() 
    return 
def myBeam_Tanjun_Sakuzu_SFDBMD(l,E,I,P,la):
    E ,I =symbols('E ,I ')
    b3=Beam(l,E,I)
    E ,I,R1,R2 =symbols('E ,I,R1,R2 ')
    b3.apply_load   (R1, 0,-1)
    b3.apply_load   (P ,la,-1) 
    b3.apply_load   (R2, l,-1)
    # b3.bc_deflection.append([(0,0),(l,0)]) # あってはいけない
    # b3.apply_support(0,'pin'   )           # あってはいけない
    # b3.apply_support(l,'roller')           # あってはいけない
    b3.bc_deflection = [(0, 0), (l, 0)]
    b3.solve_for_reaction_loads(R1,R2)
    # 
    print("#",b3.slope())       # ←←←　???　Eの復活
    # b3.plot_slope()
    # plot(b3.slope)
    b3.plot_bending_moment()
    plot(MYshear_stress(b3.shear_stress()),(x,0,l))
    return 
myl,myP,myla=1.0,1000,0.4
R1,R2=myBeam_Tanjun_reaction_loads(myl,myE,myI,myP,myla) ;print("#",myRound1(R1,1),myRound1(R2,1))
myl,myE,myI,myP,myla=1.0,200*(10**9),400*(10**6),1000,0.4
myBeam_Tanjun_Sakuzu_draw  (myl,myE,myI,myP,myla)
myBeam_Tanjun_Sakuzu_SFDBMD(myl,myE,myI,myP,myla)
# -600.0 -400.0
# (-300.0*SingularityFunction(x, 0, 2) + 500*SingularityFunction(x, 0.4, 2) - 200.0*SingularityFunction(x, 1.0, 2) + 64.0)/(E*I)

000.png
111.png
2222.png
https://docs.sympy.org/latest/modules/physics/continuum_mechanics/beam.html

いつもの? sympyの実行環境と 参考のおすすめです。

(テンプレート)

いつもと違うおすすめです。

sympy.functions.elementary.piecewise.ExprCondPair(expr, cond)

sympy.calculus.util.maximum(f, symbol, domain=Reals)

参考文献

＞一部に等分布荷重を受ける単純支持はり
＞JSME p89
　
＞片持ちばり 1)
＞構造力学公式集 表5.1 片持ちばりの公式　p128

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?