06
Jun
07

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

02
Mar
07

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

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

21
Apr
04

### 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.
25
Oct
00

### Koch’s Snowflake

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

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