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