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) 1.3. Downwind Method (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-momentum Equation 6.2. Discretizing the 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. ut+cux=f(x,t)with the following initial condition:u(x,0)=u0(x)where c is a positive real number. Here, we want to to find the solution to (1). When f(x,t)=0, its analytical solution is given by: u(x,t)=u0(x-ct). Simply put, u0 propagates along x with a constant speed of 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, we use the following approximation.utu(x,t+𝛥t)-u(x,t)𝛥tWhile to approach ux, we use the following approximation.uxu(x+𝛥x,t)-u(x-𝛥x,t)2𝛥xNext, substituting (3) and (4) back to (1) gives usut-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𝛥t2𝛥x[u(x+𝛥x,t)-u(x-𝛥x,t)]+f(x,t)𝛥tu(x,t+𝛥t)u(x,t)-c𝛥t2𝛥x[u(x+𝛥x,t)-u(x-𝛥x,t)]+f(x,t)𝛥t 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) By using finite difference method, we can approximate both ut and ux. To approximate ut, we use the following approximation.utu(x,t+𝛥t)-u(x,t)𝛥tWhile to approach ux, we use the following approximation.uxu(x,t)-u(x-𝛥x,t)𝛥xBy substituting (5) and (6) to (1), we then have: 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𝛥x[u(x,t)-u(x-𝛥x,t)]+f(x,t)u(x,t+𝛥t)u(x,t)-c𝛥t𝛥x[u(x,t)-u(x-𝛥x,t)]+f(x,t)𝛥tLet us consider c𝛥t𝛥x=r, Uj,k=u(xj,tk), Fj,k=f(xj,tk), xj=j𝛥x and tk=k𝛥t, we then have the following shorter terms:Uj,k+1=Uj,k-r(Uj,k-Uj-1,k)+Fj,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=0u(x,0)=e-(x-2)20x10The implementation is provided in the listing below. The length of the rod is divided into 100 segments (dx=0.1) and the time is divided into 400 segments (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. 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
1.3. Downwind Method (c<0) By using finite difference method, we can approximate both ut and ux. To approximate ut, we use the following approximation.utu(x,t+𝛥t)-u(x,t)𝛥tWhile to approach ux, we use the following approximation.uxu(x+𝛥x,t)-u(x,t)𝛥xSubstituting (9) and (10) back to (1) yields the following.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𝛥x[u(x+𝛥x,t)-u(x,t)]+f(x,t)𝛥tu(x,t+𝛥t)u(x.t)-c𝛥t𝛥x[u(x+𝛥x,t)-u(x,t)]+f(x,t)𝛥tLet us consider c𝛥t𝛥x=r, Uj,k=u(xj,tk), Fj,k=f(xj,tk), xj=j𝛥x and tk=k𝛥t, we can then have the following shorter terms.Uj,k+1=Uj,k-r(Uj+1,k-Uj,k)+Fj,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=0u(x,0)=e-(x-8)20x10The implementation is provided in the listing below. The length of the rod is divided into 100 segments (dx=0.1) and the time is divided into 400 segments (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. 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
1.4. The Lax MethodThe Lax method is based on the FTCS method. utu(x,t+𝛥t)-u(x,t)𝛥tHowever, u(x,t) is replaced by an average of its two neighbors.u(x,t)=u(x+𝛥x,t)+u(x-𝛥x,t)2 Thus, by substituting (14) to (13), the approximation of ut can be rewritten asutu(x,t+𝛥t)-[u(x+𝛥x,t)+u(x-𝛥x,t)2]𝛥tOn the other hand, ux can be approximated asuxu(x+𝛥x,t)-u(x-𝛥x,t)2𝛥xSubstituting (15) and (16) back to (1) yieldsut-cux+f(x,t)u(x,t+𝛥t)-[u(x+𝛥x,t)+u(x-𝛥x,t)2]𝛥t-cu(x+𝛥x,t)-u(x-𝛥x,t)2𝛥x+f(x,t)u(x,t+𝛥t)-[u(x+𝛥x,t)+u(x-𝛥x,t)2]-c𝛥t2𝛥x[u(x+𝛥x,t)-u(x-𝛥x,t)]+f(x,t)𝛥tu(x,t+𝛥t)[u(x+𝛥x,t)+u(x-𝛥x,t)2]-c𝛥t2𝛥x[u(x+𝛥x,t)-u(x-𝛥x,t)]+f(x,t)𝛥tNow, let us consider r=c𝛥t𝛥x, Uj,k=u(xj,tk), Fj,k=f(xj,tk), xj=j𝛥x, tk=k𝛥t, Uj,k+1=u(xj,tk+𝛥t) and Uj-1,k=u(xj-𝛥x,tk), we can then have the following shorter terms.Uj,k+1=Uj+1,k+Uj-1,k2-r2(Uj+1,k-Uj,k)+Fj,k𝛥t The Lax method can handle both upwind and downwind direction. There is no need for us to check the sign of 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) and the time is divided into 400 segments (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
Figure 4:Lax method for ut+0.5uxx=0
1.5. The Lax-Wendroff Method2 We start from the Taylor series definition off(x+𝛥x)f(x+𝛥x)=f(x)+𝛥xf(x)+𝛥x22!f(x)+𝛥x33!f(x)+Thus, taking only up to the second order, we can writeu(x+𝛥x,t)=u(x,t)+𝛥xux+𝛥x22uxxRemember thatut=u(x,t+𝛥t)-u(x,t)𝛥tand ux=u(x+𝛥x,t)-u(x,t)𝛥xThus, ut=-cux+f(x,t) can be rewritten as u(x,t+𝛥t)-u(x,t)𝛥t=-cu(x+𝛥x,t)-u(x,t)𝛥x+f(x,t)By substituting (19) to (22)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)𝛥tu(x,t+𝛥t)u(x,t)-c𝛥t𝛥x[u(x,t)+𝛥xux+𝛥x22uxx-u(x,t)]+f(x,t)𝛥tu(x,t+𝛥t)u(x,t)-c𝛥tux-c𝛥t𝛥x2uxx+f(x,t)𝛥tSince c is the propagation speed of u(x,t)=u0(x-ct), thus, let us take 𝛥x=-c𝛥t. As the result, now we have the following equation.u(x,t+𝛥t)=u(x,t)-c𝛥tux+c2𝛥t22uxx+f(x,t)𝛥t Next, we need to approximate ux and uxx. For ux, we use central difference method as follows.uxu(x+𝛥x,t)-u(x-𝛥x,t)2𝛥xAs for uxx, we do not use the central difference method, instead, we take uxu(x+𝛥x,t)-u(x,t)𝛥x. Therefore, we can express uxx as follows.uxxux(x+𝛥x,t)-ux(x,t)𝛥xu(x+𝛥x,t)-u(x,t)𝛥x-u(x,t)-u(x-𝛥x,t)𝛥x𝛥xu(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)𝛥x2Now that we have (24) and (25), we can substitute them back to (23): u(x,t+𝛥t)u(x,t)-c𝛥tux+c2𝛥t22uxx+f(x,t)𝛥tu(x,t)-c𝛥t[u(x+𝛥x,t)-u(x-𝛥x,t)2𝛥x]+c2𝛥t22[u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)𝛥x2]+f(x,t)𝛥tu(x,t)-c𝛥t2𝛥x[u(x+𝛥x,t)-u(x-𝛥x,t)]+c2𝛥t22𝛥x2[u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)]+f(x,t)𝛥tNow let us consider r=c𝛥t𝛥x, Uj,k=u(xj,tk), Fj,k=f(xj,tk), xj=j𝛥x, tk=k𝛥t, Uj,k+1=u(xj,tk+𝛥t) and Uj-1,k=u(xj-𝛥x,tk), we can then have the following shorter terms.Uj,k+1=Uj,k-r2(Uj+1,k-Uj-1,k)+r22(Uj+1,k-2Uj,k+Uj-1,k)+Fj,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) and the time into 200 segments as well (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
Figure 6:Lax method for ut+0.5uxx=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.𝜕T(x,t)𝜕t+v𝜕T(x,t)𝜕zAdvection=-UAS𝜌ACp(T(x,t)-TS(x,t))Fourier's lawwhereCp is the liquid specific heat capacityA is cross sectional are of the tube where the liquid is flowingAS is the heat transfer area (blue surface)U is the overall heat transfer coefficient between the tube to its environment𝜌 is the density of the liquid
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. While source temperature of the cold flow is kept at TC0. Those two temperatures never change. This condition is called Dirichlet condition.
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. 𝜕𝜕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))or𝜕𝜕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))where: d1 and d2 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 and c2 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 with 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
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 and d2 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 and d4 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 and v2 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 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.𝜕𝜕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)=v32x2TW(x,t)+d3(TH(x,t)-TW(x,t))-d4(TW(x,t)-TC(x,t))Where5:
v1=0.6366197724v2=0.1657863990v3=0.00003585526318d1=0.9792d2=0.2986d3=2.3923d4=2.8708
The initial and boundary conditions are:TH(0,t)=360 and TH(x,0)=360TC(L,t)=300 and TC(x,0)=300TW(x,0)=330, TW(0,t)x=0, and TW(L,t)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.
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) in a one-dimensional rod (see Drawing 5) over a certain distance (x) and time (t) is given by (31). 𝛼 is the diffusivity coefficient. When g(x,t)=0, the system is homogenous.𝛼²uxx-ut=g(x,t)t0and0xLwhen at t=0, the following initial condition is applied.u(x,0)=f(x)
Drawing 5:Temperature distribution (u) in a very thin rod
2.1. Explicit MethodWith finite difference method, we can approach ut and uxx. To approximate ut, we use the following approximation.utu(x,t+𝛥t)-u(x,t)𝛥tuxu(x+𝛥x,t)-u(x,t)𝛥xWhile to approach ux, we use the following approximation.uxxux(x+𝛥x,t)-ux(x,t)𝛥xu(x+𝛥x,t)-u(x,t)𝛥x-u(x,t)-u(x-𝛥x,t)𝛥x𝛥xu(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)(𝛥x)2Substituting (32) and (33) back to (31) gives us the following equations.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)𝛥tu(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)𝛥tru(x+𝛥x,t)+(1-2r)u(x,t)+ru(x-𝛥x,t)+g(x,t)𝛥tu(x,t)+r[u(x+𝛥x,t)-2u(x,t)+u(x-𝛥x,t)]where r=𝛼2𝛥t(𝛥x)2. Now, let us consider Uj,k=u(xj,tk), Gj,k=g(xj,tk), xj=j𝛥x, tk=k𝛥t, Uj,k+1=u(xj,tk+𝛥t) and Uj-1,k=u(xj-𝛥x,tk), we can then have the following shorter terms.Uj,k+1=Uj,k+r[Uj+1,k-2Uj,k+Uj-1,k]+Gj,k𝛥t 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)=Bwhere A and 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)/xoru(x,t)/t) is kept constant. However, here we will only implement the derivative of the solution with respect to x.u(x,t)x=Cwhere C is a real number. At the boundary when x=0, we can write u(𝛥x,t)-u(-𝛥x,t)2𝛥x=Cthereforeu(-𝛥x,t)=u(𝛥x,t)-2𝛥xCAt another boundary when x=Lu(L+𝛥x,t)-u(L-𝛥x,t)2𝛥x=Dwhere D is a real number. Thereforeu(L+𝛥x,t)=u(L-𝛥x,t)+2𝛥xDThe implementation of Neumann boundary condition will be done by adding additional phantom nodes as shown in Drawing 6.
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:1.14uxx-ut=0,0x10with the following boundary and initial conditions:u(0,t)=u(10,t)=0u(x,0)=f(x)={60,0x50,5<x10MATLAB 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.
Drawing 7:Simulation result: Dirichlet-Dirichlet
2.1.5. Example 2Given the following heat transfer equation:uxx-ut=0,0x𝜋with the following boundary and initial conditions:u(0,t)=20ux(𝜋,t)=3u(x,0)=0MATLAB 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.
Drawing 8:Simulation result: Dirichlet-Neumann
2.1.6. Example 3Given the following heat equation:uxx-ut=0,0x1with the following initial and boundary conditions:ux(0,t)=0u(1,t)=100u(x,0)=f(x)={60,0x1/20,1/2<x1MATLAB 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.
Drawing 9:Simulation result: Neumann-Dirichlet
2.2. Implicit Method The main idea of an implicit method is to calculate uxx not only at the current time, but also at the time one step ahead. Let us consider:uxx(x,t)=(1-𝜃)uxx(x,t)+𝜃uxx(x,t+𝛥t),0𝜃1where 𝜃 acts as a weighting factor. When 𝜃=0, (47) becomes similar to (35). What we currently have in this case is an explicit method. On the other hand, when 𝜃=1, in this case, what we have is an implicit method. However, when 𝜃=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 𝜃=12 gives us:uxx=12ux(x+𝛥x,t)-ux(x,t)𝛥x+12ux(x+𝛥x,t+𝛥t)-ux(x,t+𝛥t)𝛥x=u(x+𝛥x,t)-u(x,t)𝛥x-u(x,t)-u(x-𝛥x,t)𝛥x2𝛥x+u(x+𝛥x,t+𝛥t)-u(x,t+𝛥t)𝛥x-u(x,t+𝛥t)-u(x-𝛥x,t+𝛥t)𝛥x2𝛥x=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𝛥x2which can be shortened asuxx=Uj+1,k-2Uj,k+Uj-1,k+Uj+1,k+1-2Uj,k+1+Uj-1,k+12𝛥x2Substituting (48) back to: ut=𝛼²uxx gives us:Uj,k+1-Uj,k𝛥tut=𝛼2Uj+1,k-2Uj,k+Uj-1,k+Uj+1,k+1-2Uj,k+1+Uj-1,k+12𝛥x2uxxUj,k+1-Uj,k=12𝛼2𝛥t𝛥x2(Uj+1,k-2Uj,k+Uj-1,k+Uj+1,k+1-2Uj,k+1+Uj-1,k+1) Since r=𝛼2𝛥t𝛥x2, (49) can then be simplified as:Uj,k+1-Uj,k=r2(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=r2Uj+1,k+(1-r)Uj,k+r2Uj-1,k+r2Uj+1,k+1+r2Uj-1,k+12(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 Thus, the well-known Crank-Nicolson scheme can be expressed as follows.-rUj-1,k+1+2(1+r)Uj,k+1-rUj+1,k+1=rUj-1,k+2(1-r)Uj,k+rUj+1,kor in matrix form:[-r2(1+r)-r00-r2(1+r)-r-r2(1+r)-r-r2(1+r)00-r]A[p((k+1)𝛥t)U1,k+1U2,k+1UN-2,k+1q((k+1)𝛥t)]U=[r2(1-r)r00r2(1-r)rr2(1-r)rr2(1-r)00r][p(k𝛥t)U1,kU2,kUN-2,kq(k𝛥t)]BRed means the values are known. In (52), we can notice that some known values reside in U. We must remove these known values from U and put them to the right side. The goal is to have a perfect AU=B system where all the unknowns are all in U. [2(1+r)-r0-r2(1+r)-r-r2(1+r)-r-r2(1+r)]A[U1,k+1U2,k+1UN-2,k+1UN-1,k+1]Unew=[2(1-r)rr2(1-r)rr2(1-r)rr2(1-r)]b[U1,kU2,kUN-2,kUN-1,k]Unow+[rp((k+1)𝛥t)+rp(k𝛥t)00rq((k+1)𝛥t)+rq(k𝛥t)]cB=(b·Unow)+c Now, we can find Unew=A-1B. In real implementation, this might not be plausible and we have to perform an iteration-based algorithm to find the optimal Unew. In Maple, we can use the Unew=LinearSolve(A,B) command. In MATLAB, we can use the popular Unew=A\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
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:uxx+uyy=f(x,y)0xa,0ybThe 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 and uyy can be approximated as:uxx=𝜕2u𝜕x2u(x-𝛥x,y)-2u(x,y)+u(x+𝛥x,y)𝛥x2uyy=𝜕2u𝜕y2u(x,y-𝛥y)-2u(x,y)+u(x,y+𝛥y)𝛥y2Hence: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)For simplification purpose, let us take 𝛥x=𝛥x=h, therefore: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)Let us take Uj,k=u(xj,yk), Fj,k=f(xj,yk), xj=jh, yk=kh, Uj,k+1=u(xj,yk+h) and Uj-1,k=u(xj-h,yk), hence: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,kUj-1,k+Uj+1,k+Uj,k-1+Uj,k+1-4Uj,k=h2Fj,kTherefore, the resulting discrete formula is as follows:Uj-1,k+Uj+1,k+Uj,k-1+Uj,k+1-4Uj,k=h2Fj,kIn order to solve (59), we need to reformulate it into a matrix form. In order to do so, we first define the following girds:
Drawing 11:Grid discretization of the 2D plate.
Next, based on the Drawing 11, we can then get the following equation:[-410101-410101-4001101-4101101-4101101-4101100-410101-410101-4]A[U1,1U2,1UNy-1,1U1,2U2,2UNy-1,2U1,Nx-1U2,Nx-1UNy-1,Nx-1]U~=[-(U1,0+U0,1)+h2F1,1-U2,0+h2F2,1-(UNy-1,0+UNy,1)+h2FNy,1-U0,2+h2F1,2h2F2,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]cOnce we have (60), we can solve it to acquire U~. In MATLAB, this can be easily done by using: U~=A\c. Let us now take a more simple example. Given the following diffusion PDE:uxx+uyy=f(x,y)where:0x1,0y1Nx=Ny=4 p(y)=0q(x)=0r(y)=100sin(𝜋y)s(x)=0f(x,y)=0h=0.25Find u(x,y) at its steady state condition. First, we create the matrix U as follows:Red colored numbers represents values at the boundaries. These values are calculated from p(y), q(x), r(y), and s(x). Next, we generate an equation similar (60) and calculate the U~.[-4101000001-4101000001-4001000100-4101000101-4101000101-4001000100-4100000101-4100000101-4][U1,1U2,1U3,1U1,2U2,2U3,2U1,3U2,3U3,3]=[-U1,0-U0,1-U2,0-U3,0-U4,1-U0,20-U4,2-U1,4-U0,3-U2,4-U3,4-U4,3] 3.2. Neumann Boundary 4. Two-Dimensional Unsteady Advection-DiffusionLet us consider the flowing PDE:𝜕T𝜕t+u𝜕T𝜕x+v𝜕T𝜕y2D Advection=𝛼(𝜕2T𝜕x2+𝜕2T𝜕y22D Diffusion)+f(61) is the general transient heat equation in 2D space. It is composed of advection process and diffusion process. 𝛼 is defines as:𝛼=K𝜌cpwhere:cp is the specific heat at constant (J(kgK)𝜌 is the fluid density (kg/m3)K is the thermal conductivity (W(mK))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+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)Rearranging (63) give us: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 Taking: Ti,j,k=T(xi,yj,tk), xj=j𝛥x, yk=k𝛥y, and tk=k𝛥t, we can rewrite (64) as follows:Ti,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+𝛥tfi,j,k 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+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)Rearranging (66) give us: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 Taking: Ti,j,k=T(xi,yj,tk), xj=j𝛥x, yk=k𝛥y, and tk=k𝛥t, we can rewrite (67) as follows:Ti,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+𝛥tfi,j,k4.3. DOWNWIND-UPWIND MethodThe downwind-upwind method is used when the x-velocity is negative and the y-velocity is positive.Ti,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+𝛥tfi,j,k4.4. UPWIND-DOWNWIND MethodThe downwind-upwind method is used when the x-velocity is positive and the y-velocity is negative.Ti,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+𝛥tfi,j,k4.5. MATLAB ImplementationFirst, we apply a grid discretization to the the two-dimensional rectangle as follows:
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
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. thus we have:T(i,2)-T(i,0)2𝛥x=AHence:T(i,0)=T(i,2)-2𝛥xA4.6.2. Right EdgeAssuming that boundary condition for the left edge is 𝜕T(x,h,t)𝜕x=B. thus we have:T(i,Nx+1)-T(i,Nx-1)2𝛥x=BHenceT(i,Nx+1)=T(i,Nx-1)+2𝛥xB4.6.3. Top EdgeAssuming that boundary condition for the left edge is 𝜕T(0,y,t)𝜕y=C. thus we have:T(2,j)-T(0,j)2𝛥y=CHenceT(0,j)=T(2,j)-2𝛥yC4.6.4. Bottom EdgeAssuming that boundary condition for the left edge is 𝜕T(,y,t)𝜕y=D. thus we have:T(Nx+1,j)-T(Nx-1,j)2𝛥y=DHenceT(Ny+1,j)=T(Ny-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 (1m×1m, and thickness of 0.001m) is heated only at the left edge. The temperature applied heat is 800K. 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:dz𝜌cp𝜕T𝜕t=dzK(𝜕2T𝜕x2+𝜕2T𝜕y2)-hu(T-T)-hd(T-T)-𝜀u𝜎(T4-T4)-𝜀d𝜎(T4-T4)wheredz=0.01mhu=hd=h=10W/(m2K)𝜀u=𝜀d=𝜀=0.5𝜎=5.67×10-8W/(m2K4)The material is copper where:cp=385J/(kgK)K=400W/(mK)𝜌=8960kg/m3 We can simplify (79) as follows:𝜕T𝜕t=K𝜌cp(𝜕2T𝜕x2+𝜕2T𝜕y2)-2hdz𝜌cp(T-T)-2𝜀𝜎dz𝜌cp(T4-T4)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:Ks=W/(mK)kg/m3J/(kgK)(Km2+Km2)-W/(m2K)mkg/m3J/(kgK)K-W/(m2K4)mkg/m3J/(kgK)K4Ks=Ks-Ks-KsThe 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. 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:0=K𝜌cp(𝜕2T𝜕x2+𝜕2T𝜕y2)-hudz𝜌cp(T-T)-hddz𝜌cp(T-T)-𝜀u𝜎dz𝜌cp(T4-T4)-𝜀d𝜎dz𝜌cp(T4-T4)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).
Drawing 14:Heat transfer in a two-dimensional channel flow
Conservation of mass:𝜕u𝜕x+𝜕v𝜕yDivergence=0Conservation of linear momentum:𝜕u𝜕t+u𝜕u𝜕x+v𝜕u𝜕y=-1𝜌𝜕p𝜕x+𝜈(𝜕2u𝜕x2+𝜕2u𝜕y2)𝜕v𝜕t+u𝜕v𝜕x+v𝜕v𝜕y=-1𝜌𝜕p𝜕yGradient+𝜈(𝜕2v𝜕x2+𝜕2v𝜕y2)LaplacianConservation of energy:𝜌cp(𝜕T𝜕t+u𝜕T𝜕x+v𝜕T𝜕y)=k(𝜕2T𝜕x2+𝜕2T𝜕y2)where:u(x,y,t) is the x-velocity at location (x,y) at time tv(x,y,t) is the y-velocity at location (x,y) at time tp(x,y,t) is the pressure at location (x,y) at time tT(x,y,t) is the temperature at location (x,y) at time t𝜈 is the fluid viscosity coefficient.cp is the specific heat at constant𝜌 is the fluid densityk is the thermal diffusivity rate In planar Poiseuille flow, the flow is only in x-direction. u becomes a function of only y. Therefore: v=0, 𝜕v𝜕t=0, 𝜕v𝜕y=0, 𝜕2v𝜕y2=0, 𝜕u𝜕t=0, 𝜕u𝜕x=0, and 𝜕p𝜕x=0. As a result, (82) to (84) can be simplified as follows:1𝜌(𝜕p𝜕x)=𝜈𝜕2u𝜕y2Since 𝜕p𝜕x is known, we are left with only finding 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𝜌(dpdx)=𝜈d2udy2 as follows:u(y)=1𝜌𝜈(dpdx)y22+C1y+C2At no-slip condition, u(0)=0 and u(ly)=0. By using these conditions, we can find C1 and C2. Let us the first boundary condition:u(0)=0=C2Next, let us take the second boundary:u(ly)=0=1𝜌𝜈(dpdx)ly22+C1lyC1=-1𝜌𝜈(dpdx)ly2Hence:u(y)=1𝜌𝜈(dpdx)y22+(-1𝜌𝜈(dpdx)ly2)C1y=1𝜌𝜈(dpdx)y2(y-ly) 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.
Drawing 15:Heat transfer in a Two-dimensional model of a cross-current heat exchanger
Conservation of mass:𝜕u𝜕x+𝜕v𝜕y=0Conservation of linear momentum:𝜕u𝜕t+u𝜕u𝜕x+v𝜕u𝜕y=-1𝜌𝜕p𝜕x+𝜈(𝜕2u𝜕x2+𝜕2u𝜕y2)𝜕v𝜕t+u𝜕v𝜕x+v𝜕v𝜕y=-1𝜌𝜕p𝜕y+𝜈(𝜕2v𝜕x2+𝜕2v𝜕y2)Conservation of energy:𝜌cp(𝜕TH𝜕t+u𝜕TH𝜕x+v𝜕TH𝜕y)=𝛼1(𝜕2TH𝜕x2+𝜕2TH𝜕y2)𝜌cp(𝜕TC𝜕t+u𝜕TC𝜕x+v𝜕TC𝜕y)=𝛼2(𝜕2TC𝜕x2+𝜕2TC𝜕y2)
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-axis (u) and the velocity along the y-axis (v). Additionally, within the fluid, there is pressure (p). (93) and (94) are the Navier Stokes equations. However, there are only two equations with three unknowns (u, v, and 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):𝜕u𝜕x+𝜕v𝜕y=0Conservation of linear momentum:𝜕u𝜕t+u𝜕u𝜕x+v𝜕u𝜕y=-1𝜌𝜕p𝜕x+𝜈(𝜕2u𝜕x2+𝜕2u𝜕y2)𝜕v𝜕t+u𝜕v𝜕x+v𝜕v𝜕y=-1𝜌𝜕p𝜕y+𝜈(𝜕2v𝜕x2+𝜕2v𝜕y2)Poisson equation for the pressure:𝜕2p𝜕x2+𝜕2p𝜕y2𝜌𝜕t(𝜕u𝜕x+𝜕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.u(x,y,t+𝛥t)=u(x,y,t)-𝛥t(u𝜕u𝜕x+v𝜕u𝜕y+1𝜌𝜕p𝜕x+𝜈(𝜕2u𝜕x2+𝜕2u𝜕y2))v(x,y,t+𝛥t)=v(x,y,t)-𝛥t(u𝜕v𝜕x+v𝜕v𝜕y+1𝜌𝜕p𝜕y+𝜈(𝜕2u𝜕x2+𝜕2u𝜕y2))Next, we take the derivative of u(x,y,t+𝛥t) with respect to x and v(x,y,t+𝛥t) with respect to y:𝜕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)𝜕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)We then add (98) to (99), giving us the following:𝜕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𝛥tWe set the left hand side of (100) to zero since the continuity equation (92) must be fulfilled. Rearranging (100) then gives us:𝜕2p𝜕x2+𝜕2p𝜕y2=𝜌𝛥t[𝜕u(x,y,t)𝜕x+𝜕v(x,y,t)𝜕y]-𝜌[A+B]ignore this partSince 𝜌𝛥t𝜌, we can then ignore the rightmost term of (101). 6.1. Discretizing the x-momentum Equation(93) is discretized as follows: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) Rearranging (102) gives us:u(x,y,t+𝛥t)=-𝛥t2𝜌𝛥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)-𝛥t2𝛥xu(x,y,t)[u(x+𝛥x,y,t)-u(x-𝛥x,y,t)]-𝛥t2𝛥yv(x,y,t)[u(x,y+2𝛥y,t)-u(x,y-𝛥y,t)]+u(x,y,t) 6.2. Discretizing the y-momentum Equation(94) is discretized as follows: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)Rearranging (103) gives us:v(x,y,t+𝛥t)=-𝛥t2𝜌𝛥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)-𝛥t2𝛥xu(x,y,t)[v(x+𝛥x,y,t)-v(x-𝛥x,y,t)]-𝛥t2𝛥yv(x,y,t)[v(x,y+𝛥y,t)-v(x,y-𝛥y,t)]+v(x,y,t) 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𝜌𝛥t(u(x+𝛥x,y)-u(x-𝛥x,y)2𝛥x+v(x,y+𝛥y)-v(x,y-𝛥y)2𝛥y)b(x,y) Rearranging gives us: 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) where:b(x,y)=𝜌𝛥t(u(x+𝛥x,y)-u(x-𝛥x,y)2𝛥x+v(x,y+𝛥y)-v(x,y-𝛥y)2𝛥y) 6.4. Cavity Flow
Drawing 16:Typical boundary conditions for a cavity flow.
6.5. Channel Flow6.5.1. Example 1
Drawing 17:Typical boundary conditions for a channel flow.
6.5.2. Example 2
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.
Drawing 19:An Euler-Bernoulli Beam
For a one-dimensional elastic beam with a length of L, its deflection is governed by the Euler-Bernoulli equation as follows:d2dx2[EId2w(x)dx2]=f(x)In (107), w(x) and f(x) are the beam's deflection and the external load at location x, respectively. E is the Young's modulus, and I is the area moment of inertia. We start by assuming that E and I are constants. Hence, (107) can be rewritten as:d4w(x)dx4=1EIf(x)The initial conditions to (108) are made by:d2w(x)dx2=1EIM(x)Internal moment at xd3w(x)dx3=1EIF(x)Internal shear force at xIn a central difference method we know that:dw(x)dx|xiwi-1-wi+12𝛥xd2w(x)dx2|xiwi-1-2wi+wi+1𝛥x2d3w(x)x3|xi-wi-2+2wi-1-2wi+1+wi+22𝛥x3d4w(x)dx|xiwi-2-4wi-1+6wi-4wi+1+wi+2𝛥x4where wi=w(i𝛥x) and 𝛥x is the segment length. Hence, we can rewrite (108) as follows"d4w(x)dx4=1EIf(x)i=0nwi-2-4wi-1+6wi-4wi+1+wi+2(𝛥x)4=1EIfiwhere fi=f(i𝛥x).We can also express in its matrix form as follows:[10001001-46-410001-46-410001-46-410010001]A[w-2w-1w0w1wnwn+1wn+2]W=𝛥x4EI[L2L1f0f1fnR1R2]Fwhere w-1, w-2, wn, wn+1, L1, L2, R1, and R2 are all the properties of additional nodes that we must add . We call them as the phantom nodes (see Drawing 20).
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 with respect to x and shear force is the third derivative of w with respect to x. The combinations of those three boundaries create nine possible configurations. However, we will focus only on 6 configuration as shown in Drawing 21.
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. This then translates into:w0=0Zero horizontal tangent means that dw(0)dx=0. which then translates into:w-1-w12𝛥x=0w-1=w1 b) Pinned-end BoundaryNo position displacement means that w(0)=0, which translates to:w0=0No bending moment means thatd2w(0)dx2=0, which then translates to:w-1-2w0+w1(𝛥x)2=0w-1=-w1 c) Free-end BoundaryNo bending moment means that dw2(0)dx2=0, which then translates into:w-1-2w0+w1(𝛥x)2=0w-1=-w1+2w0 No shear force means that d3w(0)dx3=0, which then translate into:-w-2+2w-1-2w1+w22(𝛥x)3=0w-2=2w-1-2w1+w2Substituting (116) to (117) gives us:w-2=2w-1-2w1+w2=2(2w0-w1)-2w1+w2=4w0-2w1-2w1+w2=4w0-4w1+w2 7.2.2. Boundary Conditions at the Right Enda) Fixed-end Boundary (Clamped)No position displacement means that w(L)=0, which then translates into:wn=0Zero horizontal tangent means that dw(L)dx=0, which then translates into:wn-1-wn+12𝛥x=0wn+1=wn-1 b) Pinned-end Boundary No position displacement means that w(L)=0, which then translates into:wn=0No bending moment means that d2w(L)dx2=0, which then translates into:wn-1-2wn+wn+1(𝛥x)2=0wn+1=-wn-1c) Free-end BoundaryNo bending moment means that dw2(L)dx2=0, which translates into:wn-1-2wn+wn+1(𝛥x)2=0wn+1=-wn-1+2wnNo shear force means that d3w(L)dx3=0, which then translates into:-wn-2+2wn-1-2wn+1+wn+22(𝛥x)3=0wn+2=wn-2-2wn-1+2wn+1By substituting (123) to (124), we can then acquire: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 To summarize, we put (112) to (125) into a table as follows.
LeftRight
Fixed-endL2=×L1=w1w0=0R1=wn-1R2=×wn=0
Pinned-endL2=×L1=-w1w0=0R1=-wn-1R2=×wn=0
Free-endL2=4w0-4w1+w2L1=-w1+2w0R1=-wn-1+2wnR2=wn-2-4wn-1+4wn
Table 1:All possible boundary conditions
Now, we can update the L1, L2. R1, and R2 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 and F contain only numbers and not variables. Currently, L1, L2, R1, R2 still contain some variables (w). This is not a desirable format since we can not apply the linear solve command to solve Aw=F. Free-freeFree-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.[10001001-46-410001-46-410001-46-410010001]A[w-2w-1w0w1wnwn+1wn+2]W=𝛥x4EI[4w0-4w1+w2-w1+2w0f0f1fn-wn-1+2wnwn-2-4wn-1+4wn]FStarting from the left end, we write down the equations made by Row 1 to Row 4:Row 1:w-2=4w0-4w1+w2Row 2 : w-1=-w1+2w0Row 3 : w-2-4w-1+6w0-4w1+w2=f0Row 4 : w-1-4w0+6w1-4w2+w3=f0By substituting Row 1 and Row 2 to Row 3, we get the following:4w0-4w1+w2-4(-w1+2w0)+6w0-4w1+w2=fo2w0-4w1+2w2=f0Next, by substituting Row 1 and Row 2 to Row 4, we get the following:-w1+2w0-4wo+6w1-4w2+w3=f1-2w0+5w1-4w2+w3=f1As for the right end, we can do the similar steps. By substituting Row n+1 and Row n+2 to Row n, we get the following:wn-2-4wn-1+6wn-4(-wn-1+2wn)+wn-2-4wn-1+4wn=fn2wn-2-4wn-1+2wn=fnNext, by substituting Row n+1 and Row n+2 to Row n-1, we get the following:wn-3-4wn-2+6wn-1-4wn-wn-1+2wn=fn-1wn-3-4wn-2+5wn-1-2wn=fn-1Now, we can remove Row 1, Row 2, Row n+1, and Row n+2. These rows are all made by the phantom nodes and must first be removed. As a result, we can update (127) into:[2-4200-25-41001-46-410001-46-4100001-46-41001-45-2002-42]A[w0w1w2wn-2wn-1wn]W=𝛥x4EI[f0fn]F 7.3.1. Fixed-fixed (fixed ended)For a case where both ends are fixed, we use the following equations.[7-4100-46-41001-46-410001-46-4100001-46-4001-47]A[w1w2wn-2wn-1]W=𝛥x4EI[f1fn-1]F 7.3.2. Pinned-pinned (simply supported)For a case where both ends are pinned, we use the following equations.[5-4100-46-41001-46-410001-46-4100001-46-4001-45]A[w1w2wn-2wn-1]W=𝛥x4EI[f1fn-1]F 7.3.3. Fixed-pinned (cantilever, simply supported)For a case where the left end is fixed while the right end is pinned, we use the following equations.[7-4100-46-41001-46-410001-46-4100001-46-4001-45]A[w1w2wn-2wn-1]W=𝛥x4EI[f1fn-1]F 7.3.4. Pinned-fixed (simply supported, cantilever)For a case where the left end is pinned while the right end is fixed, we use the following equations. [5-4100-46-41001-46-410001-46-4100001-46-4001-47]A[w1w2wn-2wn-1]W=𝛥x4EI[f1fn-1]F 7.3.5. Fixed-free (cantilever)For a case where the left end is fixed while the right end is free, we use the following equations.[7-4100-46-41001-46-410001-46-410001-46-41001-46-4101-45-2002-42]A[w1w2wn-1wn]W=𝛥x4EI[f1fn]F7.3.6. Fixed-explicit momentAn 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 that causes bending is applied at the right end:A manually-defined constant bending moment of M or dw2(L)dx2=M. This then translates into:wn-1-2wn+wn+1𝛥x2=Mwn+1=M𝛥x2-wn-1+2wnNo shear force means that d3w(L)dx3=0, which then translates into:-wn-2+2wn-1-2wn+1+wn+22𝛥x3=0wn+2=wn-2-2wn-1+2wn+1 By substituting to (139), to (141) we can then acquire: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+4wnWe can summarize the boundary as follows:
Left Boundary(Fixed-end)Right Boundary(Applied pure torque)
L1=×L2=w1w0=0R1=M𝛥x2-wn-1+2wnR2=wn-2-4wn-1+4wn
Applying these boundaries to the general system matrix gives us:[10001001-46-410001-46-410001-46-410010001]A[w-2w-1w0=0w1wnwn+1wn+2]W=𝛥x4EI[×w100fnM𝛥x2-wn-1+2wn2M𝛥x2+wn-2-4wn-1+4wn]FSince 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 From the third-to-the last row:wn-3-4wn-2+6wn-1-4wn+M𝛥x2-wn-1+2wn=fn-1wn-3-4wn-2+5wn-1-2wn=fn-1-M𝛥x2[7-4100-46-41001-46-4100001-46-41001-47]AThus, (142) cab be rewritten as:[7-4100-46-41001-46-410001-46-4100001-46-41001-451-2002-42]A[w1w2w3w4wn-2wn-1wn]W=1EI[00000-M𝛥x22M𝛥x2]FWe can then fine 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:𝜇𝜕2w𝜕t2(x,t)=-EI𝜕4w𝜕x4(x,t)+f(x,t)where 𝜇 is the mass per unit length, E is the Young's modulus, and I is the area moment of inertia. 𝜕2w𝜕t2(x,t)=-EI𝜇i=0nwi-2,j-4wi-1,j+6wi,j-4wi+1,j+wi+2,j(𝛥x)4+fi,j𝜕2w𝜕t2(x,t)=-EI𝜇𝛥x4[10001001-46-410001-46-410001-46-410010001]A[w-2w-1w0w1wnwn+1wn+2]W+[L2L1f0f1fnR1R2]FDiscretizing 𝜕2w𝜕t2(x,t) in (147) gives us: w(x,t-𝛥t)-2w(x,t)+w(x,t+𝛥t)𝛥t2=-EI𝜇𝛥x4[10001001-46-410001-46-410001-46-410010001]A[w-2w-1w0w1wnwn+1wn+2]W+[L2L1f0f1fnR1R2]F w(x,t+𝛥t)=-EI𝛥t2𝜇𝛥x4[10001001-46-410001-46-410001-46-410010001]A[w-2w-1w0w1wnwn+1wn+2]W+𝛥t2[L2L1f0f1fnR1R2]F+2w(x,t)-w(x,t-𝛥t) Another way is to use a MATLAB built-in ODE solver. In order to do so, let us take v=𝜕w𝜕t and 𝛼=EI𝜇𝛥x4, hence:𝜕v𝜕t(x,t)=-EI𝜇𝛥x4[10001001-46-410001-46-410001-46-410010001]A[w-2w-1w0w1wnwn+1wn+2]W+[L2L1f0f1fnR1R2]FNow, we can write:[𝜕w𝜕t𝜕v𝜕t]=[0I-𝛼A0][WV]+[0F]where V=[v-2vn+2]. and 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