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

Finite Element Method (FEM) Solution to Poisson’s equation on Triangular Mesh

Finite Element Method (FEM) Solution to Poisson’s equation on Triangular Mesh
solved in Mathematica 4.2, 5/24/06; mesh generated with Gmsh; old version: Matlab and DistMesh

Here is some code to solve Poisson’s equation on an unstructured grid of triangular elements using the Finite Element Method (FEM):
(* runtime: 0.1 second *)
n = 10; nodes = Flatten[Table[{x, y}, {y, -1, 1, 2.0/(n - 1)}, {x, -1, 1, 2.0/(n - 1)}], 1]; elements = Flatten[Table[{{i, i + 1, i + n}, {i + 1, i + n + 1, i + n}} + (j - 1) n, {j, 1, n - 1}, {i, 1, n - 1}], 2]; nnode = Length[nodes]; 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[(xy =Join[Transpose[nodes[[#]]], {{1, 1, 1}}];deta = Inverse[xy][[All, {1, 2}]]; Kglobal[[#, #]] +=Det[xy]deta.Transpose[deta]/2; F[[#]] += Det[xy]/6) &, 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]];

Here is some code for an interpolated plot:
(* runtime: 10 seconds *)
n = 275; x1 = y1 = -1; x2 = y2 = 1; image = Table[0, {n}, {n}];
xIntersect[{{x1_, y1_}, {x2_, y2_}}] := If[y1 == y2, Infinity, x1 + (y - y1)(x2 - x1)/(y2 - y1)];
Scan[(plist = nodes[[#]]; xlist = nodes[[#, 1]]; ylist = nodes[[#, 2]]; p1 = plist[[1]]; {dx2, dy2} = plist[[2]] - p1; {dx3, dy3} = plist[[3]] - p1; Do[y = y1 + (y2 - y1)(i - 1)/(n - 1); {xa, xb,xc} = Map[xIntersect, Partition[Sort[plist, #1[[2]] < #2[[2]] &], 2, 1, 1]]; jlist = Round[n({xc, If[Abs[xb - xc] < Abs[xa - xc], xb, xa]} - x1)/(x2 - x1)]; If[jlist =!= {}, Do[x = x1 + (x2 - x1)(j - 1)/(n - 1); {xi, eta} = ((x - p1[[1]]){dy3, -dy2} + (y - p1[[2]]){-dx3, dx2})/(dx2 dy3 - dx3 dy2); image[[i, j]] = phi[[#]].{1 - xi - eta, xi, eta}, {j, Max[1, Min[jlist]], Min[n, Max[jlist]]}]], {i, Max[1, Round[n(Min[ylist] - y1)/(y2 - y1)]], Min[n, Round[n(Max[ylist] - y1)/(y2 - y1)]]}]) &, elements];
ListDensityPlot[image, Mesh -> False, Frame -> False, ColorFunction -> PlotColor, Epilog -> Map[{Line[Map[({x, y} = nodes[[#, {1, 2}]]; n{(x - x1)/(x2 -x1), (y - y1)/(y2 - y1)}) &, #]]} &, elements]];

Crackle

Here is a technique for generating periodic cellular textures using Voronoi diagrams.
(* runtime: 8.5 minutes *)
n = 275; SeedRandom[0]; nodes = Table[Random[], {100}, {2}];
dist[p1_, p2_] := Module[{d = Map[If[# > n/2, n - #, #] &, Abs[p2 - p1]]}, d.d];
DensityPlot[Module[{dlist = Sort[Map[dist[{j, i}, #] &, nodes]]}, dlist[[2]] -dlist[[1]]], {i, 1, n}, {j, 1, n}, PlotPoints -> n, Mesh -> False, Frame -> False];

Mathematica’s ComputationalGeometry package can also be used to generate Voronoi diagrams, Delaunay triangulation, and convex hulls. This is useful for mesh generation:
(* runtime: 5 seconds *)
<< DiscreteMath`ComputationalGeometry`; SeedRandom[0]; nodes = Table[Random[], {100}, {2}];
DiagramPlot[nodes, PlotRange -> {{0, 1}, {0, 1}}];
PlanarGraphPlot[nodes];
PlanarGraphPlot[nodes, ConvexHull[nodes]];

POV-Ray also has a built-in function for this:
// runtime: 0.5 second
camera{orthographic location <0,0,-4> look_at 0 angle 90}
plane{z,0 pigment{crackle color_map{[0 rgb 0][0.5 rgb <0,1,0>][1 rgb 1]}} finish{ambient 1}}

Link: Cellular Texture Tutorial – by Jim Scott