Archive for the 'FEM/FEA' Category

26
May
06

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:

24
May
06

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

17
Nov
04

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




Welcome !

You will find here some of my favorite hobbies and interests, especially science and art.

I hope you enjoy it!

Subscribe to the RSS feed to stay informed when I publish something new here.

I would love to hear from you! Please feel free to send me an email : bugman123-at-gmail-dot-com

Archives

Blog Stats

  • 543,881 hits