博文
利用Mathematica进行有限元编程(一)：杆单元分析
|||
首发：qixinbo.info
本文是对Mathematica有限元分析与工程应用一书的学习笔记。
这里的杆单元是总体坐标与局部坐标一致的一维有限元，因此不需要坐标变换，更易分析。
线性杆单元的形函数线性杆单元的位移函数是坐标x的一次函数，可用以下代码求解：
| 1 2 3 4 | Clear[u1, u2]; u = a1 + a2 x; solut = Solve[{u1 == (u /. x -> 0), u2 == (u /. x -> 1)}, {a1, a2}]; u = u /. solut[[1]] |
其中u1和u2是杆端点x=0和x=1处的位移，结果为：
| 1 | u1 + (-u1 + u2) x |
形函数则可以通过取系数函数Coefficient获得：
| 1 | NShape = {Coefficient[u, u1], Coefficient[u, u2]} |
结果为：
| 1 | {1 - x, x} |
二次杆单元的位移函数是坐标x的二次函数，求解代码为：
| 1 2 3 4 5 | Clear[u1, u2, u3, a1, a2, a3]; u = a1 + a2 x + a3 x^2; solut = Solve[{u1 == (u /. x -> 0), u2 == (u /. x -> 1/2), u3 == {u /. x -> 1}}, {a1, a2, a3}]; u = u /. solut[[1]] |
结果为：
| 1 | u1 + (-3 u1 + 4 u2 - u3) x + 2 (u1 - 2 u2 + u3) x^2 |
其中u1、u2和u3分别是杆端点x=0、杆中点x=0.5、杆端点x=1处的位移。
那么形函数就是：
| 1 | NShape = Coefficient[u, {u1, u2, u3}] |
结果为：
| 1 | {1 - 3 x + 2 x^2, 4 x - 4 x^2, -x + 2 x^2} |
一次杆单元的刚度矩阵可以类比弹簧元的刚度矩阵。
设弹簧的刚度系数为k，节点i，j的位移分别为ui和<span class="MathJax" id="MathJax-Element-2-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">uj" role="presentation" style="display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;padding:0px;margin:0px;position:relative;">uj，受到的力分别为<span class="MathJax" id="MathJax-Element-3-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">Fi" role="presentation" style="display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;padding:0px;margin:0px;position:relative;">Fi，<span class="MathJax" id="MathJax-Element-4-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML">Fj" role="presentation" style="display:inline;line-height:normal;text-align:left;word-spacing:normal;word-wrap:normal;float:none;direction:ltr;max-width:none;max-height:none;min-width:0px;min-height:0px;border:0px;padding:0px;margin:0px;position:relative;">Fj。易知：
写成矩阵形式为：
所以单元刚度矩阵为：
同理，一次杆单元的刚度矩阵为：
二次杆单元的刚度矩阵可通过能量法求解。应用能量法的示意图：
应变能为：
其中：
所以单元刚度矩阵为：
代入上面的形函数，那么求解代码为：
| 1 2 3 | dNx = L D[NShape, x]; K = EA Integrate[Transpose[{dNx}].{dNx}, {x, 0, 1}]/L; K // MatrixForm |
注意：代码中的L这个系数是因为使用了之前位于0到1上的形函数，如果直接采用0到L，则形函数的具体形式会发生改变。
输出结果为：
以下分别建立了弹簧元、线性杆单元、二次杆单元的单元刚度矩阵：
| 1 2 3 4 5 6 | GenerateSpringKm[k_] := Module[{y}, y = k*{{1, -1}, {-1, 1}}; y] GenerateLinearRodKm[k_] := Module[{y}, y = EE*AA/L*{{1, -1}, {-1, 1}}; y] GenerateQuadraticRodKm[k_] := Module[{y}, y = EE*AA/L*{{7, 1, -8}, {1, 7, -8}, {-8, -8, 16}}; y] |
以弹簧元为例，单独的一个弹簧元有两个节点，其单元刚度矩阵是2行2列的矩阵。对于含有n个节点的弹簧元系统，则是n行n列的整体刚度矩阵。组装过程为：假设某个弹簧元的节点的整体编号为i和j，其单刚中的四个元素要分别找到其对应的整体编号位置，然后叠加上去，即11对应ii，12对应ij，21对应ji，22对应jj。代码为：
| 1 2 3 4 5 6 7 | GlobalK = 0 IdentityMatrix[3]; AssembleSpringKm[p1_, p2_, m_] := Module[{j, k}, f = {p1, p2}; For[j = 1, j <= 2, j++, For[k = 1, k <= 2, k++, GlobalK[[f[[j]], f[[k]]]] += m[[j, k]]; ]]; GlobalK] |
代码中p1和p2是节点的整体编号(输入时要注意节点号的顺序)，m是该单元的单刚，GlobalK是整刚。
二次杆单元的整体刚度矩阵则为：
| 1 2 3 4 5 6 7 8 | AssembleQuadraticRod[p1_, p2_, p3_, m_] := Module[{i, j, k}, f = {p1, p2}; For[i = 1, i <= 3, i++, For[j = 1, j <= 3, j++, For[k = 1, k <= 3, k++, GlobalK[[f[[j]], f[[k]]]] += m[[j, k]]; ]]]; GlobalK] |
建立好整体刚度矩阵后，代入边界条件，即可求出位移量，比如：
| 1 2 3 4 5 6 7 | U = {u1x, u2x, u3x}; u1x = u3x = 0; F2x = 17.5; Load = {F1x, F2x, F3x}; tmp = Sort[Join[U, Load]]; unknown = Drop[tmp, 3]; solut = Solve[GlobalK.U == Load, unknown] |
这里使用的边界条件是编号1和3上的位移为0,2上的力为17.5，然后利用Solve求解。
模块3:得到单元应力已知单元的位移向量u(在整体坐标系下的值)，求单元上的力：
| 1 2 3 4 5 6 7 8 9 10 | SpringElementForce[k_, u_, i_] := Module[{force}, force = k.u; Print["The spring force of the ele", i, " is ", force, "kN"]; ]; RodElementForce[k_, u_, AA_, i_] := Module[{force}, force = k.u; sigma = force/AA; Print["The force of the ele", i, " is ", force, "kN, ", "the stress is ", sigma, " Pa."]; ] |
i是单元编号，u是位移向量，k是单刚。
https://blog.sciencenet.cn/blog-441611-997011.html
上一篇：Mathematica版GetData——用Mathematica提取图片中的数据点
下一篇：求解偏微分方程开源有限元软件deal.II学习--Step 3
