Juliaで、主応力の方向を図示する方法。力の流れが見えて面白いです

解説記事

プログラミング言語のJuliaで主応力の方向を図示してみたので、方法をレポートします。
FEAソフトは使わないので、比較的容易に実現できます。 道具として丸暗記しがちな材料力学の公式ですが、視覚化すると印象に残ります。

スポンサーリンク

はじめに

この記事ではプログラミング言語のJuliaを使って、主応力の方向を図示する方法を紹介します。
FEAソフトの使用が方法として思い浮かびますが、今回は、材料力学の公式を使って計算しました。
公式は便利ですが、手計算で一度に複数の位置の計算を行うのは大変です。
そこでJuliaを使って計算を行い、それらをグラフ上にプロットし視覚化すれば、簡便に全体像が見えるのではないかと思い試しました。
課題は先端に集中荷重を受ける片持ち梁とし、断面の位置ごとに主応力と法線角度を計算して「場」をプロットします。

完成図

先端に集中荷重を受ける片持ち梁を課題にしました。
計算してプロットした結果は以下の通りです。課題の条件など詳細は、後程説明します。

片持ち梁の例題
例題の仕様図

主応力の方向を図示

主応力の方向を流線として表示
主応力の方向を矢線として表示
主応力の方向を矢線として表示

曲げモーメント図とせん断力図

曲げ応力、せん断応力の分布図

使用ツール

Julia 1.5.3 を、Jupyter Notebook(IJulia)上で使用しました。

以下は、使用したツールの一覧です。

  • Juliaは無償のプログラミング言語。 
  • Jupyter Notebook(IJulia)は、Juliaで計算したりドキュメントを作成するための計算環境。ブラウザで動く。
  • Pyplot
    Python向けの視覚化ツールであるMatplotlibを、Juliaで使うためのモジュール。

IJuliaでPyplotを使っている様子

Ijuliaで作業している様子
Ijuliaで作業している様子

▼ 細かい関数の確認はJulia REPL(一行ずつ打ち込む端末)を使いました。

環境

OS:Windows 10 Pro 64-bit 20H2
ブラウザ:Chrome

導入方法

この記事の作業では、以下の3つのインストールが必要です。
まずは1.のJuliaを公式サイトからダウンロードしてインストールします。
2.と3.はJuliaのREPLという端末で操作して、インストールします。(決まったコマンドを打ち込むだけです)。

▼ 詳しくは、下記1.~3.のリンク先をご覧ください。

  1. Juliaは無償のプログラミング言語です。 
  2. Jupyter Notebook(IJulia)は、Juliaで計算したりドキュメントを作成するための計算環境です。ブラウザで動くので、書類を作る感覚で、対話的に計算したりメモをかいたりできます。
  3. Pyplot
    Python向けの視覚化ツールであるMatplotlibを、Juliaで使うためのモジュールです。


▼1.と2.の基本環境の導入については、別記事にまとめています。

環境を整える方法は何通りかあるので、ややこしいですね。
個人的な例ですが、私の場合は以下の順番で導入しました。

  1. 以前、Python向けの科学計算環境であるAnacondaをインストールしていた
    ※ この時点で、Jupyter NotebookやMatplotlibがインストールされている
  2. Julia言語をインストール
  3. REPL端末で Pkg.add("PyPlot") を実行し、JuliaでプロットするためのPyplotをインストール

「色々しているうちに、いつの間にかパッケージが入ってる」なんてこともあります。
だめもとでusing PyPlotとテストプロットを実行して、エラーが出るか確認すると手っ取り早いかもしれません。

課題の仕様

長さL 200 mm × 高さh 50 mm × 幅b 10 mmの板材の先端に、W = 100 N の荷重がかかる片持ち梁です。

片持ち梁の例題

計算式

梁の正面に計25個(=5×5)の点を設定し、それらの各点について以下4項目の数値を材料力学の公式を使って計算し、結果をプロットしました。

  1. 曲げモーメントとせん断力
  2. 曲げ応力
  3. せん断応力
  4. 主応力の大きさと法線の角度

最終的に求めたいのは4.の主応力とその方向ですが、その準備として1~3.を順番に計算します。

設定した格子点
今回の例題に設定した格子点の図

1.曲げモーメントとせん断力

曲げモーメントMとせん断力Fは、今回の課題ではそれぞれ以下の式で計算ができます。

\[\text{曲げモーメント}\quad M=-W(L-x) \mathrm{\: [Nmm]}\]

\[\text{せん断力}\quad F=W \mathrm{\: [N]} \]

材料力学でおなじみのBMD,SFDを描くとこんなかんじです。

曲げモーメント図とせん断力図
BMDとSFD

2.曲げ応力

曲げ応力 \(\sigma_{x}\) は、以下の式で計算します。

\[\sigma_{x}=\frac{My}{I} \mathrm{\: [N/mm^2]}\]

ここで、

  • y:y軸上の位置(図心からの距離)。例えば、中央ならy=0です。
  • I:断面二次モーメント。今回の長方形断面の場合、式は \(bh^3/12\) です。
片持ち梁の例題の模式図
片持ち梁の例題の模式図

3.せん断応力

せん断応力 \(\tau_{xy}\) は、以下の式で計算します。

\[\tau_{xy}=\frac{FQ}{bI} \mathrm{\: [N/mm^2]}\]

ここで、

  • F:せん断力 [N]
  • b:断面幅 [mm]
  • I:断面二次モーメント
    今回の長方形断面の場合、式は \(bh^3/12\) です。
  • Q:中立軸(z軸)に関する断面一次モーメント(式は以下参照)

Qはz軸に関する断面一次モーメント(面積モーメントの総和)で、今回の長方形断面の場合、以下の式で計算できます。

\[Q = \int_{e_1}^\frac{h}{2} y\mathrm {b}\,\mathrm{d}y\tag{3.1}\]

ここで、

  • b:断面幅 [mm]
  • h:断面の高さ [mm]
  • e1:y軸上の位置(図心からの距離) [mm]

手計算の場合は(式3.1)を以下のように展開すると計算できますが、

\[\Bigg[b\frac{1}{2}y^2\Bigg]^{\frac{h}{2}}_{e_1}=\: \frac{b}{2} \left((^h/_2)^2-{e_1}^2\right)\]

Juliaでは、QuadGKという追加パッケージを入れると定積分ができます。

GitHub - JuliaMath/QuadGK.jl: adaptive 1d numerical Gauss–Kronrod integration in Julia
adaptive 1d numerical Gauss–Kronrod integration in Julia - GitHub - JuliaMath/QuadGK.jl: adaptive 1d numerical Gauss–Kronrod integration in Julia

試しに、以下の式に今回の条件(幅b=10、高さh=50)を代入してJulia REPLで計算してみます。
面積モーメントの総和なので、0になるはずです。

\[Q = \int_{-\frac{h}{2}}^\frac{h}{2} \mathrm {10}y\,\mathrm{d}y\]

julia> using QuadGK

julia> h=50
50

julia> b=10
10

julia> quadgk(y -> b*y,-h/2,h/2, rtol=1e-8)
(0.0, 0.0)

無事0が返ってきました。(0.0, 0.0)の内の二つ目の要素は、誤差推定値(error estimate)とのことです(QuadGKのReadmeより)。

4.主応力の大きさと法線の角度

これまでに求めた曲げ応力とせん断応力から、主応力の大きさと法線(主軸)の角度を計算します。

主応力とは?

主応力とは、せん断応力が0になる面での垂直応力のことです。
例えば、チョークをねじると斜めに亀裂が入って割れますが、これは軸から45°傾いた線を法線とする面に生じた主応力(引張応力)によるものだと考えられます。当記事では計算で求めていますが、モールの応力円という図を描いて把握することもできます。

▼ ご参考に、別記事に簡単な説明を書いています。
参考 モールの応力円とは?意味と書き方を、計算をすっとばして説明するよ【超初心者向け】

以下の手順で進めます。

  1. 応力の9成分を行列表示する
  2. 1.から固有値(主応力)と固有ベクトル(方向余弦)を求める
  3. 試しにJuliaで計算してみる

1.行列として表示する

まず、応力の9成分を行列として表示します。

\[[\sigma] = \begin{bmatrix}
\sigma_{x} & \tau_{yx} & \tau_{zx} \\[0.3em]
\tau_{xy} & \sigma_{y} & \tau_{zy} \\[0.3em]
\tau_{xz} & \tau_{yz} & \sigma_{z}
\end{bmatrix}
\]

2. 固有値と固有ベクトルを求める

主応力は行列 \([\sigma]\) の固有値、その軸の向きは固有ベクトルを計算すると求まります。
結果は3軸分が出てきます。

次式だと、σ1が固有値で主応力の値です。[l,m,n]はσ1に対応する固有ベクトルで、各要素は方向余弦です。

\[\begin{bmatrix}
\sigma_{x} & \tau_{yx} & \tau_{zx} \\[0.3em]
\tau_{xy} & \sigma_{y} & \tau_{zy} \\[0.3em]
\tau_{xz} & \tau_{yz} & \sigma_{z}
\end{bmatrix}
\begin{bmatrix}
\ l_{1} \\[0.3em]
\ m_{1} \\[0.3em]
\ n_{1}
\end{bmatrix}
= \sigma_{1}
\begin{bmatrix}
\ l_{1} \\[0.3em]
\ m_{1} \\[0.3em]
\ n_{1}
\end{bmatrix}
\]

3. 試しにJuliaで計算してみると?

ここで、試しにJuliaで具体的な数値を入れて、簡単な条件で計算してみます。

例題

せん断応力 \(\tau_{xy}=100 \, Mpa\) のみが発生する状態から、主応力とその向きを計算する。

せん断応力のみが発生する状態の図

この場合、計算しなくても(モールの応力円を思い出すと)「最大主応力は100 Mpa、向きは元のx軸から45°傾いた軸を法線とする面」だと分かります。

その通りの結果になるか、確認してみます。

まず、9成分を行列表示するとこうなります。Juliaでの計算の都合上、Aという名前にしました。

\[A = \begin{bmatrix}
0 & 100 & 0 \\[0.3em]
100 & 0 & 0 \\[0.3em]
0 & 0 & 0
\end{bmatrix}
\]

Julia Replを立ち上げ、まずは線形代数の機能を使うので using LinearAlgebra と打ち込みます。
そして、固有値は eigvals(A) 、固有ベクトルはeigvecs(A)を実行すると求まります。

REPLで実行すると、以下のような結果になります。

julia> using LinearAlgebra
julia> A=[0 100 0;100 0 0;0 0 0];

julia> eigvals(A) #固有値を求める
 3-element Array{Float64,1}:
 -100.0
    0.0
  100.0

julia> eigvecs(A) #固有ベクトルを求める
 3×3 Array{Float64,2}:
  0.707107   0.0  0.707107
  -0.707107  0.0  0.707107
  0.0        1.0  0.0

この結果から「主応力は最小-100と最大100 Mpa。主軸角度は最大主応力の場合X軸から45°」だと分かります。

結果には xyzの3軸分が含まれていて、固有値が3つの要素の配列、固有ベクトルが3×3の配列として出力されています。

例えば eigvals(A) で求めた固有値のうち、正の要素の100は、3行目に存在します。
そして、それに対応する固有ベクトルは eigvecs(A) の3列目です。
[0.707107, 0.707107 ,0.0]とあります。(要素は方向余弦なので、cosθの値です。)
度に換算すると、[45°,45°,90°]となります。(90°はz軸からの角度です。)

ここで、 eigvecs(A) の1列目 [0.707107, -0.707107 ,0.0] と3列目 [0.707107, 0.707107 ,0.0] の違いに注意する必要があります。符号の組み合わせによって、どの象限に軸があるかが分かります。

Juliaでの計算結果と主軸の象限の対応
Juliaでの計算結果と主軸の象限の対応

最初は扱いを間違えてプロット結果がアッチコッチ向いてしまい、大混乱でした。

また、モールの応力円を描いて解いた図と見比べると、結果が一致しています。
これで、個々の点の主応力とその方向をJuliaで計算できることが確認できました。
あとはコードを書いて、この計算を梁全体に適用します。

今回の課題で、モールの応力円を描いてみると?

コードを書く

最終工程です。あとは格子状に点を設けてその各点で計算を実行し、結果をプロットします。

コードは、以下の通りです。※ 表示されていない場合は、Githubで見られます。

▼ Ijuliaで作業している様子

Ijuliaで作業している様子
Ijuliaで作業している様子

プロットツールについて

プロットには、Pyplotを使いました。これはPython向けのデータ可視化ライブラリのMatplotlibを、juliaで使うためのモジュールです。

使う際は、using PyPlot という文言で呼び出します。
インストール方法は、下記のGithubのReadmeや当記事の冒頭をご覧ください。

GitHub - JuliaPy/PyPlot.jl: Plotting for Julia based on matplotlib.pyplot
Plotting for Julia based on matplotlib.pyplot. Contribute to JuliaPy/PyPlot.jl development by creating an account on GitHub.

ベクトル場の図示機能を備えており、そして流線まで描けるので今回の目的をスムーズに実現できました。
使い方は、Matplotlibのドキュメントが詳しいです。

ドキュメントのサンプルをそのままJuliaで実行するとエラーがでることがあります。
例えばパラメータの囲みだと、python向けではシングルクォーテーション(”)ですが、Juliaではダブル(“”)にしないとエラーが出ました。
そこで、初めて使われる方はチュートリアルでまずはシンプルなグラフを描いて、動作を確認すると良いと思います。

Plots経由では、流線が描けませんでした

少しややこしいですが、JuliaにはPlotsというツールセットがあります。
これはグラフや視覚化のためのもので、GRやPyplotなど様々なプロットツールを共通の文法で簡単に使えるようになります。
しかしPlots経由でPyplotを呼び出すと、Streamplotが使えませんでした。Matplotlibとしての機能は制限されてしまう様です。

結果のまとめ

計算した結果を、以下にまとめます。

条件

  • 長さ200 mm × 高さ50 mm × 幅10 mmの板材の先端に、W = 100 N の荷重がかかる片持ち梁
  • 使用した計算式とコードは、前項に記載
片持ち梁の例題

主応力の方向を図示

最大主応力(引張)、最小主応力(圧縮)の向きを流線と矢線でプロットしました。
応力の流れや最大最小主応力の直交が確認できます。

主応力の方向を流線として表示
主応力の方向を矢線として表示
主応力の方向を矢線として表示

曲げモーメント図とせん断力図

主応力を求める過程で計算したので、グラフにしました。

曲げ応力、せん断応力の分布図

こちらも主応力を求める過程で計算したので、等高線機能を使ってプロットしました。
計算した点の数が少ないので粗いですが、曲げ応力は梁の上側が正の値で引張応力、下側は圧縮応力が発生しているのを確認できます。

計算値は以下の通りです。

julia> display(sigma) #曲げ応力
5×5 Array{Float64,2}:
  4.8   3.6   2.4   1.2   0.0
  2.4   1.8   1.2   0.6   0.0
 -0.0  -0.0  -0.0  -0.0  -0.0
 -2.4  -1.8  -1.2  -0.6  -0.0
 -4.8  -3.6  -2.4  -1.2  -0.0

julia> display(tau) #せん断応力
5×5 Array{Float64,2}:
 0.0    0.0    0.0    0.0    0.0
 0.225  0.225  0.225  0.225  0.225
 0.3    0.3    0.3    0.3    0.3
 0.225  0.225  0.225  0.225  0.225
 0.0    0.0    0.0    0.0    0.0

参考文献

書籍

  1. 竹園茂男, ほか. 「2.1.4 主応力と主せん断応力」. 弾性力学入門: 基礎理論から数値解法まで, 森北出版, 2007年, pp. 19–26.
  2. 日本機械学会. 「3・8 せん断によるはりの応力とたわみ」. 機械工学便覧 基礎編 α3 材料力学, 初版, 日本機械学会, 2005年, pp. 30–31.
  3. 中島正貴. 「6.5 はりのせん断応力」. 機械系教科書シリーズ 19 材料力学, コロナ社, 2005年, pp. 102–106.
  4. 石村園子. 「3.3 固有値と固有ベクトル」. やさしく学べる線形代数, 2000年, pp. 132–38.

本について補足

参考文献のうち、特に気に入っている本の感想です。

▲ 材料力学から更に先へ勉強を進めたい際に、足掛かりとなる本です。
弾性力学になると大量の数式や記号がわんさかでてきて戦意喪失しがちなのですが、この本は式の展開が丁寧且つ図も豊富で、心が折れずに読み進められています。
具体的な値を使った易しい例題が随所にあるのも、良かったです。一般式だけだと実際にどう式を扱うのか確証を持てず不安になるので、都度理解を確認出来て助かりました。

▲ 私が持っているのは初版で、資格取得時にとてもお世話になった本です。今も折に触れて読んでいます。数学の勉強が猛烈に不足している頃に読むと辛かったのですが、微積分を勉強して読みなおすと抵抗なく読めるようになりました。
今回の記事の関連だと、はりに生じる応力(曲げ、せん断、断面一次、二次モーメント等)の解説が、導出が段階的で丁寧でとても理解しやすく有難いです。

線形代数については、放送大学で講座を履修した際に使った市販本を別記事で紹介しています。初めて勉強される方は参考になるかもしれません。

Webサイト

固有値と固有ベクトル

5.1: Eigenvalues and Eigenvectors
In this section, we define eigenvalues and eigenvectors. These form the most important facet of the structure theory of square matrices. As such, eigenvalues an...

Webでは上記のLibretextsという教科書サイトをよく読んでいます。LibreTextsは、オンラインのオープンアクセス教科書プロジェクトで、様々な単元の解説を無料で閲覧できます。
レイアウトもきれいで見てて疲れにくいので、スムーズに読み進められます。

Juliaの資料

プロット 関連

Matplotlibは本来Python向けですが、JuliaではPyplotというモジュールを呼び出せば使用できます。導入にあたっては、上記のGithubのReadmeを読みました。

Matlotlibの簡潔でわかりやすいチュートリアルです。Matlotlibは機能が豊富なので、まず何をどうすれば良いのか最初はとまどいました。
そこで、上記のチュートリアルでまずはシンプルなグラフを描くところから始めました。まずは抵抗を減らして、馴染むことができます。

さいごに

本格的な解析はFEAソフトを使うのが普通だと思いますが、簡単な形状で明確に公式が存在していればJuliaでさっと計算して、手軽に視覚化できるのが便利で面白いです。

今回は直接呼び出す形での Pyplotを初めて使いました。Plots経由のバックエンドとして使うと制限があったので少々不便に感じていましたが、Matplotlibとしての機能をJuliaから使えるのはとても快適でした。

タイトルとURLをコピーしました