INSTRUCTIONS
Please use only lowercase letters in your project file names.
o The project will be submitted online only.
o Only one code in an M-file form (i.e. ending with .m) should be submitted and an
accompanying PDF file should be uploaded with additional information asked in
various tasks described in the project.
o The code must be named in all lower-case letters as below: lastname_project.m (e.g.,
Wang_project.m)
o Also upload a PDF file which must contain additional information asked in various
tasks of the project statement.
o The PDF file name must be the same as the code file name, but just with a pdf
extension, i.e., Wang_project.pdf
o All submissions will be thoroughly checked for plagiarism and awarded a nil score
and/or disciplinary action initiated, if found plagiarized. No two project submissions, if
carried out independently, can be similar or identical. You are strongly advised to build
your own independent problem-solving thinking and skills by applying them to this
project problem.
Flow Reactor
Chemical Engineers make extensive use of ideal/steady state reactors in design work
that often ignore the start-up (unsteady state) portion of the reaction. It can be important
for many reasons to consider these start-up conditions and to know how long it takes
until a system reaches a steady state. For this project, you will be looking at the time
and position dependent concentration in a flow reactor that varies along the length of
the reactor and does not vary within the cross-section. You will solve for both the steady
state and unsteady state concentrations in the system.
PROBLEM STATEMENT: Figure 1 below depicts an elongated reactor that varies
along the length of the reactor and is constant in the cross-sectional area (plug flow
reactor/distributed-parameter system). The conditions in the reactor are maintained
such that the reaction 𝐴 → 𝐵 occurs via a first order decay in the concentration of A
(−𝑟𝐴 = 𝑘𝐶). There is also diffusion/ dispersion in the system.
Figure 1. Flow reactor that varies along the “x” direction with constant cross-sectional
area. A mass balance can be performed on a finite segment of length Δx (Equation 1):
where V is volume (m3 ), Q is volumetric flow rate (m3 /hr), C is the concentration of
reactant “A” at specific locations along the reactor (mol/m3 ), D is a diffusion
coefficient (m2 /hr), A is the reactor’s cross-sectional area (m2 ), and k is the rate
constant for the reaction term (hr-1 ).
When we take the limit as Δx and Δt approach zero, we get the equation 2
where U=Q/A, the average velocity of the fluid flowing through the reactor (m/s).
The reactor is subject to the Boundary Condition at the inlet (at x=0)
The reactor is also subject to the Boundary Condition that the concentration at the
outlet (x=L) is not affected by the diffusion rate :
The Initial Condition is
The parameters for the problem to be solved are:
(dx and dt were chosen to satisfy conditions of dynamic stability)
TASKS TO BE CARRIED OUT
TASK: A) If we take the Equation 2 to be steady state, 𝜕𝐶 / 𝜕𝑡 = 0, then it reduces to
Equation 3
This equation is a second order ODE, rather than a PDE. Solve this ODE analytically
and upload the analytical solution in your PDF. Solution MUST BE typed only (not
hand-written).
TASK: B) Now develop a code to solve Equation 3 numerically for distinct values of
concentration at each x position (x=0 to x=L). The code should produce a plot of the
concentration data against the position as a blue line with axes labels “concentration”
and “position” and a legend with the title “steady state data”.
Hints: You will need to build matrices and solve using a linear equation solving method
(e.g. Naïve-Gauss Elimination or Gauss-Jordan). Since the matrix size will change
depending on the chosen value of Δx, it could be useful to consider applying the
boundary conditions for the first and last rows and then implementing a loop for the
remaining rows. This method means you don’t have to manually change the matrices
for each new Δx. Consult problems in lab sessions to further think through your
implementation.
CAUTION: No built-in functions should be used for matrix equation solutions. Usage
of Gauss-type methods is permitted.
TASK: C) You need to now develop a code to numerically solve the time-dependent
PDE (Equation 2). This code must satisfy all of the below:
o The code should include the code from TASK: B in addition to the new code that is
written to solve TASK: C
o The code should run for a time until the reactor reaches the steady state
o The code should produce a new plot of the concentration data against position shown
via several lines, where each line is the concentration data for a specific time. The plot
should include the axes labels “concentration” and “position” and a plot title “nonsteady state data.”
o Within the code, the plotting commands from TASK: B should be retained and
modified such that the plot from TASK: B also shows the final time-step data from
TASK:C as empty circles, i.e., a blue line for TASK: B and empty circles from the final
time step of TASK: C.
o Overall, combining TASK: B and: C, only one code should be developed, and this
code should generate two plots, as described above.
Hints: You can use the finite-difference on the time-dependent term as:
where Ci is the concentration at the position being solved for at the current time, Cpi is
the previous time value concentration, and ∆𝑡 is the change in time.
The matrices can be set up similar to TASK: B. The steady state can be defined as the
condition when the concentrations from the previous time step are very close to the
current time step (i.e. 𝑎𝑏𝑠(𝐶𝑖 − 𝐶𝑝𝑖) < 𝑡𝑜𝑙𝑒𝑟𝑎𝑛𝑐𝑒). You can think of this as the situation
in which the program will run WHILE the condition is not satisfied. However, be
cautious that if your code is incorrect, then the loop may run indefinitely so consider
adding safety checks such as the BREAK statement, if needed.
The reactor is initially filled with water so the first 𝐶𝑝𝑖 values are zero. The previous
values become the current values at the next time step (i.e. at the next time iteration).
TASK: D) Finally, investigate the effect of changing each of the following parameter
values:
Based on your observations:
o Write a paragraph (or longer) and describe in your PDF about the most important/
interesting things you observed by changing these parameters and/or learned in this
project activity.
755
BENG603
CHE
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
MATLAB Implementation:
Solving ODEs and Simple PDEs
11/4/2020 (LAB8)
1
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
In today’s lab session, we will solve some
ordinary differential equations and simple partial
differential equations using finite difference
methods.
2
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Organization of Files
It is important for you to save all the files you create
during the lab sections, as they will be useful resources
when solving more sophisticated problems.
•
•
Find the folder you created on your Desktop for this course. You
must have named it CHE603 or che603 or some other name.
Inside this folder, create a new subfolder titled lab8 to keep all files
organized at a central location.
755
BENG603
CHE
NOTE: If you are not using your personal computer, DO NOT FORGET to upload your lab8 folder to your
UNH Box account, so that you have easy access to your files later. Even if you are using your own
computer, it is always a good practice to sync the MATLAB code folders for the course and keep them
updated under your free UNH Box account.
CHE603
lab8
3
755
BENG603
CHE
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
SECTION 1: Practice Problem Based on
Solving an ODE from Fluid Mechanics
4
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving an ODE – Flow Through a Pipe
Problem 1—Flow through a pipe:
In Fluid Mechanics, one learns about flow through pipes and about the Navier-Stokes’ equations. The model can be
described well in cylindrical coordinate systems (r, z, θ) with constant density, viscosity, and cross-sectional area. The
reduced equation for the steady state flow along the length of the pipe that depends on radial position is:
0=−
𝛥𝑃 µ 𝑑
𝑑𝑣
+
𝑟
𝐿
𝑟 𝑑𝑟
𝑑𝑟
where ΔP is the pressure drop in the pipe, r is the radial position in the pipe, µ is the viscosity of the fluid, and L is the
length of the pipe.
It is also subject to the Boundary Conditions: at the edge of the pipe (r = R), the velocity is zero (v=0) and at the center of the
pipe (r = 0), the velocity is a maximum (dv/dr=0).
The following is given: ΔP = -100 Pa, L = 2 m, R = 0.05 m, µ = 0.001 Pa*s.
755
BENG603
CHE
Solve the given ODE model numerically using finite difference methods using a Δr = 0.001.
Consider using centered difference for derivatives.
5
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
Solving an ODE – Flow Through a Pipe
How to think about this problem?
1) We need to evaluate the ODE for this problem using finite-difference methods.
2) In finite difference methods, we divide the total range of the independent variable for the problem (in this case it will
be the radius “r” of the pipe from 0 to R = 0.05 m) into finite pieces.
3) We divide it into finite pieces by dividing the radius by an increment “Δr” (for this problem it is Δr = 0.001 m). By
doing that, we will get a number N of equations to solve.
4) After applying the centered difference method, we get a set of equations to solve. Just like we did in previous labs, we
can convert those to matrix form to get a matrix A of coefficients that multiply the “unknown” velocities at finite
difference values, and a vector B with the values on the right-hand side of the equations (see below).
𝑎11 𝑎12 … 𝑎1𝑛
𝑣1
𝑏1
𝑎21 𝑎22 … 𝑎2𝑛
𝑣2
𝑏2
∗
=
:
∶
∶
:
:
𝑎𝑛1 𝑎𝑛2 … 𝑎𝑛𝑛
𝑣𝑛
𝑏𝑛
5) Finally, we use an algorithm, such as the Gauss-Jordan method, to solve the vector with the values of the velocities.
6
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving an ODE – Flow Through a Pipe
The first step for solving this problem is to realize that this flow is axisymmetric, therefore we only need to consider half
of the pipe for our calculations (see figure below). Also, based on the reduced equation provided in the problem, the
velocity is a function of the radius only.
The figure above also shows the boundary conditions for the problem. Now, we will derive the numerical expressions for
these boundary conditions:
B.C. 1: The problem states that at r = 0, dv/dr = 0. Therefore:
755
BENG603
CHE
𝑑𝑣(𝑧, 𝑟)
ቤ
=0
𝑑𝑟 𝑟=0
𝑣 2,1 − 𝑣 1,1 = 0
where v(2,1) and v(1,1) are two values for velocities at r = 0 but different points of the z-axis.
Note that (2,1) and (1,1) are simply indices where values are stored instead of actual values.
7
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving an ODE – Flow Through a Pipe
When converting this operation to the matricial form, we can write:
−1 1
𝑣(1,1)
= [0]
𝑣(2,1)
Therefore, we get the values for A(1,1) = -1, A(1,2) = 1 and B(1) = 0.
For this problem, the matrix A will be the matrix of the coefficients that multiply the
velocities in the equations we get after applying the centered difference method, and
B is the vector of the values that appear on the right-hand side of such equations.
755
BENG603
CHE
(Think of how we did the same procedure to solve other problems in previous labs!)
8
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving an ODE – Flow Through a Pipe
B.C. 2: The problem also states that at r = R, v = 0. Therefore:
𝑣(𝑟)ቚ
𝑟=𝑅
=0
𝑣 𝑁 =0
where N is number of nodes in the problem
Note that the number of values “N” of finite-differences along “r” in the problem is based on “dr”.
For example, if “dr=0.001” and R=0.05, then we get 50 divisions
When converting this operation to the matrix form, we get:
755
BENG603
CHE
1 𝑣(𝑁) = [0]
Therefore, we get the values for A(N,N) = 1 and B(N) = 0.
9
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
Solving an ODE – Flow Through a Pipe
Now that we obtained the matrix values for our boundary conditions, we will use the centered difference method to get
the expressions for the velocity as a function of r only. We start with our given equation:
𝛥𝑃 µ 𝑑
𝑑𝑣
+
𝑟
𝐿
𝑟 𝑑𝑟
𝑑𝑟
which can be rewritten as:
1 𝑑
𝑑𝑣
1 𝑑𝑟 𝑑𝑣
𝑑2𝑣
𝛥𝑃
𝑟(𝑖)
=
+𝑟 𝑖
=
𝑟(𝑖) 𝑑𝑟
𝑑𝑟
𝑟(𝑖) 𝑑𝑟 𝑑𝑟
𝑑𝑟 2
µ𝐿
0=−
Applying the centered difference method to this equation, we get:
1 𝑣 𝑖+1 −𝑣 𝑖−1
𝑣 𝑖 + 1 − 2𝑣 𝑖 + 𝑣 𝑖 − 1
+𝑟 𝑖
𝑟(𝑖)
2𝛥𝑟
𝛥𝑟 2
=
𝛥𝑃
µ𝐿
Refer to previous lab sessions if you need to refresh your
knowledge on finite-difference methods!
10
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving an ODE – Flow Through a Pipe
Multiplying both sides of the previous equation by 𝛥𝑟 2 , and reorganizing it we get:
𝛥𝑟
𝛥𝑟
𝛥𝑃 𝛥𝑟
1+
𝑣 𝑖 + 1 − 2𝑣 𝑖 + 1 −
𝑣 𝑖−1 =
2𝑟(𝑖)
2𝑟 𝑖
µ𝐿
2
This corresponds to:
• Elements of matrix A:
𝛥𝑟
𝐴 𝑖, 𝑖 + 1 = 1 +
2𝑟(𝑖)
𝐴 𝑖, 𝑖 = −2
𝐴 𝑖, 𝑖 − 1 = 1 −
𝛥𝑟
2𝑟(𝑖)
755
BENG603
CHE
• Elements of vector B:
𝛥𝑃 𝛥𝑟
𝐵 𝑖 =
µ𝐿
2
11
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
Solving an ODE – Flow Through a Pipe
Finally, we implement these expressions in MATLAB:
EXERCISE 1:
1) Create a new M-File, call it fluid_pipeflow.m,
and save it in your lab8 folder.
2) Type the commands shown on the right to your
M-File, and save it.
Do not run the code yet!
12
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
Solving an ODE – Flow Through a Pipe
Finally, we implement these expressions in MATLAB:
EXERCISE 1:
3) Type the commands shown on the right to your
M-File, and save it.
Do not run the code yet!
13
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
Solving an ODE – Flow Through a Pipe
Finally, we implement these expressions in MATLAB:
EXERCISE 1:
4) Type the commands shown on the right to your
M-File, save it, then run it!
Velocity profile as a function of radius of pipe is
shown in the plot.
Ask Instructor/TA if you did not
understand any step in this
implementation!
14
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
Solving an ODE – Flow Through a Pipe
We will now compare the numerical solution with the analytical solution given by the equation below. The solution
incorporates Boundary Conditions describing the no-slip condition at the wall (r=R), and the maximum velocity at the
center of the pipe (r=0). We will also add plotting commands to our code to plot velocity vs. radius for both numerical and
analytical solutions.
𝑣 𝑟 = −𝑣𝑚𝑎𝑥
𝑟
1−
𝑅
2
,
𝑣𝑚𝑎𝑥
𝛥𝑃 𝑅 2
=
4µ𝐿
EXERCISE 1:
5) Adjust the end of your code according to the commands shown below,
save it, then run it again!
Now the plot is updated with the analytical solution which perfectly matches
the numerical solution!
15
755
BENG603
CHE
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
SECTION 2: Practice Problem Based on
Solving a simple PDE—1D Heat PDE
16
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving a PDE – One-Dimensional Heat Equation
Problem 2—One-dimensional heat equation:
We have learned that a simple PDE in Chemical engineering is the One-Dimensional Heat Equation:
𝜕𝑇
𝜕2𝑇
=𝛼 2
𝜕𝑡
𝜕𝑥
The initial condition is for the “t” dimension, while the boundary conditions are for the length dimension “x”.
For example, the initial condition can be: T(x, t=0)=10
And the boundary conditions can be T(x=0, t) = 50 and T(x=L, t) = 100, where x is the length of the solid.
It is given that: L = 1 , t = 50 and α = 8.4*10-5.
755
BENG603
CHE
Plot the temperature as a function of time and length, using dt = 0.5 and dx = 0.01.
17
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving a PDE – One-Dimensional Heat Equation
For solving this PDE by a numerical method, the first step is to convert the PDE to Finite-Difference Equations.
For example, in our last homework assignment, we solved an equation for temperature at various nodes of a twodimensional heated plate. Here we have similar types of two dimensions, t and x.
Since we have two independent variables (t and x), then we just need an extra index for our finite-difference expressions
(for example, i = t and j = x).
Based on the Finite-Difference methods for 1st and 2nd Order ODEs, we get the following expressions, respectively:
𝝏𝑻 𝑻𝒊+𝟏,𝒋 − 𝑻𝒊,𝒋
=
𝝏𝒕
∆𝒕
and
755
BENG603
CHE
𝝏𝟐 𝑻 𝑻𝒊,𝒋+𝟏 − 𝟐𝑻𝒊,𝒋 +𝑻𝒊,𝒋−𝟏
=
𝝏𝒙𝟐
∆𝒙𝟐
Please, refer to previous lab sessions if you need to refresh
your knowledge on finite-difference methods!
18
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving a PDE – One-Dimensional Heat Equation
Therefore, the PDE for the problem will now look like this:
𝑻𝒊+𝟏,𝒋 − 𝑻𝒊,𝒋
𝑻𝒊,𝒋+𝟏 − 𝟐𝑻𝒊,𝒋 +𝑻𝒊,𝒋−𝟏
=𝜶
∆𝒕
∆𝒙𝟐
This equation can be rearranged to provide the temperature at the next time increment:
𝑻𝒊+𝟏,𝒋
𝑻𝒊,𝒋+𝟏 − 𝟐𝑻𝒊,𝒋 +𝑻𝒊,𝒋−𝟏
= 𝑻𝒊,𝒋 + 𝜶∆𝒕
∆𝒙𝟐
Note that both ∆𝒕 and ∆𝒙 need to be chosen for solving this.
755
BENG603
CHE
Try to obtain this same expression on your own!
Ask Instructor/TA if you have any questions!
19
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving a PDE – One-Dimensional Heat Equation
Finally, we implement this expression in MATLAB:
EXERCISE 2:
1) Create a new M-File, call it heatPDE.m, and save
it in your lab8 folder.
2) Type the commands shown on the right to your
M-File, and save it.
Do not run the code yet!
755
BENG603
CHE
Take notice of the nested for loop!
In its first iteration for time (i=1), it will calculate all the temperature values for the given time at all
different length nodes.
After the script finishes these calculations, the time is incremented (i=2), then the temperatures are
calculated at all length nodes again.
And so on…
20
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Solving a PDE – One-Dimensional Heat Equation
Since for this problem, temperature is a function of both t and x, the data we will get is in the form of a
matrix. Therefore, its plot will form a surface, not a curve. The commands below will plot this surface!
EXERCISE 2:
3) Type the commands shown on the right to your
M-File, and save it, then run it.
Your plot should look similar to one shown on the right side.
Warmer temperatures are in yellow and colder in blue.
755
BENG603
CHE
You can rotate the plot in MATLAB to see all dimensions
You can put a % sign in front of code line:
s.EdgeColor = ‘none’ and rerun to see a bit
different looking plot.
21
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
Lab Submission
Paste all work and codes in an MS-Word file, save it as a PDF, and upload the PDF online
when /lab8/ submission opens.
22
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
Lab Submission Problem – Newton’s law of cooling
Problem 1-Newton’s law of cooling
In a homicide or accidental death investigation, it is critical to estimate the time of death. From experimental
observations, it is known that Newton’s law of cooling describes how the surface temperature of an object changes
proportional to temperature difference between the object and surroundings, as given by the following ODE:
𝒅𝑻
= −𝒌 (𝑻 − 𝑻𝒂 )
𝒅𝒕
Where Ta is ambient temperature and k is a constant. Suppose at t = 0, a corpse is discovered by detectives and its
temperature is measured to be 29.5°C. It is assumed that the body temperature at the time of death was 37°C. After
two hours, the temperature of the corpse was recorded as 23.5°C. The ambient temperature is given as Ta = 20 °C.
755
BENG603
CHE
a) Determine k and the time of death.
b) Solve the ODE numerically and plot the temperature variation as a function of time.
23
Acid Structure
& Nucleic
to VMD
1: Intro
LAB 8:
PDEs
and Simple
ODEs
Solving
LAB
755
BENG603
CHE
Lab Submission Problem – Reaction Diffusion
Problem 2: Reaction Diffusion Problem
Consider the reaction diffusion ordinary differential equation given below:
𝒅𝟐 𝒄
= 𝝋𝟐 𝒄
𝟐
𝒅𝒙
The boundary conditions are: c(0)=0 and c(1)=1.
In chemical engineering lingo, the quantity φ is called the Thiele modulus for this problem. It is given that φ2 = 2.5
a) Solve this problem analytically using the methods learned during lectures.
b) Then, solve numerically using finite differences (consider 11 nodes including 2 for boundary conditions) and plot
numerical and analytical solutions.
24
755
BENG603
CHE
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
MATLAB Implementation:
Solving PDEs
11/18/2020 (LAB9)
1
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
755
BENG603
CHE
In today’s lab session, we will solve some partial
differential equations using finite difference
methods.
2
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
Organization of Files
It is important for you to save all the files you create
during the lab sections, as they will be useful resources
when solving more sophisticated problems.
•
•
Find the folder you created on your Desktop for this course. You
must have named it CHE603 or che603 or some other name.
Inside this folder, create a new subfolder titled lab9 to keep all files
organized at a central location.
755
BENG603
CHE
NOTE: If you are not using your personal computer, DO NOT FORGET to upload your lab9 folder to your
UNH Box account, so that you have easy access to your files later. Even if you are using your own
computer, it is always a good practice to sync the MATLAB code folders for the course and keep them
updated under your free UNH Box account.
CHE603
lab9
3
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
Solving a PDE – Heat Flow in a Two-Dimensional Plate
Problem 1—Heat Flow in a Two-Dimensional Plate:
The temperature distribution of a two-dimensional plate (see figure below) is described by the 2D Laplace PDE:
𝜕2𝑇 𝜕2𝑇
+
=0
𝜕𝑥 2 𝜕𝑦 2
755
BENG603
CHE
Consider a plate with the top and bottom sides at temperatures of T3=100 °C and T1=0 °C and left and right vertical
sides at temperatures of T2=75°C and T4=50°C, respectively.
Considering 100 equal divisions on the horizontal and vertical sides of the 10 cm x 10 cm plate, find the temperature
distribution of plate at all nodes using finite difference methods and plot the temperature distribution via a 2D plot and
highlighting higher temperature using a “warmer” color and lower temperature using a “cooler” color.
4
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
Solving a PDE – Heat Flow in a Two-Dimensional Plate
For solving this PDE by a numerical method, the first step is to convert the PDE to Finite-Difference Equations.
For example, in one of our past homework assignments, we solved an equation for temperature at various nodes of a
two-dimensional heated plate. Back then, we used finite-difference methods to get matrices, which we then solved using
Gauss-Jordan method to find the temperature at the nodes.
Using this method, however, becomes tiresome when the number of nodes defined is too large! This time, we will use for
loops to solve this problem instead. But first, we need to apply the Centered-Difference Method to obtain the needed
expression for the temperature at a given node, T(i,j), (see figure below).
Based on the Finite-Difference methods for 2nd Order ODEs, we get the following expressions:
Along x-direction:
𝜕2 𝑇
𝜕𝑥 2
=
𝑇𝑖+1,𝑗 −2𝑇𝑖,𝑗 +𝑇𝑖−1,𝑗
∆𝑥 2
755
BENG603
CHE
and
Along y-direction:
𝜕2 𝑇
𝜕𝑦 2
=
𝑇𝑖,𝑗+1 −2𝑇𝑖,𝑗 +𝑇𝑖,𝑗−1
∆𝑦 2
5
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
755
BENG603
CHE
Solving a PDE – Heat Flow in a Two-Dimensional Plate
Therefore, the PDE for the problem will now look like this:
𝑇𝑖+1,𝑗 − 2𝑇𝑖,𝑗 +𝑇𝑖−1,𝑗 𝑇𝑖,𝑗+1 − 2𝑇𝑖,𝑗 +𝑇𝑖,𝑗−1
+
=0
∆𝑥 2
∆𝑦 2
This equation can be rearranged to provide the temperature at a given node (see Figure on previous slide):
𝑻𝒊,𝒋
𝑻𝒊+𝟏,𝒋 +𝑻𝒊−𝟏,𝒋 +𝑻𝒊,𝒋+𝟏 +𝑻𝒊,𝒋−𝟏
=
𝟒
Try to obtain this same expression on your own!
Ask Instructor/TA if you have any questions!
6
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
755
BENG603
CHE
Solving a PDE – Heat Flow in a Two-Dimensional Plate
Finally, we implement this expression in MATLAB:
EXERCISE 1:
1) Create a new M-File, call it problem1.m, and
save it in your lab9 folder.
2) Type the commands shown on the right to your
M-File, and save it.
Do not run the code yet!
Take notice of the nested for loop!
In its first iteration for x (i=2), it will calculate all the temperature values for the given x value at all
different values of y. After the script finishes these calculations, the x value is incremented (i=2), then
the temperatures are calculated at all values of y again. And so on…
7
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
Solving a PDE – Heat Flow in a Two-Dimensional Plate
Since for this problem, temperature is a function of both x and y, the data we will get is in the form of a
matrix. Therefore, its plot will form a surface, not a curve. We can hide the z-axis by assigning colors that
change according to their values.
EXERCISE 1:
3) Type the commands shown on the right to your
M-File, and save it, then run it.
Your plot should look similar to this:
Notice how most of the
temperature values for the nodes
are very low.
755
BENG603
CHE
Why do you think this happens?
8
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
755
BENG603
CHE
Solving a PDE – Heat Flow in a Two-Dimensional Plate
Most of the values for the temperatures at the nodes are close to zero because we are initializing their
temperatures at zero. Therefore, it is in our best interest to find those values when the temperature
distribution reaches the steady-state. The code below ensures that we achieve that!
EXERCISE 1:
4) Edit your script by adding the commands shown
on the right, save it, then run it.
Your plot should look similar to this:
The while loop ensures that the temperatures at
all nodes are recalculated until the difference
between the new and old values are lower the
specified tolerance, i.e. until the system reaches
steady-state!
9
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
755
BENG603
CHE
Solving a PDE – Heat Flow in a Two-Dimensional Plate
Now, let’s see the effect that decreasing the number of nodes has on the temperature distribution.
EXERCISE 1:
5) First, delete the line of code that says “shading interp”, so that you can see the grids again.
6) Next, change the values of dx and dy to 0.2 to get 50 equal intervals on each axis, save it, then run it.
7) Then, change the values of dx and dy to 1 to get 10 equal intervals on each axis, save it, then run it.
100 intervals
50 intervals
10 intervals
10
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
Solving a PDE – One-Dimensional Wave Equation
Problem 2—One-Dimensional Wave Equation:
Solve the following PDE, which is known as the one-dimensional wave equation:
2
𝜕2𝑢
𝜕
𝑢
2
=
𝑐
𝜕𝑡 2
𝜕𝑥 2
Boundary conditions: 𝑢(𝑡, 0) = 𝑢(𝑡, 𝐿) = 0
Initial conditions: 𝑢 𝑡 = 0, 𝑥 = sin(2π𝑥) and
𝜕𝑢
𝜕𝑡
𝑡 = 0, 𝑥 = 0
Consider L=1 and c2=1.
755
BENG603
CHE
Plot u as a function of x only.
11
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
755
BENG603
CHE
Solving a PDE – One-Dimensional Wave Equation
Just like the previous problem, the first step is to convert the PDE to Finite-Difference Equations.
Based on the centered-difference method for 2nd Order ODEs, we get the following expressions:
Along t-dimension:
𝜕2 𝑢
𝜕𝑡 2
=
𝑢𝑖+1,𝑗 −2𝑢𝑖,𝑗 +𝑢𝑖−1,𝑗
∆𝑡 2
and
Along x-dimension:
𝜕2 𝑢
𝜕𝑥 2
=
𝑢𝑖,𝑗+1 −2𝑢𝑖,𝑗 +𝑢𝑖,𝑗−1
∆𝑥 2
12
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
755
BENG603
CHE
Solving a PDE – One-Dimensional Wave Equation
Therefore, the PDE for the problem will now look like this:
𝑢𝑖+1,𝑗 − 2𝑢𝑖,𝑗 +𝑢𝑖−1,𝑗
𝑢
− 2𝑢𝑖,𝑗 +𝑢𝑖,𝑗−1
2 𝑖,𝑗+1
=
𝑐
∆𝑡 2
∆𝑥 2
This equation can be rearranged to provide the displacement (u) with an increment in the time:
𝑢𝑖+1,𝑗
∆𝑡 2 𝑐 2
= 2𝑢𝑖,𝑗 − 𝑢𝑖−1,𝑗 +
𝑢𝑖,𝑗+1 − 2𝑢𝑖,𝑗 +𝑢𝑖,𝑗−1
∆𝑥 2
Try to obtain this same expression on your own!
Ask Instructor/TA if you have any questions!
13
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
Solving a PDE – One-Dimensional Wave Equation
Now, we need to convert the initial conditions to numerical form. The first one is straight-forward:
Given form
𝑢 𝑡 = 0, 𝑥 = sin(2π𝑥)
Numerical form
𝑢 1, 𝑗 = sin(2π𝑥(𝑗))
Meanwhile, for the second initial condition, we need to use the forward-difference method:
Given form
𝜕𝑢
𝜕𝑡
𝑥, 𝑡 = 0 = 0
Numerical form
𝑢 2, 𝑗 − 𝑢 1, 𝑗
=0
∆𝑡
755
BENG603
CHE
𝑢 2, 𝑗 − 𝑢 1, 𝑗 = 0
Try to obtain this same expression on your own!
Ask Instructor/TA if you have any questions!
14
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
Solving a PDE – One-Dimensional Wave Equation
Finally, we implement this expression in MATLAB:
EXERCISE 2:
1) Create a new M-File, call it problem2.m, and
save it in your lab9 folder.
2) Type the commands shown on the right to your
M-File, and save it.
Do not run the code yet!
755
BENG603
CHE
Notice how we are introducing a vector x. We do
this because we need to calculate the initial
conditions for all positions of x!
15
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
755
BENG603
CHE
Solving a PDE – One-Dimensional Wave Equation
Finally, we implement this expression on MATLAB:
EXERCISE 2:
3) Type the commands shown on the right to your
M-File, and save it, then run it.
Your plot should look similar to this:
In this plot, each curve represents the
displacement of the wave as a function of its
position for a given value of time.
16
Acid Structure
& Nucleic
LAB 1: Intro to VMD
PDEs
9: Solving
LAB
Practice Problem (No submission) – Two-Dimensional Wave Equation
Problem 3—Two-Dimensional Wave Equation:
Solve the following PDE, which is known as the two-dimensional wave equation:
2
2
𝜕2𝑢
𝜕
𝑢
𝜕
𝑢
2
=𝑐
+
𝜕𝑡 2
𝜕𝑥 2 𝜕𝑦 2
The solution of this equation represents the shape of a biological membrane (see figures).
755
BENG603
CHE
Boundary conditions: 𝑢(0, 𝑦, 𝑡) = 𝑢(𝑎, 𝑦, 𝑡) = 𝑢(𝑥, 0, 𝑡) = 𝑢(𝑥, 𝑏, 𝑡) = 0
Initial conditions: 𝑢 𝑥, 𝑦, 𝑡 = 0 = 𝑥𝑦 2 − 𝑥 3 − 𝑦 and
𝜕𝑢
𝜕𝑡
𝑥, 𝑦, 𝑡 = 0 = 0
Consider a=2, b=3 and c2=36. Plot u as a function of x and y only.
17
Purchase answer to see full
attachment