構造力学・強度計算実践への一歩 板計算の説明(有限要素法) | 2013-04-28 Mori design office |
および特異例として梁理論値、円板理論値と比較検証します(C項)。 最後に、計算式を誘導する為の理論と手順を述べます(D項)。 → 板計算の画面(FEM) |
A. 計算式 |
1つの要素(矩形)において、 横長2a、縦長2bとして 次の剛性方程式が成立し、頂点(節点)の 変位・角度(δ,θX,θY)に対し、頂点に作用する外力が算出できる式となっている。 剛性方程式 [K]{ω} = {F}12x12 12x1 12x1 ---(A1) ここに、 |
![]() |
変位ベクトル を{ω}とし、要素の頂点 1,2,3,4 における、z方向変位δZと x,yに対する角度θX,θYのベクトルとなっている。
{ω} = [δZ1 θX1 θY1 δZ2 θX2 θY2 δZ3 θX3 θY3 δZ4 θX4 θY4 ]T ---(A2) 外力ベクトル を{F}とし、要素の頂点 1,2,3,4 に作用する、z方向力fとモーメントmのベクトルとなっている。 {F} = [ f1 mX1 mY1 f2 mX2 mY2 f3 mX3 mY3 f4 mX4 mY4 ]T ---(A3) なお、 モーメントの方向は、例えば mX は 3次元でのmY であり、Y軸中心回転(右手法則)とする。 モーメントMは、板理論で用いられている、単位長さ当たりのモーメントとする。 mX1 = bMX1 mX2 = -bMX2 mX3 = -bMX3 mX4 = bMX4 mY1 = aMY1 mY2 = aMY2 mY3 = -aMY3 mY4 = -aMY4 ---(A4) |
K01,01 = D(10a4+10b4+(7-2ν)a2b2)/(10a3b3) K01,02 = D((1+4ν)a2+10b2)/(10a2b) K01,03 = D(10a2+(1+4ν)b2)/(10ab2) | [kgf/m] [kgf] [kgf] |
K01,04 = D(5a4-10b4-(7-2ν)a2b2)/(10a3b3) K01,05 = D((1-ν)a2+10b2)/(10a2b) K01,06 = D(5a2-(1+4ν)b2)/(10ab2) | [kgf/m] [kgf] [kgf] |
K01,07 = D(-5a4-5b4+(7-2ν)a2b2)/(10a3b3) K01,08 = -D((1-ν)a2-5b2)/(10a2b) K01,9 = D(5a2-(1-ν)b2)/(10ab2) | [kgf/m] [kgf] [kgf] |
K01,10 = D(-10a4+5b4-(7-2ν)a2b2)/(10a3b3) K01,11 = D(-(1+4*ν)a2+5b2)/(10a2b) K01,12 = D(10a2+(1-ν)b2)/(10ab2) | [kgf/m] [kgf] [kgf] |
K02,01 = K01,02 K02,02 = D(4(1-ν)a2+20b2)/(15ab) K02,03 = Dν | [kgf] [kgf・m] [kgf・m] |
K02,04 = -D((1-ν)a2+10b2)/(10a2b) K02,05 = D(-(1-ν)a2+10b2)/(15ab) K02,06 = 0 | [kgf] [kgf・m] [kgf・m] |
K02,07 = -D(-(1-ν)a2+5b2)/(10a2b) K02,08 = D((1-ν)a2+5b2)/(15ab) K02,09 = 0 | [kgf] [kgf・m] [kgf・m] |
K02,10 = D(-(1+4*ν)a2+5b2)/(10a2b) K02,11 = D(-4*(1-ν)a2+10b2)/(15ab) K02,12 = 0 | [kgf] [kgf・m] [kgf・m] |
K03,01 = K01,03 K03,02 = K02,03 K03,03 = D(20a2+4*(1-ν)b2)/(15ab) | [kgf] [kgf・m] [kgf・m] |
K03,04 = D(5a2-(1+ 4ν)b2)/(10ab2) K03,05 = 0 K03,06 = D(10a2-(4-4ν)b2)/(15ab) | [kgf] [kgf・m] [kgf・m] |
K03,07 = D(-5a2+(1-ν)b2)/(10ab2) K03,08 = 0 K03,09 = D(5a2+(1-ν)b2)/(15ab) | [kgf] [kgf・m] [kgf・m] |
K03,10 = -D(10a2+(1-ν)b2)/(10ab2) K03,11 = 0 K03,12 = D(10a2-(1-ν)b2)/(15ab) | [kgf] [kgf・m] [kgf・m] |
・・・ ・・・中略 ・・・ K10,10 = D(10a4+10b4+(7-2ν)a2b2)/(10a3b3) K10,11 = D((1+4ν)a2+10b2)/(10a2b) K10,12 = D(-10a2-(1+4ν)b2)/(10ab2) | [kgf/m] [kgf] [kgf] |
K11,10 = K10,11 K11,11 = D(8(1-ν)a2+40b2)/(30ab) K11,12 = -Dν | [kgf] [kgf・m] [kgf・m] |
K12,10 = K10,12 K12,11 = K11,12 K12,12 = D(20a2+4(1-ν)b2)/(15ab) なお、 板の曲げ剛性 D = Eh3/(12(1-ν2)) | [kgf] [kgf・m] [kgf・m] [kgf・m] |
B. 計算方法 (プログラミングへの要約) |
既知外力より変位を算出する手順を述べる。ここでは、外部モーメントは板の外周のみ作用するとする。 @12元連立方程式 上記式(A1)を全要素に適用させる。連立式は、要素数の12倍の数となる。 A組立合成(全体剛性行列) 共通節点のは加算合成して、節点総数に対する式を得る。 例えば、板の区割り数を横ia、縦ib の場合 要素総数Nelt=ia・ib、節点総数Nnod =(ia+1)・(ib+1) 全要素 に式(A1)を適応させると、連立式数Nelt x 3 であり、 共通節点のは加算合成で、連立式数 Nnodx3となる。 ・外周のモーメント 板の辺a全にMy、辺b全にMxが作用する場合、 要素毎に計算、式(A4)を計算する。 節点に作用するモーメントが、2個の要素にまたがるばあいは、各々要素に対し1/2づづにする 。 この値を外力{F}に代入する。 B変位の既知処理 変位に関し、既に既知である場合、連立式から除外する。 ・例:境界条件は支持方法であり、下記の既知節点処理となる。 ・全周単純支持では、変位z=0の既知であり、既知節点数 No = 2(ia+ib)だけを除外して 連立式数は Nnodx3 - 2(ia+ib) ・全周固定支持では、θx=0,θy=0の既知も加わる故、 No = 2(ia+ib)・3 個の式を除外して 連立式数 Nnodx3 - 2(ia+ib)x3 のを解くことになる。 ・強制変位δzの場合 変位δzに相当の力 {Fz} = [K]{δ} 算出する。 外力Fと相殺して、{Fo} = {F} - {Fz} を算出する。 全剛性式 [K]{ω} = {Fo} であり、上述と同様に既知節点を除外して、連立式を立てる。 ・この既知節点除外のプログラミングは、 例えば、既知節点が1,2,3,・・のKOZ個の場合、既知節点配列をknow[0,1,2,3,・・]とすると、 for(i=1;i<=Nnod3;i++){index[i]=1;} のパラメータに for(i=1;i<=KOZ;i++){index[know[i]]=0;} m = 0; for(i=1;i<= Nnod3;i++){ if (index[i]>0){ m = m + 1; index[m] = i;} } for(i=1;i<= m;i++){ F[i] = F[index[i]]; for(j=1;j<= m;j++){K[(i-1)*index[i]+ j] = K[(index[i]-1)*Nnod3+ index[j]];} } ・これにて、求めるべき式 [K]{ω} = {F}の K,Fが得られたことになる。 C連立式の解法 ・得られた [K]は正方行列mxmで、逆行列 [K]-1 で変位ωが求められる。 ・プログラミングとしては、ガウスジョルダン法やコレスキー法などがある。 ・この方法は、処理結果が外力F[i]の場所に変位の値が出る。即ち、F[i]=変位の値 に置換される。 D変位{ω}の算出 ・Cで得られた値に対し、Bの既知処理の逆処理を行う for(i=1;i<=Nnod3;i++){ω[i] = 0;} for(i=1;i<=mm;i++){ω[index[i]] = F[i];} for(i=1;i<=Nnod3;i++){ ω[i] = ω[i] + 強制変位δ[i]; } E内力の算出 ・変位の値が算出できたので、各要素ごとに剛性方程式(A1)に代入して内力F及びモーメントMが算出できる。 ・算出した内力は、各々節点に関し単純加算して、節点の内力値を算出する。 ・モーメントは、式(A4)を参照して正負符号も含めて、単純平均でMx,Myを算出 (例えば同一節点が2個あれば、その合計値を2で割り単純平均)する。 ・応力 σx= 6Mx/h2 σy= 6My/h2 |
C.計算結果 (精度検証) |
一般荷重での収束性 下記(a)〜(d)は、収束性が良く、要素数1600で便覧値に近い。
ポアソン比ν
☆はり理論:長さa=30cm、幅b=1cm、断面2次モーメント I = b・T3/12 = 2.25cm4 断面係数 Z = b・T2/6 = 1.5cm3 θ=F・a2/(2EI)、δ=F・a2/(3EI)、σ=-F・a/Z モーメント
・このモーメントの場合は、幅(1cm)の分割を多く(100)の方が精度が高い。 参考
分割数を大きくするほど、円周部のモーメント値に乱れが大きくなる。 (g)典型的条件に条件追加によるプログラムチェック(数値およびグラフ割愛) 全周単純支持の条件に ・左右辺において、θy=0 及び Mx=全周固定での値 ・上下辺において、θx=0 及び My=全周固定での値 結果値は、全周固定の条件と同一となる。 荷重零の条件に ・変位δ=集中荷重が有りでの変位 (集中荷重の点を 変位拘束点として) 結果値は、集中荷重が既知で 荷重下位置が変位未知の場合 と同一となる。 などが考えられる。 | ![]() 使用した物理定数値 ポアソン比 ν= 0.3 ヤング率 E= 2100000 [kgf/cm2] |
D.理論 (式の誘導) |
近似多項式 四角形要素において、定数項α1〜α12を用いて、任意のたわみを近似多項式
歪は、上記の内挿関数を用いて次のようになる。
---(10)
---(12)
12(1-ν )2 要素剛性行列の誘導計算
・板の関連資料:板計算説明(級数) FEMの関連:FEM要約(四角形)、 FEM入門(三角形) Copyright(C) Mori Design office All rights Reserved. |