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

静力学と材料力学 演習問題1,2,3 P27 「材料力学演習(20221021)」をsympyでやってみた。

Last updated at Posted at 2023-05-05

オリジナル

演習問題1,2,3 P27
http://zairikiweb.starfree.jp/zai_enshuh/zai_enshuh_ver.3.3.pdf#page=32

「材料力学演習(20221021)一括(ver.3.3)」を勉強したい。＞sympy
https://qiita.com/mrrclb48z/items/c3274daf5f7a99cfe027#%E4%BD%9C%E6%A5%AD%E4%B8%AD%E7%A7%81%E3%81%AE%E8%A7%A3%E7%AD%94sympy

sympyで

？？？？？四捨五入　#問題2. dAB= 44.0000000000000

# ver0.1
from sympy import *
NAB,NBC,NCD =symbols('NAB NBC NCD',real=True)
P1,P2       =symbols('P1  P2'     ,real=True)
RA,RB       =symbols('RA  RB'     ,real=True)
A  ,EE      =symbols('A   EE'     ,real=True)
lAB,lBC,lCD =symbols('lAB lBC lCD',real=True)
RA ,RD      =symbols('RA  RD'     ,real=True)
myAsyuhen=Eq( RA +NAB   ,0)
myBsyuhen=Eq(-NAB+P1+NBC,0)
myCsyuhen=Eq(-NBC-P2+NCD,0)
myDsyuhen=Eq( NCD+RD    ,0)
dAB=NAB*lAB/(A*EE)
dBC=NBC*lBC/(A*EE)
dCD=NCD*lCD/(A*EE)
myHenkei=Eq(dAB+dBC+dCD,0)
ans=solve([myAsyuhen,
           myBsyuhen,
           myCsyuhen,
           myDsyuhen,
           myHenkei],[RA,NAB,NBC,NCD,RD])
mySubs1={P1 :15 ,P2 :10,
         lAB:0.6,lBC:1.0,lCD:1.4,
         A  :500*10**(-6),EE:200*10**(+6)}
RA__ans=ans[RA ].subs(mySubs1)
NAB_ans=ans[NAB].subs(mySubs1)
NBC_ans=ans[NBC].subs(mySubs1)
NCD_ans=ans[NCD].subs(mySubs1)
RD__ans=ans[RD ].subs(mySubs1)
print("#問題1. RA =",round(RA__ans,2))
print("#問題1. RD =",round(RD__ans,2))
mySubs2={RA :RA__ans,
         NAB:NAB_ans,NBC:NBC_ans,NCD:NCD_ans,
         RD :RD__ans}
print("#問題2. NAB=",round(NAB_ans,2))
print("#問題2. NBC=",round(NBC_ans,2))
print("#問題2. NCD=",round(NCD_ans,2))
dAB__ans=dAB.subs(mySubs1).subs(mySubs2)
dCD__ans=dCD.subs(mySubs1).subs(mySubs2)
print("#問題2. dAB=",round(dAB__ans*10**6,2))
print("#問題2. dCD=",round(dCD__ans*10**6,1))
print('#############################################')
from sympy import *
NAB,NBC,NCD =symbols('NAB NBC NCD',real=True)
P1, P2      =symbols('P1  P2'     ,real=True)
RA, RB      =symbols('RA  RB'     ,real=True)
A , EE      =symbols('A   EE'     ,real=True)
A50         =symbols('A50'        ,real=True)
lAB,lBC,lCD =symbols('lAB lBC lCD',real=True)
RA ,RD      =symbols('RA  RD'     ,real=True)
myAsyuhen=Eq( RA +NAB   ,0)
myBsyuhen=Eq(-NAB+P1+NBC,0)
myCsyuhen=Eq(-NBC-P2+NCD,0)
myDsyuhen=Eq( NCD+RD    ,0)
dAB=NAB*lAB/(A*EE)
dBC=NBC*lBC/(A50*EE)
dCD=NCD*lCD/(A*EE)
myHenkei=Eq(dAB+dBC+dCD,0)
ans=solve([myAsyuhen,
           myBsyuhen,
           myCsyuhen,
           myDsyuhen,
           myHenkei],[RA,NAB,NBC,NCD,RD])
mySubs1={P1 :15 ,P2 :10,
         lAB:0.6,lBC:1.0,lCD:1.4,
         A  :500*10**(-6),EE:200*10**(+6),
         A50:500*10**(-6)*0.5}
RA__ans=ans[RA ].subs(mySubs1)
NAB_ans=ans[NAB].subs(mySubs1)
NBC_ans=ans[NBC].subs(mySubs1)
NCD_ans=ans[NCD].subs(mySubs1)
RD__ans=ans[RD ].subs(mySubs1)
print("#問題3. RA =",round(RA__ans,2))
print("#問題3. RD =",round(RD__ans,2))
#問題1. RA = -7.33
#問題1. RD = -2.33
#問題2. NAB= 7.33
#問題2. NBC= -7.67
#問題2. NCD= 2.33
#問題2. dAB= 44.0000000000000
#問題2. dCD= 32.7
#############################################
#問題3. RA = -9.25
#問題3. RD = -4.25

・ver0.2
・2.33 RD=-4.25 でした。問題3. -9.25 2.33
・round関数が、44.0000000000000でした。
・辞書形式で、最後まで。大変かも？

# ver0.2
# ver0.2 マトリックス変位法で(直列構造)
from sympy import *
var('fx1,fx2,fx3,fx4,k1,k2,k3,u1,u2,u3,u4')
var('P1,P2,lAB,lBC,lCD,A,E', positive=True)
f=Matrix([[fx1],                   \
          [fx2],                   \
          [fx3],                   \
          [fx4]])
K=Matrix([[ k1,-k1  ,     0,  0],  \
          [-k1,k1+k2,-k2   ,  0],  \
          [  0,  -k2, k2+k3,-k3],  \
          [  0,    0,   -k3 ,k3]])
u=Matrix([[ u1],                   \
          [ u2],                   \
          [ u3],                   \
          [ u4]])
dic0=solve(Eq(f,K*u),[fx1,fx4,u2,u3])
rep1={fx2: P1,fx3:-P2,k1 :A*E/lAB,k2 :A*E/lBC,k3 :A*E/lCD,u1:0,u4:0}
rep2={ P1: 15, P2: 10,lAB:0.6    ,lBC:1.0    ,lCD:1.4}
rep3={ A :500, E :200*10**(-6)}
dic1={k:v.subs(rep1) for k,v in dic0.items()}      #;print("#",dic1)
dic2={k:v.subs(rep2) for k,v in dic1.items()}      #;print("#",dic2)
dic3={k:v.subs(rep3) for k,v in dic2.items()}      #;print("#",dic3)
print("# 問題1.",round( dic2[fx1],2)                              ,round( dic2[fx4],2))  
print("# 問題2.",round(-dic2[fx1],2),round(-dic2[fx1]-rep2[ P1],2),round(-dic2[fx4],2))  
print("# 　〃  ",round( dic3[ u2],2),round( dic3[ u3]          ,2))  
print() 
rep4={fx2: P1,fx3:-P2,k1 :A*E/lAB,k2 :Rational(1,2)*A*E/lBC,k3 :A*E/lCD,u1:0,u4:0}
dic4={k:v.subs(rep4).subs(rep2) for k,v in dic0.items()}      #;print("#",dic1)
print("# ??? 2.33 RD=-4.25 になるそうです。問題3.",round( dic4[fx1],2)                              ,round( dic2[fx4],2))  
# 問題1. -7.33 2.33
# 問題2. 7.33 -7.67 -2.33
# 　〃   44.0000000000000 -32.67
# ??? 2.33 RD=-4.25 でした。問題3. -9.25 2.33

・ver0.3

# ver0.3
from sympy import *

# 記号定義
var('P l AE')
u_syms    = symbols('u1:5')  # u1, u2, u3, u4
f_syms    = symbols('f1:5')  # f1, f2, f3, f4
k1, k2,k3 = symbols('k1 k2 k3')

def calc_member_forces(f_nodes):
    """
    f_nodes: 各節点にかかる外力（スカラー）のリスト（長さ n）

    戻り値:
    forces: 各部材に発生する内力（長さ n-1 のリスト）
            引張: 正、圧縮: 負
    """
    n = len(f_nodes)  # 節点数から自動判定
    if n < 2:
        raise ValueError("節点数は2以上である必要があります")

    # 部材数 = n - 1、本数分の記号を定義 M1, M2, ..., M_{n-1}
    M = symbols(f'M1:{n}')  # インデックスは 0〜n-2

    equations = []
    for i in range(n):
        eq = 0
        if i > 0:
            eq += M[i - 1]  # 左の部材の力
        if i < n - 1:
            eq -= M[i]      # 右の部材の力
        eq -= f_nodes[i]    # 外力
        equations.append(Eq(eq, 0))  # ∑F = 0

    # 解を求める
    sol = solve(equations, M)

    # 結果を順番に整形
    forces = [sol[m] for m in M]
    return forces

# 行列埋め込み関数
def myEmbed_matrix(big, small, row, col):
    r, c = small.shape
    big[row:row + r, col:col + c] = small
    return big

f  = Matrix(f_syms)
K1 = Matrix([[k1,-k1],[-k1,k1]])
K2 = Matrix([[k2,-k2],[-k2,k2]])
K3 = Matrix([[k3,-k3],[-k3,k3]])
K  = myEmbed_matrix(zeros(4, 4), K1, 0, 0)  \
    +myEmbed_matrix(zeros(4, 4), K2, 1, 1)  \
    +myEmbed_matrix(zeros(4, 4), K3, 2, 2)
u  = Matrix(u_syms)
# ↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
var('A,E')
var('l1,l2,l3')
var('P1,P2,P3,P4')
var('u1,u2,u3,u4')
var('f1,f2,f3,f4')
rep = {k1: A*E/l1,k2: A*E/l2,k3:A*E/l3,                 \
       f1: -P1   ,f2: P2    ,f3:-P3   ,f4:-(-P1+P2-P3)  \
      }
sol = solve([Eq(f,K*u).subs(rep),Eq(u1,0)],u)
# ----------------------------------------------------------------------------------------------
rep_val={l1: 0.6,l2: 1.0,l3: 1.4, \
         P1:20.0,P2:15.0,P3:10.0, \
         A:500*10**(-3),E:200
        }
f1=(K*u).subs(rep).subs(sol).subs(rep_val)
u1=u              .subs(sol).subs(rep_val)
n1=calc_member_forces(f1)
print("# 問題1",format(f1[3],'.0f'),format(n1[0],'.0f'),format(n1[1],'.0f'),format(n1[2],'.0f'))
print("# 問題2",format(u1[3],'.2f'))
# ----------------------------------------------------------------------------------------------
rep_val=rep_val | {P2:35.0}
f3=((K*u).subs(rep).subs(sol)).subs(rep_val)
u3=(u              .subs(sol)).subs(rep_val)
n3=calc_member_forces(f3)
print("# 問題3",format(f3[3],'.0f'),format(n3[0],'.0f'),format(n3[1],'.0f'),format(n3[2],'.0f'), \
               format(u3[3],'.1f'))
# ↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑↑ 
# 問題1 15 20 5 15
# 問題2 0.38
# 問題3 -5 20 -15 -5 -0.1

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

(テンプレート)

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

ChatGPT先生より

Python の round() 関数は、数値の丸め処理は行いますが、表示桁数を保証するものではありません。
formatted = format(value, '.3f')
print(f"{value:.3f}") # 15.000

参考文献

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

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