University of New Hampshire Applied Mathematics for Chemical Engineers Questions

User Generated

pbare753

Engineering

University of New Hampshire

Description

You need the skill of MATLAB.

This homework is based on the PowerPoint I post and you can only use the method on the PowerPoint to solve the problem; otherwise, it would be defined as the wrong solutions and I will not accept it.

The homework is about the applied mathematics for chemical engineers and using the MATLAB to solve two related problem which is about the PDE and second order ODE.

PS: There are four questions in this homework and I have done half of them but I am not sure if they are correct.

Unformatted Attachment Preview

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
User generated content is uploaded by users for the purposes of learning and should be used following Studypool's honor code & terms of service.

Explanation & Answer

here it is up to part B, I'm working on part C now...


Anonymous
I was struggling with this subject, and this helped me a ton!

Studypool
4.7
Trustpilot
4.5
Sitejabber
4.4

Related Tags