Python Exercise 4

Fredrik Eikeland Fossan

Feb 26, 2019

Elliptic and parabolic partial differential equations

Problem 1: Stationary temperature field in beam cross-section

In this problem we shall solve solve problem 1 of last weeks theory exercise, using python.

Figure 1: Beam cross-section.

a) Create a python script that solves for the stationary temperature field in the beam cross-section. Choose $ \Delta x = \Delta y = \frac{1}{N} $, and plot the solution for different values of $ N $. See problem 1 from the theory exercise 4 for details on discretization and treatment of boundary conditions.

Hint.

Figure 2: Sketch of discretization with $ N=3 $

Figure 2 shows the discretized grid for $ N=3 $. The nodes with unknown temperature are colored black ($ T_1, \dots T_9 $), whereas ghost nodes are colored grey and nodes where the temperature is known are colored blue.

If the components of $ \boldsymbol{T} $ is numbered by starting in the lower left part of the grid and the index is increased in the $ x $-direction first, the function plotSurfaceNeumannDirichlet could be used to plot the solution.

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
def plotSurfaceNeumannDirichlet(Temp, Ttop, Tleft, l, N, nxTicks=4, nyTicks=4):
   """ Surface plot of the stationary temperature in quadratic beam cross-section.
       Note that the components of T has to be started in the
       lower left part of the grid with increasing indexes in the
       x-direction first.


        Args:
            Temp(array):  the unknown temperatures, i.e. [T_1 .... T_(NxN)]
            Ttop(float):  temperature at the top boundary
            Tleft(float): temperature at the left boundary
            l(float):     height/width of the sides
            N(int):       number of nodes with unknown temperature in x/y direction
            nxTicks(int): number of ticks on x-label (default=4)
            nyTicks(int): number of ticks on y-label (default=4)
    """
    x = np.linspace(0, l, N + 1)
    y = np.linspace(0, l, N + 1)
    
    X, Y = np.meshgrid(x, y)
    
    T = np.zeros_like(X)
    
    T[-1,:] = Ttop
    T[:,0] = Tleft
    k = 1
    for j in range(N):
        T[j,1:] = Temp[N*(k-1):N*k]
        k+=1

    fig = plt.figure()
    ax = fig.gca(projection='3d')
    surf = ax.plot_surface(X, Y, T, rstride=1, cstride=1, cmap=cm.coolwarm,
                           linewidth=0, antialiased=False)
    ax.set_zlim(0, Ttop)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('T [$^o$C]')
    
    xticks=np.linspace(0.0, l, nxTicks+1)
    ax.set_xticks(xticks)
    
    yticks=np.linspace(0.0, l, nyTicks+1)
    ax.set_yticks(yticks)

Problem 2: Transient temperature field in one-dimensional beam

In this problem we shall solve a parabolic PDE. We are going to find the temperature field along a one-dimensional beam as function of time (transient temperature field). The temperature field is known in the initial state $ t=0 $, and the temperature is known at the boundary.

This problem is goverened by the one-dimensional transient heat conduction equation, $$ \begin{equation} \frac{\partial T}{\partial t}=k \frac{\partial^2 T}{\partial x^2} \label{eq:pde} \end{equation} $$ where $ k $ is the heat conduction coefficient. In this exercise we set $ k=1 $. The beam has length $ L=1 $.

Boundary conditions: $ T(0,t)=100 $ and $ T(1,t)=0 $.

Initial conditions: $ T(x,0)=0 $ Find the temperature field along the beam as function of time. Discretize the length with a suitable number of nodes.

a) Discretize \eqref{eq:pde} with the FTCS scheme and and write a script that solves the problem with given boundary conditions. Explore with different values for $ r $ (the numerical Fourier number), e.g. $ r=0.2 $, $ r=0.5 $ and $ r=1.0 $. Plot the temperature field for $ t=0.001 $, $ t=0.025 $ and $ t=0.4 $.

b) Same as a), but now using the $ \theta $-scheme instead of the FTCS scheme. Use both $ \theta = \frac{1}{2} $ (Crank) and $ \theta = 1 $ (Laasonen). (Due to the canceled lecture on 28th of February, the curriculum for the $ \theta $-scheme was not covered before the deadline of this exercise. Because of this you may still get this exercise approved even though you are not able to complete this last sub-exercise.)