qlearner的个人博客分享 http://blog.sciencenet.cn/u/qlearner

博文

利用Mathematica进行有限元编程(一):杆单元分析

已有 4958 次阅读 2016-8-17 16:19 |个人分类:编程学习|系统分类:科研笔记| 有限元, 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的位移分别为uiui<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;">ujuj,受到的力分别为<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;">FiFi<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;">FjFj。易知:

<span class="MathJax" id="MathJax-Element-5-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML" display="block"><mtable columnalign="right left" rowspacing="3pt" columnspacing="0em" displaystyle="true">Fi=kui&#x2212;kujFj=kuj&#x2212;kui" role="presentation" style="display:inline;font-size:18px;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;">FiFj=kuikuj=kujkuiFi=kui−kujFj=kuj−kui

写成矩阵形式为:
<span class="MathJax" id="MathJax-Element-6-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">[<mtable rowspacing="4pt" columnspacing="1em">k&#x2212;k&#x2212;kk][<mtable rowspacing="4pt" columnspacing="1em">uiuj]=[<mtable rowspacing="4pt" columnspacing="1em">FiFj]" role="presentation" style="display:inline;font-size:18px;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;">[kkkk][uiuj]=[FiFj][k−k−kk][uiuj]=[FiFj]

所以单元刚度矩阵为:
<span class="MathJax" id="MathJax-Element-7-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">k<mrow class="MJX-TeXAtom-ORD">ij=[<mtable rowspacing="4pt" columnspacing="1em">k&#x2212;k&#x2212;kk]" role="presentation" style="display:inline;font-size:18px;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;">kij=[kkkk]kij=[k−k−kk]

同理,一次杆单元的刚度矩阵为:
<span class="MathJax" id="MathJax-Element-8-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">k<mrow class="MJX-TeXAtom-ORD">ij=EAL[<mtable rowspacing="4pt" columnspacing="1em">1&#x2212;1&#x2212;11]" role="presentation" style="display:inline;font-size:18px;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;">kij=EAL[1111]kij=EAL[1−1−11]

二次杆单元的刚度矩阵可通过能量法求解。应用能量法的示意图:

应变能为:
<span class="MathJax" id="MathJax-Element-9-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">U=&#x222B;<mi mathvariant="normal">&#x03A9;12&#x03C3;&#x03F5;dV=&#x222B;0L12EA&#x03F5;2dx=&#x222B;0L12EAaTdNTdxdNdxadx=12aTKa" role="presentation" style="display:inline;font-size:18px;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;">U=Ω12σϵdV=L012EAϵ2dx=L012EAaTdNTdxdNdxadx=12aTKaU=∫Ω12σϵdV=∫0L12EAϵ2dx=∫0L12EAaTdNTdxdNdxadx=12aTKa

其中:
<span class="MathJax" id="MathJax-Element-10-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">&#x03C3;=E&#x03F5;;&#x03F5;=dudx=dNdxa,a=<mo stretchy="false">[u1,u2,u3<mo stretchy="false">]T" role="presentation" style="display:inline;font-size:18px;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;">σ=Eϵ;ϵ=dudx=dNdxa,a=[u1,u2,u3]Tσ=Eϵ;ϵ=dudx=dNdxa,a=[u1,u2,u3]T

所以单元刚度矩阵为:
<span class="MathJax" id="MathJax-Element-11-Frame" tabindex="0" data-mathml="<math xmlns="http://www.w3.org/1998/Math/MathML" display="block">K=EA&#x222B;<mrow class="MJX-TeXAtom-ORD">0L<mo stretchy="false">[dNdx<mo stretchy="false">]T<mo stretchy="false">[dNdx<mo stretchy="false">]dx" role="presentation" style="display:inline;font-size:18px;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;">K=EAL0[dNdx]T[dNdx]dxK=EA∫0L[dNdx]T[dNdx]dx

代入上面的形函数,那么求解代码为:


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:建立单元刚度矩阵

以下分别建立了弹簧元、线性杆单元、二次杆单元的单元刚度矩阵:

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行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是单刚。




http://blog.sciencenet.cn/blog-441611-997011.html

上一篇:Mathematica版GetData——用Mathematica提取图片中的数据点
下一篇:求解偏微分方程开源有限元软件deal.II学习--Step 3

1 yangb919

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-15 04:38

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部