1import three;
2import palette;
3
4int N = 26;
5real[] C = array(N,0);
6real[][] A = new real[N][N];
7for(int i = 0; i < N; ++i)
8  for(int j = 0; j < N; ++j)
9    A[i][j] = 0;
10
11real Tb = 100; // deg C
12real h = 240; // 240 W/m^2 K
13real k = 240; // W/m K
14real Tinf = 20; // deg C
15real L = 12; // cm
16real t = 2; // cm
17
18real delta = 0.01; // 1 cm = 0.01 m
19
20// (1,2)-(2,2)-(3,2)-...-(13,2)
21//   |     |     |          |
22// (1,1)-(2,1)-(3,1)-...-(13,1)
23//
24//        |
25//       \ /
26//        V
27//
28// 13-14-15-...-24-25
29//  |  |  | ...  |  |
30//  0- 1- 2-...-11-12
31
32// but, note zero-based array indexing, so counting starts at 0
33int indexof(int m, int n)
34{
35  return 13(n-1)+m-1;
36}
37
38int i = 0;
39
40// fixed temperature bottom left
41A[i][indexof(1,1)] = 1; C[i] = Tb;
42++i;
43// fixed temperature middle left
44A[i][indexof(1,2)] = 1; C[i] = Tb;
45++i;
46
47// interior nodes
48for(int m = 2; m<13; ++m)
49{
50  A[i][indexof(m,2)] = -4;
51  A[i][indexof(m-1,2)] = A[i][indexof(m+1,2)] = 1;
52  A[i][indexof(m,1)] = 2;
53  C[i] = 0;
54  ++i;
55}
56
57// convective bottom side nodes
58for(int m = 2; m<13; ++m)
59{
60  A[i][indexof(m,1)] = -(2+h*delta/k);
61  A[i][indexof(m-1,1)] = A[i][indexof(m+1,1)] = 0.5;
62  A[i][indexof(m,2)] = 1;
63  C[i] = -h*delta*Tinf/k;
64  ++i;
65}
66
67// convective bottom right corner node
68A[i][indexof(13,2)] = A[i][indexof(12,1)] = 0.5;
69A[i][indexof(13,1)] = -(1+h*delta/k);
70C[i] = -h*delta*Tinf/k;
71++i;
72
73// convective middle right side node
74A[i][indexof(13,2)] = -(2+h*delta/k);
75A[i][indexof(13,1)] = 1;
76A[i][indexof(12,2)] = 1;
77C[i] = -h*delta*Tinf/k;
78++i;
79
80real[] T = solve(A,C);
81
82pen[] Palette = Gradient(256,blue,cyan,yellow,orange,red);
83
84real[][] T = {T[0:13],T[13:26],T[0:13]};
85T = transpose(T);
86
87size3(15cm);
88real w = 10;
89real h = 5;
90currentprojection = orthographic(2*(L,h,w),Y);
91draw((L,t,0)--(L,0,0)--(L,0,w)--(0,0,w)--(0,-h,w));
92draw((0,t,w)--(0,t+h,w)--(0,t+h,0)--(0,t,0));
93draw((L,0,w)--(L,t,w)--(0,t,w)--(0,t,0)--(L,t,0)--(L,t,w));
94
95real wo2 = 0.5*w;
96draw((0,0,wo2)--(0,t,wo2)--(L,t,wo2)--(L,0,wo2)--cycle);
97
98// centre points
99surface square = surface(shift(-0.5,-0.5,wo2)*unitsquare3);
100surface bottomsquare = surface(shift(-0.5,-0.5,wo2)*scale(1,0.5,1)*unitsquare3);
101surface topsquare = surface(shift(-0.5,0,wo2)*scale(1,0.5,1)*unitsquare3);
102surface leftsquare = surface(shift(-0.5,-0.5,wo2)*scale(0.5,1,1)*unitsquare3);
103surface rightsquare = surface(shift(0,-0.5,wo2)*scale(0.5,1,1)*unitsquare3);
104surface NEcorner = surface(shift(0,0,wo2)*scale(0.5,0.5,1)*unitsquare3);
105surface SEcorner = surface(shift(0,-0.5,wo2)*scale(0.5,0.5,1)*unitsquare3);
106surface SWcorner = surface(shift(-0.5,-0.5,wo2)*scale(0.5,0.5,1)*unitsquare3);
107surface NWcorner = surface(shift(-0.5,0,wo2)*scale(0.5,0.5,1)*unitsquare3);
108
109material lookupColour(int m,int n)
110{
111  int index = round(Palette.length*(T[m-1][n-1]-60)/(100-60));
112  if(index >= Palette.length) index = Palette.length-1;
113  return emissive(Palette[index]);
114}
115
116draw(shift(0,1,0)*rightsquare,lookupColour(1,2));
117for(int i = 2; i < 13; ++i)
118{
119  draw(shift(i-1,1,0)*square,lookupColour(i,2));
120}
121draw(shift(12,1,0)*leftsquare,lookupColour(13,2));
122
123draw(shift(0,2,0)*SEcorner,lookupColour(1,3));
124draw(shift(0,0,0)*NEcorner,lookupColour(1,1));
125for(int i = 2; i < 13; ++i)
126{
127  draw(shift(i-1,0,0)*topsquare,lookupColour(i,1));
128  draw(shift(i-1,2,0)*bottomsquare,lookupColour(i,3));
129}
130draw(shift(12,2,0)*SWcorner,lookupColour(13,3));
131draw(shift(12,0,0)*NWcorner,lookupColour(13,1));
132
133// annotations
134draw("$x$",(0,-h/2,w)--(L/4,-h/2,w),Y,Arrow3(HookHead2(normal=Z)),BeginBar3(Y));
135draw("$y$",(0,0,1.05*w)--(0,2t,1.05*w),Z,Arrow3(HookHead2(normal=X)),
136     BeginBar3(Z));
137draw("$z$",(L,-h/2,0)--(L,-h/2,w/4),Y,Arrow3(HookHead2(normal=X)),BeginBar3(Y));
138
139draw("$L$",(0,-h/4,w)--(L,-h/4,w),-Y,Arrows3(HookHead2(normal=Z)),
140     Bars3(Y),PenMargins2);
141draw("$w$",(L,-h/4,0)--(L,-h/4,w),-Y,Arrows3(HookHead2(normal=X)),
142     Bars3(Y),PenMargins2);
143draw("$t$",(1.05*L,0,0)--(1.05*L,t,0),-2Z,Arrows3(HookHead2(normal=Z)),
144     Bars3(X),PenMargins2);
145
146label(ZY()*"$T_b$",(0,t+h/2,wo2));
147
148label("$h$,$T_\infty$",(L/2,t+h/2,0),Y);
149path3 air = (L/2,t+h/3,w/3.5)--(1.5*L/2,t+2*h/3,w/8);
150draw(air,EndArrow3(TeXHead2));
151draw(shift(0.5,0,0)*air,EndArrow3(TeXHead2));
152draw(shift(1.0,0,0)*air,EndArrow3(TeXHead2));
153