材料力学_(4)_モールの応力円（英:Mohr's stress circle）「傾斜断面の応力」をsympyでやってみたい。

・記事タイトル変更予定です。 材料力学_(３)は、作業中です。
・間違いをご指摘いただけると幸いです。勉強中です。

ネット(youtube)

福山大学工学部機械システム工学科 様
・＞2階のテンソル(2:48)
・私は、行列による座標変換を理解していません。申し訳ありません。
途中3:44から (0:00〜6:06)

ネット(pdf)

・「材料力学演習(20221021)一括(ver.3.3).pdf」第5章　組み合わせ応力(129頁〜154頁,P134〜159) 基本事項３(モールの応力円) 134頁,P139
？？？タブレットだと、開く事ができないかもしれません。

「材料力学演習(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

文献

・材料力学 (JSMEテキストシリーズ) 大型本 – 2007/8/30 135頁

・JSMEテキストシリーズ　材料力学

発行日：2007年8月
体裁：A4判/並製/2色刷り/208頁
目次(p2)

・演習材料力学 (JSMEテキストシリーズ) 大型本 – 2010/11/1

Mathematicaで

勉強中

sympyで(classCircleで)

・できました。できると思いませんでした。
・記号の式の順番まで、対応していません。
　自作関数で、文字列に変換したらできるような気もします。
class sympy.geometry.ellipse.Circle(*args, **kwargs)
https://docs.sympy.org/latest/modules/geometry/ellipses.html#sympy.geometry.ellipse.Circle

from sympy import *
var('σ,τ,A,A0,θ,τxy,σx,σy',real=True)
eqσ=Eq(σ, τxy*sin(2*θ) + (σx - σy)*cos(2*θ)/2 + (σx + σy)/2) ;print("# ",eqσ)
eqτ=Eq(τ,-τxy*cos(2*θ) + (σx - σy)*sin(2*θ)/2)               ;print("# ",eqτ)  ;print()

eqσ=Eq(eqσ.lhs-(σx + σy)/2,eqσ.rhs-(σx + σy)/2)              #;print("# ",eqσ)
eqτ=Eq(eqτ.lhs,eqτ.rhs)                                      #;print("# ",eqτ)  ;print()

var('x,y',real=True)
f=(eqσ.lhs**2+eqτ.lhs**2)-(eqσ.rhs**2+eqτ.rhs**2)            #; print("# ",f)
f=f.subs({σ:x,τ:y})                                          #; print("# ",f)
ci=Circle(f)                                                 #; print(c)
cc=Point(ci.center.x.factor(),ci.center.y)                   #; print(cc)
po=poly((ci.radius**2).simplify() ,τxy)                      #; print("#",p)
ci=Circle(cc,sqrt( factor(po.nth(0))+τxy**2) )                ; print("# ",ci)
eq=Eq((σ-ci.center.x)**2+τxy**2,ci.radius**2)                 ; print("# ",eq)
#  Eq(σ, σx/2 + σy/2 + τxy*sin(2*θ) + (σx - σy)*cos(2*θ)/2)
#  Eq(τ, -τxy*cos(2*θ) + (σx - σy)*sin(2*θ)/2)

#  Circle(Point2D((σx + σy)/2, 0), sqrt(τxy**2 + (σx - σy)**2/4))
#  Eq(τxy**2 + (σ - (σx + σy)/2)**2, τxy**2 + (σx - σy)**2/4)

sympyで(テキストより)

・多項式の分割と約分のため、polyにしました。
nth(*N)
Returns the n-th coefficient of f where N are the exponents of the generators in the term of interest.
https://docs.sympy.org/latest/modules/polys/reference.html#sympy.polys.polytools.Poly.nth

# ver0.1
from sympy import *
var('σ,τ,A,A0,θ,τxy,σx,σy',real=True)
eqσ=Eq(σ,τxy*sin(2*θ) + (σx - σy)*cos(2*θ)/2 + (σx + σy)/2)  ;print("# ",eqσ)
eqτ=Eq(τ,-τxy*cos(2*θ) + (σx - σy)*sin(2*θ)/2)               ;print("# ",eqτ)  ;print()

eqσ=Eq(eqσ.lhs-(σx + σy)/2,eqσ.rhs-(σx + σy)/2)              #;print("# ",eqσ)
eqτ=Eq(eqτ.lhs,eqτ.rhs)                                      #;print("# ",eqτ)  ;print()

var('σxMσy',real=True)
rep={σx - σy:σxMσy}
eqσ=eqσ.subs(rep)                                              # ;print("# ",eqσ)
eqτ=eqτ.subs(rep)                                              # ;print("# ",eqτ)  ;print()

g=(eqσ.rhs**2+eqτ.rhs**2).expand()                             #;print("# ",f)
p=poly(g,τxy,σxMσy)                                            #;print("#",p)
g=p.nth(2,0).simplify()*τxy**2+p.nth(0,2).simplify()*σxMσy**2  #;print("# ",g)
g=g.subs({σxMσy:σx - σy})

var('h,k',real=True)
f=(eqσ.lhs**2+eqτ.lhs**2).expand()                             #;print("# ",f)
p=poly(f,τ)                                                    #;print("# ",p)
f=p.nth(0)                                                     #;print("# ",f)
q=(σ-h)**2-k
ans=solve(q-f,[h,k])
f=q.subs({h:factor(ans[0][0]),k:ans[0][1]})                    #;print("# ",f)

eq=Eq(f,g)  ;print("# ",eq)  
#  Eq(σ, σx/2 + σy/2 + τxy*sin(2*θ) + (σx - σy)*cos(2*θ)/2)
#  Eq(τ, -τxy*cos(2*θ) + (σx - σy)*sin(2*θ)/2)
#  Eq((σ - (σx + σy)/2)**2, τxy**2 + (σx - σy)**2/4)

参考