Mathematica有限元解平面三角形划分习题5.4.pdf

(*许韬 硕士1402*) (*u [x_ ,y_] :=a0+a1*x+a2*y; v [x_ ,y_] :=b0 +b1 *x+b2 *y*) (*结构离散化与编号*) 30 0 nodeMatrix = 30 20 ;(*通过矩阵赋值是为了数据的导入*) 0 20 0 0 n = Length [nodeMatrix [[All, 1 ]]]; node [i_] := nodeMatrix [[i]]; Array[node, n ]; ListPlot [Array [node, n ], PlotMarkers → {Automatic, Large }, ImageSize → Small] 20 ● ● 15 10 5 ● ● 5 10 15 20 25 30 (*节点位移阵列*) nodeDisplacementMatrix = Array [u, {2 n }]; {%} // MatrixForm (*添加位移边界条件*) nodeDisplacementMatrix [[2 * 1 - 0]] = 0; nodeDisplacementMatrix [[2 * 3 - 1]] = 0; nodeDisplacementMatrix [[2 * 3 - 0]] = 0; nodeDisplacementMatrix [[2 * 4 - 1]] = 0; nodeDisplacementMatrix [[2 * 4 - 0]] = 0; (*节点外荷载阵列*) nodeForceMatrix = Array [r, {2 n }]; {%} // MatrixForm (*添加外力边界条件*) nodeForceMatrix [[2 * 2]] = -1000; (*节点约束反力阵列*) nodeForceMatrix [[2 * 1 - 1]] = 0; nodeForceMatrix [[2 * 1 - 0]] = ry1; nodeForceMatrix [[2 * 2 - 1]] = 0; nodeForceMatrix [[2 * 3 - 1]] = rx3; nodeForceMatrix [[2 * 3 - 0]] = ry3; nodeForceMatrix [[2 * 4 - 1]] = rx4; nodeForceMatrix [[2 * 4 - 0]] = ry4; {nodeForceMatrix } // MatrixForm (*添加约束反力*) (u [1] u [2] u [3] u [4] u [5] u [6] u [7] u [8] ) (r [1] r [2] r [3] r [4] r [5] r [6] r [7] r [8] ) (0 ry1 0 - 1000 rx3 ry3 rx4 ry4 ) 2 习题5.4.nb (*各个单元描述*) elasticity = 3 * 105 ; thickness = 5; µ = 0.25; (*均匀板就先把材料参数定义成常量吧*) elementNodeMatrix = 1 2 4 ; 3 4 2 element [i_] := elementNodeMatrix [[i]]; x [i_] := {nodeMatrix [[elementNodeMatrix [[i, 1]], 1]], nodeMatrix [[elementNodeMatrix [[i, 2]], 1]], nodeMatrix [[elementNodeMatrix [[i, 3]], 1]]}; y [i_] := {nodeMatrix [[elementNodeMatrix [[i, 1]], 2]], nodeMatrix [[elementNodeMatrix [[i, 2]], 2]],