Go to Qiita Advent Calendar 2024 Top
LoginSignup
1
1

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

sympyのBeamで「有名な問題の利用(3)-梁×梁【材料力学・構造力学11】（基礎編）」作図も

Last updated at Posted at 2024-11-03

・タイトルを変更予定です。
・作図は、sympyとFreeCADの2通りです。

パイソニスタの方へ
・計算と作図を別々にしました。
　ひとつにする方法のアドバイスをいただけると幸いです。

オリジナル Youtube

「有名な問題の利用(3)-梁×梁【材料力学・構造力学11】（基礎編）」　　

(0:00~9:36)
筑波大学山本研究室 OPEN COURSEWARE 様

wolframalphaで(オリジナル 様の方法を参考に)

解
R = (16 P)/17, y = -(i l^3 P)/(51 e)

sympyで(オリジナル 様の方法を参考に)

実行環境　https://live.sympy.org/

# 有名な問題の利用(3)-梁×梁【材料力学・構造力学11】（基礎編）
# ver 0.1 
from sympy import *
var('E,I,P,l',positive=True)
var('R'      ,positive=True)
eq   =Eq((P-R)*l**3*Rational(1, 3)/(E*I),R*l**3*Rational(1,48)/(E*I))
sol_R=solve(eq,R)[0]
print("#",               sol_R  )
print("#",eq.rhs.subs({R:sol_R}))
# 16*P/17
# P*l**3/(51*E*I)

sympyのintegrate(積分)で

・ChatGPT-3.5先生へ「sympyの積分で積分定数を表示する方法を教えて下さい。」
https://qiita.com/mrrclb48z/items/71baeb8ea795843632b1

# 有名な問題の利用(3)-梁×梁【材料力学・構造力学11】（基礎編）
# ver 0.2 integrate(積分)で
from sympy import *
var('E,I,P,l',positive=True)
var('R,x'    ,positive=True)
C1 = symbols('C1') 
C2 = symbols('C2') 
def integrate_with_constant_C2(expression, variable):
    return integrate(expression, variable) + C2
def integrate_with_constant_C1(expression, variable):
    return integrate(expression, variable) + C1
def myTawami_KatamotiHari(E,I,myP,myl):
    w=integrate_with_constant_C2(integrate_with_constant_C1((l-x)*myP,x),x)
    sol=solve([Eq(     w   .subs({x:0}),0),
               Eq(diff(w,x).subs({x:0}),0)
          ],[C1,C2])
    return w.subs({C1:sol[C1],C2:sol[C2]})/(E*I)
def myTawami_TanjunHari(E,I,myP,myl):
    w=integrate_with_constant_C2(integrate_with_constant_C1(x*myP,x),x)
    sol=solve([Eq(     w   .subs({x:0}),0),
               Eq(diff(w,x).subs({x:0}),0)
          ],[C1,C2])
    return w.subs({C1:sol[C1],C2:sol[C2]})/(E*I)
eq   =Eq( myTawami_KatamotiHari(E,I,(P-R),x).subs({x:l}),
          myTawami_TanjunHari  (E,I,   R ,x).subs({x:l*Rational(1,2)}))
sol_R=solve(eq,R)[0]
print("#",               sol_R  )
print("#",eq.rhs.subs({R:sol_R}))
# 16*P/17
# P*l**3/(51*E*I)

sympyで、エネルギー法で

・??? 私は、エネルギーの符号の意味を理解していません。???偶然かも

# 有名な問題の利用(3)-梁×梁【材料力学・構造力学11】（基礎編）
# ver 0.3
from sympy import *
var('E,I,P,l',positive=True)
var('x,R')
U_DC=  1/(2*E*I)*integrate(((P-R)           *x)**2,(x,0,l))
U_AB=2*1/(2*E*I)*integrate( (R*Rational(1,2)*x)**2,(x,0,l*Rational(1,2)))
sol_R=solve(Eq(diff(U_DC,R),-diff(U_AB,R)),R)[0]   # ????　わかりませんでした。
print("#",sol_R)
print("#",diff(U_AB,R).subs({R:sol_R}))
# # 16*P/17
# # P*l**3/(51*E*I)
print()
print("#",diff(U_AB+U_AB,R))  # ??? できませんでした
print("#",diff(U_AB-U_AB,R))  # ??? できませんでした

sympyで、片持ちはりのたわみ性マトリックス　より

・??? 勉強中

sympyのBeamで

・Beamで作図をしました。

# 有名な問題の利用(3)-梁×梁【材料力学・構造力学11】（基礎編）
# ver 0.4
from sympy import *
from sympy.physics.continuum_mechanics.beam import Beam
from sympy import symbols
E,I,P= symbols('E,I,P')
l=symbols('l',positive=True)
var("R,x,M,A1,A2,M1,M2")
var("p0,m0")
R1,M1=symbols('R1,M1')
R2,M2=symbols('R2,M2')
def myBeam_Tanjun___Sakuzu(l,E,I,P):
    b=Beam(l,E,I)
    b.apply_load(P ,l*Rational(1,2),-1) 
    b.apply_support(0,'pin')
    b.apply_support(l,'roller')
    p = b.draw()  
    p.show() 
    return 
def myBeam_Katamoti_Sakuzu(l,E,I,P):
    b=Beam(l,E,I)
    b.apply_load(P ,l,-1) 
    b.apply_support(0, 'fixed')
    p = b.draw()  
    p.show() 
    return 
def myBeam_Katamoti_Tawami(l,E,I,P):
    b=Beam(l,E,I)
    b.apply_load(P,l,-1) 
    b.apply_load(R,0,-1)   
    b.apply_load(M,0,-2)
    b.bc_deflection.append((0,0))
    b.bc_slope     .append((0,0))
    b.solve_for_reaction_loads(R,M)
    return b.deflection()
def myBeam_Tanjun___Tawami(l,E,I,P):
    b=Beam(l,E,I)
    b.apply_load(P,l*Rational(1,2),-1) 
    b.apply_load(R1,0,-1)
    b.apply_load(R2,l,-1)
    b.bc_deflection.append((0,0))
    b.bc_deflection.append((l,0))
    b.solve_for_reaction_loads(R1,R2)
    return b.deflection()
sol_R=solve(Eq(myBeam_Katamoti_Tawami(l,E,I,P).subs({x:l})              .subs({P:P-R}),    
               myBeam_Tanjun___Tawami(l,E,I,P).subs({x:l*Rational(1,2)}).subs({P:  R})
              ),R)[0]
print("#",                                                                  sol_R  )
print("#",myBeam_Tanjun___Tawami(l,E,I,P).subs({x:l*Rational(1,2)}).subs({P:sol_R}))
# 16*P/17
# -P*l**3/(51*E*I)
# 作図
myBeam_Katamoti_Sakuzu(1,1,1,1)
myBeam_Tanjun___Sakuzu(1,1,1,1)

DC間 片持ち梁

111.png

AB間 単純梁

222.png

FreeCADで作図

import FreeCAD
import Part
import Draft
import Mesh
#########################################################################################################
# 計算
##############################################################
from sympy import *
l=1
D=  Point(0,0,0)
C=D+Point(l,0,0)
A=C+Point(0,-l*Rational(1,2),0)
B=C+Point(0, l*Rational(1,2),0)
P=C+Point(0,               0,2/10)
##########################################################################################################
### 作図用 私の以下の３D作図バージョン古いです。2024-11-03
##########################################################################################################
### 3D作図
def myXYZ2Txt(A):
     return '(' + str(A.x) + ',' + str(A.y) + ',' + str(A.z) + ')'
    #return ''
def myTxtXYZ(A,myWedgei):
    P5x=float(A.x)
    P5y=float(A.y)
    P5z=float(A.z)
    p5 = FreeCAD.Vector(P5x, P5y, P5z)
    myText = Draft.makeText(myWedgei, p5)
    myText.Label = myWedgei
    #FreeCADGui.ActiveDocument.ActiveObject.FontSize = '0.2 mm'
    #FreeCADGui.ActiveDocument.ActiveObject.FontSize = '0.2 mm'
    FreeCADGui.ActiveDocument.ActiveObject.FontSize = '0.05 mm'
    return
def myTxtXYZ_S(*xy_tx):
    for i in range(1,int(len(xy_tx)/2)+1):
        myTxtXYZ(xy_tx[2*i-2],xy_tx[2*i-1]+myXYZ2Txt(xy_tx[2*i-2]) )
    return
def myLine(A,B):
    Ax,Ay,Az=float(A.x),float(A.y),float(A.z)
    Bx,By,Bz=float(B.x),float(B.y),float(B.z)
    pl = FreeCAD.Placement()
    pl.Rotation.Q = (0.4247081540122249, 0.17592004639554645, 0.33985110062924484, 0.8204732460821097)
    pl.Base = FreeCAD.Vector(-3.9166066876399563, -2.1670824762243774, 1.7495260956243028)
    points = [FreeCAD.Vector(Ax,Ay,Az), FreeCAD.Vector(Bx,By,Bz)]
    line = Draft.make_wire(points, placement=pl, closed=False, face=True, support=None)
    Draft.autogroup(line)
    return
def myLine_S(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    return
def myLine_C(*args):
    for i in range(1,len(args)):
        myLine(args[i-1],args[i])
    myLine(args[i],args[0])
    return
def myLine_H(*args):
    for i in range(1,len(args)):
        myLine(args[0],args[i])
    return
#######################################################
myLine_S  (A,B)
myLine_S  (C,D)
myLine_S  (C,P)
myTxtXYZ_S(A,"A",B,"B",C,"C",D,"D",P,"P")
#
doc = App.activeDocument()
App.ActiveDocument.addObject("App::Origin", "Origin")
App.ActiveDocument.getObject('Origin').Visibility = True
App.ActiveDocument.recompute()
Gui.activeDocument().activeView().viewAxonometric()
Gui.SendMsgToActiveView("ViewFit")

??? 全体図　isometric方向

111.png

拡大図

・矢印の書き方を勉強中。
・黄緑色は、CAD操作で着色です。勉強中。
・座標軸は、CAD操作で非表示です。
222.png

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

(テンプレート) 　開発環境

・以下ができたら、助かります。指定と全部です

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

・以下と同じに見えてきた。

1
1
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
1
1

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?