3次元の主応力とその方向を求めてみた

どうも。最近ものづくりできてないせいで研究で書いたソースをネタにブログを書くハメになっているだんです。

今回は3次元各軸の垂直応力、せん断応力から主応力とその方向(方向余弦)を求めるプログラムを書いたので紹介します。研究の内容には全く関係ないしコードは全部オリジナルなので好きに使ってください。

こちらを参考にさせていただきました。

構造力学メモ（その１） 3次元の主応力とミーゼス応力｜Kou｜note

求め方

基本的には3×3行列の固有値と固有ベクトルを求める操作になります。

応力テンソルを以下のように表します。

このとき以下になる固有値 $λ$ が主応力で

以下のようになる固有ベクトル $x$ が主応力の方向余弦です

ただし参考にさせていただいたサイトによると偏差応力を使った方が簡単らしいので平均応力 $σ_m$ を垂直応力から引いた偏差応力テンソル $s$ の固有値を求めます。

このとき

$|s-λI|=-λ^3+J_2λ+J_3=0$

となります。 $J_2$ と $J_3$ は応力の不変量というやつで

$J_2=s_{xy}s_{yx}+s_{yz}s_{zy}+s_{zx}s_{xz}-s_{xx}s_{yy}-s_{yy}s_{zz}-s_{zz}s_{xx}$

$J_3=s_{xx}s_{yy}s_{zz}+s_{xy}s_{yz}s_{zx}+s_{xz}s_{zy}s_{yx}-s_{yz}s_{zy}s_{xx}-s_{zx}s_{xz}s_{yy}-s_{xy}s_{yx}s_{zz}$

で求められます(前述のサイト様が超参考になりました)。

三次方程式の解はいくつか求める方法があるらしいですが、コンピューターに解かせるのでビエトの解を使います。これは

$x^3-Ax-B=0$

の解が

$θ=\frac{1}{3} \arccos{(\frac{3\sqrt{3}B}{2\sqrt{A^3}})}$ として

$x=\frac{2}{3}\sqrt{3A}\cos{θ},\frac{2}{3}\sqrt{3A}\cos{(θ+240°)},\frac{2}{3}\sqrt{3A}\cos{(θ+120°)}$

となるというものです。

これで求めた偏差応力テンソルの固有値 $λ_1,λ_2,λ_3$ に平均応力を足して応力テンソルの固有値つまり主応力を得ます。

$σ_1=σ_m+λ_1\\σ_2=σ_m+λ_2\\σ_3=σ_m+λ_3$

続いて各固有値に対応する固有ベクトルを求めます。

これもいくつか方法がありますがここではガウス・ザイデル法を用いました。

ガウス＝ザイデル法 - Wikipedia

ざっくりいうと $A\vec{x}=\vec{b}$ を解くための方法です。こんな簡単な操作で解が求まるの素晴らしいですね。さすがガウス先生です(？)

ここでは固有値を求めたいので定義式どおり $A=σ-σ_nI,b=\vec{0}$ として3回計算します。(2回計算した時点で直交するベクトルを求めた方が早い気もする)(誤差だよ誤差)

注意点として最後に正規化しないと単位ベクトルになりません。

プログラム

C++で書きました。というか例によって最初C#で書いてブログ用に書き直しました。

他でも使えそうなガウス・ザイデル法のところだけ関数化してあります。引数の配列とdimを書き換えればn×n行列で解けるはずです(たぶん)。

#include "iostream"
#include "cmath"

void GaussSeidelMethod(double A[3][3],double b[3],double x[3]){
    const int dim = 3;

    double x_[dim];
    for(int i = 0;i < dim;i++)x_[i] = 1;
    
    for (int i = 0; i < 1000; i++){
        double _x[dim];
        for(int j = 0;j < dim;j++){
            _x[j] = b[j];
            for(int k = 0;k < j;k++){
                _x[j] -= A[j][k]*_x[k];
            }
            for(int k = j + 1;k < dim;k++){
                _x[j] -= A[j][k]*x_[k];
            }
            _x[j] /= A[j][j];
        }
        
        bool stop = true;
        double th = 0.000001;
        for(int j = 0;j < dim;j++){
            if(std::abs(_x[j] - x_[j]) > th){
                stop = false;
                break;
            }
        }

        if(stop){
            for(int j = 0;j < dim;j++){
                x[j] = _x[j];
            }
        }else{
            for(int j = 0;j < dim;j++){
                x_[j] = _x[j];
            }
        }
    }
}

int main(int argc, char* argv[])
{
    double sx;
    double sy;
    double sz;
    double txy;
    double tzx;
    double tyz;

    std::cin >> sx >> sy >> sz >> txy >> tyz >> tzx;

    double sm = (sx + sy + sz) / 3.0;         //平均応力
    double s[3][3] = {                        //偏差応力テンソル
        { sx - sm, txy, tzx},
        { txy, sy - sm, tyz},
        { tzx, tyz, sz - sm},
    };
    
    //第2不変量
    double J2 = -(s[0][0] * s[1][1] + s[1][1] * s[2][2] + s[2][2] * s[0][0])
                +(s[0][1] * s[1][0] + s[2][1] * s[1][2] + s[2][0] * s[0][2]);

    //第3不変量
    double J3 = s[0][0] * s[1][1] * s[2][2] 
              + s[0][1] * s[1][2] * s[2][0] 
              + s[0][2] * s[2][1] * s[1][0] 
              - s[1][2] * s[2][1] * s[0][0] 
              - s[2][0] * s[0][2] * s[1][1] 
              - s[0][1] * s[1][0] * s[2][2];

    double alpha = acos(3.0 * sqrt(3.0) * J3 / 2.0 / sqrt(J2 * J2 * J2)) / 3.0;

    double sn[3] = {
        2.0 / 3.0 * sqrt(3.0 * J2) * cos(alpha) + sm,
        2.0 / 3.0 * sqrt(3.0 * J2) * cos(alpha + 240.0 / 180.0 * M_PI) + sm,
        2.0 / 3.0 * sqrt(3.0 * J2) * cos(alpha + 120.0 / 180.0 * M_PI) + sm
    };

    std::cout << sn[0] << "\t" << sn[1] << "\t" << sn[2] << std::endl;

    double cosx[3];
    double cosy[3];
    double cosz[3];

    for (int i = 0; i < 3; i++)
    {
        
        double A[3][3] = {
            { sx - sn[i], txy, tzx},
            { txy, sy - sn[i], tyz},
            { tzx, tyz, sz - sn[i]},
        };

        double b[3] = {0,0,0};
        double x[3];
        GaussSeidelMethod(A,b,x);
        double abs = sqrt(x[0]*x[0]+x[1]*x[1]+x[2]*x[2]);

        cosx[i] = x[0] / abs;
        cosy[i] = x[1] / abs;
        cosz[i] = x[2] / abs;
    }

    std::cout << cosx[0] << "\t" << cosx[1] << "\t" << cosx[2] << std::endl;
    std::cout << cosy[0] << "\t" << cosy[1] << "\t" << cosy[2] << std::endl;
    std::cout << cosz[0] << "\t" << cosz[1] << "\t" << cosz[2] << std::endl;
    return 0;
}