FEM Solution to Poisson’s equation on Quad Mesh

FEM Solution to Poisson’s equation on Quad Mesh – Mathematica 4.2, AutoCAD 2000, 6/10/07
Here is some code to solve Poisson’s equation on an unstructured grid of quadrilateral elements:
(* runtime: 0.05 second *)
n = 10; nodes = Flatten[Table[{x, y}, {x, -1, 1,2.0/(n - 1)}, {y, -1, 1, 2.0/(n - 1)}], 1]; nnode = Length[nodes]; elements =Flatten[Table[{i, i + 1, i + n + 1, i + n} + (j - 1) n, {j, 1, n - 1}, {i, 1, n - 1}], 1]; nelem = Length[elements];
fixed = Map[Position[nodes, #][[1, 1]] &, Select[nodes, Max[Abs[#]] == 1 &]];
Kglobal = Table[0, {nnode}, {nnode}]; F = phi = Table[0, {nnode}];
Scan[(plist = nodes[[#]]; dphi = {plist[[2]] - plist[[1]],plist[[4]] - plist[[1]]}; B = Inverse[Transpose[dphi].dphi]; C1 = {{2, -2}, {-2, 2}}B[[1, 1]] + {{3, 0}, {0, -3}}B[[1, 2]] + {{2, 1}, {1,2}}B[[2, 2]];C2 = {{-1, 1}, {1, -1}}B[[1, 1]] + {{-3, 0}, {0, 3}} B[[1, 2]] + {{-1, -2}, {-2, -1}}B[[2, 2]]; Kglobal[[#, #]] += Det[dphi]Join[Thread[Join[C1, C2]], Thread[Join[C2, C1]]]/6; F[[#]] += Det[Join[{{1, 1, 1}}, Transpose[nodes[[#[[{1, 2, 3}]], {1, 2}]]]]]/4) &, elements];
free = Complement[Range[nnode], fixed]; phi[[free]] = LinearSolve[Kglobal[[free, free]], F[[free]]];
Mean[x_] := Plus @@ x/Length[x]; PlotColor[x_] := Hue[2(1 - x)/3];
Show[Graphics[Table[{PlotColor[Mean[phi[[elements[[i]]]]]/Max[phi]], Polygon[nodes[[elements[[i]]]]]}, {i, 1, nelem}], AspectRatio -> 1]];

Links

• FEM_50 – simple Matlab Poisson solver for quads and triangles
• FEM Poisson – paper by Nasser Abbasi
• Poisson solver – simple finite volume C code by Jeroen Valk

Advertisements

CFM56-5 Turbofan Jet Engine

I saw this cool looking poster of a turbofan engine schematic and it inspired me to make a 3D rendering of it. This engine is used for the Airbus. The streamlines shown here are just for visual effect. In reality, the streamlines would be turbulent. The blue streamlines show where cold air enters and bypasses the engine (turbofans typically have a high bypass ratio). Green streamlines are shown through the compressor stages. The red streamlines show where hot gas exits the burners and turbine stage.

CFM56-5 Turbofan Jet Engine – drawn in AutoCAD 2000, AutoLisp, rendered in POV-Ray 3.6.1, 3/2/07

Finite Element Analysis (FEA)

Here we see the deformation and stress distribution between two gears where the left gear is fixed but the right gear has an applied torque. This was solved using Finite Element Analysis (FEA) to solve the linear elasticity equations for an unstructured grid of quad elements. The color represents the equivalent (Von Mises) stress.

Finite Element Analysis (FEA) – AutoCAD 2000, AutoLisp, Mathematica 4.2, 5/26/06

(* runtime: 0.8 second *)
n = 18; dtheta = 2Pi/n; nodes = Flatten[Table[r{Sin[theta], Cos[theta]}, {r, 0.4, 1, 0.2}, {theta, 0, 2Pi - dtheta, dtheta}], 1];
elements = Flatten[Table[n j + Append[Table[i + {0,1, n + 1, n}, {i, 1, n - 1}], {1, n + 1, 2n, n}], {j, 0, 2}], 1]; nnode = Length[nodes]; nelem = Length[elements];
i = 3n/2 + 1; fixed = {{i - 1, "x"}, {i, "x"}, {i, "y"}, {i + 1, "x"}}; loads = {{n + 1, {0, -200}}};
Dkl := 0.25{{-1 + y, 1 - y, 1 + y, -1 - y}, {-1 + x, -1 - x, 1 + x, 1 - x}}; Dkm = 0.25{{-1, 1, 1, -1}, {-1, -1, 1, 1}};
e = 1500; nu = 0.3; Epss = {{1, nu, 0}, {nu, 1, 0}, {0, 0, (1 - nu)/2}}e/(1 - nu^2); thk = 1;
Kglobal = Table[0, {2nnode}, {2nnode}];
Scan[(element = #; enodes = nodes[[element]]; Klocal = Table[0, {8}, {8}]; Kee = Table[0, {4}, {4}]; Kre = Table[0, {8}, {4}]; Scan[({x, y} = #/Sqrt[3]; t = Inverse[Dkl.enodes].Dkl; strain = {{t[[1, 1]], 0, t[[1, 2]], 0, t[[1, 3]], 0, t[[1, 4]], 0}, {0, t[[2, 1]], 0, t[[2, 2]], 0, t[[2, 3]], 0, t[[2, 4]]}, {t[[2, 1]], t[[1, 1]], t[[2, 2]], t[[1, 2]], t[[2, 3]], t[[1, 3]], t[[2, 4]], t[[1, 4]]}}; ttm = -2Inverse[Dkm.enodes].{{x, 0}, {0, y}}; ct = {{ttm[[1,1]], 0, ttm[[1, 2]], 0}, {0, ttm[[2, 1]], 0, ttm[[2, 2]]}, {ttm[[2, 1]], ttm[[1, 1]], ttm[[2, 2]], ttm[[1, 2]]}}Det[Dkm.enodes]/Det[Dkl.enodes]; Dj = thk Abs[Det[Dkl.enodes]]; Klocal += Transpose[strain].(Epss.strain) Dj; Kee += Transpose[ct].(Epss.ct) Dj; Kre += Transpose[strain].(Epss.ct) Dj) &, {{-1, -1}, {1, -1}, {1, 1}, {-1, 1}}];Klocal -= Kre.Inverse[Kee].Transpose[Kre];Do[Kglobal[[2element[[i]] - di, 2element[[j]] - dj]] += Klocal[[2i -di, 2j - dj]], {di, 0, 1}, {dj, 0, 1}, {i, 1, 4}, {j, 1,4}]) &, elements];
Scan[(i = 2#[[1]] - If[#[[2]] == "x", 1, 0]; Kglobal[[i, i]] *= 10^6) &, fixed];
paload = Table[0, {2nnode}]; Scan[(paload[[2#[[1]] - {1, 0}]] = #[[2]]) &, loads];
del = Partition[LinearSolve[Kglobal, paload], 2];
Do[nodes2 = nodes + scale del; Show[Graphics[Map[Polygon[nodes2[[#]]] &, elements], AspectRatio -> Automatic]], {scale, 0, 1, 0.1}];

Links:

• 
  

  Gears – nice FEA animation by Noran Engineering, Inc. (where I currently work)

• 
  
• Mathematica FEA – code for 2D beam using quad elements, by Chee Kiang Yew
• Mini-FEA – simple Java applet
• IMTEK Mathematica Supplement – free FEA and meshing packages
• Rubber Membrane – interesting animation of buckling behavior
• Physica – has the ability to do “MultiPhysics” (for example, fluid-structure interaction)
• Real-Time FEA – cool Java applet for 3D beam by Matthias Müller-Fischer
• AceFEM – FEA environment for Mathematica

Fundamental Modes of Vibration for a Violin-Shaped Membrane

Mode displacement plots of the 19th, 28th, and 45th modes, respectively, from left to right, computed with GeoStar 2.8.Click here to see a holographic interferogram of a real violin.

Koch’s Snowflake

Koch’s Snowflake  is another type of L-System.

Links

• The self-similar rippling of leaf edges and torn plastic sheets – hyperbolic surfaces resembling Koch’s snowflake
• L-Systems – Mathematica code by unknown author, hosted by Jan Poland