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

 1234 Clear[u1, u2];u = a1 + a2 x;solut = Solve[{u1 == (u /. x -> 0), u2 == (u /. x -> 1)}, {a1, a2}];u = u /. solut[[1]]

 1 u1 + (-u1 + u2) x

 1 NShape = {Coefficient[u, u1], Coefficient[u, u2]}

 1 {1 - x, x}

 12345 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

 1 NShape = Coefficient[u, {u1, u2, u3}]

 1 {1 - 3 x + 2 x^2, 4 x - 4 x^2, -x + 2 x^2}

<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-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

 123 dNx = L D[NShape, x];K = EA Integrate[Transpose[{dNx}].{dNx}, {x, 0, 1}]/L;K // MatrixForm

 123456 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]

 1234567 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]

 12345678 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]

 1234567 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]

 12345678910 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

亓欣波

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