片持ち梁の主応力や方向をプロットするコード

cantilever-stress-plot.jl

# 片持ち梁について、以下の計算を行う。
# 曲げモーメント、せん断力、主応力、主応力の角度
# 最終的に、各数値をプロットした図を出力する。
# x=0を固定端とする

using PyPlot
using LinearAlgebra
using QuadGK

W = 100 #荷重
b = 10 #断面幅
h = 50  #高さ
L = 200 #長さ

x = range(0,L,length=5) # 梁の長さ方向に、要素5個分の配列を作る
y = range(-h/2,h/2,length=5) # 断面の高さ方向に、要素5個分の配列を作る
X = repeat(reshape(x, 1, :), length(y), 1) #xを1行に並べて、yの要素の数だけ繰り返す配列
Y = repeat(y, 1, length(x))                #yを1列に並べて、xの要素の数だけ繰り返す配列

#自由端からの距離を求める
l(x)=L-x
#モーメントを求める関数
M(x)=-W*l(x)
#長方形の断面二次モーメント
I = (b*h^3)/12 
#z軸に関する断面一次モーメント
Q = quadgk.(Y -> b*Y,Y,h/2, rtol=1e-8)
#曲げ応力を計算する
sigma = @. (M(X)*Y)/I
#せん断応力を計算する
g = [(W*Q[i][1])/(I*b) for i in 1:length(Y)]
#せん断応力の配列のカタチを整える
tau = reshape(g, length(y), length(x)) 
#sigmaとtauの値を、行列の要素に入れる
B=[[sigma[i] tau[i] 0;tau[i] 0 0;0 0 0] for i=1:length(sigma)]
#各部位の3軸分の主応力の値を取得する
C=eigvals.(B)
#各部位の3軸分の方向余弦を取得する
D=eigvecs.(B)

# 最大主応力のIndexを取得する
E=argmax.(C)
E1=argmin.(C)

#最大主応力に対応する方向余弦を取得する
F=[[D[i][1,E[i]]] for i=1:length(E)]

#最小主応力に対応する方向余弦を取得する
F1=[[D[i][1,E1[i]]] for i=1:length(E)]
F2=[[D[i][2,E1[i]]] for i=1:length(E)]

#矢線の大きさ
r = 3

#最大主応力のx,y要素を作成
u = r*F
v = r*broadcast.(sin,(broadcast.(acos,F)))

#最大主応力が0の場合は、要素をゼロにする。
for i in 1:length(E)
    u[i]=ifelse(C[i][E[i]]==0, u[i] .* 0,u[i].*1)
    v[i]=ifelse(C[i][E[i]]==0, v[i] .* 0,v[i].*1)
end

#最小主応力のx,y要素を作成
umin = r*F1
vmin = r*F2

#最小主応力のx,y要素が第4象限に存在したら、第二象限に移動
#最小主応力が0だったら、x,y要素を0にする
for i in 1:length(E)
    umin[i]=ifelse(umin[i][1]>0, umin[i] * -1, umin[i]*1)
    vmin[i]=ifelse(vmin[i][1]<0, vmin[i] * -1, vmin[i]*1)
    umin[i]=ifelse(C[i][E1[i]]==0, umin[i] .* 0,umin[i].*1)
    vmin[i]=ifelse(C[i][E1[i]]==0, vmin[i] .* 0,vmin[i].*1)
end


#Pyplotへ渡すために成形する（最大主応力）
#Matrixにしないと、エラーが出てしまうので
u1 = vcat(u...)
v1 = vcat(v...)

u2 =reshape(u1, length(y), length(x))
v2 =reshape(v1, length(y), length(x))

#Pyplotへ渡すために成形する（最小主応力）
#Matrixにしないと、エラーが出てしまうので
u3 = vcat(umin...)
v3 = vcat(vmin...)

u4 =reshape(u3, length(y), length(x))
v4 =reshape(v3, length(y), length(x))


X1=reshape(X, length(y), length(x))
Y1=reshape(Y, length(y), length(x))

## デバッグ用
display(sigma)
display(tau)
display(map(M,x))


## Plotする
## Pyplotで表示するコード。

fig = figure("pyplot_streamplot",figsize=(14.4,19.2),dpi=100)
font= "BIZ UDGothic"
tc = "#0000CC" #Tensionのcolor
cc = "#CC0066" #Compressionのcolor

#主応力の方向を流線として表示する
subplot(411)
streamplot(X1,Y1,u2,v2, color=tc,linewidth=2)
streamplot(X1,Y1,u4,v4, color=cc,linewidth=2)
gca().invert_yaxis()
gca().set_xlabel("x (mm)", fontsize=15)
gca().set_ylabel("y (mm)", fontsize=15)
gca().set_title("主応力の方向(流線)", fontname = font, fontsize=20,y=1.03)
xticks(fontsize=15)
yticks(fontsize=15)


#主応力の方向を場の矢線として表示する(引張/圧縮)
subplot(412)
quiver(X1, Y1, u2, -v2,units="xy", headwidth=5,headlength=5,width=0.5, color=tc, label="Tension")
quiver(X1, Y1, u4, -v4,units="xy", headwidth=5,headlength=5,width=0.5, color=cc,label="Compression")
gca().invert_yaxis()
gca().set_xlabel("x (mm)", fontsize=15)
gca().set_ylabel("y (mm)", fontsize=15)
xticks(fontsize=15)
yticks(fontsize=15)
gca().set_title("主応力の方向(引張/圧縮)", fontname = font, fontsize=20,y=1.03)
gca().legend(fontsize=15,framealpha=0.5)

#主応力の方向を場の矢線として表示する(引張)
subplot(413)
quiver(X1, Y1, u2, -v2,units="xy",headwidth=5,headlength=5,width=0.5, color=tc,label="Tension")
gca().invert_yaxis()
gca().set_xlabel("x (mm)", fontsize=15)
gca().set_ylabel("y (mm)", fontsize=15)
xticks(fontsize=15)
yticks(fontsize=15)
gca().set_title("主応力の方向(引張)", fontname = font, fontsize=20,y=1.03)
gca().legend(fontsize=15,framealpha=0.5) #凡例を表示する

#主応力の方向を場の矢線として表示する(圧縮)
subplot(414)
quiver(X1, Y1, u4, -v4,units="xy",headwidth=5,headlength=5,width=0.5,color=cc,label="Compression")
gca().invert_yaxis()
gca().set_xlabel("x (mm)", fontsize=15)
gca().set_ylabel("y (mm)", fontsize=15)
xticks(fontsize=15)
yticks(fontsize=15)
gca().set_title("主応力の方向(圧縮)", fontname = font, fontsize=20,y=1.03)
gca().legend(fontsize=15,framealpha=0.5) #凡例を表示する

#サブプロットのレイアウトを整える
tight_layout(pad=0.4, w_pad=0.5, h_pad=5)

#plt.savefig("e:\\test1-1.png") #PNGファイルを出力する場合はコメントアウト解除

###先のプロットとは別の画像として描画
fig2 = figure("pyplot_streamplot2",figsize=(14.4,19.2),dpi=100)
font= "BIZ UDGothic"

#曲げモーメント図
subplot(411)
p=plot(x,map(M,x), lw=2, c="blue")

gca().set_title("曲げモーメント図", fontname = font, fontsize=20,y=1.1)
gca().set_xlim(0, L)

gca().xaxis.set_ticks_position("top")
gca().yaxis.set_ticks_position("left")

gca().xaxis.set_label_position("top")

gca().spines["top"].set_position(("data",0)) # Y軸位置を調整
xticks(fontsize=15)
yticks(fontsize=15)

#y軸の値を対数表示にする。linthreshはプロットが線形になる範囲を指定
gca().set_yscale("symlog", linthresh=1e5)
gca().set_ylim(-2e4,0) #y軸の範囲を指定

gca().set_xlabel("x (mm)",labelpad=-40,fontsize=15)
gca().set_ylabel("M (Nmm)", labelpad=-40,fontsize=15)


#せん断力図
subplot(412)
p=plot(x,repeat([W],length(x),1), lw=2, c="blue")

gca().set_title("せん断力図", fontname = font, fontsize=20,y=1.1)
xticks(x,fontsize=15)
yticks(fontsize=15)
gca().set_xlim(0, 200) #x軸の範囲を指定
gca().set_ylim(0, W+50) #y軸の範囲を指定

gca().set_xlabel("x (mm)", fontsize=15)
gca().set_ylabel("F (N)", labelpad=-1, fontsize=15)


#曲げ応力の分布
subplot(413)
cs=contourf(X1,Y1, sigma)
xticks(x,fontsize=15)
yticks(y,fontsize=15)
gca().invert_yaxis()
clb=plt.colorbar(fraction=0.05, pad=0.04) #カラーバーの位置を指定
clb.ax.text(0,7,"(N/mm^2)", fontsize=10, rotation=0,ha="center") #カラーバー付近に単位を記載

gca().set_xlabel("x (mm)", fontsize=15)
gca().set_ylabel("y (mm)", fontsize=15)
gca().set_title("曲げ応力の分布", fontname = font, fontsize=20,y=1.03)


#せん断応力の分布
subplot(414)
cs2=contourf(X1,Y1, tau)
xticks(x,fontsize=15)
yticks(y,fontsize=15)
gca().invert_yaxis()
gca().set_xlabel("x (mm)", fontsize=15)
gca().set_ylabel("y (mm)", fontsize=15)
gca().set_title("せん断応力の分布", fontname = font ,fontsize=20,y=1.03)
clb2=plt.colorbar(fraction=0.05, pad=0.04) #カラーバーの位置を指定
clb2.ax.text(0,0.34,"(N/mm^2)", fontsize=10, rotation=0,ha="center") #カラーバー付近に単位を記載
#サブプロットのレイアウトを整える
tight_layout(pad=0.4, w_pad=0.5, h_pad=3)
#plt.savefig("e:\\test2.png") #PNGファイルを出力する場合はコメントアウト解除