I would like to use the Newton Raphson method for the found values of T
Department of Mechanical Engineering ME 2173: Numerical Methods Using MATLAB Fall 2019 MATLAB Project 1 Publish your script and upload to Blackboard. Include your MATLAB files as well. NO OTHER FORMATS WILL BE ALLOWED Project Introduction and Theory The ideal gas equation has been widely used to relate Temperature, Pressure and Volume of Gases. The equation can be written as1, ?? = ?? ? Where, ?? = ????? ?????? ?? ???, ? = ??? ???????? , ? = ???????????, and ? = ??? ????????. There are several other equations that characterize gas better than the Ideal Gas equation, like the van der Waals equation. This equation, however, is not as accurate for temperatures above the critical temperature, so equations like the Redlich–Kwong (RK) equation of state would generally be used at higher temperatures. The dimensionless RK equation is shown below2 ? = 3? ?? − 3? − 9? √???(?? + 3?) (1) Where, ? = ? 100 (2) ? = ? 1000 (3) In Equations 2 and 3, ? is the sum of the numbers in your abc123. Project Experiment and Objectives Consider a piston cylinder as shown in Figure 1. This cylinder has a piston on top, attached to a pressure gauge, a heat source at the bottom, a non-specific gas trapped inside. Figure 1. Gas expansion in a piston-cylinder device In this example, pressure and volume are known. However, temperature is unknown! The goal of this project to determine temperature ? for different pressure ?. Tasks Task 1: Finding approximate root and calculating isothermal expansion work Let’s assume the gas is being heated. Increasing the temperature will increase the volume (i.e. as heat increases, gas expands). As no external pressure is being added, consider it as 50 for now. As gas expands, one can measure the volume, therefore, molar volume of gas can be calculated at any moment, but the temperature is unknown. For this task, you need to determine temperature ? at ? = 50 and ?? = 0.8 using the brute force method and the isothermal expansion work between two values of ?? using the symbolic toolbox. The steps to follow are listed below: i. Rewrite the equation 1 for root finding, ?(?) = 3? ?? − 3? − 9? √???(?? + 3?) − ? ii. Use Equations 2 and 3 to determine ? and ?. Remember, Q is the sum of the numbers in your abc123 iii. Determine the value of ? for which ?(?) ≈ 0 using Brute force (for loop) approach and a tolerance of 5%. Vector of temperatures 0 ≤ ? ≤ 100 iv. Use the symbolic toolbox to calculate the isothermal expansion work at the temperature T found in the previous step. In this case the volume is changing from ?? = 0.7 to ?? = 0.9 and the work is calculated as the integral of the pressure function as ? = ∫ ???? 0.9 0.7 Task 2: Determine suitable root finding method From Task 1, you have found the root of the function using the brute force method. In this section, you need to determine a better suited method to find the root of the equation, (i.e. determine ? so that ?(?) ≈ 0) considering ? = 50 and ?? = 0.8 while using different methods. The steps to follow are listed below: i. Compare the two bracketing methods with two open methods from the list below, a. Bracketing methods: i. Bisection method. ii. False position method. b. Open methods: i. Newton Raphson method. ii. Modified Newton Raphson method. iii. Secant method. ii. Make a comparison between results obtained with both open and closed methods. a. Provide a comparison of accuracy between the 4 methods. b. Provide a comparison of computational efficiency between the 4 methods. iii. Discuss which of the methods is the better option for root finding based on the comparisons. Task 3: Final task! Pressure is not constant anymore Now that the best root finding method has been established, it is time for the final task! To determine Temperature when Pressure is no longer constant, that is, find ? so that ?(?) ≈ 0 while the pressure ? is varying. The list of steps to achieve this task is shown below: i. Similarly to Task 2, rewrite the equation for root finding (keep in mind that now ? and ? are variables) ?(?, ?) = 3? ?? − 3? − 9? √???(?? + 3?) − ? ii. Use the script from task 2 and convert it into a FUNCTION FILE that a. Takes a, b, ??, ?, and an initial value of T as inputs. b. Determines the value for ? at which ?(?, ?) is approximately 0 using the root finding method you chose as best in Task 2. c. Outputs that specific value of ?. iii. In a different script file, use the “For Loop” set up so that each iteration employs a different value of ? and calls the function file, from the previous step, to determine ?. Create an array to store the resulting ? after each iteration. iv. Run the script file with ? = 1: 20 and create a matrix where the first column is the vector of pressures, P, and the second column is the vector of temperatures, T. v. Plot ? vs. ?. References 1. Clapeyron, E. (1834). "Mémoire sur la puissance motrice de la chaleur". Journal de l'École Polytechnique (in French). XIV: 153–90. Facsimile at the Bibliothèque nationale de France (pp. 153–90) 2. Murdock, James W. (1993), Fundamental fluid mechanics for the practicing engineer, CRC Press, pp. 25–27, ISBN 0- 8247-8808-7 %1-D Heat equation %example 1 at page 782 %lambdat=c.k/h^2 %T(x,t)=temperature along the rod %by finite difference method %(T(x,t+dt)-T(x,t))/dt = (T(x+dx,t)-2T(x,t)+T(x-dx,t))/dx^2 %solve for T(x,t+dt) by iteration %heat constant close all clear all c=1; L=1; %length of the rod N=5; %# of elements %dicretize Xspace h=L/N; x_vec=0:h:L; %discretize time T=0.5; M=50; k=T/M; t_vec=0:k:T; %temperature matrix T_mat=zeros(length(x_vec),length(t_vec)); %boundary conditions T_mat(1:end/2)=200; T_mat(end,:)=0; %initial conditions T_mat(:,1)= sin(pi.*x_vec); [tt,xx]=meshgrid(t_vec,x_vec); subplot(2,1,1); mesh(xx,tt,T_mat); title('Analytical vs numerical in 3D'); T_mat_calc=zeros(length(x_vec),length(t_vec)); T_mat_calc(:,1)= sin(pi.*x_vec); lamda=(c*k/(h^2)); for tdx=1:length(t_vec)-1 for idx=2:length(x_vec)-1 T_mat(idx,tdx+1)=T_mat(idx,tdx)+lamda*((T_mat(idx +1,tdx)-2*T_mat(idx,tdx)+T_mat(idx-1,tdx))); T_mat_calc(idx,tdx+1)= (exp((-pi^2)*(t_vec(tdx +1))))*sin(pi*(x_vec(idx))); end end %plot subplot(2,1,2) [tt,xx]=meshgrid(t_vec,x_vec); mesh(xx,tt,T_mat); figure plot(x_vec,T_mat(:,1),x_vec,T_mat(:,11),x_vec,T_mat(:,21),x_vec,T_mat(:,31),x_vec,T_mat(:,41),x_vec,T_mat(:,51)); 1 xlabel('Rod length (m)'); ylabel('Temperature (C)'); legend('Initially','At t=0.1sec','At t=0.2 secs','At t=0.3 secs',' At t=0.4 secs',' At t=0.5 secs'); title('Numerical Solution'); figure plot(x_vec,T_mat_calc(:,1),x_vec,T_mat_calc(:,11),x_vec,T_mat_calc(:,21),x_vec,T_mat_calc(:,31),x_vec,T_mat_calc(:,41),x_vec,T_mat_calc(:,51)); xlabel('Rod length (m)'); ylabel('Temperature (C)'); legend('Initially','At t=0.1sec','At t=0.2 secs','At t=0.3 secs',' At t=0.4 secs',' At t=0.5 secs'); title('Analytical Solution'); 2 3 Published with MATLAB® R2019b 4