a Numerical Methods with MATLABfor Complex Systems Auralius Manurung Mechanical Engineering, Universitas Pertamina

Table of Contents

1. Case 1: One-Dimensional Linear Advection 1.1. Forward in Time Centered in Space (FTCS) 1.2. Upwind Method (c>0)$(c>0)$ 1.3. Downwind Method (c<0)$(c<0)$ 1.4. The Lax Method 1.5. The Lax-Wendroff Method 1.6. Cross Current Heat Exchanger as An Example 1.6.1. Two-Equation Model 1.6.2. Three-Equation Model 2. Case 2: One-Dimensional Heat Equation 2.1. Explicit Method 2.1.1. Dirichlet Boundary Conditions 2.1.2. Neumann Boundary Condition 2.1.3. MATLAB Implementation 2.1.4. Example 1 2.1.5. Example 2 2.1.6. Example 3 2.2. Implicit Method 3. Two-Dimensional Poisson Equation 3.1. Dirichlet Boundary 3.2. Neumann Boundary 4. Two-Dimensional Unsteady Advection-Diffusion 4.1. UPWIND-UPWIND Method 4.2. DOWNWIND-DOWNWIND Method 4.3. DOWNWIND-UPWIND Method 4.4. UPWIND-DOWNWIND Method 4.5. MATLAB Implementation 4.6. Neumann Boundary Condition 4.6.1. Left Edge 4.6.2. Right Edge 4.6.3. Top Edge 4.6.4. Bottom Edge 4.7. Energy Balance Model: Out-of-Plane Heat Transfer for a Thin Plate 4.7.1. Transient Case 4.7.2. Stationary Case 5. Navier-Stokes with Heat Transfer 5.1. Channel Flow with Heat Equation 5.2. Two Counter-Direction Channel Flows with Heat Equations 6. Navier-Stokes in A 2D Rectangular World with Explicit Pressure Equation 6.1. Discretizing the x$x$-momentum Equation 6.2. Discretizing the y$y$-momentum Equation 6.3. Discretizing the Pressure Equation 6.4. Cavity Flow 6.5. Channel Flow 6.5.1. Example 1 6.5.2. Example 2 7. Steady Simple Euler-Bernoulli Beam 7.1. Mathematical Modeling 7.2. Engineering Interpretation 7.2.1. Boundary Conditions at the Left End 7.2.2. Boundary Conditions at the Right End 7.3. System Matrices for Different Boundaries 7.3.1. Fixed-fixed (fixed ended) 7.3.2. Pinned-pinned (simply supported) 7.3.3. Fixed-pinned (cantilever, simply supported) 7.3.4. Pinned-fixed (simply supported, cantilever) 7.3.5. Fixed-free (cantilever) 7.3.6. Fixed-explicit moment 8. Unsteady Simple Euler-Bernoulli Beam

1. Case 1: One-Dimensional Linear Advection1 Given a hyperbolic PDE as follows. (1)ut+cux=f(x,t)$$u_t+c{u}_x=f(x,t)$$with the following initial condition:(2)u(x,0)=u0(x)$$u(x,0)=u_0(x)$$where c${\displaystyle c}$ is a positive real number. Here, we want to to find the solution to (1). When f(x,t)=0${\displaystyle f(x,t)=0}$, its analytical solution is given by: u(x,t)=u0(x-ct)${\displaystyle u(x,t)=u_0(x-ct)}$. Simply put, u0${\displaystyle u_0}$ propagates along x${\displaystyle x}$ with a constant speed of c${\displaystyle c}$. This kind of PDE is widely used for real-time modeling and control applications. Therefore, it is essential to investigate its non-homogenous version to allow the implementation of a PDE system further. A PDE system is composed of more than one interacting PDE. For example, there are at least two PDEs in a heat exchanger: one for the hot system and another one for the cold system. The two PDEs interact by exchanging heat between them.

1.1. Forward in Time Centered in Space (FTCS) This method is a finite difference method but with central difference for the distance to increase the approximation. To approximate ut${\displaystyle u_t}$, we use the following approximation.(3)ut≈u(x,t+𝛥t)-u(x,t)

𝛥t$$u_t\approx \frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}$$While to approach ux${\displaystyle u_x}$, we use the following approximation.(4)ux≈u(x+𝛥x,t)-u(x-𝛥x,t)

2𝛥x$$u_x\approx \frac{u(x+𝛥x,t)-u(x-𝛥x,t)}{2𝛥x}$$Next, substituting (3) and (4) back to (1) gives usa

ut	≈-cux+f(x,t)
u(x,t+𝛥t)-u(x,t) 𝛥t	≈-cu(x+𝛥x,t)-u(x-𝛥x,t) 2𝛥x+f(x,t)
u(x,t+𝛥t)-u(x,t)	≈-c𝛥t 2𝛥xau(x+𝛥x,t)-u(x-𝛥x,t)a+f(x,t)𝛥t
u(x,t+𝛥t)	≈u(x,t)-c𝛥t 2𝛥xau(x+𝛥x,t)-u(x-𝛥x,t)a+f(x,t)𝛥t

$$\begin{array}{rl}{u}_t& \approx -c{u}_x+f(x,t)\\ \frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}& \approx -c\frac{u(x+𝛥x,t)-u(x-𝛥x,t)}{2𝛥x}+f(x,t)\\ u(x,t+𝛥t)-u(x,t)& \approx -\frac{c𝛥t}{2𝛥x}[u(x+𝛥x,t)-u(x-𝛥x,t)]+f(x,t)𝛥t\\ u(x,t+𝛥t)& \approx u(x,t)-\frac{c𝛥t}{2𝛥x}[u(x+𝛥x,t)-u(x-𝛥x,t)]+f(x,t)𝛥t\end{array}$$ FTCS method is unconditionally unstable. It should never be used. Nevertheless, FTCS is the fundamental building block for another methods, such as the Lax method.

1.2. Upwind Method (c>0)${\displaystyle (c>0)}$ By using finite difference method, we can approximate both ut${\displaystyle u_t}$ and ux${\displaystyle u_x}$. To approximate ut${\displaystyle u_t}$, we use the following approximation.(5)ut≈u(x,t+𝛥t)-u(x,t)

𝛥t$$u_t\approx \frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}$$While to approach ux${\displaystyle u_x}$, we use the following approximation.(6)ux≈u(x,t)-u(x-𝛥x,t)

𝛥x$$u_x\approx \frac{u(x,t)-u(x-𝛥x,t)}{𝛥x}$$By substituting (5) and (6) to (1), we then have: a

ut	≈-cux+f(x,t)
u(x,t+𝛥t)-u(x,t) 𝛥t	≈-cu(x,t)-u(x-𝛥x,t) 𝛥x+f(x,t)
u(x,t+𝛥t)+u(x,t)	≈-c𝛥t 𝛥xau(x,t)-u(x-𝛥x,t)a+f(x,t)
u(x,t+𝛥t)	≈u(x,t)-c𝛥t 𝛥xau(x,t)-u(x-𝛥x,t)a+f(x,t)𝛥t

$$\begin{array}{rl}{u}_t& \approx -c{u}_x+f(x,t)\\ \frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}& \approx -c\frac{u(x,t)-u(x-𝛥x,t)}{𝛥x}+f(x,t)\\ u(x,t+𝛥t)+u(x,t)& \approx -\frac{c𝛥t}{𝛥x}[u(x,t)-u(x-𝛥x,t)]+f(x,t)\\ u(x,t+𝛥t)& \approx u(x,t)-\frac{c𝛥t}{𝛥x}[u(x,t)-u(x-𝛥x,t)]+f(x,t)𝛥t\end{array}$$Let us consider c𝛥t

𝛥x=r${\displaystyle \frac{c𝛥t}{𝛥x}=r}$, Uj,k=u(xj,tk)${\displaystyle U_{j,k}=u(x_j,t_k)}$, Fj,k=f(xj,tk)${\displaystyle F_{j,k}=f(x_j,t_k)}$, xj=j𝛥x${\displaystyle x_j=j𝛥x}$ and tk=k𝛥t${\displaystyle t_k=k𝛥t}$, we then have the following shorter terms:(7)

Uj,k+1=Uj,k-r(Uj,k-Uj-1,k)+Fj,k𝛥t

$$\boxed{U_{j,k+1}=U_{j,k}-r(U_{j,k}-U_{j-1,k})+F_{j,k}𝛥t}$$ MATLAB implementation of the upwind method (7) is as follows.

1function next_u_array = upwind(u_array, f_array, c, dt, dx)2 3N = length(u_array);4next_u_array = u_array;5 6for i = 2:N7 next_u_array(i) = u_array(i)-c*dt/dx*(u_array(i)-u_array(i-1))+f_array(i)*dt;;8end9 10end

To increase the efficiency, the listing above can be vectorized as follows.

1function next_u_array = upwind(u_array, f_array, c, dt, dx)2 3N = length(u_array);4next_u_array = u_array;5 6next_u_array(2:N) = u_array(2:N) - (c*dt/dx).* ...7 (u_array(2:N)- u_array(1:N-1)) + f_array(2:N).*dt; 8 9end

Next, we will implement the PDE in (8) into the MATLAB upwind function above. ut+0.5ux=0(8)u(x,0)=e-(x-2)20⩽x⩽10$$\begin{array}{c}{u}_t+0.5u_x=0\\ u(x,0)=e^{-(x-2{)}^2}\\ 0⩽x⩽10\end{array}$$The implementation is provided in the listing below. The length of the rod is divided into 100 segments (dx=0.1${\displaystyle dx=0.1}$) and the time is divided into 400 segments (dt=0.05${\displaystyle dt=0.05}$).

1L = 10;2dx = 0.1;3x = 0:dx:L;4 5dt = 0.05;6T = 20;7t = 0:dt:T;8 9c = 0.5; % Upwind10 11u = zeros(length(x),1);12f = zeros(length(x),1);13 14for k = 1:length(x)15 u(k) = exp(-(x(k)-2).^2);16end17 18figure;19h = plot(0,0);20title('Upwind')21ylim([0 1]);22 23for k = 1:length(t)24 u = upwind(u,f, c, dt, dx);25 set(h,'XData',x,'Ydata',u);26 drawnow;27end

The analytical solution to the PDE on (8) is u(x,t)=e-((x-0.5t)-2)2${\displaystyle u(x,t)=e^{-((x-0.5t)-2{)}^2}}$. However, when executed, the code above will generate a numerical solution as in Figure 1. As we can see, the numerical solution is not exactly the same as its analytical solution. Over time, its amplitude is decreasing.

Figure 1:Upwind method for ut+0.5uxx=0${\displaystyle u_t+0.5u_{xx}=0}$

1.3. Downwind Method (c<0)${\displaystyle (c<0)}$ By using finite difference method, we can approximate both ut${\displaystyle u_t}$ and ux${\displaystyle u_x}$. To approximate ut${\displaystyle u_t}$, we use the following approximation.(9)ut≈u(x,t+𝛥t)-u(x,t)

𝛥t$$u_t\approx \frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}$$While to approach ux${\displaystyle u_x}$, we use the following approximation.(10)ux≈u(x+𝛥x,t)-u(x,t)

𝛥x$$u_x\approx \frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}$$Substituting (9) and (10) back to (1) yields the following.a

ut	=-cux+f(x,t)
u(x,t+𝛥t)-u(x,t) 𝛥t	≈-cu(x+𝛥x,t)-u(x,t) 𝛥x+f(x,t)
u(x,t+𝛥t)+u(x,t)	≈-c𝛥t 𝛥xau(x+𝛥x,t)-u(x,t)a+f(x,t)𝛥t
u(x,t+𝛥t)	≈u(x.t)-c𝛥t 𝛥xau(x+𝛥x,t)-u(x,t)a+f(x,t)𝛥t

$$\begin{array}{rl}{u}_t& =-c{u}_x+f(x,t)\\ \frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}& \approx -c\frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}+f(x,t)\\ u(x,t+𝛥t)+u(x,t)& \approx -\frac{c𝛥t}{𝛥x}[u(x+𝛥x,t)-u(x,t)]+f(x,t)𝛥t\\ u(x,t+𝛥t)& \approx u(x.t)-\frac{c𝛥t}{𝛥x}[u(x+𝛥x,t)-u(x,t)]+f(x,t)𝛥t\end{array}$$Let us consider c𝛥t

𝛥x=r${\displaystyle \frac{c𝛥t}{𝛥x}=r}$, Uj,k=u(xj,tk)${\displaystyle U_{j,k}=u(x_j,t_k)}$, Fj,k=f(xj,tk)${\displaystyle F_{j,k}=f(x_j,t_k)}$, xj=j𝛥x${\displaystyle x_j=j𝛥x}$ and tk=k𝛥t${\displaystyle t_k=k𝛥t}$, we can then have the following shorter terms.(11)

Uj,k+1=Uj,k-r(Uj+1,k-Uj,k)+Fj,k𝛥t

$$\boxed{U_{j,k+1}=U_{j,k}-r(U_{j+1,k}-U_{j,k})+F_{j,k}𝛥t}$$ MATLAB implementation of an upwind method expressed in (11) is as follows.

1function next_u_array = downwind(u_array, f_array, c, dt, dx)2 3N = length(u_array);4next_u_array = u_array;5 6for i = 1:N-17 next_u_array(i) = u_array(i)-c*dt/dx*(u_array(i+1)-u_array(i))+f_array(i)*dt;8end9 10end

Further, to increase its efficiency, we can turn the MATLAB code above in to vectorized operations as follows.

1function next_u_array = downwind(u_array, f_array, c, dt, dx)2 3N = length(u_array);4next_u_array = u_array;5 6next_u_array(1:N-1) = u_array(1:N-1) - (c*dt/dx).* ...7 (u_array(2:N) - u_array(1:N-1)) + f_array(1:N-1).*dt; 8 9end

To test the listing above, let us solve the following first order linear PDE.ut-0.5ux=0(12)u(x,0)=e-(x-8)20⩽x⩽10$$\begin{array}{c}{u}_t-0.5u_x=0\\ u(x,0)=e^{-(x-8{)}^2}\\ 0⩽x⩽10\end{array}$$The implementation is provided in the listing below. The length of the rod is divided into 100 segments (dx=0.1${\displaystyle dx=0.1}$) and the time is divided into 400 segments (dt=0.05${\displaystyle dt=0.05}$).

1L = 10;2dx = 0.1;3x = 0:dx:L;4 5dt = 0.05;6T = 20;7t = 0:dt:T;8 9c = -0.5; % Downwind10 11u = zeros(length(x), 1);12f = zeros(length(x), 1);13 14for k = 1:length(x)15 u(k) = exp(-(x(k)-8).^2);16end17 18figure;19h = plot(0,0);20title('Downwind')21ylim([0 1]);22 23for k = 1:length(t)24 u = downwind(u,f, c, dt, dx); 25 set(h,'XData',x,'Ydata',u);26 drawnow;27end

The analytical solution of (12) is u(x,t)=e-((x+0.5t)-8)2${\displaystyle u(x,t)=e^{-((x+0.5t)-8{)}^2}}$. However, when executed, the provided listing above generates the a solutoinas in Figure 2. As we can see, the numerical solution is not exactly the same as its analytical solution. Over time, its amplitude is actually decreasing, similar to the upwind cases.

Figure 2:Downwind method for ut-0.5uxx=0${\displaystyle u_t-0.5u_{xx}=0}$

1.4. The Lax MethodThe Lax method is based on the FTCS method. (13)ut≈u(x,t+𝛥t)-

u(x,t)

𝛥t$$u_t\approx \frac{u(x,t+𝛥t)-\boxed{u(x,t)}}{𝛥t}$$However, u(x,t)${\displaystyle u(x,t)}$ is replaced by an average of its two neighbors.(14)u(x,t)=u(x+𝛥x,t)+u(x-𝛥x,t)

2$$u(x,t)=\frac{u(x+𝛥x,t)+u(x-𝛥x,t)}{2}$$ Thus, by substituting (14) to (13), the approximation of ut${\displaystyle u_t}$ can be rewritten as(15)ut≈u(x,t+𝛥t)-au(x+𝛥x,t)+u(x-𝛥x,t)

2a

𝛥t$$u_t\approx \frac{u(x,t+𝛥t)-\left[\frac{u(x+𝛥x,t)+u(x-𝛥x,t)}{2}\right]}{𝛥t}$$On the other hand, ux${\displaystyle u_x}$ can be approximated as(16)ux≈u(x+𝛥x,t)-u(x-𝛥x,t)

2𝛥x$$u_x\approx \frac{u(x+𝛥x,t)-u(x-𝛥x,t)}{2𝛥x}$$Substituting (15) and (16) back to (1) yieldsa

ut	≈-cux+f(x,t)
u(x,t+𝛥t)-au(x+𝛥x,t)+u(x-𝛥x,t) 2a 𝛥t	≈-cu(x+𝛥x,t)-u(x-𝛥x,t) 2𝛥x+f(x,t)
u(x,t+𝛥t)-au(x+𝛥x,t)+u(x-𝛥x,t) 2a	≈-c𝛥t 2𝛥xau(x+𝛥x,t)-u(x-𝛥x,t)a+f(x,t)𝛥t
u(x,t+𝛥t)	≈au(x+𝛥x,t)+u(x-𝛥x,t) 2a-
	c𝛥t 2𝛥xau(x+𝛥x,t)-u(x-𝛥x,t)a+f(x,t)𝛥t

$$\begin{array}{rl}{u}_t& \approx -c{u}_x+f(x,t)\\ \frac{u(x,t+𝛥t)-\left[\frac{u(x+𝛥x,t)+u(x-𝛥x,t)}{2}\right]}{𝛥t}& \approx -c\frac{u(x+𝛥x,t)-u(x-𝛥x,t)}{2𝛥x}+f(x,t)\\ u(x,t+𝛥t)-\left[\frac{u(x+𝛥x,t)+u(x-𝛥x,t)}{2}\right]& \approx -\frac{c𝛥t}{2𝛥x}[u(x+𝛥x,t)-u(x-𝛥x,t)]+f(x,t)𝛥t\\ u(x,t+𝛥t)& \approx \left[\frac{u(x+𝛥x,t)+u(x-𝛥x,t)}{2}\right]-\\ & \frac{c𝛥t}{2𝛥x}[u(x+𝛥x,t)-u(x-𝛥x,t)]+f(x,t)𝛥t\end{array}$$Now, let us consider r=c𝛥t

𝛥x${\displaystyle r=\frac{c𝛥t}{𝛥x}}$, Uj,k=u(xj,tk)${\displaystyle U_{j,k}=u(x_j,t_k)}$, Fj,k=f(xj,tk)${\displaystyle F_{j,k}=f(x_j,t_k)}$, xj=j𝛥x${\displaystyle x_j=j𝛥x}$, tk=k𝛥t${\displaystyle t_k=k𝛥t}$, Uj,k+1=u(xj,tk+𝛥t)${\displaystyle U_{j,k+1}=u(x_j,t_k+𝛥t)}$ and Uj-1,k=u(xj-𝛥x,tk)${\displaystyle U_{j-1,k}=u(x_j-𝛥x,t_k)}$, we can then have the following shorter terms.(17)

Uj,k+1=Uj+1,k+Uj-1,k

2-r

2(Uj+1,k-Uj,k)+Fj,k𝛥t

$$\boxed{U_{j,k+1}=\frac{U_{j+1,k}+U_{j-1,k}}{2}-\frac{r}{2}(U_{j+1,k}-U_{j,k})+F_{j,k}𝛥t}$$ The Lax method can handle both upwind and downwind direction. There is no need for us to check the sign of c${\displaystyle c}$ as in the previous methods. The listing below shows MATLAB implementation of the Lax method.

1function next_u_array = lax(u_array, f_array, c, dt, dx)2 3N = length(u_array);4next_u_array = u_array;5 6r = c*dt/dx;7 8for i = 2:N-1 % i = 1 and i = N are untouched !!!9 next_u_array(i) = 0.5*(u_array(i+1)+u_array(i-1)) - ...10 0.5*r*(u_array(i+1)-u_array(i-1)) + ...11 f_array(i)*dt;12end13 14% fill up the two ends by using upwind or downwind method 15if c > 016 next_u_array(N) = u_array(N) - r* ...17 (u_array(N)- u_array(N-1)) + f_array(N).*dt; 18 19elseif c < 020 next_u_array(1) = u_array(1) - r* ...21 (u_array(2) - u_array(1)) + f_array(1).*dt; 22 23end

There is another practical issue with the Lax method. The Lax method leaves the two end segments untouched. This can be seen from the for loop statement that starts from 2 to N-1 (see the above listing, line 8). To address this issue, we can apply the upwind and downwind method for the two untouched ends. Further, to increase the efficiency of the Lax method, we can turn the MATLAB code above in to vectorized operations as follows.

1function next_u_array = lax(u_array, f_array, c, dt, dx)2 3N = length(u_array);4next_u_array = u_array;5 6r = c*dt/dx;7 8next_u_array(2:N-1) = 0.5.*(u_array(3:N)+(u_array(1:N-2))) - ...9 (0.5*r).* (u_array(3:N)- u_array(1:N-2)) + ...10 f_array(2:N-1).*dt; 11 12 13% fill up the two ends by using upwind or downwind method 14if c > 015 next_u_array(N) = u_array(N) - r* ...16 (u_array(N)- u_array(N-1)) + f_array(N).*dt; 17 18elseif c < 019 next_u_array(1) = u_array(1) - r* ...20 (u_array(2) - u_array(1)) + f_array(1).*dt; 21 22end

Next, let us test the code above for both the upwind and the downwind cases. For the upwind case, we will use the PDE in (8). As for the downwind direction, we will use the PDE in (12). For both methods, the length of the rod is divided into 100 segments (dx=0.1${\displaystyle dx=0.1}$) and the time is divided into 400 segments (dt=0.05${\displaystyle dt=0.05}$). The listing below is the implementation of (8) and (12).

1L = 10;2dx = 0.05;3x = 0:dx:L;4 5dt = 0.05;6T = 20;7t = 0:dt:T;8 9u = zeros(length(x), 1);10f = zeros(length(x), 1);11 12%% --------------------------------------------------------------13c = 0.5; % moving to right 14 15for k = 1:length(x)16 u(k) = exp(-10*(x(k)-2).^2);17end18 19figure;20h1 = plot(0,0);21title('Lax, moving to the right c = 0.5')22ylim([0 1]);23 24for k = 1:length(t)25 u = lax(u,f, c, dt, dx);26 set(h1,'XData',x,'Ydata',u);27 drawnow;28end29 30%% --------------------------------------------------------------31c = -0.5; % moving to left32 33for k = 1:length(x)34 u(k) = exp(-10*(x(k)-8).^2);35end36 37figure;38h2 = plot(0,0);39title('Lax, moving to the left, c = -0.5')40ylim([0 1]);41 42for k = 1:length(t)43 u = lax(u,f, c, dt, dx);44 set(h2,'XData',x,'Ydata',u);45 drawnow;46end

When executed, the above listing will generate Figure 3 and Figure 4 where we can see the resulting amplitude is also decreasing over time.

Figure 3:Lax method for ut-0.5uxx=0${\displaystyle u_t-0.5u_{xx}=0}$

Figure 4:Lax method for ut+0.5uxx=0${\displaystyle u_t+0.5u_{xx}=0}$

1.5. The Lax-Wendroff Method2 We start from the Taylor series definition off(x+𝛥x)${\displaystyle f(x+𝛥x)}$(18)f(x+𝛥x)=f(x)+𝛥xf′(x)+𝛥x2

2!f′′(x)+𝛥x3

3!f′′′(x)+…$$f(x+𝛥x)=f(x)+𝛥x{f}^{\prime }(x)+\frac{𝛥x^2}{2!}{f}^{\prime \prime }(x)+\frac{𝛥x^3}{3!}{f}^{\prime \prime \prime }(x)+\dots$$Thus, taking only up to the second order, we can write(19)a

u(x+𝛥x,t)

=u(x,t)+𝛥xux+𝛥x2 2uxx

$$\begin{array}{rl}u(x+𝛥x,t)& =u(x,t)+𝛥x{u}_x+\frac{𝛥x^2}{2}{u}_{xx}\end{array}$$Remember that(20)ut=u(x,t+𝛥t)-u(x,t)

𝛥t$$u_t=\frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}$$and (21)ux=u(x+𝛥x,t)-u(x,t)

𝛥x$$u_x=\frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}$$Thus, ut=-cux+f(x,t)${\displaystyle u_t=-c{u}_x+f(x,t)}$ can be rewritten as (22)u(x,t+𝛥t)-u(x,t)

𝛥t=-cu(x+𝛥x,t)-u(x,t)

𝛥x+f(x,t)$$\frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}=-c\frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}+f(x,t)$$By substituting (19) to (22)a

u(x,t+𝛥t)-u(x,t) 𝛥t	≈-cu(x+𝛥x,t)-u(x,t) 𝛥x+f(x,t)
u(x,t+𝛥t)	≈u(x,t)-c𝛥t 𝛥x[u(x+𝛥x,t)-u(x,t)]+f(x,t)𝛥t
u(x,t+𝛥t)	≈u(x,t)-c𝛥t 𝛥x[u(x,t)+𝛥xux+𝛥x2 2uxx-u(x,t)]+f(x,t)𝛥t
u(x,t+𝛥t)	≈u(x,t)-c𝛥tux-c𝛥t𝛥x 2uxx+f(x,t)𝛥t

$$\begin{array}{rl}\frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}& \approx -c\frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}+f(x,t)\\ u(x,t+𝛥t)& \approx u(x,t)-c\frac{𝛥t}{𝛥x}[u(x+𝛥x,t)-u(x,t)]+f(x,t)𝛥t\\ u(x,t+𝛥t)& \approx u(x,t)-c\frac{𝛥t}{𝛥x}[\cancel{u(x,t)}+𝛥x{u}_x+\frac{𝛥x^2}{2}{u}_{xx}-\cancel{u(x,t)}]+f(x,t)𝛥t\\ u(x,t+𝛥t)& \approx u(x,t)-c𝛥t{u}_x-c\frac{𝛥t𝛥x}{2}{u}_{xx}+f(x,t)𝛥t\end{array}$$Since c${\displaystyle c}$ is the propagation speed of u(x,t)=u0(x-ct)${\displaystyle u(x,t)=u_0(x-ct)}$, thus, let us take 𝛥x=-c𝛥t${\displaystyle 𝛥x=-c𝛥t}$. As the result, now we have the following equation.(23)

u(x,t+𝛥t)=≈u(x,t)-c𝛥tux+c2𝛥t2

2uxx+f(x,t)𝛥t

$$\boxed{u(x,t+𝛥t)=\approx u(x,t)-c𝛥t{u}_x+c^2\frac{𝛥t^2}{2}{u}_{xx}+f(x,t)𝛥t}$$ Next, we need to approximate ux${\displaystyle u_x}$ and uxx${\displaystyle u_{xx}}$. For ux${\displaystyle u_x}$, we use central difference method as follows.(24)ux≈u(x+𝛥x,t)-u(x-𝛥x,t)

2𝛥x$$u_x\approx \frac{u(x+𝛥x,t)-u(x-𝛥x,t)}{2𝛥x}$$As for uxx${\displaystyle u_{xx}}$, we do not use the central difference method, instead, we take ux≈u(x+𝛥x,t)-u(x,t)

𝛥x${\displaystyle u_x\approx \frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}}$. Therefore, we can express uxx${\displaystyle u_{xx}}$ as follows.

	uxx	≈ux(x+𝛥x,t)-ux(x,t) 𝛥x
		≈u(x+𝛥x,t)-u(x,t) 𝛥x-u(x,t)-u(x-𝛥x,t) 𝛥x 𝛥x
		≈u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t) 𝛥x2(25)

$$\begin{array}{rl}{u}_{xx}& \approx {\displaystyle \frac{u_x(x+𝛥x,t)-u_x(x,t)}{𝛥x}}\\ & \approx {\displaystyle \frac{\frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}-\frac{u(x,t)-u(x-𝛥x,t)}{𝛥x}}{𝛥x}}\\ & \approx {\displaystyle \frac{u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)}{𝛥x^2}}\end{array}$$Now that we have (24) and (25), we can substitute them back to (23):

	u(x,t+𝛥t)	≈u(x,t)-c𝛥tux+c2𝛥t2 2uxx+f(x,t)𝛥t
		a≈u(x,t)-c𝛥t[u(x+𝛥x,t)-u(x-𝛥x,t) 2𝛥x]+ c2𝛥t2 2[u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t) 𝛥x2]+f(x,t)𝛥t
		a≈u(x,t)-c𝛥t 2𝛥x[u(x+𝛥x,t)-u(x-𝛥x,t)]+ c2𝛥t2 2𝛥x2[u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)]+f(x,t)𝛥t

$$\begin{array}{rl}u(x,t+𝛥t)& \approx u(x,t)-c𝛥t{u}_x+c^2\frac{𝛥t^2}{2}{u}_{xx}+f(x,t)𝛥t\\ & \begin{array}{c}\approx u(x,t)-c𝛥t\left[\frac{u(x+𝛥x,t)-u(x-𝛥x,t)}{2𝛥x}\right]+\\ c^2\frac{𝛥t^2}{2}\left[{\displaystyle \frac{u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)}{𝛥x^2}}\right]+f(x,t)𝛥t\end{array}\\ & \begin{array}{c}\approx u(x,t)-\frac{c𝛥t}{2𝛥x}[u(x+𝛥x,t)-u(x-𝛥x,t)]+\\ \frac{c^2𝛥t^2}{2𝛥x^2}[u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)]+f(x,t)𝛥t\end{array}\end{array}$$Now let us consider r=c𝛥t

𝛥x${\displaystyle r=\frac{c𝛥t}{𝛥x}}$, Uj,k=u(xj,tk)${\displaystyle U_{j,k}=u(x_j,t_k)}$, Fj,k=f(xj,tk)${\displaystyle F_{j,k}=f(x_j,t_k)}$, xj=j𝛥x${\displaystyle x_j=j𝛥x}$, tk=k𝛥t${\displaystyle t_k=k𝛥t}$, Uj,k+1=u(xj,tk+𝛥t)${\displaystyle U_{j,k+1}=u(x_j,t_k+𝛥t)}$ and Uj-1,k=u(xj-𝛥x,tk)${\displaystyle U_{j-1,k}=u(x_j-𝛥x,t_k)}$, we can then have the following shorter terms.(26)

Uj,k+1=Uj,k-r

2(Uj+1,k-Uj-1,k)+r2

2(Uj+1,k-2Uj,k+Uj-1,k)+Fj,k𝛥t

$$\boxed{U_{j,k+1}=U_{j,k}-\frac{r}{2}(U_{j+1,k}-U_{j-1,k})+\frac{r^2}{2}(U_{j+1,k}-2U_{j,k}+U_{j-1,k})+F_{j,k}𝛥t}$$ MATLAB implementation of (26) is as follows (see line 6 to 11).

1function next_u_array = lax_wendroff(u_array, f_array, c, dt, dx)2 3N = length(u_array);4next_u_array = u_array;5 6r = c*dt/dx;7for i = 2:N-18 next_u_array(i) = u_array(i)-0.5*r*(u_array(i+1)-u_array(i-1)) ...9 + 0.5*r^2 * (u_array(i+1)-2*u_array(i)+u_array(i-1)) ... 10 + f_array(i)*dt;11end12 13% fill up the two ends by using upwind or downwind method 14if c > 015 next_u_array(N) = u_array(N) - r* ...16 (u_array(N)- u_array(N-1)) + f_array(N)*dt; 17 18elseif c < 019 next_u_array(1) = u_array(1) - r* ...20 (u_array(2) - u_array(1)) + f_array(1)*dt; 21 22end

Here, we encounter the same issue as in the previous Lax method. The Lax-Wendroff method also leaves the two end segments untouched. This can be seen clearly from the for-statement that starts from 2 to N-1 (see line 7 of the listing above). To address this issue, we apply upwind and downwind method for the two ends only (see line 14 to 20 of the listing above). Further, to increase its efficiency, we can turn the MATLAB code above in to vectorized operations as follows.

1function next_u_array = lax_wendroff(u_array, f_array, c, dt, dx)2 3N = length(u_array);4next_u_array = u_array;5 6r = c*dt/dx;7 8next_u_array(2:N-1) = u_array(2:N-1)-0.5*r*(u_array(3:N)-u_array(1:N-2)) ...9 + 0.5*r^2 * (u_array(3:N)-2*u_array(2:N-1)+u_array(1:N-2)) ... 10 + f_array(2:N-1)*dt;11 12% fill up the two ends by using upwind or downwind method 13if c > 014 next_u_array(N) = u_array(N) - r* ...15 (u_array(N)- u_array(N-1)) + f_array(N)*dt; 16 17elseif c < 018 next_u_array(1) = u_array(1) - r* ...19 (u_array(2) - u_array(1)) + f_array(1)*dt; 20 21end

Next, let us test the code listings above for both the upwind and the downwind cases. For the upwind case, we will use the same system as in (8) As for the downwind direction, we will also use the same system as in (12). We will divide the length of the rod into 200 segments (dx=0.05)${\displaystyle (dx=0.05)}$ and the time into 200 segments as well (dt=0.05)${\displaystyle (dt=0.05)}$. See the listing below.

1L = 10;2dx = 0.05;3x = 0:dx:L;4 5dt = 0.05;6T = 20;7t = 0:dt:T;8 9u = zeros(length(x), 1);10f = zeros(length(x), 1);11 12%% ------------------------------------------------------------13c = 0.5; % moving to right14 15for k = 1:length(x)16 u(k) = exp(-10*(x(k)-2).^2);17end18 19h = figure;20h1 = plot(0,0);21title('Lax-Wendroff, moving to the right c = 0.5')22ylim([0 1]);23 24for k = 1:length(t)25 u = lax_wendroff(u,f, c, dt, dx);26 set(h1,'XData',x,'Ydata',u);27 drawnow28end29 30%% ------------------------------------------------------------31c = -0.5; % moving to left32 33for k = 1:length(x)34 u(k) = exp(-10*(x(k)-8).^2);35end36 37h = figure;38h1 = plot(0,0);39title('Lax-Wendroff, moving to the left, c = -0.5')40ylim([0 1]);41 42for k = 1:length(t)43 u = lax_wendroff(u,f, c, dt, dx);44 set(h1,'XData',x,'Ydata',u);45 drawnow46end

When executed, the listing above will generate the Figure 5 and Figure 6. As shown by those two figures, the, Lax-Wendroff method does provide more stable solutions. Even tough the amplitudes are still decreasing, the decrease in the amplitude is much smaller than other methods that have been previously discussed.

Figure 5:Lax method for ut-0.5uxx=0${\displaystyle u_t-0.5u_{xx}=0}$

Figure 6:Lax method for ut+0.5uxx=0${\displaystyle u_t+0.5u_{xx}=0}$

1.6. Cross Current Heat Exchanger as An Example3A heat exchanger is modeled with two or more interacting non-homogenous PDEs. The simplest model involves two PDEs only: the hot flow and the cold flow. Yet another simple method, three PDEs are involved: the hot flow, the cold flow, and the wall between them. 1.6.1. Two-Equation ModelFlowing liquid inside a container can be modeled as an advective PDE. When external environment is involved, the PDE becomes non-homogenous (see Drawing 1). Interaction with the environment is modeled with Fourier's Law. The overall process can be modeled as follows.(27)𝜕T(x,t)

𝜕t+v𝜕T(x,t)

𝜕z⏠⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏢Advection=-UAS

𝜌ACp(T(x,t)-TS(x,t))⏠⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏢Fourier's law$$\underset{\text{Advection}}{\underbrace{\frac{𝜕T(x,t)}{𝜕t}+v\frac{𝜕T(x,t)}{𝜕z}}}=-\underset{\text{Fourier's law}}{\underbrace{\frac{U{A}_S}{𝜌A{C}_p}(T(x,t)-T_S(x,t))}}$$whereCp${\displaystyle C_p}$ is the liquid specific heat capacityA${\displaystyle A}$ is cross sectional are of the tube where the liquid is flowingAS${\displaystyle A_S}$ is the heat transfer area (blue surface)U${\displaystyle U}$ is the overall heat transfer coefficient between the tube to its environment𝜌${\displaystyle 𝜌}$ is the density of the liquid

As𝜌ACpUAs(Ts-T)TS(x,t)T(x,t)A

Drawing 1:Flowing liquid in red container while transferring heat to its environment

A diagram of a simplest heat exchanger is presented in Drawing 2. From Drawing 2, we can see that hot liquid flows from left to right and cold liquid flows from right to left. Source temperature of the hot flow is kept at TH0${\displaystyle T_{H0}}$. While source temperature of the cold flow is kept at TC0${\displaystyle T_{C0}}$. Those two temperatures never change. This condition is called Dirichlet condition.

L

TH(0,t)=TH0

TC(L,t)=TC0

TH(x,0)=TH0

TC(x,0)=TC0

c2c1d1d2TH(x,t)c1${\displaystyle c_1}$ and c2${\displaystyle c_2}$ are in m/s${\displaystyle m/s}$d1${\displaystyle d_1}$ and d2${\displaystyle d_2}$ are in 1

s${\displaystyle \frac{1}{s}}$Boundary conditionBoundary conditionInitial conditionInitial conditionTC(x,t)DirichletDirichlet

Drawing 2:Cross-current heat exchanger modeled with two PDEs

The system model of a heat exchanger based on Drawing 2 can be described as follows. (28)a

𝜕 𝜕tTH(x,t)+c1𝜕 𝜕xTH(x,t)	=-d1(TH(x,t)-TC(x,t))
𝜕 𝜕tTC(x,t)-c2𝜕 𝜕xTC(x,t)	=-d2(TC(x,t)-TH(x,t))

$$\begin{array}{rl}\frac{𝜕}{𝜕t}{T}_H(x,t)+c_1\frac{𝜕}{𝜕x}{T}_H(x,t)& =-d_1(T_H(x,t)-T_C(x,t))\\ \frac{𝜕}{𝜕t}{T}_C(x,t)-c_2\frac{𝜕}{𝜕x}{T}_C(x,t)& =-d_2(T_C(x,t)-T_H(x,t))\end{array}$$or(29)a

𝜕 𝜕tTH(x,t)	=-c1𝜕 𝜕xTH(x,t)-d1(TH(x,t)-TC(x,t))
𝜕 𝜕tTC(x,t)	=c2𝜕 𝜕xTC(x,t)+d2(TH(x,t)-TC(x,t))

$$\begin{array}{rl}\frac{𝜕}{𝜕t}{T}_H(x,t)& =-c_1\frac{𝜕}{𝜕x}{T}_H(x,t)-d_1(T_H(x,t)-T_C(x,t))\\ \frac{𝜕}{𝜕t}{T}_C(x,t)& =c_2\frac{𝜕}{𝜕x}{T}_C(x,t)+d_2(T_H(x,t)-T_C(x,t))\end{array}$$where: • d1${\displaystyle d_1}$ and d2${\displaystyle d_2}$ describe convective heat transfer between the the two flowing liquids: from hot liquid to cold liquid and from cold liquid to hot liquid, respectively. • c1${\displaystyle c_1}$ and c2${\displaystyle c_2}$ describes advective heat transfer inside the the hot and cold liquid, respectively. This is analogous to the actual movement of the liquids since the heat moves together with the liquids.• For the simulation, we will set the parameters according to the paper by Zobiri et.al.4 In their paper, it is mentioned that they use upwind-downwind method. However, we will solve the system above by using both the upwind-downwind method and the Lax-Wendroff method so that we can compare their results. The listing below shows the implementation of a cross-current heat exchanger as in (28) which runs for 20s${\displaystyle 20s}$ with 0.01s${\displaystyle 0.01s}$ step-time.

1L = 1;2dx = 0.2;3x = 0:dx:L;4 5dt = 0.01;6T = 20;7t = 0:dt:T;8 9% Define the parameters10c1 = 0.001711593407;11c2 = -0.01785371429;12d1 = -0.03802197802;13d2 = 0.3954285714;14 15% ======================================================================16% Upwind-downwind method17% ======================================================================18TH = zeros(length(x), 1);19TC = zeros(length(x), 1);20f = zeros(length(x), 1);21 22% Intial condition23TH(1:length(x)) = 273+30;24TC(1:length(x)) = 273+10;25 26figure;27hold on28h1 = plot(0,0,'r');29h2 = plot(0,0,'b');30title('Cross-current heat exchanger');31ylim([270 320]);32xlabel('x');33ylabel('Temperature');34legend('Hot flow','Cold flow');35 36for k = 1 : length(t) 37 set(h1,'XData',x,'Ydata',TH);38 set(h2,'XData',x,'Ydata',TC);39 drawnow;40 41 f = d1.*(TH-TC);42 TH = upwind(TH, f, c1, dt, dx);43 f = d2.*(TH-TC);44 TC = downwind(TC, f, c2, dt, dx);45end46 47% ======================================================================48% Lax-Wnderoff49% ======================================================================50TH = zeros(length(x), 1);51TC = zeros(length(x), 1);52f = zeros(length(x), 1);53 54% Intial conditoin55TH(1:length(x)) = 273+30;56TC(1:length(x)) = 273+10;57 58figure;59hold on60h3 = plot(0,0,'r');61h4 = plot(0,0,'b');62title('Cross-current heat exchanger');63ylim([270 320]);64xlabel('x');65ylabel('Temperature');66legend('Hot flow','Cold flow');67 68for k = 1:length(t)69 set(h3,'XData',x,'Ydata',TH);70 set(h4,'XData',x,'Ydata',TC);71 drawnow;72 73 f = d1.*(TH-TC);74 TH = lax_wendroff(TH, f, c1, dt, dx);75 f = d2.*(TH-TC);76 TC = lax_wendroff(TC, f, c2, dt, dx);77end

From line 41 to 44 of the listing above, we use the upwind method for the hot flow since the flow moves from left to right. As for the cold flow, we use the downwind method since the cold flow moves from right to left. As for the Lax-Wendroff method, there is only one function that we call (line 73 to 76). However, we still need to call that function twice, one for the cold flow, another one is for the hot flow. When executed, the provided lisiting above will generate animations of the temperature evolutions for both methods. Figure 7 and Figure 8 below show the temperature distribution of the modeled heat exchanger after 20 seconds for two different methods. From Figure 8, we can see that Lax-Wendroff method generates flat curve in its middle area. Theoretically, this makes more sense when compared to the results acquired from the upwind-downwind method. However, if we think with a practical perspective, the upwind-and downwind method might be more desirable since a theoretically perfect advection might be impossible in a real world.

Figure 7:After 20 s, with upwind-downwind method

Figure 8:After 20 s, with Lax-Wendroff method

1.6.2. Three-Equation Model

L

TH(0,t)=TH0

TC(L,t)=TC0

TH(x,0)=TH0

TC(x,0)=TC0

v2v1TH(x,t)TC(x,t)d1d2TW(x,t)TW(x,0)=TW0

∂

∂tTW(L,t)=0

∂

∂tTW(0,t)=0

v3NeumannNeumannDirichletDirichletd3d4

Drawing 3:Cross-current heat exchanger modeled with three PDEs

In three-equation model (see Drawing 3), we consider the wall between the hot flow and the cold flow. In general, we can summarize how the heat are being exchanged as follows.• d1${\displaystyle d_1}$ and d2${\displaystyle d_2}$ describe convective heat transfer from the the liquid to the wall: from the hot liquid to the wall and from the cold liquid to the wall, respectively. • d3${\displaystyle d_3}$ and d4${\displaystyle d_4}$ describe convective heat transfer from the the wall to the liquid: from the wall to the hot liquid and from the wall to the cold liquid, respectively. • v1${\displaystyle v_1}$ and v2${\displaystyle v_2}$ describes advective heat transfer inside the the hot and cold liquid, respectively. This is analogous to the actual movement of the liquids since the heat moves together with the liquids.• v3${\displaystyle v_3}$ describes diffusive/conductive heat transfer inside the tube wall.• Therefore, the three-equation model of a cross-current heat exchanger can be expressed as follows.(30)a

𝜕 𝜕tTH(x,t)	=-v1𝜕 𝜕xTH(x,t)-d1(TH(x,t)-TW(x,t))
𝜕 𝜕tTC(x,t)	=v2𝜕 𝜕xTC(x,t)+d2(TW(x,t)-TC(x,t))
𝜕 𝜕tTW(x,t)	=v3∂2 ∂x2TW(x,t)+d3(TH(x,t)-TW(x,t))-d4(TW(x,t)-TC(x,t))

$$\begin{array}{rl}\frac{𝜕}{𝜕t}{T}_H(x,t)& =-v_1\frac{𝜕}{𝜕x}{T}_H(x,t)-d_1(T_H(x,t)-T_W(x,t))\\ \frac{𝜕}{𝜕t}{T}_C(x,t)& =v_2\frac{𝜕}{𝜕x}{T}_C(x,t)+d_2(T_W(x,t)-T_C(x,t))\\ \frac{𝜕}{𝜕t}{T}_W(x,t)& =v_3\frac{{\partial }^2}{\partial x^2}{T}_W(x,t)+d_3(T_H(x,t)-T_W(x,t))-d_4(T_W(x,t)-T_C(x,t))\end{array}$$Where5:

v1=0.6366197724${\displaystyle v_1=0.6366197724}$v2=0.1657863990${\displaystyle v_2=0.1657863990}$v3=0.00003585526318${\displaystyle v_3=0.00003585526318}$

d1=0.9792${\displaystyle d_1=0.9792}$d2=0.2986${\displaystyle d_2=0.2986}$d3=2.3923${\displaystyle d_3=2.3923}$d4=2.8708${\displaystyle d_4=2.8708}$

The initial and boundary conditions are:TH(0,t)=360${\displaystyle T_H(0,t)=360}$ and TH(x,0)=360${\displaystyle T_H(x,0)=360}$TC(L,t)=300${\displaystyle T_C(L,t)=300}$ and TC(x,0)=300${\displaystyle T_C(x,0)=300}$TW(x,0)=330${\displaystyle T_W(x,0)=330}$, ∂TW(0,t)

∂x=0${\displaystyle \frac{\partial T_W(0,t)}{\partial x}=0}$, and ∂TW(L,t)

∂x=0${\displaystyle \frac{\partial T_W(L,t)}{\partial x}=0}$ In (30), the third equation is different than the first two equations. The third equation is a non-homogeneous one-dimensional heat equation (diffusion) with Neumann boundary conditions. We will discuss this in the next section. Nevertheless, we will directly show the implementation of the three-equation model of a heat exchanger in MATLAB. For the first and second equations, we will use upwind and downwind method, respectively (see line 51 to 55). The rest of the code can be understood after reading the next section.

1L = 1;2tF = 5;3Nx = 10;4Nt = 100;5dt = tF/Nt;6dx = L/Nx;7 8v1 = 0.6366197724;9v2 = 0.1657863990;10v3 = 0.00003585526318;11d1 = 0.9792;12d2 = 0.2986;13d3 = 2.3923;14d4 = 2.8708;15 16TH = zeros(1, Nx+1);17TC = zeros(1, Nx+1);18TW = zeros(1, Nx+1+2); % +2 phantom nodes19fH = zeros(1, Nx+1);20fC = zeros(1, Nx+1);21fW = zeros(1, Nx+1+2); % +2 phantom nodes22 23% Intial condition24TH(1:Nx+1) = 360;25TC(1:Nx+1) = 300;26TW(1:Nx+1+2) = 330;27 28h=figure;29hold on30h1 = plot(0,0,'r-*');31h2 = plot(0,0,'b-*');32h3 = plot(0,0,'m-*');33title('Cross-current heat exchanger')34ylim([270 400])35xlabel('x')36ylabel('Temperature')37legend('Hot flow', 'Cold flow', 'Wall')38 39xa = 1:Nx+1;40xb = 0:Nx+2;41for k = 1 : Nt42 set(h1,'XData',xa,'Ydata',TH);43 set(h2,'XData',xa,'Ydata',TC);44 set(h3,'XData',xb,'Ydata',TW);45 drawnow;46 47 TW = apply_bc(TW, L/Nx, ["Neumann", "Neumann"], [0, 0]); 48 49 % TW(2:end-1) : it means we ignore phantom nodes50 51 fH = -d1.*(TH-TW(2:end-1));52 TH = upwind(TH, fH, v1, dt, dx);53 54 fC = d2.*(TH-TW(2:end-1));55 TC = downwind(TC, fC, -v2, dt, dx);56 57 fW(2:end-1) = d3.*(TH-TW(2:end-1)) - d4.*(TW(2:end-1)-TC);58 TW = diffusion_1d(TW, fW, sqrt(v3), dx, dt);59end

When executed for 20 seconds, temperature distributions of different parts of the heat exchanger are shown in Drawing 4.

Phantom node∂TW(0,t)

∂x=0${\displaystyle \frac{\partial T_W(0,t)}{\partial x}=0}$Phantom node∂TW(L,t)

∂x=0${\displaystyle \frac{\partial T_W(L,t)}{\partial x}=0}$TH(0,t)=360TC(L,t)=300n-th segment

Drawing 4:Temperature distributions after 5 seconds

Figure 9:Animation of the temperature evolution of the heat exchanger

2. Case 2: One-Dimensional Heat Equation6 The equation of temperature distribution (u${\displaystyle u}$) in a one-dimensional rod (see Drawing 5) over a certain distance (x${\displaystyle x}$) and time (t${\displaystyle t}$) is given by (31). 𝛼${\displaystyle 𝛼}$ is the diffusivity coefficient. When g(x,t)=0${\displaystyle g(x,t)=0}$, the system is homogenous.(31)𝛼²uxx-ut=g(x,t)t⩾0 and 0⩽x⩽L$$\begin{array}{c}𝛼²u_{xx}-u_t=g(x,t)\\ t⩾0\text{and}0⩽x⩽L\end{array}$$when at t=0${\displaystyle t=0}$, the following initial condition is applied.u(x,0)=f(x)$$u(x,0)=f(x)$$

x=0x=Lu(x,0)=f(x)Boundary conditionBoundary condition

Drawing 5:Temperature distribution (u${\displaystyle u}$) in a very thin rod

2.1. Explicit MethodWith finite difference method, we can approach ut${\displaystyle u_t}$ and uxx${\displaystyle u_{xx}}$. To approximate ut${\displaystyle u_t}$, we use the following approximation.(32)ut≈u(x,t+𝛥t)-u(x,t)

𝛥t$$u_t\approx \frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}$$(33)ux≈u(x+𝛥x,t)-u(x,t)

𝛥x$$u_x\approx \frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}$$While to approach ux${\displaystyle u_x}$, we use the following approximation.

	uxx	≈ux(x+𝛥x,t)-ux(x,t) 𝛥x
		≈u(x+𝛥x,t)-u(x,t) 𝛥x-u(x,t)-u(x-𝛥x,t) 𝛥x 𝛥x
		≈u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t) (𝛥x)2(34)

$$\begin{array}{rl}{u}_{xx}& \approx {\displaystyle \frac{u_x(x+𝛥x,t)-u_x(x,t)}{𝛥x}}\\ & \approx {\displaystyle \frac{\frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}-\frac{u(x,t)-u(x-𝛥x,t)}{𝛥x}}{𝛥x}}\\ & \approx {\displaystyle \frac{u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)}{(𝛥x{)}^2}}\end{array}$$Substituting (32) and (33) back to (31) gives us the following equations.a

ut	≈𝛼2uxx+g(x,t)
u(x,t+𝛥t)-u(x,t) 𝛥t	≈𝛼2u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t) 𝛥x2+g(x,t)
u(x,t+𝛥t)-u(x,t)	≈𝛼2𝛥t 𝛥x2(u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t))+g(x,t)𝛥t
u(x,t+𝛥t)	≈𝛼2𝛥t 𝛥x2u(x+𝛥x,t)+(1-2𝛼2𝛥t 𝛥x2)u(x,t)+𝛼2𝛥t 𝛥x2u(x-𝛥x,t)+g(x,t)𝛥t
	≈ru(x+𝛥x,t)+(1-2r)u(x,t)+ru(x-𝛥x,t)+g(x,t)𝛥t
	≈u(x,t)+rau(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)a

$$\begin{array}{rl}{u}_t& \approx {𝛼}^2{u}_{xx}+g(x,t)\\ \frac{u(x,t+𝛥t)-u(x,t)}{𝛥t}& \approx {𝛼}^2{\displaystyle \frac{u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)}{𝛥x^2}}+g(x,t)\\ u(x,t+𝛥t)-u(x,t)& \approx \frac{{𝛼}^2𝛥t}{𝛥x^2}(u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t))+g(x,t)𝛥t\\ u(x,t+𝛥t)& \approx \frac{{𝛼}^2𝛥t}{𝛥x^2}u(x+𝛥x,t)+(1-2\frac{{𝛼}^2𝛥t}{𝛥x^2})u(x,t)+\frac{{𝛼}^2𝛥t}{𝛥x^2}u(x-𝛥x,t)+g(x,t)𝛥t\\ & \approx ru(x+𝛥x,t)+(1-2r)u(x,t)+ru(x-𝛥x,t)+g(x,t)𝛥t\\ & \approx u(x,t)+r[u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)]\end{array}$$where r=𝛼2𝛥t

(𝛥x)2${\displaystyle r=\frac{{𝛼}^2𝛥t}{(𝛥x{)}^2}}$. Now, let us consider Uj,k=u(xj,tk)${\displaystyle U_{j,k}=u(x_j,t_k)}$, Gj,k=g(xj,tk)${\displaystyle G_{j,k}=g(x_j,t_k)}$, xj=j𝛥x${\displaystyle x_j=j𝛥x}$, tk=k𝛥t${\displaystyle t_k=k𝛥t}$, Uj,k+1=u(xj,tk+𝛥t)${\displaystyle U_{j,k+1}=u(x_j,t_k+𝛥t)}$ and Uj-1,k=u(xj-𝛥x,tk)${\displaystyle U_{j-1,k}=u(x_j-𝛥x,t_k)}$, we can then have the following shorter terms.(35)a

Uj,k+1=Uj,k+raUj+1,k-2Uj,k+Uj-1,ka+Gj,k𝛥t

$$\begin{array}{c}\boxed{U_{j,k+1}=U_{j,k}+r[U_{j+1,k}-2U_{j,k}+U_{j-1,k}]+G_{j,k}𝛥t}\end{array}$$ 2.1.1. Dirichlet Boundary ConditionsIn (31), certain conditions are applied and kept true at all time. These conditions are called boundary conditions. There are several types of boundary conditions. One of them is Dirichlet boundary condition. A typical Dirichlet boundary condition would be u(0,t)=Au(L,t)=B$$\begin{array}{c}u(0,t)=A\\ u(L,t)=B\end{array}$$where A and B${\displaystyle B}$ are real numbers. The implementation of the Dirichlet boundary condition is very straightforward by directly controlling values at the two end nodes. 2.1.2. Neumann Boundary ConditionAnother type of boundary condition is the Neumann boundary condition. In Neumann boundary condition, the derivative of the solution function (∂u(x,t)/∂x or ∂u(x,t)/∂t)${\displaystyle (\partial u(x,t)/\partial x\text{or}\partial u(x,t)/\partial t)}$ is kept constant. However, here we will only implement the derivative of the solution with respect to x${\displaystyle x}$.(36)∂u(x,t)

∂x=C$$\frac{\partial u(x,t)}{\partial x}=C$$where C${\displaystyle C}$ is a real number. At the boundary when x=0${\displaystyle x=0}$, we can write (37)u(𝛥x,t)-u(-𝛥x,t)

2𝛥x=C $$\frac{u(𝛥x,t)-u(-𝛥x,t)}{2𝛥x}=C$$therefore(38)u(-𝛥x,t)=u(𝛥x,t)-2𝛥xC$$u(-𝛥x,t)=u(𝛥x,t)-2𝛥xC$$At another boundary when x=L${\displaystyle x=L}$(39)u(L+𝛥x,t)-u(L-𝛥x,t)

2𝛥x=D $$\frac{u(L+𝛥x,t)-u(L-𝛥x,t)}{2𝛥x}=D$$where D${\displaystyle D}$ is a real number. Therefore(40)u(L+𝛥x,t)=u(L-𝛥x,t)+2𝛥xD$$u(L+𝛥x,t)=u(L-𝛥x,t)+2𝛥xD$$The implementation of Neumann boundary condition will be done by adding additional phantom nodes as shown in Drawing 6.

u(0,tk)u(L,tk)u(𝛥x,tk)u(L+𝛥x,tk)u(-𝛥x,tk)……x=0x=L

Drawing 6:Extra phantom nodes for Neumann boundary condition

2.1.3. MATLAB ImplementationMATLAB implementation of (35) is as follows.

1% alpha^2 * u_xx = u_tt2% Divide 0<=x<=L with Delta_x increments3% Divide 0<=t<=tF with Delta_t increments4 5function u_array = diffusion_1d(u_array, g_array, alpha, Delta_x, Delta_t)6 7r = alpha^2 * Delta_t / Delta_x^2;8 9j = 2 : length(u_array)-1;10u_array(j) = ...11 u_array(j) + r*(u_array(j+1) -2 * u_array(j) + u_array(j-1)) + ...12 g_array(j)*Delta_t;13end

The implementation above is already vectorized. In line 9 of the listing above, we can see that the diffusion is not applied to the first and last nodes. At those two end-nodes, we will apply the boundary conditions. Therefore, we also create two more functions for the initial and boundary conditions (see the listing below).

1% Initial condition2function u_array = apply_ic(u_array, f, Delta_x)3 4% Apply the initial condition U(x,0)=f(x)5j = 1:length(u_array);6u_array(j) = f((j-1)*Delta_x);7 8end9 10%==========================================================11 12% Boundary condition13function u_array = apply_bc(u_array, Delta_x, types, values)14 15% left-end16if lower(types(1)) == "dirichlet"17 u_array(1) = values(1);18elseif lower(types(1)) == "neumann"19 u_array(1) = u_array(3) - 2*Delta_x * values(1);20end21 22% right-end23if lower(types(2)) == "dirichlet"24 u_array(end) = values(2); 25elseif lower(types(2)) == "neumann"26 u_array(end) = u_array(end-2) + 2 * Delta_x * values(2);27end28end

In the Neumann boundary condition, we need to add additional phantom nodes. There are some examples provided in the next section to understand how to use the aforementioned three functions. 2.1.4. Example 1Given the following heat transfer equation:(41)a

1.14uxx-ut=0, 0⩽x⩽10

$$\begin{array}{c}1.14u_{xx}-u_t=0,0⩽x⩽10\end{array}$$with the following boundary and initial conditions:(42)a

u(0,t)=u(10,t)=0

u(x,0)=f(x)={a 60, 0⩽x⩽5 0, 5<x⩽10 |

$$\begin{array}{c}u(0,t)=u(10,t)=0\\ u(x,0)=f(x)=\{\begin{array}{ll}60,& 0⩽x⩽5\\ 0,& 5<x⩽10\end{array}\end{array}$$MATLAB implementation of (41) is as follows.

1L = 10;2alpha = sqrt(1.14);3 4tF = 20;5Nx = 10;6Nt = 1000;7 8u_array = zeros(1,Nx+1); % +1 because we start from 0 to Nx9g_array = zeros(1,Nx+1); % homogenous10u_array = apply_ic(u_array,@f, L/Nx);11 12figure13hold on14h = plot(0,0,'-*');15ylim([0 70])16 17x = 0 : Nx;18 19for k = 1 : Nt20 set(h, 'XData', x, 'YData', u_array)21 drawnow;22 23 u_array = apply_bc(u_array, L/Nx, ["Dirichlet", "Dirichlet"], [0, 0]); 24 u_array = diffusion_1d(u_array, g_array, alpha, L/Nx, tF/Nt);25end26xlabel('n-th segment')27ylabel('u(x,t)')28 29%% Intial Condition30function u=f(x)31L = 10;32u = zeros(1,length(x));33 34for k = 1 : length(x)35 if x(k) >= 0 && x(k) <=L/236 u(k) = 60;37 elseif x(k) > L/2 && x(k) <= L38 u(k) = 0;39 end40end41 42end

The result is shown in Drawing 7 below.

t=0→u(x,0)=f(x)={a

60,	0⩽x⩽5
0,	5<x⩽10

|t=4t=8t=12t=16u(0,t)=0u(10,t)=0

Drawing 7:Simulation result: Dirichlet-Dirichlet

2.1.5. Example 2Given the following heat transfer equation:(43)a

uxx-ut=0, 0⩽x⩽𝜋

$$\begin{array}{c}{u}_{xx}-u_t=0,0⩽x⩽𝜋\end{array}$$with the following boundary and initial conditions:(44)a

u(0,t)=20

ux(𝜋,t)=3

u(x,0)=0

$$\begin{array}{c}u(0,t)=20\\ u_x(𝜋,t)=3\\ u(x,0)=0\end{array}$$MATLAB implementation of (43) is as follows.

1L = pi;2alpha = 1;3 4tF = 10;5Nx = 10;6Nt = 1000;7 8u_array = zeros(1,Nx+1+1); % +1 -> start from 0 to Nx9 % another +1 -> add one phantom nodes10g_array = zeros(1,length(u_array));11u_array = apply_ic(u_array,@f, L/Nx);12 13figure14hold on15h = plot(0,0,'-*');16ylim([0 32])17 18x = 0 : Nx+1;19 20for k = 1 : Nt21 set(h, 'XData', x, 'YData', u_array)22 drawnow;23 24 u_array = apply_bc(u_array, L/Nx, ["Dirichlet", "Neumann"], [20, 3]);25 u_array = diffusion_1d(u_array, g_array, alpha, L/Nx, tF/Nt);26end27xlabel('n-th segment')28ylabel('u(x,t)')29 30%% Intial Condition31function u=f(x)32 u = zeros(1,length(x));33end

The result is shown in Drawing 8 below.

u(0,t)=20Phantom nodesux(𝜋,t)=3t=0→u(x,0)=f(x)=0t=4t=2t=6t=8

Drawing 8:Simulation result: Dirichlet-Neumann

2.1.6. Example 3Given the following heat equation:(45)a

uxx-ut=0, 0⩽x⩽1

$$\begin{array}{c}{u}_{xx}-u_t=0,0⩽x⩽1\end{array}$$with the following initial and boundary conditions:(46)a

ux(0,t)=0

u(1,t)=100

u(x,0)=f(x)={a 60, 0⩽x⩽1/2 0, 1/2<x⩽1 |

$$\begin{array}{c}{u}_x(0,t)=0\\ u(1,t)=100\\ u(x,0)=f(x)=\{\begin{array}{ll}60,& 0⩽x⩽1/2\\ 0,& 1/2<x⩽1\end{array}\end{array}$$MATLAB implementation of (45) is as follows.

1L = 1;2alpha = 1;3 4tF = 1;5Nx = 10;6Nt = 1000;7 8u_array = zeros(1,Nx+1+1); % +1 because we start from 0 to Nx9 % another +1 -> add one phantom nodes10g_array = zeros(1,length(u_array));11u_array = apply_ic(u_array,@f, L/Nx);12 13hf = figure;14hold on15h = plot(0,0,'-*');16ylim([0 105])17xlabel('n-th segment')18ylabel('u(x,t)')19 20x = 0 : Nx+1;21 22for k = 1 : Nt23 set(h, 'XData', x, 'YData', u_array)24 drawnow;25 26 u_array = apply_bc(u_array, L/Nx, ["Neumann", "Dirichlet"], [0, 100]); 27 u_array = diffusion_1d(u_array, g_array, alpha, L/Nx, tF/Nt);28end29 30%% Initial condition31function u=f(x)32L = 1;33u = zeros(1,length(x));34 35for k = 1 : length(x)36 if x(k) >= 0 && x(k) <=L/237 u(k) = 60;38 elseif x(k) > L/2 && x(k) <= L39 u(k) = 0;40 end41end42 43end

The result is shown in Drawing 9 below.

t=0→u(x,0)=f(x)={a

60,	0⩽x⩽1
0,	1<x⩽10

|t=0.4t=0.2t=0.6t=0.8u(1,t)=100ux(0,t)=0Phantom nodes

Drawing 9:Simulation result: Neumann-Dirichlet

2.2. Implicit Method The main idea of an implicit method is to calculate uxx${\displaystyle u_{xx}}$ not only at the current time, but also at the time one step ahead. Let us consider:(47)uxx(x,t)=(1-𝜃)uxx(x,t)+𝜃uxx(x,t+𝛥t), 0⩽𝜃⩽1$$u_{xx}(x,t)=(1-𝜃)u_{xx}(x,t)+𝜃u_{xx}(x,t+𝛥t),0⩽𝜃⩽1$$where 𝜃${\displaystyle 𝜃}$ acts as a weighting factor. When 𝜃=0${\displaystyle 𝜃=0}$, (47) becomes similar to (35). What we currently have in this case is an explicit method. On the other hand, when 𝜃=1${\displaystyle 𝜃=1}$, in this case, what we have is an implicit method. However, when 𝜃=0.5${\displaystyle 𝜃=0.5}$, the method is called Crank-Nicolson method. In this section, let us focus on mathematical derivations of the Crank-Nicolson method.Expanding (47) bu first setting 𝜃=1

2${\displaystyle 𝜃=\frac{1}{2}}$ gives us:

	uxx	=1 2ux(x+𝛥x,t)-ux(x,t) 𝛥x+1 2ux(x+𝛥x,t+𝛥t)-ux(x,t+𝛥t) 𝛥x
		=u(x+𝛥x,t)-u(x,t) 𝛥x-u(x,t)-u(x-𝛥x,t) 𝛥x 2𝛥x+u(x+𝛥x,t+𝛥t)-u(x,t+𝛥t) 𝛥x-u(x,t+𝛥t)-u(x-𝛥x,t+𝛥t) 𝛥x 2𝛥x
		a=u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t) 2𝛥x2 +u(x+𝛥x,t+𝛥t)-2u(x,t+𝛥t)+u(x-𝛥x,t+𝛥t) 2𝛥x2

$$\begin{array}{rl}{u}_{xx}& ={\displaystyle \frac{1}{2}}{\displaystyle \frac{u_x(x+𝛥x,t)-u_x(x,t)}{𝛥x}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{u_x(x+𝛥x,t+𝛥t)-u_x(x,t+𝛥t)}{𝛥x}}\\ & ={\displaystyle \frac{\frac{u(x+𝛥x,t)-u(x,t)}{𝛥x}-\frac{u(x,t)-u(x-𝛥x,t)}{𝛥x}}{2𝛥x}}+{\displaystyle \frac{\frac{u(x+𝛥x,t+𝛥t)-u(x,t+𝛥t)}{𝛥x}-\frac{u(x,t+𝛥t)-u(x-𝛥x,t+𝛥t)}{𝛥x}}{2𝛥x}}\\ & \begin{array}{c}={\displaystyle \frac{u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)}{2𝛥x^2}}\\ +{\displaystyle \frac{u(x+𝛥x,t+𝛥t)-2u(x,t+𝛥t)+u(x-𝛥x,t+𝛥t)}{2𝛥x^2}}\end{array}\end{array}$$which can be shortened as(48)uxx=Uj+1,k-2Uj,k+Uj-1,k+Uj+1,k+1-2Uj,k+1+Uj-1,k+1

2𝛥x2$$u_{xx}=\frac{U_{j+1,k}-2U_{j,k}+U_{j-1,k}+U_{j+1,k+1}-2U_{j,k+1}+U_{j-1,k+1}}{2𝛥x^2}$$Substituting (48) back to: ut=𝛼²uxx${\displaystyle u_t=𝛼²u_{xx}}$ gives us:(49)a

Uj,k+1-Uj,k 𝛥t⏠⏣⏣⏣⏡⏣⏣⏣⏢ut	=𝛼2Uj+1,k-2Uj,k+Uj-1,k+Uj+1,k+1-2Uj,k+1+Uj-1,k+1 2𝛥x2⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢uxx
Uj,k+1-Uj,k	=1 2𝛼2𝛥t 𝛥x2(Uj+1,k-2Uj,k+Uj-1,k+Uj+1,k+1-2Uj,k+1+Uj-1,k+1)

$$\begin{array}{rl}\underset{u_t}{\underbrace{\frac{U_{j,k+1}-U_{j,k}}{𝛥t}}}& ={𝛼}^2\underset{u_{xx}}{\underbrace{\frac{U_{j+1,k}-2U_{j,k}+U_{j-1,k}+U_{j+1,k+1}-2U_{j,k+1}+U_{j-1,k+1}}{2𝛥x^2}}}\\ U_{j,k+1}-U_{j,k}& =\frac{1}{2}\frac{{𝛼}^2𝛥t}{𝛥x^2}(U_{j+1,k}-2U_{j,k}+U_{j-1,k}+U_{j+1,k+1}-2U_{j,k+1}+U_{j-1,k+1})\end{array}$$ Since r=𝛼2𝛥t

𝛥x2${\displaystyle r=\frac{{𝛼}^2𝛥t}{𝛥x^2}}$, (49) can then be simplified as:(50)a

Uj,k+1-Uj,k	=r 2(Uj+1,k-2Uj,k+Uj-1,k+Uj+1,k+1-2Uj,k+1+Uj-1,k+1)
Uj,k+1+rUj,k+1	=r 2Uj+1,k+(1-r)Uj,k+r 2Uj-1,k+r 2Uj+1,k+1+r 2Uj-1,k+1
2(1+r)Uj,k+1	=rUj+1,k+2(1-r)Uj,k+rUj-1,k+rUj+1,k+1+rUj-1,k+1
-rUj-1,k+1+2(1+r)Uj,k+1-rUj+1,k+1	=rUj-1,k+2(1-r)Uj,k+rUj+1,k

$$\begin{array}{rl}{U}_{j,k+1}-U_{j,k}& =\frac{r}{2}(U_{j+1,k}-2U_{j,k}+U_{j-1,k}+U_{j+1,k+1}-2U_{j,k+1}+U_{j-1,k+1})\\ U_{j,k+1}+r{U}_{j,k+1}& =\frac{r}{2}{U}_{j+1,k}+(1-r)U_{j,k}+\frac{r}{2}{U}_{j-1,k}+\frac{r}{2}{U}_{j+1,k+1}+\frac{r}{2}{U}_{j-1,k+1}\\ 2(1+r)U_{j,k+1}& =r{U}_{j+1,k}+2(1-r)U_{j,k}+r{U}_{j-1,k}+r{U}_{j+1,k+1}+r{U}_{j-1,k+1}\\ -r{U}_{j-1,k+1}+2(1+r)U_{j,k+1}-r{U}_{j+1,k+1}& =r{U}_{j-1,k}+2(1-r)U_{j,k}+r{U}_{j+1,k}\end{array}$$ Thus, the well-known Crank-Nicolson scheme can be expressed as follows.(51)

-rUj-1,k+1+2(1+r)Uj,k+1-rUj+1,k+1=rUj-1,k+2(1-r)Uj,k+rUj+1,k

$$\boxed{-r{U}_{j-1,k+1}+2(1+r)U_{j,k+1}-r{U}_{j+1,k+1}=r{U}_{j-1,k}+2(1-r)U_{j,k}+r{U}_{j+1,k}}$$or in matrix form:[a

-r	2(1+r)	-r		⋯	0
0	-r	2(1+r)	-r		⋮
		-r	2(1+r)	-r
			⋱
⋮				-r	2(1+r)
0	⋯			0	-r

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

p((k+1)𝛥t)

U1,k+1

U2,k+1

⋮

UN-2,k+1

q((k+1)𝛥t)

]⏠⏣⏣⏣⏣⏡⏣⏣⏣⏣⏢U=(52)[a

r	2(1-r)	r		⋯	0
0	r	2(1-r)	r		⋮
		r	2(1-r)	r
			⋱
⋮				r	2(1-r)
0	⋯			0	r

][a

p(k𝛥t)

U1,k

U2,k

⋮

UN-2,k

q(k𝛥t)

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢B$$\begin{array}{c}\underset{A}{\underbrace{\left[\begin{array}{cccccc}-r& 2(1+r)& -r& & \cdots & 0\\ 0& -r& 2(1+r)& -r& & ⋮\\ & & -r& 2(1+r)& -r& \\ & & & \ddots & & \\ ⋮& & & & -r& 2(1+r)\\ 0& \cdots & & & 0& -r\end{array}\right]}}\underset{U}{\underbrace{\left[\begin{array}{c}p((k+1)𝛥t)\\ U_{1,k+1}\\ U_{2,k+1}\\ ⋮\\ U_{N-2,k+1}\\ q((k+1)𝛥t)\end{array}\right]}}=\\ \underset{B}{\underbrace{\left[\begin{array}{cccccc}r& 2(1-r)& r& & \cdots & 0\\ 0& r& 2(1-r)& r& & ⋮\\ & & r& 2(1-r)& r& \\ & & & \ddots & & \\ ⋮& & & & r& 2(1-r)\\ 0& \cdots & & & 0& r\end{array}\right]\left[\begin{array}{c}p(k𝛥t)\\ U_{1,k}\\ U_{2,k}\\ ⋮\\ U_{N-2,k}\\ q(k𝛥t)\end{array}\right]}}\end{array}$$Red means the values are known. In (52), we can notice that some known values reside in U${\displaystyle U}$. We must remove these known values from U${\displaystyle U}$ and put them to the right side. The goal is to have a perfect AU=B${\displaystyle AU=B}$ system where all the unknowns are all in U${\displaystyle U}$. [a

2(1+r)	-r		⋯	0
-r	2(1+r)	-r		⋮
		⋱
		-r	2(1+r)	-r
⋯			-r	2(1+r)

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

U1,k+1

U2,k+1

⋮

UN-2,k+1

UN-1,k+1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢Unew=(53)[a

2(1-r)	r		⋯	⋯
r	2(1-r)	r
	⋱
		r	2(1-r)	r
⋯			r	2(1-r)

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢b[a

U1,k

U2,k

⋮

UN-2,k

UN-1,k

]⏠⏣⏣⏣⏡⏣⏣⏣⏢Unow+[a

rp((k+1)𝛥t)+rp(k𝛥t)

0

⋮

0

rq((k+1)𝛥t)+rq(k𝛥t)

]⏠⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏢c⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢B=(b·Unow)+c$$\begin{array}{c}\underset{A}{\underbrace{\left[\begin{array}{ccccc}2(1+r)& -r& & \cdots & 0\\ -r& 2(1+r)& -r& & ⋮\\ & & \ddots & & \\ & & -r& 2(1+r)& -r\\ \cdots & & & -r& 2(1+r)\end{array}\right]}}\underset{U_{new}}{\underbrace{\left[\begin{array}{c}{U}_{1,k+1}\\ U_{2,k+1}\\ ⋮\\ U_{N-2,k+1}\\ U_{N-1,k+1}\end{array}\right]}}=\\ \underset{B={\displaystyle (}b·U_{now}{\displaystyle )}+c}{\underbrace{\underset{b}{\underbrace{\left[\begin{array}{ccccc}2(1-r)& r& & \cdots & \cdots \\ r& 2(1-r)& r& & \\ & \ddots & & & \\ & & r& 2(1-r)& r\\ \cdots & & & r& 2(1-r)\end{array}\right]}}\underset{U_{now}}{\underbrace{\left[\begin{array}{c}{U}_{1,k}\\ U_{2,k}\\ ⋮\\ U_{N-2,k}\\ U_{N-1,k}\end{array}\right]}}+\underset{c}{\underbrace{\left[\begin{array}{c}rp((k+1)𝛥t)+rp(k𝛥t)\\ 0\\ ⋮\\ 0\\ rq((k+1)𝛥t)+rq(k𝛥t)\end{array}\right]}}}}\end{array}$$ Now, we can find Unew=A-1B${\displaystyle U_{new}=A^{-1}B}$. In real implementation, this might not be plausible and we have to perform an iteration-based algorithm to find the optimal Unew${\displaystyle U_{new}}$. In Maple, we can use the Unew=LinearSolve(A,B)${\displaystyle U_{new}=\text{LinearSolve(A,B)}}$ command. In MATLAB, we can use the popular Unew=A\B${\displaystyle U_{new}=A\backslash B}$ command. Other alternatives are using the Jacobi method and the Gauss-Siedel method as explained the Greenberg's book.3. Two-Dimensional Poisson Equation73.1. Dirichlet Boundary

u(0,y)=p(y)u(x,b)=s(x)u(a,y)=r(y)bau(x,0)=q(x)xyu(x,y)DyxD

Drawing 10:Heat distribution on a two-dimensional plane with Dirichlet boundary conditions

Heat transfer on a two-dimensional plane (see Drawing 10) is described as follows:(54)uxx+uyy=f(x,y)0≤x≤a,0≤y≤b$$\begin{array}{c}{u}_{xx}+u_{yy}=f(x,y)\\ 0\le x\le a,0\le y\le b\end{array}$$The PDE in (54) is known as the Poisson equation. Its homogenous form is known as the Laplace equation. Finite difference method for the system above is as follows. First, uxx${\displaystyle u_{xx}}$ and uyy${\displaystyle u_{yy}}$ can be approximated as:(55)a

uxx=𝜕2u 𝜕x2≈u(x-𝛥x,y)-2u(x,y)+u(x+𝛥x,y) 𝛥x2

uyy=𝜕2u 𝜕y2≈u(x,y-𝛥y)-2u(x,y)+u(x,y+𝛥y) 𝛥y2

$$\begin{array}{l}{u}_{xx}=\frac{{𝜕}^{2}u}{𝜕x^2}\approx \frac{u(x-𝛥x,y)-2u(x,y)+u(x+𝛥x,y)}{𝛥x^2}\\ u_{yy}=\frac{{𝜕}^{2}u}{𝜕y^2}\approx \frac{u(x,y-𝛥y)-2u(x,y)+u(x,y+𝛥y)}{𝛥y^2}\end{array}$$Hence:(56)a

uxx+uyy	=f(x,y)
u(x-𝛥x,y)-2u(x,y)+u(x+𝛥x,y) 𝛥x2+u(x,y-𝛥y)-2u(x,y)+u(x,y+𝛥y) 𝛥y2	=f(x,y)

$$\begin{array}{rl}{u}_{xx}+u_{yy}& =f(x,y)\\ \frac{u(x-𝛥x,y)-2u(x,y)+u(x+𝛥x,y)}{𝛥x^2}+\frac{u(x,y-𝛥y)-2u(x,y)+u(x,y+𝛥y)}{𝛥y^2}& =f(x,y)\end{array}$$For simplification purpose, let us take 𝛥x=𝛥x=h${\displaystyle 𝛥x=𝛥x=h}$, therefore:(57)a

uxx+uyy	=f(x,y)
u(x-h,y)-2u(x,y)+u(x+h,y) h2+u(x,y-h)-2u(x,y)+u(x,y+h) h2	=f(x,y)

$$\begin{array}{rl}{u}_{xx}+u_{yy}& =f(x,y)\\ \frac{u(x-h,y)-2u(x,y)+u(x+h,y)}{h^2}+\frac{u(x,y-h)-2u(x,y)+u(x,y+h)}{h^2}& =f(x,y)\end{array}$$Let us take Uj,k=u(xj,yk)${\displaystyle U_{j,k}=u(x_j,y_k)}$, Fj,k=f(xj,yk)${\displaystyle F_{j,k}=f(x_j,y_k)}$, xj=jh${\displaystyle x_j=jh}$, yk=kh${\displaystyle y_k=kh}$, Uj,k+1=u(xj,yk+h)${\displaystyle U_{j,k+1}=u(x_j,y_k+h)}$ and Uj-1,k=u(xj-h,yk)${\displaystyle U_{j-1,k}=u(x_j-h,y_k)}$, hence:(58)a

u(x-h,y)-2u(x,y)+u(x+h,y) h2+u(x,y-h)-2u(x,y)+u(x,y+h) h2	=f(x,y)
Uj-1,k-2Uj.k+Uj+1,k+Uj,k-1-2Uj.k+Uj,k+1	=h2Fj,k
Uj-1,k+Uj+1,k+Uj,k-1+Uj,k+1-4Uj,k	=h2Fj,k

$$\begin{array}{rl}\frac{u(x-h,y)-2u(x,y)+u(x+h,y)}{h^2}+\frac{u(x,y-h)-2u(x,y)+u(x,y+h)}{h^2}& =f(x,y)\\ U_{j-1,k}-2U_{j.k}+U_{j+1,k}+U_{j,k-1}-2U_{j.k}+U_{j,k+1}& =h^2{F}_{j,k}\\ U_{j-1,k}+U_{j+1,k}+U_{j,k-1}+U_{j,k+1}-4U_{j,k}& =h^2{F}_{j,k}\end{array}$$Therefore, the resulting discrete formula is as follows:(59)

a

Uj-1,k+Uj+1,k+Uj,k-1+Uj,k+1-4Uj,k

=h2Fj,k

$$\boxed{\begin{array}{rl}{U}_{j-1,k}+U_{j+1,k}+U_{j,k-1}+U_{j,k+1}-4U_{j,k}& =h^2{F}_{j,k}\end{array}}$$In order to solve (59), we need to reformulate it into a matrix form. In order to do so, we first define the following girds:

Left edgeRight edgeTop edgeBottom edgek=0k=1… …j=0j=1⋮⋮k=Nxj=Nyhhab

Drawing 11:Grid discretization of the 2D plate.

Next, based on the Drawing 11, we can then get the following equation:[a

-4	1	0	1				⋯	0
1	-4	1	0	1				⋮
	⋱	⋱	⋱	⋱	⋱
0	1	-4	0	0	1
1	0	1	-4	1	0	1
	1	0	1	-4	1	0	1
		⋱	⋱	⋱	⋱	⋱	⋱	⋱
		1	0	1	-4	1	0	1
			⋱	⋱	⋱	⋱	⋱
			1	0	0	-4	1	0
				1	0	1	-4	1
⋮					⋱	⋱	⋱	⋱
0	…				1	0	1	-4

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

U1,1

U2,1

⋮

UNy-1,1

U1,2

U2,2

⋮

UNy-1,2

⋮

U1,Nx-1

U2,Nx-1

⋮

UNy-1,Nx-1

]⏠⏣⏣⏣⏣⏡⏣⏣⏣⏣⏢U=[a

-(U1,0+U0,1)+h2F1,1

-U2,0+h2F2,1

⋮

-(UNy-1,0+UNy,1)+h2FNy,1

-U0,2+h2F1,2

h2F2,2

⋮

-UNy,2+h2FNy-1,2

⋮

-(U1,Nx+U0,Nx-1)+h2F1,Nx-1

-U2,Nxh2F2,Nx-1

⋮

-(UNy-1,Nx+UNy,Nx-1)+h2FNy-1,Nx-1

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢c (60)$$\begin{array}{c}\underset{A}{\underbrace{\left[\begin{array}{rrrrrrrrr}-4& 1& 0& 1& & & & \cdots & 0\\ 1& -4& 1& 0& 1& & & & ⋮\\ & \ddots & \ddots & \ddots & \ddots & \ddots & & & \\ 0& 1& -4& 0& 0& 1& & & \\ 1& 0& 1& -4& 1& 0& 1& & \\ & 1& 0& 1& -4& 1& 0& 1& \\ & & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots \\ & & 1& 0& 1& -4& 1& 0& 1\\ & & & \ddots & \ddots & \ddots & \ddots & \ddots & \\ & & & 1& 0& 0& -4& 1& 0\\ & & & & 1& 0& 1& -4& 1\\ ⋮& & & & & \ddots & \ddots & \ddots & \ddots \\ 0& \dots & & & & 1& 0& 1& -4\end{array}\right]}}\underset{\tilde{U}}{\underbrace{\left[\begin{array}{c}{U}_{1,1}\\ U_{2,1}\\ ⋮\\ U_{N_y-1,1}\\ U_{1,2}\\ U_{2,2}\\ ⋮\\ U_{N_y-1,2}\\ ⋮\\ U_{1,N_x-1}\\ U_{2,N_x-1}\\ ⋮\\ U_{N_y-1,N_x-1}\end{array}\right]}}=\underset{c}{\underbrace{\left[\begin{array}{c}-(U_{1,0}+U_{0,1})+h^2{F}_{1,1}\\ -U_{2,0}+h^2{F}_{2,1}\\ ⋮\\ -(U_{N_y-1,0}+U_{N_y,1})+h^2{F}_{N_y,1}\\ -U_{0,2}+h^2{F}_{1,2}\\ h^2{F}_{2,2}\\ ⋮\\ -U_{N_y,2}+h^2{F}_{N_y-1,2}\\ ⋮\\ -(U_{1,N_x}+U_{0,N_x-1})+h^2{F}_{1,N_x-1}\\ -U_{2,N_x}{h}^2{F}_{2,N_x-1}\\ ⋮\\ -(U_{N_y-1,N_x}+U_{N_y,N_x-1})+h^2{F}_{N_y-1,N_x-1}\end{array}\right]}}\\ \end{array}$$Once we have (60), we can solve it to acquire U${\displaystyle \tilde{U}}$. In MATLAB, this can be easily done by using: U=A\c${\displaystyle \tilde{U}=A\backslash c}$. Let us now take a more simple example. Given the following diffusion PDE:uxx+uyy=f(x,y)$$u_{xx}+u_{yy}=f(x,y)$$where:0⩽x⩽1, 0⩽y⩽1${\displaystyle 0⩽x⩽1,0⩽y⩽1}$Nx=Ny=4${\displaystyle N_x=N_y=4}$ p(y)=0${\displaystyle p(y)=0}$q(x)=0${\displaystyle q(x)=0}$r(y)=100sin(𝜋y)${\displaystyle r(y)=100\sin (𝜋y)}$s(x)=0${\displaystyle s(x)=0}$f(x,y)=0${\displaystyle f(x,y)=0}$h=0.25${\displaystyle h=0.25}$Find u(x,y)${\displaystyle u(x,y)}$ at its steady state condition. First, we create the matrix U${\displaystyle U}$ as follows:

U=a

a 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 70.71 100 70.71 0

a

j=0

j=1

j=2

j=3

j=4

k=0k=1k=2k=3k=4Red colored numbers represents values at the boundaries. These values are calculated from p(y)${\displaystyle p(y)}$, q(x)${\displaystyle q(x)}$, r(y)${\displaystyle r(y)}$, and s(x)${\displaystyle s(x)}$. Next, we generate an equation similar (60) and calculate the U${\displaystyle \tilde{U}}$.a

a -4 1 0 1 0 0 0 0 0 1 -4 1 0 1 0 0 0 0 0 1 -4 0 0 1 0 0 0 1 0 0 -4 1 0 1 0 0 0 1 0 1 -4 1 0 1 0 0 0 1 0 1 -4 0 0 1 0 0 0 1 0 0 -4 1 0 0 0 0 0 1 0 1 -4 1 0 0 0 0 0 1 0 1 -4

a

a U1,1 U2,1 U3,1 U1,2 U2,2 U3,2 U1,3 U2,3 U3,3

=a

a -U1,0-U0,1 -U2,0 -U3,0-U4,1 -U0,2 0 -U4,2 -U1,4-U0,3 -U2,4 -U3,4-U4,3

$$\left[\begin{array}{c}\begin{array}{rrrrrrrrr}-4& 1& 0& 1& 0& 0& 0& 0& 0\\ 1& -4& 1& 0& 1& 0& 0& 0& 0\\ 0& 1& -4& 0& 0& 1& 0& 0& 0\\ 1& 0& 0& -4& 1& 0& 1& 0& 0\\ 0& 1& 0& 1& -4& 1& 0& 1& 0\\ 0& 0& 1& 0& 1& -4& 0& 0& 1\\ 0& 0& 0& 1& 0& 0& -4& 1& 0\\ 0& 0& 0& 0& 1& 0& 1& -4& 1\\ 0& 0& 0& 0& 0& 1& 0& 1& -4\end{array}\end{array}\right]\left[\begin{array}{c}\begin{array}{c}{U}_{1,1}\\ U_{2,1}\\ U_{3,1}\\ U_{1,2}\\ U_{2,2}\\ U_{3,2}\\ U_{1,3}\\ U_{2,3}\\ U_{3,3}\end{array}\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}-U_{1,0}-U_{0,1}\\ -U_{2,0}\\ -U_{3,0}-U_{4,1}\\ -U_{0,2}\\ 0\\ -U_{4,2}\\ -U_{1,4}-U_{0,3}\\ -U_{2,4}\\ -U_{3,4}-U_{4,3}\end{array}\end{array}\right]$$ 3.2. Neumann Boundary 4. Two-Dimensional Unsteady Advection-DiffusionLet us consider the flowing PDE:(61)𝜕T

𝜕t+u𝜕T

𝜕x+v𝜕T

𝜕y⏠⏣⏣⏣⏡⏣⏣⏣⏢2D Advection=𝛼(𝜕2T

𝜕x2+𝜕2T

𝜕y2⏠⏣⏣⏣⏡⏣⏣⏣⏢2D Diffusion)+f$$\frac{𝜕T}{𝜕t}+\underset{\text{2D Advection}}{\underbrace{u\frac{𝜕T}{𝜕x}+v\frac{𝜕T}{𝜕y}}}=𝛼\left(\underset{\text{2D Diffusion}}{\underbrace{\frac{{𝜕}^{2}T}{𝜕x^2}+\frac{{𝜕}^{2}T}{𝜕y^2}}}\right)+f$$(61) is the general transient heat equation in 2D space. It is composed of advection process and diffusion process. 𝛼${\displaystyle 𝛼}$ is defines as:(62)𝛼=K

𝜌cp$$𝛼=\frac{\mathcal{K}}{𝜌c_p}$$where:cp${\displaystyle c_p}$ is the specific heat at constant (J(kg K)${\displaystyle (\text{J}(\text{kg}\text{K})}$𝜌${\displaystyle 𝜌}$ is the fluid density (kg/m3)${\displaystyle (\text{kg}/{\text{m}}^3)}$K${\displaystyle \mathcal{K}}$ is the thermal conductivity (W(m K))${\displaystyle \left(\text{W}\right(\text{m}\text{K}\left)\right)}$f${\displaystyle f}$ is the additional forcing component In order to discretize (61), we will use the simplest methods, which are the upwind/downwind method for the advection and the explicit FTCS method for the diffusion. 4.1. UPWIND-UPWIND MethodThe upwind-upwind method is used when both velocities are positive or zero. We start by discretizing (61) as follows: T(x,y,t+𝛥t)-T(x,y,t)

𝛥t+u(x,y,t)T(x,y,t)-T(x-𝛥x,y,t)

𝛥x+ …(63)v(x,y,t)T(x,y,t)-T(x,y-𝛥y,t)

𝛥y=𝛼(T(x+𝛥x,y,t)-2T(x,y,z)+T(x-𝛥x,y,t)

𝛥x2)+ …𝛼(T(x,y+𝛥y,t)-2T(x,y,z)+T(x,y-𝛥y,t)

𝛥y2)+f(x,y,t)$$\begin{array}{c}\frac{T(x,y,t+𝛥t)-T(x,y,t)}{𝛥t}+u(x,y,t)\frac{T(x,y,t)-T(x-𝛥x,y,t)}{𝛥x}+\dots \\ v(x,y,t)\frac{T(x,y,t)-T(x,y-𝛥y,t)}{𝛥y}=𝛼\left(\frac{T(x+𝛥x,y,t)-2T(x,y,z)+T(x-𝛥x,y,t)}{𝛥x^2}\right)+\dots \\ 𝛼\left(\frac{T(x,y+𝛥y,t)-2T(x,y,z)+T(x,y-𝛥y,t)}{𝛥y^2}\right)+f(x,y,t)\end{array}$$Rearranging (63) give us:(64)a

T(x,y,t+𝛥t)=	𝛼𝛥t 𝛥x2(T(x+𝛥x,y,t)-2T(x,y,z)+T(x-𝛥x,y,t))+
	𝛼𝛥t 𝛥y2(T(x,y+𝛥y,t)-2T(x,y,z)+T(x,y-𝛥y,t))-
	u(x,y,t)𝛥t 𝛥x(T(x,y,t)-T(x-𝛥x,y,t))-
	v(x,y,t)𝛥t 𝛥y(T(x,y,t)-T(x,y-𝛥y,t)) + T(x,y,t)+f(x,y,t)𝛥t

$$\begin{array}{rl}T(x,y,t+𝛥t)=& \frac{𝛼𝛥t}{𝛥x^2}(T(x+𝛥x,y,t)-2T(x,y,z)+T(x-𝛥x,y,t))+\\ & \frac{𝛼𝛥t}{𝛥y^2}(T(x,y+𝛥y,t)-2T(x,y,z)+T(x,y-𝛥y,t))-\\ & \frac{u(x,y,t)𝛥t}{𝛥x}(T(x,y,t)-T(x-𝛥x,y,t))-\\ & \frac{v(x,y,t)𝛥t}{𝛥y}(T(x,y,t)-T(x,y-𝛥y,t))+T(x,y,t)+f(x,y,t)𝛥t\end{array}$$ Taking: Ti,j,k=T(xi,yj,tk)${\displaystyle T_{i,j,k}=T(x_i,y_j,t_k)}$, xj=j𝛥x${\displaystyle x_j=j𝛥x}$, yk=k𝛥y${\displaystyle y_k=k𝛥y}$, and tk=k𝛥t${\displaystyle t_k=k𝛥t}$, we can rewrite (64) as follows:(65)

aTi,j,k+1=𝛼𝛥t((Ti+1,j,k-2Ti,j,k+Ti-1,j,k

𝛥x2)+(Ti,j+1,k-2Ti,j,k+Ti,j-1,k

𝛥y2))- 𝛥t(ui,j,k

𝛥x(Ti,j,k-Ti-1,j,k)+vi,j,k

𝛥y(Ti,j,k-Ti,j-1,k))+ Ti,j,k + 𝛥t fi,j,k

$$\boxed{\begin{array}{c}{T}_{i,j,k+1}=𝛼𝛥t\left(\right(\frac{T_{i+1,j,k}-2T_{i,j,k}+T_{i-1,j,k}}{𝛥x^2})+(\frac{T_{i,j+1,k}-2T_{i,j,k}+T_{i,j-1,k}}{𝛥y^2}\left)\right)-\\ 𝛥t(\frac{u_{i,j,k}}{𝛥x}(T_{i,j,k}-T_{i-1,j,k})+\frac{v_{i,j,k}}{𝛥y}(T_{i,j,k}-T_{i,j-1,k}))+T_{i,j,k}+\\ 𝛥t{f}_{i,j,k}\end{array}}$$ 4.2. DOWNWIND-DOWNWIND MethodThe downwind-downwind method is used when both velocities are negative or zero. The discretization process is similar to the upwind method. The only difference is in the advection part. The discretization is done as follows: T(x,y,t+𝛥t)-T(x,y,t)

𝛥t+u(x,y,t)T(x+𝛥x,y,t)-T(x,y,t)

𝛥x+ …(66)v(x,y,t)T(x,y+𝛥y,t)-T(x,y,t)

𝛥y=𝛼(T(x+𝛥x,y,t)-2T(x,y,z)+T(x-𝛥x,y,t)

𝛥x2)+ …𝛼(T(x,y+𝛥y,t)-2T(x,y,z)+T(x,y-𝛥y,t)

𝛥y2)+f(x,y,t)$$\begin{array}{c}\frac{T(x,y,t+𝛥t)-T(x,y,t)}{𝛥t}+u(x,y,t)\frac{T(x+𝛥x,y,t)-T(x,y,t)}{𝛥x}+\dots \\ v(x,y,t)\frac{T(x,y+𝛥y,t)-T(x,y,t)}{𝛥y}=𝛼\left(\frac{T(x+𝛥x,y,t)-2T(x,y,z)+T(x-𝛥x,y,t)}{𝛥x^2}\right)+\dots \\ 𝛼\left(\frac{T(x,y+𝛥y,t)-2T(x,y,z)+T(x,y-𝛥y,t)}{𝛥y^2}\right)+f(x,y,t)\end{array}$$Rearranging (66) give us:(67)a

T(x,y,t+𝛥t)=	𝛼𝛥t 𝛥x2(T(x+𝛥x,y,t)-2T(x,y,z)+T(x-𝛥x,y,t))+
	𝛼𝛥t 𝛥y2(T(x,y+𝛥y,t)-2T(x,y,z)+T(x,y-𝛥y,t))-
	u(x,y,t)𝛥t 𝛥x(T(x+𝛥x,y,t)-T(x,y,t))-
	v(x,y,t)𝛥t 𝛥y(T(x,y+𝛥y,t)-T(x,y,t)) + T(x,y,t)+f(x,y,t)𝛥t

$$\begin{array}{rl}T(x,y,t+𝛥t)=& \frac{𝛼𝛥t}{𝛥x^2}(T(x+𝛥x,y,t)-2T(x,y,z)+T(x-𝛥x,y,t))+\\ & \frac{𝛼𝛥t}{𝛥y^2}(T(x,y+𝛥y,t)-2T(x,y,z)+T(x,y-𝛥y,t))-\\ & \frac{u(x,y,t)𝛥t}{𝛥x}(T(x+𝛥x,y,t)-T(x,y,t))-\\ & \frac{v(x,y,t)𝛥t}{𝛥y}(T(x,y+𝛥y,t)-T(x,y,t))+T(x,y,t)+f(x,y,t)𝛥t\end{array}$$ Taking: Ti,j,k=T(xi,yj,tk)${\displaystyle T_{i,j,k}=T(x_i,y_j,t_k)}$, xj=j𝛥x${\displaystyle x_j=j𝛥x}$, yk=k𝛥y${\displaystyle y_k=k𝛥y}$, and tk=k𝛥t${\displaystyle t_k=k𝛥t}$, we can rewrite (67) as follows:(68)

aTi,j,k+1=𝛼𝛥t((Ti+1,j,k-2Ti,j,k+Ti-1,j,k

𝛥x2)+(Ti,j+1,k-2Ti,j,k+Ti,j-1,k

𝛥y2))- 𝛥t(ui,j,k

𝛥x(Ti+1,j,k-Ti,j,k)+vi,j,k

𝛥y(Ti,j+1,k-Ti,j,k))+ Ti,j,k+ 𝛥t fi,j,k

$$\boxed{\begin{array}{c}{T}_{i,j,k+1}=𝛼𝛥t\left(\right(\frac{T_{i+1,j,k}-2T_{i,j,k}+T_{i-1,j,k}}{𝛥x^2})+(\frac{T_{i,j+1,k}-2T_{i,j,k}+T_{i,j-1,k}}{𝛥y^2}\left)\right)-\\ 𝛥t(\frac{u_{i,j,k}}{𝛥x}(T_{i+1,j,k}-T_{i,j,k})+\frac{v_{i,j,k}}{𝛥y}(T_{i,j+1,k}-T_{i,j,k}))+T_{i,j,k}+\\ 𝛥t{f}_{i,j,k}\end{array}}$$4.3. DOWNWIND-UPWIND MethodThe downwind-upwind method is used when the x${\displaystyle x}$-velocity is negative and the y${\displaystyle y}$-velocity is positive.(69)

aTi,j,k+1=𝛼𝛥t((Ti+1,j,k-2Ti,j,k+Ti-1,j,k

𝛥x2)+(Ti,j+1,k-2Ti,j,k+Ti,j-1,k

𝛥y2))- 𝛥t(ui,j,k

𝛥x(Ti+1,j,k-Ti,j,k)+vi,j,k

𝛥y(Ti,j,k-Ti,j-1,k))+ Ti,j,k+ 𝛥t fi,j,k

$$\boxed{\begin{array}{c}{T}_{i,j,k+1}=𝛼𝛥t\left(\right(\frac{T_{i+1,j,k}-2T_{i,j,k}+T_{i-1,j,k}}{𝛥x^2})+(\frac{T_{i,j+1,k}-2T_{i,j,k}+T_{i,j-1,k}}{𝛥y^2}\left)\right)-\\ 𝛥t(\frac{u_{i,j,k}}{𝛥x}(T_{i+1,j,k}-T_{i,j,k})+\frac{v_{i,j,k}}{𝛥y}(T_{i,j,k}-T_{i,j-1,k}))+T_{i,j,k}+\\ 𝛥t{f}_{i,j,k}\end{array}}$$4.4. UPWIND-DOWNWIND MethodThe downwind-upwind method is used when the x${\displaystyle x}$-velocity is positive and the y${\displaystyle y}$-velocity is negative.(70)

aTi,j,k+1=𝛼𝛥t((Ti+1,j,k-2Ti,j,k+Ti-1,j,k

𝛥x2)+(Ti,j+1,k-2Ti,j,k+Ti,j-1,k

𝛥y2))- 𝛥t(ui,j,k

𝛥x(Ti,j,k-Ti-1,j,k)+vi,j,k

𝛥y(Ti,j+1,k-Ti,j,k))+ Ti,j,k+ 𝛥t fi,j,k

$$\boxed{\begin{array}{c}{T}_{i,j,k+1}=𝛼𝛥t\left(\right(\frac{T_{i+1,j,k}-2T_{i,j,k}+T_{i-1,j,k}}{𝛥x^2})+(\frac{T_{i,j+1,k}-2T_{i,j,k}+T_{i,j-1,k}}{𝛥y^2}\left)\right)-\\ 𝛥t(\frac{u_{i,j,k}}{𝛥x}(T_{i,j,k}-T_{i-1,j,k})+\frac{v_{i,j,k}}{𝛥y}(T_{i,j+1,k}-T_{i,j,k}))+T_{i,j,k}+\\ 𝛥t{f}_{i,j,k}\end{array}}$$4.5. MATLAB ImplementationFirst, we apply a grid discretization to the the two-dimensional rectangle as follows:

ROWCOLUMNi${\displaystyle i}$ or y${\displaystyle y}$j${\displaystyle j}$ or x${\displaystyle x}$𝛥y𝛥x(1,1)(Nx,Ny)vu

Drawing 12:A discretized two-dimensional plane

Next, we can then implement (65) and (68) in MATLAB, as follows:

1function T = advection_diffusion(T, ... % Temperature matrix2 u, ... % x-velocity matrix3 v, ... % y-velocity matrix 4 dt, ... % time step5 dx, ... % delta x6 dy,... % delta y7 alpha, ... % thermal diffusivity8 f) % external forcing term9r = size(T,1);10c = size(T,2);11 12for i = 2:r-113 for j = 2:c-114 if (u(i,j)>=0 && v(i,j)>=0) % upwind-upwind15 T(i,j) = alpha*dt * ( (T(i+1,j)-2*T(i,j)+T(i-1,j))/dx^2 + ...16 (T(i,j+1)-2*T(i,j)+T(i,j-1))/dy^2 ) - ...17 dt * (u(i,j)/dx * (T(i,j)-T(i-1,j)) + ... 18 v(i,j)/dy * (T(i,j)-T(i,j-1)) ) + ...19 T(i,j) + f(i,j)*dt;20 21 elseif (u(i,j)<0 && v(i,j)<0) % downwind-downwind22 T(i,j) = alpha*dt * ( (T(i+1,j)-2*T(i,j)+T(i-1,j))/dx^2 + ...23 (T(i,j+1)-2*T(i,j)+T(i,j-1))/dy^2 ) - ...24 dt * (u(i,j)/dx * (T(i+1,j)-T(i,j)) + ... 25 v(i,j)/dy * (T(i,j+1)-T(i,j)) ) + ...26 T(i,j) + f(i,j)*dt;27 28 elseif (u(i,j)>=0 && v(i,j)<0) % upwind-downwind29 T(i,j) = alpha*dt * ( (T(i+1,j)-2*T(i,j)+T(i-1,j))/dx^2 + ...30 (T(i,j+1)-2*T(i,j)+T(i,j-1))/dy^2 ) - ...31 dt * (u(i,j)/dx * (T(i,j)-T(i-1,j)) + ... 32 v(i,j)/dy * (T(i,j+1)-T(i,j)) ) + ...33 T(i,j) + f(i,j)*dt;34 35 elseif (u(i,j)<0 && v(i,j)>=0) % downwind-upwind36 T(i,j) = alpha*dt * ( (T(i+1,j)-2*T(i,j)+T(i-1,j))/dx^2 + ...37 (T(i,j+1)-2*T(i,j)+T(i,j-1))/dy^2 ) - ...38 dt * (u(i,j)/dx * (T(i+1,j)-T(i,j)) + ... 39 v(i,j)/dy * (T(i,j)-T(i,j-1)) ) + ...40 T(i,j) + f(i,j)*dt;41 end42 end43end44end

In the code listing above, we can see that the implemented functions never update the boundaries. Thus, we can conveniently set the boundary condition only once for the Dirichlet boundary condition. However, for he Neumann boundary condition, we will need to create phantom nodes. This topic will be further discussed in the next section. Now that we have covered all possible directions for the velocity vector and implemented them in MATLAB, we will then test our MATLAB implementation in the following scenario. First, we will define a circular vector field on a thin plate. This vector field represents the velocity vector. Next, we constrained the temperature at a certain location in that plate (heat source) and watch how the heat propagates.

Figure 10:Heat advection on a plate with circular velocity vector field

4.6. Neumann Boundary Condition

Phantom nodesLeft edgeRight edgeTop edgeBottom edgej=1j=2… …i=1i=2⋮⋮yxj=Nxi=NyActual nodesℓh

Drawing 13:A discretized two-dimensional plane

4.6.1. Left EdgeAssuming that boundary condition for the left edge is 𝜕T(x,0,t)

𝜕x=A${\displaystyle \frac{𝜕T(x,0,t)}{𝜕x}=A}$. thus we have:(71)T(i,2)-T(i,0)

2𝛥x=A$$\frac{T(i,2)-T(i,0)}{2𝛥x}=A$$Hence:(72)

T(i,0)=T(i,2)-2𝛥xA

$$\boxed{T(i,0)=T(i,2)-2𝛥xA}$$4.6.2. Right EdgeAssuming that boundary condition for the left edge is 𝜕T(x,h,t)

𝜕x=B${\displaystyle \frac{𝜕T(x,h,t)}{𝜕x}=B}$. thus we have:(73)T(i,Nx+1)-T(i,Nx-1)

2𝛥x=B$$\frac{T(i,N_x+1)-T(i,N_x-1)}{2𝛥x}=B$$Hence(74)

T(i,Nx+1)=T(i,Nx-1)+2𝛥xB

$$\boxed{T(i,N_x+1)=T(i,N_x-1)+2𝛥xB}$$4.6.3. Top EdgeAssuming that boundary condition for the left edge is 𝜕T(0,y,t)

𝜕y=C${\displaystyle \frac{𝜕T(0,y,t)}{𝜕y}=C}$. thus we have:(75)T(2,j)-T(0,j)

2𝛥y=C$$\frac{T(2,j)-T(0,j)}{2𝛥y}=C$$Hence(76)

T(0,j)=T(2,j)-2𝛥yC

$$\boxed{T(0,j)=T(2,j)-2𝛥yC}$$4.6.4. Bottom EdgeAssuming that boundary condition for the left edge is 𝜕T(ℓ,y,t)

𝜕y=D${\displaystyle \frac{𝜕T(\ell ,y,t)}{𝜕y}=D}$. thus we have:(77)T(Nx+1,j)-T(Nx-1,j)

2𝛥y=D$$\frac{T(N_x+1,j)-T(N_x-1,j)}{2𝛥y}=D$$Hence(78)

T(Ny+1,j)=T(Ny-1,j)+2𝛥yD

$$\boxed{T(N_y+1,j)=T(N_y-1,j)+2𝛥yD}$$ 4.7. Energy Balance Model: Out-of-Plane Heat Transfer for a Thin PlateUp to now, we only deal with homogenous equation. However, when the system interacts with its environment or with other system, the governing equations become non non-homogenous. To demonstrate this, we will implement the following scenario in MATLAB and compare the result with COMSOL. This scenario is taken from the COMSOL Application Gallery8. A very thin plate (1 m×1 m${\displaystyle 1\text{m}\times 1\text{m}}$, and thickness of 0.001 m${\displaystyle 0.001\text{m}}$) is heated only at the left edge. The temperature applied heat is 800 K${\displaystyle 800\text{K}}$. Other edges are thermally insulated. Heat transfer to the external environment only happens from the two surfaces of the plate. 4.7.1. Transient CaseFor the transient case, the governing heat equation of the system above is as follows:(79)a

dz𝜌cp𝜕T 𝜕t	=dzK(𝜕2T 𝜕x2+𝜕2T 𝜕y2)-hu(T-T∞)-hd(T-T∞)
	-𝜀u𝜎(T4-T4∞)-𝜀d𝜎(T4-T4∞)

$$\begin{array}{rl}{d}_z𝜌c_p\frac{𝜕T}{𝜕t}& =d_z\mathcal{K}(\frac{{𝜕}^{2}T}{𝜕x^2}+\frac{{𝜕}^{2}T}{𝜕y^2})-h_{\mathrm{u}}(T-T_{\infty }^{})-h_{\mathrm{d}}(T-T_{\infty }^{})\\ & -{𝜀}_{\mathrm{u}}𝜎(T^4-T_{\infty }^4)-{𝜀}_{\mathrm{d}}𝜎(T^4-T_{\infty }^4)\end{array}$$wheredz=0.01 m${\displaystyle d_z=0.01\text{m}}$hu=hd=h=10 W/(m2K)${\displaystyle h_{\text{u}}=h_{\text{d}}=h=10\text{W}/\left({\text{m}}^2\text{K}\right)}$𝜀u=𝜀d=𝜀=0.5${\displaystyle {𝜀}_{\text{u}}={𝜀}_{\text{d}}=𝜀=0.5}$𝜎=5.67×10-8 W/(m2K4)${\displaystyle 𝜎=5.67\times {10}^{-8}\text{W}/\left({\text{m}}^2{\text{K}}^4\right)}$The material is copper where:cp=385 J/(kg K)${\displaystyle c_p=385\text{J}/\left(\text{kg}\text{K}\right)}$K=400 W/(m K)${\displaystyle \mathcal{K}=400\text{W}/\left(\text{m}\text{K}\right)}$𝜌=8960 kg/m3${\displaystyle 𝜌=8960\text{kg}/{\text{m}}^3}$ We can simplify (79) as follows:(80)a

𝜕T 𝜕t

=K 𝜌cp(𝜕2T 𝜕x2+𝜕2T 𝜕y2)-2h dz𝜌cp(T-T∞)-2𝜀𝜎 dz𝜌cp(T4-T4∞)

$$\begin{array}{rl}\frac{𝜕T}{𝜕t}& =\frac{\mathcal{K}}{𝜌c_p}(\frac{{𝜕}^{2}T}{𝜕x^2}+\frac{{𝜕}^{2}T}{𝜕y^2})-\frac{2h}{d_z𝜌c_p}(T-T_{\infty }^{})-\frac{2𝜀𝜎}{d_z𝜌c_p}(T^4-T_{\infty }^4)\end{array}$$The first thing we can do to see of the formula in (80) is correct is by checking the correctness of the unis. To check the units, we do it as follows:K

s= W/(m K)

kg/m3 J/(kg K)(K

m2+K

m2)-W/(m2K)

m kg/m3 J/(kg K)K-W/(m2K4)

m kg/m3 J/(kg K)K4K

s=K

s-K

s-K

s$$\begin{array}{c}\frac{\text{K}}{\text{s}}=\frac{\text{W}/\left(\text{m}\text{K}\right)}{\text{kg}/{\text{m}}^3\text{J}/\left(\text{kg}\text{K}\right)}(\frac{\text{K}}{{\text{m}}^2}+\frac{\text{K}}{{\text{m}}^2})-\frac{\text{W}/\left({\text{m}}^2\text{K}\right)}{\text{m}\text{kg}/{\text{m}}^3\text{J}/\left(\text{kg}\text{K}\right)}\text{K}-\frac{\text{W}/\left({\text{m}}^2{\text{K}}^4\right)}{\text{m}\text{kg}/{\text{m}}^3\text{J}/\left(\text{kg}\text{K}\right)}{\text{K}}^4\\ \frac{\text{K}}{\text{s}}=\frac{\text{K}}{\text{s}}-\frac{\text{K}}{\text{s}}-\frac{\text{K}}{\text{s}}\end{array}$$The units are indeed correct. MATLAB implementation of (80) is as follows:

1dz = 0.01;2cp = 385;3K = 400;4rho = 8960;5alpha = K/(rho*cp);6beta = dz*rho*cp;7 8dx = 0.01;9dy = 0.01;10dt = 0.001; % Small enough to not cause instability11 12ly = 1; % 1x1 m^213lx = 1;14Nx = lx/dx;15Ny = ly/dy;16Nt = tF/dt;17 18tF = 10000; % Simulation time, modify accordingly19 20h = 10;21sigma = 5.67e-8;22e = 0.5;23T_inf = 300;24 25clear T_array u_array v_array f_array;26 27% Some extra nodes for phantom nodes28u_array = zeros(Nx+2,Ny+1); % Disable advection29v_array = zeros(Nx+2,Ny+1); % Disable advection30T_array = 297.15*ones(Nx+2,Ny+1); 31f_array = zeros(Nx+2,Ny+1);32 33% Preparation for plotting34hfig = figure;35axis off;36hold on;37title('Diffusion')38caxis([300 800])39 40for k = 1 : Nt 41 % Heat source at the center of the plate, keep it true at all time42 T_array = apply_bc(T_array, dx, dy, ...43 ["Dirichlet", "Neumann", "Neumann", "Neumann"], ...44 [800, 0, 0, 0]);45 46 % Calculate the forcing component47 f_array(2:end-1,2:end-1) = -(2*h/beta*(T_array(2:end-1,2:end-1) - T_inf) ...48 + 2*e*sigma/beta*(T_array(2:end-1,2:end-1).^4 - T_inf^4)); 49 50 % Solve the PDE 51 T_array = advection_diffusion_upwind(T_array, ...52 u_array, ...53 v_array, ...54 dt, ...55 dx, ...56 dy, ...57 alpha, ...58 f_array);59end60 61heatmap(T_array, false);

In the listing above, we disable the advection by setting u=v=0${\displaystyle u=v=0}$. The main iteration starts from line 40 to line 59. There are three things we put in the iteration: boundary condition update, external forcing component calculation, and PDE solving. After 10000 second, we will acquire the following result:

Figure 11:Temperature distribution of the plate

We can compare the result with COMSOL. However, the provided COMSOL file only simulates the stationary case.4.7.2. Stationary CaseFor the stationary case, the governing heat equation of the system above is as follows:(81)a

0	=K 𝜌cp(𝜕2T 𝜕x2+𝜕2T 𝜕y2)-hu dz𝜌cp(T-T∞)-hd dz𝜌cp(T-T∞)
	-𝜀u𝜎 dz𝜌cp(T4-T4∞)-𝜀d𝜎 dz𝜌cp(T4-T4∞)

$$\begin{array}{rl}0& =\frac{\mathcal{K}}{𝜌c_p}(\frac{{𝜕}^{2}T}{𝜕x^2}+\frac{{𝜕}^{2}T}{𝜕y^2})-\frac{h_{\mathrm{u}}}{d_z𝜌c_p}(T-T_{\infty }^{})-\frac{h_{\mathrm{d}}}{d_z𝜌c_p}(T-T_{\infty }^{})\\ & -\frac{{𝜀}_{\mathrm{u}}𝜎}{d_z𝜌c_p}(T^4-T_{\infty }^4)-\frac{{𝜀}_{\mathrm{d}}𝜎}{d_z𝜌c_p}(T^4-T_{\infty }^4)\end{array}$$The system in (81) is a 2D Poisson PDE, which we can solve it by using the five-point difference method.5. Navier-Stokes with Heat Transfer In this section, we will model a cross-current heat exchanger by using the famous Navier-Stokes equation for fluid flows with heat transfer. We will start with one flow first. After that, we will implement both flows and the occurring heat transfer between them.5.1. Channel Flow with Heat EquationWe will model the flow of the fluid as a two-dimensional flow. The flow is driven by pressure gradient (planar Poiseuille flow).

Flow InNo-slip wallNo-slip wall 𝜕p

𝜕x=constant${\displaystyle \frac{𝜕p}{𝜕x}=\text{constant}}$ v(x,y,t)=0u(x,y,t)=u(y)kdiffusionu(y)advectionT(x,y,t)T(0,y,t)𝜕T(x,ly,t)

𝜕y=0T(x,y,t)𝜕T(x,o,t)

𝜕y=0𝜕T(lx,y,t)

𝜕x=0Flow Outxyℓxℓyvu

Drawing 14:Heat transfer in a two-dimensional channel flow

Conservation of mass:(82)𝜕u

𝜕x+𝜕v

𝜕y⏠⏣⏣⏣⏡⏣⏣⏣⏢Divergence=0$$\underset{\text{Divergence}}{\underbrace{\frac{𝜕u}{𝜕x}+\frac{𝜕v}{𝜕y}}}=0$$Conservation of linear momentum:(83)𝜕u

𝜕t+u𝜕u

𝜕x+v𝜕u

𝜕y=-1

𝜌𝜕p

𝜕x+𝜈(𝜕2u

𝜕x2+𝜕2u

𝜕y2)(84)𝜕v

𝜕t+u𝜕v

𝜕x+v𝜕v

𝜕y=-1

𝜌𝜕p

𝜕y⏠⏣⏣⏡⏣⏣⏢Gradient+𝜈(𝜕2v

𝜕x2+𝜕2v

𝜕y2)⏠⏣⏣⏣⏡⏣⏣⏣⏢Laplacian$$\begin{array}{c}\frac{𝜕u}{𝜕t}+u\frac{𝜕u}{𝜕x}+v\frac{𝜕u}{𝜕y}=-\frac{1}{𝜌}\frac{𝜕p}{𝜕x}+\mathcal{𝜈}(\frac{{𝜕}^{2}u}{𝜕x^2}+\frac{{𝜕}^{2}u}{𝜕y^2})\\ \frac{𝜕v}{𝜕t}+u\frac{𝜕v}{𝜕x}+v\frac{𝜕v}{𝜕y}=-\underset{\text{Gradient}}{\underbrace{\frac{1}{𝜌}\frac{𝜕p}{𝜕y}}}+\underset{\text{Laplacian}}{\underbrace{𝜈(\frac{{𝜕}^{2}v}{𝜕x^2}+\frac{{𝜕}^{2}v}{𝜕y^2})}}\end{array}$$Conservation of energy:(85)𝜌cp(𝜕T

𝜕t+u𝜕T

𝜕x+v𝜕T

𝜕y)=k(𝜕2T

𝜕x2+𝜕2T

𝜕y2)$$𝜌c_p(\frac{𝜕T}{𝜕t}+u\frac{𝜕T}{𝜕x}+v\frac{𝜕T}{𝜕y})=k(\frac{{𝜕}^{2}T}{𝜕x^2}+\frac{{𝜕}^{2}T}{𝜕y^2})$$where:u(x,y,t)${\displaystyle u(x,y,t)}$ is the x${\displaystyle x}$-velocity at location (x,y)${\displaystyle (x,y)}$ at time t${\displaystyle t}$v(x,y,t)${\displaystyle v(x,y,t)}$ is the y${\displaystyle y}$-velocity at location (x,y)${\displaystyle (x,y)}$ at time t${\displaystyle t}$p(x,y,t)${\displaystyle p(x,y,t)}$ is the pressure at location (x,y)${\displaystyle (x,y)}$ at time t${\displaystyle t}$T(x,y,t)${\displaystyle T(x,y,t)}$ is the temperature at location (x,y)${\displaystyle (x,y)}$ at time t${\displaystyle t}$𝜈${\displaystyle 𝜈}$ is the fluid viscosity coefficient.cp${\displaystyle c_p}$ is the specific heat at constant𝜌${\displaystyle 𝜌}$ is the fluid densityk${\displaystyle k}$ is the thermal diffusivity rate In planar Poiseuille flow, the flow is only in x${\displaystyle x}$-direction. u${\displaystyle u}$ becomes a function of only y${\displaystyle y}$. Therefore: v=0${\displaystyle v=0}$, 𝜕v

𝜕t=0${\displaystyle \frac{𝜕v}{𝜕t}=0}$, 𝜕v

𝜕y=0${\displaystyle \frac{𝜕v}{𝜕y}=0}$, 𝜕2v

𝜕y2=0${\displaystyle \frac{{𝜕}^{2}v}{𝜕y^2}=0}$, 𝜕u

𝜕t=0${\displaystyle \frac{𝜕u}{𝜕t}=0}$, 𝜕u

𝜕x=0${\displaystyle \frac{𝜕u}{𝜕x}=0}$, and 𝜕p

𝜕x=0${\displaystyle \frac{𝜕p}{𝜕x}=0}$. As a result, (82) to (84) can be simplified as follows:(86)1

𝜌(𝜕p

𝜕x)=𝜈𝜕2u

𝜕y2$$\frac{1}{𝜌}\left(\frac{𝜕p}{𝜕x}\right)=\mathcal{𝜈}\frac{{𝜕}^{2}u}{𝜕y^2}$$Since 𝜕p

𝜕x${\displaystyle \frac{𝜕p}{𝜕x}}$ is known, we are left with only finding u(y)${\displaystyle u(y)}$ as the solution to (86). Thus, now we have an ODE instead of a PDE. As an ODE, we can directly solve 1

𝜌(dp

dx)=𝜈d2u

dy2${\displaystyle \frac{1}{𝜌}\left(\frac{dp}{dx}\right)=\mathcal{𝜈}\frac{d^{2}u}{d{y}^2}}$ as follows:u(y)=1

𝜌𝜈(dp

dx)y2

2+C1y+C2$$u(y)=\frac{1}{𝜌𝜈}\left(\frac{dp}{dx}\right)\frac{y^2}{2}+C_{1}y+C_2$$At no-slip condition, u(0)=0${\displaystyle u(0)=0}$ and u(ly)=0${\displaystyle u(l_y)=0}$. By using these conditions, we can find C1${\displaystyle C_1}$ and C2${\displaystyle C_2}$. Let us the first boundary condition:u(0)=0=C2$$u(0)=0=C_2$$Next, let us take the second boundary:u(ly)=0=1

𝜌𝜈(dp

dx)l2y

2+C1lyC1=-1

𝜌𝜈(dp

dx)ly

2$$\begin{array}{c}u(l_y)=0=\frac{1}{𝜌𝜈}\left(\frac{dp}{dx}\right)\frac{l_y^2}{2}+C_1{l}_y\\ C_1=-\frac{1}{𝜌𝜈}\left(\frac{dp}{dx}\right)\frac{l_y}{2}\end{array}$$Hence:a

u(y)	=1 𝜌𝜈(dp dx)y2 2+(-1 𝜌𝜈(dp dx)ly 2)⏠⏣⏣⏣⏡⏣⏣⏣⏢C1y
	=1 𝜌𝜈(dp dx)y 2(y-ly)

$$\begin{array}{rl}u(y)& =\frac{1}{𝜌𝜈}\left(\frac{dp}{dx}\right)\frac{y^2}{2}+\underset{C_1}{\underbrace{(-\frac{1}{𝜌𝜈}(\frac{dp}{dx}\left)\frac{l_y}{2}\right)}}y\\ & =\frac{1}{𝜌𝜈}\left(\frac{dp}{dx}\right)\frac{y}{2}(y-l_y)\end{array}$$ 5.2. Two Counter-Direction Channel Flows with Heat EquationsGiven below a simple two-dimensional model of a cross-current heat exchanger. The model is compose of two channel flows. Between the two channel, heat transfer occurs.

Flow InFlow OutlyNo-slip wallℓx𝜕p(x,y,t)

𝜕x=P1(x)${\displaystyle \frac{𝜕p(x,y,t)}{𝜕x}=P_1(x)}$ Flow In𝜕p(x,y,t)

𝜕x=P2(x)${\displaystyle \frac{𝜕p(x,y,t)}{𝜕x}=P_2(x)}$ No-slip wallNo-slip wallNo-slip wall𝛼1𝛼2hChHu1(y)u2(y)diffusiondiffusionadvectionadvectionconvectionconvectionTH(x,y,t)TC(x,y,t)TH(0,y,0)TC(ℓx,y,0)Flow Outly

Drawing 15:Heat transfer in a Two-dimensional model of a cross-current heat exchanger

Conservation of mass:(87)𝜕u

𝜕x+𝜕v

𝜕y=0$$\frac{𝜕u}{𝜕x}+\frac{𝜕v}{𝜕y}=0$$Conservation of linear momentum:(88)𝜕u

𝜕t+u𝜕u

𝜕x+v𝜕u

𝜕y=-1

𝜌𝜕p

𝜕x+𝜈(𝜕2u

𝜕x2+𝜕2u

𝜕y2)(89)𝜕v

𝜕t+u𝜕v

𝜕x+v𝜕v

𝜕y=-1

𝜌𝜕p

𝜕y+𝜈(𝜕2v

𝜕x2+𝜕2v

𝜕y2)$$\begin{array}{c}\frac{𝜕u}{𝜕t}+u\frac{𝜕u}{𝜕x}+v\frac{𝜕u}{𝜕y}=-\frac{1}{𝜌}\frac{𝜕p}{𝜕x}+\mathcal{𝜈}(\frac{{𝜕}^{2}u}{𝜕x^2}+\frac{{𝜕}^{2}u}{𝜕y^2})\\ \frac{𝜕v}{𝜕t}+u\frac{𝜕v}{𝜕x}+v\frac{𝜕v}{𝜕y}=-\frac{1}{𝜌}\frac{𝜕p}{𝜕y}+𝜈(\frac{{𝜕}^{2}v}{𝜕x^2}+\frac{{𝜕}^{2}v}{𝜕y^2})\end{array}$$Conservation of energy:(90)𝜌cp(𝜕TH

𝜕t+u𝜕TH

𝜕x+v𝜕TH

𝜕y)=𝛼1(𝜕2TH

𝜕x2+𝜕2TH

𝜕y2)$$𝜌c_p(\frac{𝜕T_H}{𝜕t}+u\frac{𝜕T_H}{𝜕x}+v\frac{𝜕T_H}{𝜕y})={𝛼}_1(\frac{{𝜕}^2{T}_H}{𝜕x^2}+\frac{{𝜕}^2{T}_H}{𝜕y^2})$$(91)𝜌cp(𝜕TC

𝜕t+u𝜕TC

𝜕x+v𝜕TC

𝜕y)=𝛼2(𝜕2TC

𝜕x2+𝜕2TC

𝜕y2)$$𝜌c_p(\frac{𝜕T_C}{𝜕t}+u\frac{𝜕T_C}{𝜕x}+v\frac{𝜕T_C}{𝜕y})={𝛼}_2(\frac{{𝜕}^2{T}_C}{𝜕x^2}+\frac{{𝜕}^2{T}_C}{𝜕y^2})$$

6. Navier-Stokes in A 2D Rectangular World with Explicit Pressure Equation9 The Navier-Stokes equations are used to govern dynamics of moving fluid. In a two dimensional word, the fluid has two velocity components: the velocity along the x${\displaystyle x}$-axis (u${\displaystyle u}$) and the velocity along the y${\displaystyle y}$-axis (v${\displaystyle v}$). Additionally, within the fluid, there is pressure (p${\displaystyle p}$). (93) and (94) are the Navier Stokes equations. However, there are only two equations with three unknowns (u${\displaystyle u}$, v${\displaystyle v}$, and p${\displaystyle p}$). Therefore, since the fluid in this case is assumed to be incompressible, the mass conservation law must be fulfilled as well. With a help from (92), we can generate the equation for the pressure as in (95). Conservation of mass (continuity equation):(92)𝜕u

𝜕x+𝜕v

𝜕y=0$$\frac{𝜕u}{𝜕x}+\frac{𝜕v}{𝜕y}=0$$Conservation of linear momentum:(93)𝜕u

𝜕t+u𝜕u

𝜕x+v𝜕u

𝜕y=-1

𝜌𝜕p

𝜕x+𝜈(𝜕2u

𝜕x2+𝜕2u

𝜕y2)(94)𝜕v

𝜕t+u𝜕v

𝜕x+v𝜕v

𝜕y=-1

𝜌𝜕p

𝜕y+𝜈(𝜕2v

𝜕x2+𝜕2v

𝜕y2)$$\begin{array}{c}\frac{𝜕u}{𝜕t}+u\frac{𝜕u}{𝜕x}+v\frac{𝜕u}{𝜕y}=-\frac{1}{𝜌}\frac{𝜕p}{𝜕x}+\mathcal{𝜈}(\frac{{𝜕}^{2}u}{𝜕x^2}+\frac{{𝜕}^{2}u}{𝜕y^2})\\ \frac{𝜕v}{𝜕t}+u\frac{𝜕v}{𝜕x}+v\frac{𝜕v}{𝜕y}=-\frac{1}{𝜌}\frac{𝜕p}{𝜕y}+𝜈(\frac{{𝜕}^{2}v}{𝜕x^2}+\frac{{𝜕}^{2}v}{𝜕y^2})\end{array}$$Poisson equation for the pressure:(95)𝜕2p

𝜕x2+𝜕2p

𝜕y2≈𝜌

𝜕t(𝜕u

𝜕x+𝜕v

𝜕y)$$\frac{{𝜕}^{2}p}{𝜕x^2}+\frac{{𝜕}^{2}p}{𝜕y^2}\approx \frac{𝜌}{𝜕t}(\frac{𝜕u}{𝜕x}+\frac{𝜕v}{𝜕y})$$ In order to derive the Poisson equation for the pressure, first, let us expand (93) and (94) by using the forward Euler method.(96)a

u(x,y,t+𝛥t)

=u(x,y,t)-𝛥t(u𝜕u 𝜕x+v𝜕u 𝜕y+1 𝜌𝜕p 𝜕x+𝜈(𝜕2u 𝜕x2+𝜕2u 𝜕y2))

$$\begin{array}{rl}u(x,y,t+𝛥t)& =u(x,y,t)-𝛥t(u\frac{𝜕u}{𝜕x}+v\frac{𝜕u}{𝜕y}+\frac{1}{𝜌}\frac{𝜕p}{𝜕x}+𝜈(\frac{{𝜕}^{2}u}{𝜕x^2}+\frac{{𝜕}^{2}u}{𝜕y^2}\left)\right)\end{array}$$(97)a

v(x,y,t+𝛥t)

=v(x,y,t)-𝛥t(u𝜕v 𝜕x+v𝜕v 𝜕y+1 𝜌𝜕p 𝜕y+𝜈(𝜕2u 𝜕x2+𝜕2u 𝜕y2))

$$\begin{array}{rl}v(x,y,t+𝛥t)& =v(x,y,t)-𝛥t(u\frac{𝜕v}{𝜕x}+v\frac{𝜕v}{𝜕y}+\frac{1}{𝜌}\frac{𝜕p}{𝜕y}+𝜈(\frac{{𝜕}^{2}u}{𝜕x^2}+\frac{{𝜕}^{2}u}{𝜕y^2}\left)\right)\end{array}$$Next, we take the derivative of u(x,y,t+𝛥t)${\displaystyle u(x,y,t+𝛥t)}$ with respect to x${\displaystyle x}$ and v(x,y,t+𝛥t)${\displaystyle v(x,y,t+𝛥t)}$ with respect to y${\displaystyle y}$:(98)𝜕u(x,y,t+𝛥t)

𝜕x=𝜕u(x,y,t)

𝜕x-𝛥t(1

𝜌𝜕2p

𝜕x2+𝜕

𝜕x(u𝜕u

𝜕x+v𝜕u

𝜕y+𝜈(𝜕2u

𝜕x2+𝜕2u

𝜕y2))⏠⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏢A)(99)𝜕v(x,y,t+𝛥t)

𝜕y=𝜕v(x,y,t)

𝜕y-𝛥t(1

𝜌𝜕2p

𝜕y2+𝜕

𝜕y(u𝜕v

𝜕x+v𝜕v

𝜕y+𝜈(𝜕2v

𝜕x2+𝜕2v

𝜕y2))⏠⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏢B)$$\begin{array}{c}\frac{𝜕u(x,y,t+𝛥t)}{𝜕x}=\frac{𝜕u(x,y,t)}{𝜕x}-𝛥t(\frac{1}{𝜌}\frac{{𝜕}^{2}p}{𝜕x^2}+\underset{A}{\underbrace{\frac{𝜕}{𝜕x}(u\frac{𝜕u}{𝜕x}+v\frac{𝜕u}{𝜕y}+𝜈(\frac{{𝜕}^{2}u}{𝜕x^2}+\frac{{𝜕}^{2}u}{𝜕y^2}\left)\right)}})\\ \frac{𝜕v(x,y,t+𝛥t)}{𝜕y}=\frac{𝜕v(x,y,t)}{𝜕y}-𝛥t(\frac{1}{𝜌}\frac{{𝜕}^{2}p}{𝜕y^2}+\underset{B}{\underbrace{\frac{𝜕}{𝜕y}(u\frac{𝜕v}{𝜕x}+v\frac{𝜕v}{𝜕y}+𝜈(\frac{{𝜕}^{2}v}{𝜕x^2}+\frac{{𝜕}^{2}v}{𝜕y^2}\left)\right)}})\end{array}$$We then add (98) to (99), giving us the following:(100)𝜕u(x,y,t+𝛥t)

𝜕x+𝜕v(x,y,t+𝛥t)

𝜕y⏠⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏢=0=𝜕u(x,y,t)

𝜕x+𝜕v(x,y,t)

𝜕y-𝛥t

𝜌(𝜕2p

𝜕x2+𝜕2p

𝜕y2)-A𝛥t-B𝛥t$$\underset{=0}{\underbrace{\frac{𝜕u(x,y,t+𝛥t)}{𝜕x}+\frac{𝜕v(x,y,t+𝛥t)}{𝜕y}}}=\frac{𝜕u(x,y,t)}{𝜕x}+\frac{𝜕v(x,y,t)}{𝜕y}-\frac{𝛥t}{𝜌}(\frac{{𝜕}^{2}p}{𝜕x^2}+\frac{{𝜕}^{2}p}{𝜕y^2})-A𝛥t-B𝛥t$$We set the left hand side of (100) to zero since the continuity equation (92) must be fulfilled. Rearranging (100) then gives us:(101)

𝜕2p

𝜕x2+𝜕2p

𝜕y2=𝜌

𝛥t[𝜕u(x,y,t)

𝜕x+𝜕v(x,y,t)

𝜕y]-𝜌[A+B]⏠⏣⏣⏣⏡⏣⏣⏣⏢ignore this part

$$\boxed{\frac{{𝜕}^{2}p}{𝜕x^2}+\frac{{𝜕}^{2}p}{𝜕y^2}=\frac{𝜌}{𝛥t}[\frac{𝜕u(x,y,t)}{𝜕x}+\frac{𝜕v(x,y,t)}{𝜕y}]-\underset{\text{ignore this part}}{\underbrace{𝜌[A+B]}}}$$Since 𝜌

𝛥t≫𝜌${\displaystyle \frac{𝜌}{𝛥t}\gg 𝜌}$, we can then ignore the rightmost term of (101). 6.1. Discretizing the x${\displaystyle x}$-momentum Equation(93) is discretized as follows:(102)a

u(x,y,t+𝛥t)-u(x,y,t) 𝛥t+

u(x,y,t)u(x+𝛥,y,t)-u(x-𝛥x,y,t) 2𝛥x+v(x,y,t)u(x,y+𝛥y,t)-u(x,y-𝛥y,t) 2𝛥y=

-1 𝜌p(x+𝛥x,y)-p(x-𝛥x,y) 2𝛥x

+𝜈(u(x+𝛥x,y,t)-2u(x,y,t)+u(x-𝛥x,y,t) 𝛥x2+u(x,y+𝛥y,t)-2u(x,y,t)+u(x,y-𝛥y,t) 𝛥y2)

$$\begin{array}{c}\frac{u(x,y,t+𝛥t)-u(x,y,t)}{𝛥t}+\\ u(x,y,t)\frac{u(x+𝛥,y,t)-u(x-𝛥x,y,t)}{2𝛥x}+v(x,y,t)\frac{u(x,y+𝛥y,t)-u(x,y-𝛥y,t)}{2𝛥y}=\\ -\frac{1}{𝜌}\frac{p(x+𝛥x,y)-p(x-𝛥x,y)}{2𝛥x}\\ +𝜈(\frac{u(x+𝛥x,y,t)-2u(x,y,t)+u(x-𝛥x,y,t)}{𝛥x^2}+\frac{u(x,y+𝛥y,t)-2u(x,y,t)+u(x,y-𝛥y,t)}{𝛥y^2})\end{array}$$ Rearranging (102) gives us:

a

u(x,y,t+𝛥t)

=-𝛥t 2𝜌𝛥x[p(x+𝛥x,y)-p(x-𝛥x,y)]

+𝛥t𝜈(u(x+𝛥x,y,t)-2u(x,y,t)+u(x-𝛥x,y,t) 𝛥x2+u(x,y+𝛥y,t)-2u(x,y,t)+u(x,y-𝛥y,t) 𝛥y2)

-𝛥t 2𝛥xu(x,y,t)[u(x+𝛥x,y,t)-u(x-𝛥x,y,t)]-𝛥t 2𝛥yv(x,y,t)[u(x,y+2𝛥y,t)-u(x,y-𝛥y,t)]+u(x,y,t)

$$\boxed{\begin{array}{c}u(x,y,t+𝛥t)\\ =-\frac{𝛥t}{2𝜌𝛥x}[p(x+𝛥x,y)-p(x-𝛥x,y)]\\ +𝛥t𝜈(\frac{u(x+𝛥x,y,t)-2u(x,y,t)+u(x-𝛥x,y,t)}{𝛥x^2}+\frac{u(x,y+𝛥y,t)-2u(x,y,t)+u(x,y-𝛥y,t)}{𝛥y^2})\\ -\frac{𝛥t}{2𝛥x}u(x,y,t)[u(x+𝛥x,y,t)-u(x-𝛥x,y,t)]-\frac{𝛥t}{2𝛥y}v(x,y,t)[u(x,y+2𝛥y,t)-u(x,y-𝛥y,t)]+u(x,y,t)\end{array}}$$ 6.2. Discretizing the y${\displaystyle y}$-momentum Equation(94) is discretized as follows:(103)a

v(x,y,t+𝛥t)-v(x,y,t) 𝛥t+

u(x,y,t)v(x+𝛥x,y,t)-v(x-𝛥x,y,t) 2𝛥x+v(x,y,t)v(x+2𝛥y,y,t)-v(x,y-𝛥y,t) 2𝛥y=

-1 𝜌p(x,y+𝛥y)-p(x,y-𝛥y) 2𝛥y

+𝜈(v(x+𝛥x,y,t)-2v(x,y,t)+v(x-𝛥x,y,t) 𝛥x2+v(x,y+𝛥y,t)-2v(x,y,t)+v(x,y-𝛥y,t) 𝛥y2)

$$\begin{array}{c}\frac{v(x,y,t+𝛥t)-v(x,y,t)}{𝛥t}+\\ u(x,y,t)\frac{v(x+𝛥x,y,t)-v(x-𝛥x,y,t)}{2𝛥x}+v(x,y,t)\frac{v(x+2𝛥y,y,t)-v(x,y-𝛥y,t)}{2𝛥y}=\\ -\frac{1}{𝜌}\frac{p(x,y+𝛥y)-p(x,y-𝛥y)}{2𝛥y}\\ +𝜈(\frac{v(x+𝛥x,y,t)-2v(x,y,t)+v(x-𝛥x,y,t)}{𝛥x^2}+\frac{v(x,y+𝛥y,t)-2v(x,y,t)+v(x,y-𝛥y,t)}{𝛥y^2})\end{array}$$Rearranging (103) gives us:

a

v(x,y,t+𝛥t)

=-𝛥t 2𝜌𝛥y[p(x,y+𝛥y)-p(x,y-𝛥y)]

+𝛥t𝜈(v(x+𝛥x,y,t)-2v(x,y,t)+v(x-𝛥x,y,t) 𝛥x2+v(x,y+𝛥y,t)-2v(x,y,t)+v(x,y-𝛥y,t) 𝛥y2)

-𝛥t 2𝛥xu(x,y,t)[v(x+𝛥x,y,t)-v(x-𝛥x,y,t)]-𝛥t 2𝛥yv(x,y,t)[v(x,y+𝛥y,t)-v(x,y-𝛥y,t)]+v(x,y,t)

$$\boxed{\begin{array}{c}v(x,y,t+𝛥t)\\ =-\frac{𝛥t}{2𝜌𝛥y}[p(x,y+𝛥y)-p(x,y-𝛥y)]\\ +𝛥t𝜈(\frac{v(x+𝛥x,y,t)-2v(x,y,t)+v(x-𝛥x,y,t)}{𝛥x^2}+\frac{v(x,y+𝛥y,t)-2v(x,y,t)+v(x,y-𝛥y,t)}{𝛥y^2})\\ -\frac{𝛥t}{2𝛥x}u(x,y,t)[v(x+𝛥x,y,t)-v(x-𝛥x,y,t)]-\frac{𝛥t}{2𝛥y}v(x,y,t)[v(x,y+𝛥y,t)-v(x,y-𝛥y,t)]+v(x,y,t)\end{array}}$$ 6.3. Discretizing the Pressure Equation(102) is discretized as follows: p(x+𝛥x,y)-2p(x,y)+p(x-𝛥x,y)

𝛥x2+p(x,y+𝛥y)-2p(x,y)+p(x,y-𝛥y)

𝛥y2≈(104)𝜌

𝛥t(u(x+𝛥x,y)-u(x-𝛥x,y)

2𝛥x+v(x,y+𝛥y)-v(x,y-𝛥y)

2𝛥y)|⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢b(x,y)$$\begin{array}{c}\frac{p(x+𝛥x,y)-2p(x,y)+p(x-𝛥x,y)}{𝛥x^2}+\frac{p(x,y+𝛥y)-2p(x,y)+p(x,y-𝛥y)}{𝛥y^2}\approx \\ \underset{b(x,y)}{\underbrace{\frac{𝜌}{𝛥t}(\frac{u(x+𝛥x,y)-u(x-𝛥x,y)}{2𝛥x}+\frac{v(x,y+𝛥y)-v(x,y-𝛥y)}{2𝛥y})}}\end{array}$$ Rearranging gives us: (105)

a

p(x,y)

=-𝛥x2𝛥y2b(x,y)+𝛥y2[p(x-𝛥x,y)+p(x+𝛥x,y)]+𝛥x2[p(x,y+𝛥y)+p(x,y-𝛥y)] 2(𝛥x2+𝛥y2)

$$\boxed{\begin{array}{rl}p(x,y)& =\frac{-𝛥x^2𝛥y^{2}b(x,y)+𝛥y^2[p(x-𝛥x,y)+p(x+𝛥x,y)]+𝛥x^2[p(x,y+𝛥y)+p(x,y-𝛥y)]}{2(𝛥x^2+𝛥y^2)}\end{array}}$$ where:(106)

b(x,y)=𝜌

𝛥t(u(x+𝛥x,y)-u(x-𝛥x,y)

2𝛥x+v(x,y+𝛥y)-v(x,y-𝛥y)

2𝛥y)

$$\boxed{b(x,y)=\frac{𝜌}{𝛥t}(\frac{u(x+𝛥x,y)-u(x-𝛥x,y)}{2𝛥x}+\frac{v(x,y+𝛥y)-v(x,y-𝛥y)}{2𝛥y})}$$ $$ 6.4. Cavity Flow

u=Re 𝜈

𝜌 Lxyv=0u=0v=0u=0v=0u=0v=0𝜕p

𝜕y=0𝜕p

𝜕y=0𝜕p

𝜕y=0p=0lx=1ly=1

Drawing 16:Typical boundary conditions for a cavity flow.

6.5. Channel Flow6.5.1. Example 1

xyu=0v=0u=0v=0u=0v=0u=0v=0𝜕p

𝜕y=0𝜕p

𝜕y=0p=5lx=3ly=0.1p=2InletOutet

Drawing 17:Typical boundary conditions for a channel flow.

6.5.2. Example 2

xyu=0v=0u=0v=0u=0v=0u=0v=0𝜕p

𝜕y=0𝜕p

𝜕y=0p=5lx=3ly=0.1p=2InletOutet

Drawing 18:Typical boundary conditions for a channel flow.

7. Steady Simple Euler-Bernoulli Beam In this section, we will discuss a simple cantilever problem. Simple means, it single span beam. Its maximum number of supports is two, located at both ends of the beam. Multispan cantilever beam will requires a more sophisticated boundary value problem solver. 7.1. Mathematical ModelingThe Euler–Bernoulli equation describes the relationship between the beam's deflection and the applied load.

f(x)x=0x=Lx

Drawing 19:An Euler-Bernoulli Beam

For a one-dimensional elastic beam with a length of L${\displaystyle L}$, its deflection is governed by the Euler-Bernoulli equation as follows:(107)d2

dx2[EI d2w(x)

dx2]=f(x)$$\frac{d^2}{d{x}^2}\left[EI\frac{d^{2}w(x)}{d{x}^2}\right]=f(x)$$In (107), w(x)${\displaystyle w(x)}$ and f(x)${\displaystyle f(x)}$ are the beam's deflection and the external load at location x${\displaystyle x}$, respectively. E${\displaystyle E}$ is the Young's modulus, and I${\displaystyle I}$ is the area moment of inertia. We start by assuming that E${\displaystyle E}$ and I${\displaystyle I}$ are constants. Hence, (107) can be rewritten as:(108)d4w(x)

dx4=1

EI f(x)$$\frac{d^{4}w(x)}{d{x}^4}=\frac{1}{EI}f(x)$$The initial conditions to (108) are made by:a

d2w(x) dx2	=1 EI M(x)	Internal moment at x
d3w(x) dx3	=1 EI F(x)	Internal shear force at x

$$\begin{array}{rlr}\frac{d^{2}w(x)}{d{x}^2}& =\frac{1}{EI}M(x)& \text{Internal moment at }x\\ \frac{d^{3}w(x)}{d{x}^3}& =\frac{1}{EI}F(x)& \text{Internal shear force at }x\end{array}$$In a central difference method we know that:(109)a

	dw(x) dxaxi≈wi-1-wi+1 2𝛥x
	d2w(x) dx2axi≈wi-1-2wi+wi+1 𝛥x2
	d3w(x) x3axi≈-wi-2+2wi-1-2wi+1+wi+2 2𝛥x3
	d4w(x) dxaxi≈wi-2-4wi-1+6wi-4wi+1+wi+2 𝛥x4

$$\begin{array}{rl}& \frac{dw(x)}{dx}{|}_{x_i}\approx \frac{w_{i-1}-w_{i+1}}{2𝛥x}\\ & \frac{d^{2}w(x)}{d{x}^2}{|}_{x_i}\approx \frac{w_{i-1}-2w_i+w_{i+1}}{𝛥x^2}\\ & \frac{d^{3}w(x)}{x^3}{|}_{x_i}\approx \frac{-w_{\mathrm{i}-2}+2w_{i-1}-2w_{i+1}+w_{i+2}}{2𝛥x^3}\\ & \frac{d^{4}w(x)}{dx}{|}_{x_i}\approx \frac{w_{i-2}-4w_{i-1}+6w_{\mathrm{i}}-4w_{i+1}+w_{i+2}}{𝛥x^4}\end{array}$$where wi=w(i𝛥x)${\displaystyle w_i=w(i𝛥x)}$ and 𝛥x${\displaystyle 𝛥x}$ is the segment length. Hence, we can rewrite (108) as follows"(110)a

d4w(x) dx4	=1 EI f(x)
n∑i=0wi-2-4wi-1+6wi-4wi+1+wi+2 (𝛥x)4	=1 EI fi

$$\begin{array}{rl}\frac{d^{4}w(x)}{d{x}^4}& =\frac{1}{EI}f(x)\\ \sum _{i=0}^n\frac{w_{i-2}-4w_{i-1}+6w_{\mathrm{i}}-4w_{i+1}+w_{i+2}}{(𝛥x{)}^4}& =\frac{1}{EI}{f}_i\end{array}$$where fi=f(i𝛥x)${\displaystyle f_i=f(i𝛥x)}$.We can also express in its matrix form as follows:(111)a

[a 1 0 … 0 0 1 0 … 0 1 -4 6 -4 1 0 … 0 0 1 -4 6 -4 1 … 0 ⋮ ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ 0 … 0 1 -4 6 -4 1 0 … 0 1 0 0 … 0 1 ]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a w-2 w-1 w0 w1 ⋮ wn wn+1 wn+2 ]⏠⏣⏣⏣⏡⏣⏣⏣⏢W

=𝛥x4 EI[a L2 L1 f0 f1 ⋮ fn R1 R2 ]⏠⏣⏣⏣⏡⏣⏣⏣⏢F

$$\begin{array}{cc}\underset{A}{\underbrace{\left[\begin{array}{cccccccc}1& 0& & & & & \dots & 0\\ 0& 1& 0& & & & \dots & 0\\ 1& -4& 6& -4& 1& 0& \dots & 0\\ 0& 1& -4& 6& -4& 1& \dots & 0\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \dots & 0& 1& -4& 6& -4& 1\\ 0& & & & \dots & 0& 1& 0\\ 0& & & & & \dots & 0& 1\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_{-2}\\ w_{-1}\\ w_0\\ w_1\\ ⋮\\ w_n\\ w_{n+1}\\ w_{n+2}\end{array}\right]}}& =\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\begin{array}{c}{L}_2\\ L_1\\ f_0\\ f_1\\ ⋮\\ f_n\\ R_1\\ R_2\end{array}\right]}}\end{array}$$where w-1${\displaystyle w_{-1}}$, w-2${\displaystyle w_{-2}}$, wn${\displaystyle w_n}$, wn+1${\displaystyle w_{n+1}}$, L1${\displaystyle L_1}$, L2${\displaystyle L_2}$, R1${\displaystyle R_1}$, and R2${\displaystyle R_2}$ are all the properties of additional nodes that we must add . We call them as the phantom nodes (see Drawing 20).

w0… wnw1wn+1w-1w2w3w-2wn+2Phantom nodesPhantom nodesf0f1f2L1L2R1R2f3… fn

Drawing 20:. Finite difference model of a one-dimensional problem with extra phantom nodes for the boundary conditions

7.2. Engineering InterpretationThe solutions to (108) are subjected to the given initial conditions (an initial value problem). However, in engineering, it is more common to call these initial values as boundary conditions since the independent variables are spatial and not temporal. Also, (108) is actually originated from a partial differential equation. There are three possible boundary conditions for each end, which are:• Fixed-end (clamped), which means that there is no displacement at the current end. Additionally, at the current end, its horizontal tangent is also zero.• Pinned-end, which means that there is no displacement at the current end. Also, at the current end, its bending moment is also zero.• Free-end, which means that both bending moment and shear force are zero at the current end.As additional information, bending moment is the second derivative of w${\displaystyle w}$ with respect to x${\displaystyle x}$ and shear force is the third derivative of w${\displaystyle w}$ with respect to x${\displaystyle x}$. The combinations of those three boundaries create nine possible configurations. However, we will focus only on 6 configuration as shown in Drawing 21.

f(x)f(x)f(x)f(x)f(x)Both ends are fixed/supportedOne end is fixed andanother end is freeMOne end is supported andan explicit torque is applied at another end

Drawing 21:Engineering interpretation of the Euler-Bernoulli equation

Now, let us explore the boundary conditions together with their discrete forms in details. 7.2.1. Boundary Conditions at the Left Enda) Fixed-end Boundary (Clamped)• No position displacement means that w(0)=0${\displaystyle w(0)=0}$. This then translates into:(112)w0=0$$w_0=0$$• Zero horizontal tangent means that dw(0)

dx=0${\displaystyle \frac{dw(0)}{dx}=0}$. which then translates into:(113)w-1-w1

2𝛥x=0⟹w-1=w1$$\frac{w_{-1}-w_1}{2𝛥x}=0⟹w_{-1}=w_1$$ b) Pinned-end Boundary• No position displacement means that w(0)=0${\displaystyle w(0)=0}$, which translates to:(114)w0=0$$w_0=0$$• No bending moment means thatd2w(0)

dx2=0${\displaystyle \frac{d^{2}w(0)}{d{x}^2}=0}$, which then translates to:(115)w-1-2w0+w1

(𝛥x)2=0⟹w-1=-w1$$\frac{w_{-1}-2w_0+w_1}{(𝛥x{)}^2}=0⟹w_{-1}=-w_1$$ c) Free-end Boundary• No bending moment means that dw2(0)

dx2=0${\displaystyle \frac{d{w}^2(0)}{d{x}^2}=0}$, which then translates into:(116)w-1-2w0+w1

(𝛥x)2=0⟹w-1=-w1+2w0 $$\begin{array}{c}\frac{w_{-1}-2w_0+w_1}{(𝛥x{)}^2}=0⟹w_{-1}=-w_1+2w_0\\ \end{array}$$ • No shear force means that d3w(0)

dx3=0${\displaystyle \frac{d^{3}w(0)}{d{x}^3}=0}$, which then translate into:(117)-w-2+2w-1-2w1+w2

2(𝛥x)3=0⟹w-2=2w-1-2w1+w2$$\frac{-w_{-2}+2w_{-1}-2w_1+w_2}{2(𝛥x{)}^3}=0⟹w_{-2}=2w_{-1}-2w_1+w_2$$Substituting (116) to (117) gives us:(118)a

w-2	=2w-1-2w1+w2
	=2(2w0-w1)-2w1+w2
	=4w0-2w1-2w1+w2
	=4w0-4w1+w2

$$\begin{array}{rl}{w}_{-2}& =2w_{-1}-2w_1+w_2\\ & =2(2w_0-w_1)-2w_1+w_2\\ & =4w_0-2w_1-2w_1+w_2\\ & =4w_0-4w_1+w_2\end{array}$$ 7.2.2. Boundary Conditions at the Right Enda) Fixed-end Boundary (Clamped)• No position displacement means that w(L)=0${\displaystyle w(L)=0}$, which then translates into:(119)wn=0$$w_n=0$$• Zero horizontal tangent means that dw(L)

dx=0${\displaystyle \frac{dw(L)}{dx}=0}$, which then translates into:(120)wn-1-wn+1

2𝛥x=0⟹wn+1=wn-1$$\frac{w_{n-1}-w_{n+1}}{2𝛥x}=0⟹w_{n+1}=w_{n-1}$$ b) Pinned-end Boundary • No position displacement means that w(L)=0${\displaystyle w(L)=0}$, which then translates into:(121)wn=0$$w_n=0$$• No bending moment means that d2w(L)

dx2=0${\displaystyle \frac{d^{2}w(L)}{d{x}^2}=0}$, which then translates into:(122)wn-1-2wn+wn+1

(𝛥x)2=0⟹wn+1=-wn-1$$\frac{w_{n-1}-2w_n+w_{n+1}}{(𝛥x{)}^2}=0⟹w_{n+1}=-w_{n-1}$$c) Free-end Boundary• No bending moment means that dw2(L)

dx2=0${\displaystyle \frac{d{w}^2(L)}{d{x}^2}=0}$, which translates into:(123)wn-1-2wn+wn+1

(𝛥x)2=0⟹wn+1=-wn-1+2wn $$\begin{array}{c}\frac{w_{n-1}-2w_n+w_{n+1}}{(𝛥x{)}^2}=0⟹w_{n+1}=-w_{n-1}+2w_n\\ \end{array}$$• No shear force means that d3w(L)

dx3=0${\displaystyle \frac{d^{3}w(L)}{d{x}^3}=0}$, which then translates into:(124)-wn-2+2wn-1-2wn+1+wn+2

2(𝛥x)3=0⟹wn+2=wn-2-2wn-1+2wn+1$$\frac{-w_{n-2}+2w_{n-1}-2w_{n+1}+w_{n+2}}{2(𝛥x{)}^3}=0⟹w_{n+2}=w_{n-2}-2w_{n-1}+2w_{n+1}$$By substituting (123) to (124), we can then acquire:(125)a

wn+2	=wn-2-2wn-1+2wn+1
	=wn-2-2wn-1+2(-wn-1+2wn)
	=wn-2-2wn-1-2wn-1+4wn
	=wn-2-4wn-1+4wn

$$\begin{array}{rl}{w}_{n+2}& =w_{n-2}-2w_{n-1}+2w_{n+1}\\ & =w_{n-2}-2w_{n-1}+2(-w_{n-1}+2w_n)\\ & =w_{n-2}-2w_{n-1}-2w_{n-1}+4w_n\\ & =w_{n-2}-4w_{n-1}+4w_n\end{array}$$ To summarize, we put (112) to (125) into a table as follows.

	Left	Right
Fixed-end	aL2=×L1=w1w0=0${\displaystyle \begin{array}{l}{L}_2=\times \\ L_1=w_1\\ w_0=0\end{array}}$	aR1=wn-1R2=×wn=0${\displaystyle \begin{array}{l}{R}_1=w_{n-1}\\ R_2=\times \\ w_n=0\end{array}}$
Pinned-end	aL2=×L1=-w1w0=0${\displaystyle \begin{array}{l}{L}_2=\times \\ L_1=-w_1\\ w_0=0\end{array}}$	R1=-wn-1R2=×wn=0$$\begin{array}{c}{R}_1=-w_{n-1}\\ R_2=\times \\ w_n=0\end{array}$$
Free-end	aL2=4w0-4w1+w2L1=-w1+2w0${\displaystyle \begin{array}{l}{L}_2=4w_0-4w_1+w_2\\ L_1=-w_1+2w_0\end{array}}$	aR1=-wn-1+2wnR2=wn-2-4wn-1+4wn${\displaystyle \begin{array}{l}{R}_1=-w_{n-1}+2w_n\\ R_2=w_{n-2}-4w_{n-1}+4w_n\end{array}}$

Table 1:All possible boundary conditions

Now, we can update the L1${\displaystyle L_1}$, L2${\displaystyle L_2}$. R1${\displaystyle R_1}$, and R2${\displaystyle R_2}$ in (111) depending on the current boundary conditions. 7.3. System Matrices for Different BoundariesIn order to solve (111) , we must first make sure that matrix A${\displaystyle A}$ and F${\displaystyle F}$ contain only numbers and not variables. Currently, L1${\displaystyle L_1}$, L2${\displaystyle L_2}$, R1${\displaystyle R_1}$, R2${\displaystyle R_2}$ still contain some variables (w${\displaystyle w}$). This is not a desirable format since we can not apply the linear solve command to solve Aw=F${\displaystyle Aw=F}$. ➔ Free-free

f(x)Free-free boundary condition is not going to work for the Euler Bernoulli equation since the matrix that we get is not a full rank matrix. Hence, its system equation can not be solved. However, this free-free case is a a good example on how we should adapt (111) to include some boundary conditions. By including the boundaries into (111) we get the following equation.a

[a 1 0 … 0 0 1 0 … 0 1 -4 6 -4 1 0 … 0 0 1 -4 6 -4 1 … 0 ⋮ ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ 0 … 0 1 -4 6 -4 1 0 … 0 1 0 0 … 0 1 ]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a w-2 w-1 w0 w1 ⋮ wn wn+1 wn+2 ]⏠⏣⏣⏣⏡⏣⏣⏣⏢W

=𝛥x4 EI[a 4w0-4w1+w2 -w1+2w0 f0 f1 ⋮ fn -wn-1+2wn wn-2-4wn-1+4wn ]⏠⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏢F

(127)$$\begin{array}{c}\begin{array}{rl}\underset{A}{\underbrace{\left[\begin{array}{cccccccc}1& 0& & & & & \dots & 0\\ 0& 1& 0& & & & \dots & 0\\ 1& -4& 6& -4& 1& 0& \dots & 0\\ 0& 1& -4& 6& -4& 1& \dots & 0\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \dots & 0& 1& -4& 6& -4& 1\\ 0& & & & \dots & 0& 1& 0\\ 0& & & & & \dots & 0& 1\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_{-2}\\ w_{-1}\\ w_0\\ w_1\\ ⋮\\ w_n\\ w_{n+1}\\ w_{n+2}\end{array}\right]}}& =\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\begin{array}{c}4w_0-4w_1+w_2\\ -w_1+2w_0\\ f_0\\ f_1\\ ⋮\\ f_n\\ -w_{n-1}+2w_n\\ w_{n-2}-4w_{n-1}+4w_n\end{array}\right]}}\end{array}\\ \end{array}$$Starting from the left end, we write down the equations made by Row 1 to Row 4:(128)a

Row 1:	w-2=4w0-4w1+w2
Row 2 :	w-1=-w1+2w0
Row 3 :	w-2-4w-1+6w0-4w1+w2=f0
Row 4 :	w-1-4w0+6w1-4w2+w3=f0

$$\begin{array}{rl}\text{Row 1}:& w_{-2}=4w_0-4w_1+w_2\\ \text{Row 2 : }& w_{-1}=-w_1+2w_0\\ \text{Row 3 : }& w_{-2}-4w_{-1}+6w_0-4w_1+w_2=f_0\\ \text{Row 4 : }& w_{-1}-4w_0+6w_1-4w_2+w_3=f_0\end{array}$$By substituting Row 1 and Row 2 to Row 3, we get the following:4w0-4w1+w2-4(-w1+2w0)+6w0-4w1+w2=fo(129)2w0-4w1+2w2=f0$$\begin{array}{c}4w_0-4w_1+w_2-4(-w_1+2w_0)+6w_0-4w_1+w_2=f_o\\ 2w_0-4w_1+2w_2=f_0\end{array}$$Next, by substituting Row 1 and Row 2 to Row 4, we get the following:-w1+2w0-4wo+6w1-4w2+w3=f1(130)-2w0+5w1-4w2+w3=f1$$\begin{array}{c}-w_1+2w_0-4w_o+6w_1-4w_2+w_3=f_1\\ -2w_0+5w_1-4w_2+w_3=f_1\end{array}$$As for the right end, we can do the similar steps. By substituting Row n+1${\displaystyle n+1}$ and Row n+2${\displaystyle n+2}$ to Row n${\displaystyle n}$, we get the following:wn-2-4wn-1+6wn-4(-wn-1+2wn)+wn-2-4wn-1+4wn=fn(131)2wn-2-4wn-1+2wn=fn$$\begin{array}{c}{w}_{n-2}-4w_{n-1}+6w_n-4(-w_{n-1}+2w_n)+w_{n-2}-4w_{n-1}+4w_n=f_n\\ 2w_{n-2}-4w_{n-1}+2w_n=f_n\end{array}$$Next, by substituting Row n+1${\displaystyle n+1}$ and Row n+2${\displaystyle n+2}$ to Row n-1${\displaystyle n-1}$, we get the following:wn-3-4wn-2+6wn-1-4wn-wn-1+2wn=fn-1(132)wn-3-4wn-2+5wn-1-2wn=fn-1$$\begin{array}{c}{w}_{n-3}-4w_{n-2}+6w_{n-1}-4w_n-w_{n-1}+2w_n=f_{n-1}\\ w_{n-3}-4w_{n-2}+5w_{n-1}-2w_n=f_{n-1}\end{array}$$Now, we can remove Row 1, Row 2, Row n+1${\displaystyle n+1}$, and Row n+2${\displaystyle n+2}$. These rows are all made by the phantom nodes and must first be removed. As a result, we can update (127) into:(133)[a

2	-4	2	0				…	0
-2	5	-4	1	0			…	0
1	-4	6	-4	1	0		…	0
0	1	-4	6	-4	1	0	…	0
⋮	⋱	⋱	⋱	⋱	⋱	⋱		⋮
0		…	0	1	-4	6	-4	1
0			…	0	1	-4	5	-2
0				…	0	2	-4	2

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w0

w1

w2

⋮

wn-2

wn-1

wn

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W=𝛥x4

EI[a

f0

⋮

fn

]⏠⏣⏣⏣⏡⏣⏣⏣⏢F$$\underset{A}{\underbrace{\left[\begin{array}{ccccccccc}2& -4& 2& 0& & & & \dots & 0\\ -2& 5& -4& 1& 0& & & \dots & 0\\ 1& -4& 6& -4& 1& 0& & \dots & 0\\ 0& 1& -4& 6& -4& 1& 0& \dots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ 0& & \dots & 0& 1& -4& 6& -4& 1\\ 0& & & \dots & 0& 1& -4& 5& -2\\ 0& & & & \dots & 0& 2& -4& 2\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_0\\ w_1\\ w_2\\ ⋮\\ w_{n-2}\\ w_{n-1}\\ w_n\end{array}\right]}}=\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\begin{array}{c}{f}_0\\ ⋮\\ f_n\end{array}\right]}}$$ 7.3.1. Fixed-fixed (fixed ended)

f(x)For a case where both ends are fixed, we use the following equations.(134)a

7	-4	1	0				…	0
-4	6	-4	1	0			…	0
1	-4	6	-4	1	0		…	0
0	1	-4	6	-4	1	0	…	0
		⋱	⋱	⋱	⋱	⋱		⋮
⋮			⋱	⋱	⋱	⋱	⋱
0			…	0	1	-4	6	-4
0				…	0	1	-4	7

⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w1

w2

⋮

wn-2

wn-1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W=𝛥x4

EI[a

f1

⋮

fn-1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢F$$\underset{A}{\underbrace{\left[\begin{array}{ccccccccc}7& -4& 1& 0& & & & \dots & 0\\ -4& 6& -4& 1& 0& & & \dots & 0\\ 1& -4& 6& -4& 1& 0& & \dots & 0\\ 0& 1& -4& 6& -4& 1& 0& \dots & 0\\ & & \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & & \ddots & \ddots & \ddots & \ddots & \ddots & \\ 0& & & \dots & 0& 1& -4& 6& -4\\ 0& & & & \dots & 0& 1& -4& 7\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_{n-2}\\ w_{n-1}\end{array}\right]}}=\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\begin{array}{c}{f}_1\\ ⋮\\ f_{n-1}\end{array}\right]}}$$ 7.3.2. Pinned-pinned (simply supported)

f(x)For a case where both ends are pinned, we use the following equations.(135)a

5	-4	1	0				…	0
-4	6	-4	1	0			…	0
1	-4	6	-4	1	0		…	0
0	1	-4	6	-4	1	0	…	0
		⋱	⋱	⋱	⋱	⋱		⋮
⋮			⋱	⋱	⋱	⋱	⋱
0			…	0	1	-4	6	-4
0				…	0	1	-4	5

⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w1

w2

⋮

wn-2

wn-1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W=𝛥x4

EI[a

f1

⋮

fn-1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢F$$\underset{A}{\underbrace{\left[\begin{array}{ccccccccc}5& -4& 1& 0& & & & \dots & 0\\ -4& 6& -4& 1& 0& & & \dots & 0\\ 1& -4& 6& -4& 1& 0& & \dots & 0\\ 0& 1& -4& 6& -4& 1& 0& \dots & 0\\ & & \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & & \ddots & \ddots & \ddots & \ddots & \ddots & \\ 0& & & \dots & 0& 1& -4& 6& -4\\ 0& & & & \dots & 0& 1& -4& 5\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_{n-2}\\ w_{n-1}\end{array}\right]}}=\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\overline{\begin{array}{c}{f}_1\\ ⋮\\ f_{n-1}\end{array}}\right]}}$$ 7.3.3. Fixed-pinned (cantilever, simply supported)

f(x)For a case where the left end is fixed while the right end is pinned, we use the following equations.(136)a

7	-4	1	0				…	0
-4	6	-4	1	0			…	0
1	-4	6	-4	1	0		…	0
0	1	-4	6	-4	1	0	…	0
		⋱	⋱	⋱	⋱	⋱		⋮
⋮			⋱	⋱	⋱	⋱	⋱
0			…	0	1	-4	6	-4
0				…	0	1	-4	5

⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w1

w2

⋮

wn-2

wn-1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W=𝛥x4

EI[a

f1

⋮

fn-1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢F$$\underset{A}{\underbrace{\left[\begin{array}{ccccccccc}7& -4& 1& 0& & & & \dots & 0\\ -4& 6& -4& 1& 0& & & \dots & 0\\ 1& -4& 6& -4& 1& 0& & \dots & 0\\ 0& 1& -4& 6& -4& 1& 0& \dots & 0\\ & & \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & & \ddots & \ddots & \ddots & \ddots & \ddots & \\ 0& & & \dots & 0& 1& -4& 6& -4\\ 0& & & & \dots & 0& 1& -4& 5\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_{n-2}\\ w_{n-1}\end{array}\right]}}=\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\overline{\begin{array}{c}{f}_1\\ ⋮\\ f_{n-1}\end{array}}\right]}}$$ 7.3.4. Pinned-fixed (simply supported, cantilever)

f(x)For a case where the left end is pinned while the right end is fixed, we use the following equations. (137)a

5	-4	1	0				…	0
-4	6	-4	1	0			…	0
1	-4	6	-4	1	0		…	0
0	1	-4	6	-4	1	0	…	0
		⋱	⋱	⋱	⋱	⋱		⋮
⋮			⋱	⋱	⋱	⋱	⋱
0			…	0	1	-4	6	-4
0				…	0	1	-4	7

⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w1

w2

⋮

wn-2

wn-1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W=𝛥x4

EI[a

f1

⋮

fn-1

]⏠⏣⏣⏣⏡⏣⏣⏣⏢F$$\underset{A}{\underbrace{\left[\begin{array}{ccccccccc}5& -4& 1& 0& & & & \dots & 0\\ -4& 6& -4& 1& 0& & & \dots & 0\\ 1& -4& 6& -4& 1& 0& & \dots & 0\\ 0& 1& -4& 6& -4& 1& 0& \dots & 0\\ & & \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ ⋮& & & \ddots & \ddots & \ddots & \ddots & \ddots & \\ 0& & & \dots & 0& 1& -4& 6& -4\\ 0& & & & \dots & 0& 1& -4& 7\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_{n-2}\\ w_{n-1}\end{array}\right]}}=\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\overline{\begin{array}{c}{f}_1\\ ⋮\\ f_{n-1}\end{array}}\right]}}$$ 7.3.5. Fixed-free (cantilever)

f(x)For a case where the left end is fixed while the right end is free, we use the following equations.(138)[a

7	-4	1	0				…	0
-4	6	-4	1	0			…	0
1	-4	6	-4	1	0		…	0
0	1	-4	6	-4	1	0	…	0
	⋱	⋱	⋱	⋱	⋱	⋱		⋮
	…	0	1	-4	6	-4	1	0
		…	0	1	-4	6	-4	1
⋮			…	0	1	-4	5	-2
0				…	0	2	-4	2

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w1

w2

⋮

wn-1

wn

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W=𝛥x4

EI[a

f1

⋮

fn

]⏠⏣⏣⏣⏡⏣⏣⏣⏢F$$\underset{A}{\underbrace{\left[\begin{array}{ccccccccc}7& -4& 1& 0& & & & \dots & 0\\ -4& 6& -4& 1& 0& & & \dots & 0\\ 1& -4& 6& -4& 1& 0& & \dots & 0\\ 0& 1& -4& 6& -4& 1& 0& \dots & 0\\ & \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ & \dots & 0& 1& -4& 6& -4& 1& 0\\ & & \dots & 0& 1& -4& 6& -4& 1\\ ⋮& & & \dots & 0& 1& -4& 5& -2\\ 0& & & & \dots & 0& 2& -4& 2\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_{n-1}\\ w_n\end{array}\right]}}=\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\begin{array}{c}{f}_1\\ ⋮\\ f_n\end{array}\right]}}$$7.3.6. Fixed-explicit moment

MAn explicit moment in a beam is another type of boundary condition. Although it may look an applied input to the beam, it is not. Euler-Bernoulli equation only takes force/pressure as its input. It does not accept torque. A pure torque M${\displaystyle M}$ that causes bending is applied at the right end:• A manually-defined constant bending moment of M${\displaystyle M}$ or dw2(L)

dx2=M${\displaystyle \frac{d{w}^2(L)}{d{x}^2}=M}$. This then translates into:(139)wn-1-2wn+wn+1

𝛥x2=M⟹wn+1=M𝛥x2-wn-1+2wn $$\begin{array}{c}\frac{w_{n-1}-2w_n+w_{n+1}}{𝛥x^2}=M⟹w_{n+1}=M𝛥x^2-w_{n-1}+2w_n\\ \end{array}$$• No shear force means that d3w(L)

dx3=0${\displaystyle \frac{d^{3}w(L)}{d{x}^3}=0}$, which then translates into:(140)-wn-2+2wn-1-2wn+1+wn+2

2 𝛥x3=0⟹wn+2=wn-2-2wn-1+2wn+1$$\frac{-w_{n-2}+2w_{n-1}-2w_{n+1}+w_{n+2}}{2𝛥x^3}=0⟹w_{n+2}=w_{n-2}-2w_{n-1}+2w_{n+1}$$ By substituting to (139), to (141) we can then acquire:(141)a

wn+2	=wn-2-2wn-1+2wn+1
	=wn-2-2wn-1+2(M𝛥x2-wn-1+2wn)
	=wn-2-2wn-1-2wn-1+4wn+2M𝛥x2
	=wn-2-4wn-1+4wn

$$\begin{array}{rl}{w}_{n+2}& =w_{n-2}-2w_{n-1}+2w_{n+1}\\ & =w_{n-2}-2w_{n-1}+2(M𝛥x^2-w_{n-1}+2w_n)\\ & =w_{n-2}-2w_{n-1}-2w_{n-1}+4w_n+2M𝛥x^2\\ & =w_{n-2}-4w_{n-1}+4w_n\end{array}$$We can summarize the boundary as follows:

Left Boundary(Fixed-end)	Right Boundary(Applied pure torque)
aL1=×L2=w1w0=0${\displaystyle \begin{array}{l}{L}_1=\times \\ L_2=w_1\\ w_0=0\end{array}}$	aR1=M𝛥x2-wn-1+2wnR2=wn-2-4wn-1+4wn${\displaystyle \begin{array}{l}{R}_1=M𝛥x^2-w_{n-1}+2w_n\\ R_2=w_{n-2}-4w_{n-1}+4w_n\end{array}}$

Applying these boundaries to the general system matrix gives us:(142)[a

1	0					…	0
0	1	0				…	0
1	-4	6	-4	1	0	…	0
0	1	-4	6	-4	1	…	0
⋮		⋱	⋱	⋱	⋱	⋱	⋮
0	…	0	1	-4	6	-4	1
0				…	0	1	0
0					…	0	1

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w-2

w-1

w0=0

w1

⋮

wn

wn+1

wn+2

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W=𝛥x4

EI[a

×

w1

0

0

⋮

fn

M𝛥x2-wn-1+2wn

2M𝛥x2+wn-2-4wn-1+4wn

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢F$$\underset{A}{\underbrace{\left[\begin{array}{cccccccc}1& 0& & & & & \dots & 0\\ 0& 1& 0& & & & \dots & 0\\ 1& -4& 6& -4& 1& 0& \dots & 0\\ 0& 1& -4& 6& -4& 1& \dots & 0\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \dots & 0& 1& -4& 6& -4& 1\\ 0& & & & \dots & 0& 1& 0\\ 0& & & & & \dots & 0& 1\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_{-2}\\ w_{-1}\\ w_0=0\\ w_1\\ ⋮\\ w_n\\ w_{n+1}\\ w_{n+2}\end{array}\right]}}=\underset{F}{\underbrace{\frac{𝛥x^4}{EI}\left[\begin{array}{c}\times \\ w_1\\ 0\\ 0\\ ⋮\\ f_n\\ M𝛥x^2-w_{n-1}+2w_n\\ 2M𝛥x^2+w_{n-2}-4w_{n-1}+4w_n\end{array}\right]}}$$Since we have already addressed the left boundary, hence, we only need to address the right boundary.From the second-to-the last row:wn-2-4wn-1+6wn-4wn+1+wn+2=fnwn-2-4wn-1+6wn-4(M𝛥x2-wn-1+2wn)+(2M𝛥x2+wn-2-4wn-1+4wn)=fnwn-2-4wn-1+6wn-2M𝛥x2+4wn-1-8wn+wn-2-4wn-1+4wn=fn2wn-2-4wn-1+2wn=fn+2M𝛥x2$$\begin{array}{c}{w}_{n-2}-4w_{n-1}+6w_n-4w_{n+1}+w_{n+2}=f_n\\ w_{n-2}-4w_{n-1}+6w_n-4(M𝛥x^2-w_{n-1}+2w_n)+(2M𝛥x^2+w_{n-2}-4w_{n-1}+4w_n)=f_n\\ w_{n-2}-4w_{n-1}+6w_n-2M𝛥x^2+4w_{n-1}-8w_n+w_{n-2}-4w_{n-1}+4w_n=f_n\\ 2w_{n-2}-4w_{n-1}+2w_n=f_n+2M𝛥x^2\end{array}$$ From the third-to-the last row:wn-3-4wn-2+6wn-1-4wn+M𝛥x2-wn-1+2wn=fn-1(143)wn-3-4wn-2+5wn-1-2wn=fn-1-M𝛥x2$$\begin{array}{c}{w}_{n-3}-4w_{n-2}+6w_{n-1}-4w_n+M𝛥x^2-w_{n-1}+2w_n=f_{n-1}\\ w_{n-3}-4w_{n-2}+5w_{n-1}-2w_n=f_{n-1}-M𝛥x^2\end{array}$$(144)a

7	-4	1	0			…	0
-4	6	-4	1	0		…	0
1	-4	6	-4	1	0	…	0
	⋱	⋱	⋱	⋱
		⋱	⋱	⋱	⋱
0	…	0	1	-4	6	-4	1
0			…	0	1	-4	7

⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢A$$\underset{A}{\underbrace{\left[\begin{array}{cccccccc}7& -4& 1& 0& & & \dots & 0\\ -4& 6& -4& 1& 0& & \dots & 0\\ 1& -4& 6& -4& 1& 0& \dots & 0\\ & \ddots & \ddots & \ddots & \ddots & & & \\ & & \ddots & \ddots & \ddots & \ddots & & \\ 0& \dots & 0& 1& -4& 6& -4& 1\\ 0& & & \dots & 0& 1& -4& 7\end{array}\right]}}$$Thus, (142) cab be rewritten as:(145)[a

7	-4	1	0				…	0
-4	6	-4	1	0			…	0
1	-4	6	-4	1	0		…	0
0	1	-4	6	-4	1	0	…	0
⋮	⋱	⋱	⋱	⋱	⋱	⋱		⋮
0		…	0	1	-4	6	-4	1
0		…	0	1	-4	5	1	-2
0				…	0	2	-4	2

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w1

w2

w3

w4

⋮

wn-2

wn-1

wn

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W=1

EI[a

0
0
0
0
⋮
0
-M𝛥x2
2M𝛥x2

]⏠⏣⏣⏣⏣⏡⏣⏣⏣⏣⏢F$$\underset{A}{\underbrace{\left[\begin{array}{ccccccccc}7& -4& 1& 0& & & & \dots & 0\\ -4& 6& -4& 1& 0& & & \dots & 0\\ 1& -4& 6& -4& 1& 0& & \dots & 0\\ 0& 1& -4& 6& -4& 1& 0& \dots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & \ddots & \ddots & & ⋮\\ 0& & \dots & 0& 1& -4& 6& -4& 1\\ 0& & \dots & 0& 1& -4& 5& 1& -2\\ 0& & & & \dots & 0& 2& -4& 2\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_1\\ w_2\\ w_3\\ w_4\\ ⋮\\ w_{n-2}\\ w_{n-1}\\ w_n\end{array}\right]}}=\underset{F}{\underbrace{\frac{1}{EI}\left[\begin{array}{c}0\\ 0\\ 0\\ 0\\ ⋮\\ 0\\ -M𝛥x^2\\ 2M𝛥x^2\end{array}\right]}}$$We can then fine W${\displaystyle W}$ by using the linear solve function in MATLAB. 8. Unsteady Simple Euler-Bernoulli Beam Similar to the previous section, in this section, we still discuss a simple cantilever problem. The dynamic model of an Euler-Bernoulli beam is governed by the following equation:(146)𝜇𝜕2w

𝜕t2(x,t)=-EI 𝜕4w

𝜕x4(x,t)+f(x,t)$$𝜇\frac{{𝜕}^{2}w}{𝜕t^2}(x,t)=-EI\frac{{𝜕}^{4}w}{𝜕x^4}(x,t)+f(x,t)$$where 𝜇${\displaystyle 𝜇}$ is the mass per unit length, E${\displaystyle E}$ is the Young's modulus, and I${\displaystyle I}$ is the area moment of inertia. ${\displaystyle }$(147)a

𝜕2w 𝜕t2(x,t)	=-EI 𝜇n∑i=0wi-2,j-4wi-1,j+6wi,j-4wi+1,j+wi+2,j (𝛥x)4+ fi,j
𝜕2w 𝜕t2(x,t)	=-EI 𝜇𝛥x4[a 1 0 … 0 0 1 0 … 0 1 -4 6 -4 1 0 … 0 0 1 -4 6 -4 1 … 0 ⋮ ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ 0 … 0 1 -4 6 -4 1 0 … 0 1 0 0 … 0 1 ]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a w-2 w-1 w0 w1 ⋮ wn wn+1 wn+2 ]⏠⏣⏣⏣⏡⏣⏣⏣⏢W+[a L2 L1 f0 f1 ⋮ fn R1 R2 ]⏠⏣⏣⏣⏡⏣⏣⏣⏢F

$$\begin{array}{rl}\frac{{𝜕}^{2}w}{𝜕t^2}(x,t)& =-\frac{EI}{𝜇}\sum _{i=0}^n\frac{w_{i-2,j}-4w_{i-1,j}+6w_{\mathrm{i},j}-4w_{i+1,j}+w_{i+2,j}}{(𝛥x{)}^4}+f_{i,j}\\ \frac{{𝜕}^{2}w}{𝜕t^2}(x,t)& =-\frac{EI}{𝜇𝛥x^4}\underset{A}{\underbrace{\left[\begin{array}{cccccccc}1& 0& & & & & \dots & 0\\ 0& 1& 0& & & & \dots & 0\\ 1& -4& 6& -4& 1& 0& \dots & 0\\ 0& 1& -4& 6& -4& 1& \dots & 0\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \dots & 0& 1& -4& 6& -4& 1\\ 0& & & & \dots & 0& 1& 0\\ 0& & & & & \dots & 0& 1\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_{-2}\\ w_{-1}\\ w_0\\ w_1\\ ⋮\\ w_n\\ w_{n+1}\\ w_{n+2}\end{array}\right]}}+\underset{F}{\underbrace{\left[\begin{array}{c}{L}_2\\ L_1\\ f_0\\ f_1\\ ⋮\\ f_n\\ R_1\\ R_2\end{array}\right]}}\end{array}$$Discretizing 𝜕2w

𝜕t2(x,t)${\displaystyle \frac{{𝜕}^{2}w}{𝜕t^2}(x,t)}$ in (147) gives us: w(x,t-𝛥t)-2w(x,t)+w(x,t+𝛥t)

𝛥t2=-EI

𝜇𝛥x4[a

1	0					…	0
0	1	0				…	0
1	-4	6	-4	1	0	…	0
0	1	-4	6	-4	1	…	0
⋮		⋱	⋱	⋱	⋱	⋱	⋮
0	…	0	1	-4	6	-4	1
0				…	0	1	0
0					…	0	1

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w-2

w-1

w0

w1

⋮

wn

wn+1

wn+2

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W+[a

L2

L1

f0

f1

⋮

fn

R1

R2

]⏠⏣⏣⏣⏡⏣⏣⏣⏢F$$\begin{array}{c}\frac{w(x,t-𝛥t)-2w(x,t)+w(x,t+𝛥t)}{𝛥t^2}=\\ -\frac{EI}{𝜇𝛥x^4}\underset{A}{\underbrace{\left[\begin{array}{cccccccc}1& 0& & & & & \dots & 0\\ 0& 1& 0& & & & \dots & 0\\ 1& -4& 6& -4& 1& 0& \dots & 0\\ 0& 1& -4& 6& -4& 1& \dots & 0\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \dots & 0& 1& -4& 6& -4& 1\\ 0& & & & \dots & 0& 1& 0\\ 0& & & & & \dots & 0& 1\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_{-2}\\ w_{-1}\\ w_0\\ w_1\\ ⋮\\ w_n\\ w_{n+1}\\ w_{n+2}\end{array}\right]}}+\underset{F}{\underbrace{\left[\begin{array}{c}{L}_2\\ L_1\\ f_0\\ f_1\\ ⋮\\ f_n\\ R_1\\ R_2\end{array}\right]}}\end{array}$$ (148)w(x,t+𝛥t)=-EI𝛥t2

𝜇𝛥x4[a

1	0					…	0
0	1	0				…	0
1	-4	6	-4	1	0	…	0
0	1	-4	6	-4	1	…	0
⋮		⋱	⋱	⋱	⋱	⋱	⋮
0	…	0	1	-4	6	-4	1
0				…	0	1	0
0					…	0	1

]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a

w-2

w-1

w0

w1

⋮

wn

wn+1

wn+2

]⏠⏣⏣⏣⏡⏣⏣⏣⏢W+𝛥t2[a

L2

L1

f0

f1

⋮

fn

R1

R2

]⏠⏣⏣⏣⏡⏣⏣⏣⏢F+2w(x,t)-w(x,t-𝛥t)$$\begin{array}{c}w(x,t+𝛥t)=-\frac{EI𝛥t^2}{𝜇𝛥x^4}\underset{A}{\underbrace{\left[\begin{array}{cccccccc}1& 0& & & & & \dots & 0\\ 0& 1& 0& & & & \dots & 0\\ 1& -4& 6& -4& 1& 0& \dots & 0\\ 0& 1& -4& 6& -4& 1& \dots & 0\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \dots & 0& 1& -4& 6& -4& 1\\ 0& & & & \dots & 0& 1& 0\\ 0& & & & & \dots & 0& 1\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_{-2}\\ w_{-1}\\ w_0\\ w_1\\ ⋮\\ w_n\\ w_{n+1}\\ w_{n+2}\end{array}\right]}}+𝛥t^2\underset{F}{\underbrace{\left[\begin{array}{c}{L}_2\\ L_1\\ f_0\\ f_1\\ ⋮\\ f_n\\ R_1\\ R_2\end{array}\right]}}+\\ 2w(x,t)-w(x,t-𝛥t)\end{array}$$ Another way is to use a MATLAB built-in ODE solver. In order to do so, let us take v=𝜕w

𝜕t${\displaystyle v=\frac{𝜕w}{𝜕t}}$ and 𝛼=EI

𝜇𝛥x4${\displaystyle 𝛼=\frac{EI}{𝜇𝛥x^4}}$, hence:(149)a

𝜕v 𝜕t(x,t)

=-EI 𝜇𝛥x4[a 1 0 … 0 0 1 0 … 0 1 -4 6 -4 1 0 … 0 0 1 -4 6 -4 1 … 0 ⋮ ⋱ ⋱ ⋱ ⋱ ⋱ ⋮ 0 … 0 1 -4 6 -4 1 0 … 0 1 0 0 … 0 1 ]⏠⏣⏣⏣⏣⏣⏣⏣⏣⏡⏣⏣⏣⏣⏣⏣⏣⏣⏢A[a w-2 w-1 w0 w1 ⋮ wn wn+1 wn+2 ]⏠⏣⏣⏣⏡⏣⏣⏣⏢W+[a L2 L1 f0 f1 ⋮ fn R1 R2 ]⏠⏣⏣⏣⏡⏣⏣⏣⏢F

$$\begin{array}{rl}\frac{𝜕v}{𝜕t}(x,t)& =-\frac{EI}{𝜇𝛥x^4}\underset{A}{\underbrace{\left[\begin{array}{cccccccc}1& 0& & & & & \dots & 0\\ 0& 1& 0& & & & \dots & 0\\ 1& -4& 6& -4& 1& 0& \dots & 0\\ 0& 1& -4& 6& -4& 1& \dots & 0\\ ⋮& & \ddots & \ddots & \ddots & \ddots & \ddots & ⋮\\ 0& \dots & 0& 1& -4& 6& -4& 1\\ 0& & & & \dots & 0& 1& 0\\ 0& & & & & \dots & 0& 1\end{array}\right]}}\underset{W}{\underbrace{\left[\begin{array}{c}{w}_{-2}\\ w_{-1}\\ w_0\\ w_1\\ ⋮\\ w_n\\ w_{n+1}\\ w_{n+2}\end{array}\right]}}+\underset{F}{\underbrace{\left[\begin{array}{c}{L}_2\\ L_1\\ f_0\\ f_1\\ ⋮\\ f_n\\ R_1\\ R_2\end{array}\right]}}\end{array}$$Now, we can write:(150)a

a 𝜕w 𝜕t 𝜕v 𝜕t

=a 0 I -𝛼A 0 a W V +a 0 F

$$\begin{array}{rl}\left[\begin{array}{c}\frac{𝜕w}{𝜕t}\\ \frac{𝜕v}{𝜕t}\end{array}\right]& =\left[\begin{array}{cc}0& I\\ -𝛼A& 0\end{array}\right]\left[\begin{array}{c}W\\ V\end{array}\right]+\left[\begin{array}{c}0\\ F\end{array}\right]\end{array}$$where V=[v-2 … vn+2]${\displaystyle V=[v_{-2}\dots v_{n+2}]}$. and I${\displaystyle I}$ is an identity matrix. Now, we can use (150) for MATLAB a built-in ODE solver.

1https://www.uni-muenster.de/Physik.TP/en/teaching/courses/numerical_methods_for_complex_systems_i_ws2016-2017.html

2Lecture in Mathematics ETH Zurich: Numerical Methods for Conservations Laws by Randall J. LeVeque

3Taken from: Process Dynamics, Modeling, and Control, Ogunnaike et. al., page 106

4PDE Observer Design for Counter-Current Heat Flows in a Heat-Exchange by F. Zobiri, et. al., published at IFAC-PapersOnLine

5https://www.maplesoft.com/applications/view.aspx?SID=119402&view=html

6Summarized from the Greenberg's Book

7Summarized from the Greenberg's Book

8https://www.comsol.com/model/out-of-plane-heat-transfer-for-a-thin-plate-493

9https://www.youtube.com/watch?v=63Ik8uZMxCQ&ab_channel=wdflannery