# -*- coding: utf-8 -*- # 3.0 # # Discrete Heat Transfer # from sympy import * init_printing() # # Wire-Mesh Model # # We consider a 2D mesh of nodes connected by thin wires. # In equilibrium, the temperature of a node is equal to the mean of its neighbours. # At boundaries, we assume the temperature to be 0. # \$N\$ and \$M\$ denote the number of inner nodes in \$x\$- and \$y\$-direction: # N = 4; M = 4 # # x and y will hold the coordinates of the mesh nodes: # x = np.linspace(0,1,N+1) y = np.linspace(0,1,M+1) # # And the \$(N + 1) (M + 1)\$ unknowns \$u_{n, m}\$ are saved here: # u = Matrix(N+1,M+1,lambda i,j: Symbol('u_%d%d'%(i,j))) pprint(u) # # The function \$f\$ is used to compute the right-hand sides (heat sources/drains): # f = lambda x,y: 0 # # We set up the respective system of linear equations: # inner = [Eq(-u[n-1,m]-u[n,m-1]+4*u[n,m]-u[n+1,m]-u[n,m+1],f(x[n],y[n])) for n in range(1,N) for m in range(1,M)] btbnd = [Eq(u[n,0],sin(n*pi/N)) for n in range(0,N+1)] tpbnd = [Eq(u[n,M],0) for n in range(0,N+1)] lfbnd = [Eq(u[0,m],0) for m in range(1,M)] rgbnd = [Eq(u[N,m],0) for m in range(1,M)] print "equations for inner domain: " for eInner in inner: pprint(eInner) print "bottom boundary: " pprint(btbnd) print "top boundary: " pprint(tpbnd) print "left boundary: " pprint(lfbnd) print "right boundary: " pprint(rgbnd) gls = inner + btbnd + tpbnd + lfbnd + rgbnd # # ... and use SymPy to solve the system: # from itertools import chain # convert matrix with unknowns to nested list and then "flatten" it through chain unbek = list(chain.from_iterable(u.tolist())) lsg = solve(gls,unbek) # substitute the solution dictionary to matrix with unknowns u=u.subs(lsg) print "solution: " pprint(u) # # Plotting the solution # import numpy as np from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib import cm # evaluate the solution Q = np.array(u.evalf()).astype(float) X, Y = np.meshgrid(x, y) fig = plt.figure() ax = fig.gca(projection='3d') ax.scatter(X,Y,Q,) fig = plt.figure() ax = fig.gca(projection='3d') ax.plot_wireframe(X,Y,Q, rstride=1, cstride=1, antialiased=True) fig = plt.figure() ax = fig.gca(projection='3d') ax.plot_surface(X,Y,Q, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True) plt.show() # # Finite Volume Model # # Set up the unknowns # N = 4; M = 4 x = np.linspace(0,1,N+1) y = np.linspace(0,1,M+1) u = Matrix(N+1,M+1,lambda i,j: Symbol('u_%d%d'%(i,j))) # # Right-hand side: # f = lambda x,y: 0 # # System of linear equations: # part1 = [Eq(-u[n-1,m]-u[n,m-1]+4*u[n,m]-u[n+1,m]-u[n,m+1],f(x[n],y[n])) for n in range(1,N) for m in range(1,M)]+\ [Eq(u[n,0],sin(n*pi/N)) for n in range(0,N+1)] neumpart = [Eq(-u[n-1,M]-u[n,M-1]+3*u[n,M]-u[n+1,M],0) for n in range(1,N)]+\ [Eq(-u[0,m-1]+3*u[0,m]-u[1,m]-u[0,m+1],0) for m in range(1,M)]+\ [Eq(-u[N,m-1]+3*u[N,m]-u[N-1,m]-u[N,m+1],0) for m in range(1,M)]+\ [Eq(-u[N,M-1]+2*u[N,M]-u[N-1,M],0)]+\ [Eq(-u[0,M-1]+2*u[0,M]-u[1,m],0)] gls = part1+neumpart # # and solve the system # unbek = list(chain.from_iterable(u.tolist())) lsg = solve(gls,unbek) u=u.subs(lsg) print "solution: " print(u.evalf()) # Q = np.array(u.evalf()).astype(float) X, Y = np.meshgrid(x, y) fig = plt.figure() ax = fig.gca(projection='3d') ax.scatter(X,Y,Q,) fig = plt.figure() ax = fig.gca(projection='3d') ax.plot_wireframe(X,Y,Q, rstride=1, cstride=1, antialiased=True) fig = plt.figure() ax = fig.gca(projection='3d') ax.plot_surface(X,Y,Q, rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=True) plt.show()