A Mathcad‐based educational experience to address the design of nonisothermal plug flow reactors

Mathcad is a simple‐to‐use and intuitive mathematical software that helps students to minimize the mathematical difficulties involved in solving engineering problems. The design of nonisothermal plug flow reactors (PFR) is a fundamental issue within the field of chemical reaction engineering; however, its teaching–learning process is hindered by students' mathematical difficulties in solving ordinary differential equations. In this paper, the software Mathcad was conveniently integrated into an educational experience through the resolution of two real case studies. In the first one, a simple liquid‐phase reaction is considered in a PFR working at different operating conditions, whereas the second case evaluates a PFR taking place multiple reactions (parallel reactions) with a heat exchanger attached. The assessment of this experience, which was held into two 5‐h Mathcad workshops, revealed that Mathcad made the design of non‐isothermal PFR more appealing, facilitated the understanding of the design process, and brought another dimension to the way the students perform complex calculations.


| INTRODUCTION
Chemical reaction engineering (CRE) is the area that studies the interface between two pillars in Chemical Engineering, such as reaction analysis and reactor design. The first one requires combining the principles of stoichiometry, thermodynamics, and kinetics to understand the rate laws, mechanisms, and the conversions of chemical reactors. The second one covers reactor modeling and the identification of its optimal conditions to carry out the reactions [23]. There is no doubt that CRE is a fundamental course of any chemical engineering program and, therefore, a solid knowledge of this matter is expected from these graduate students.
Unfortunately, as in many other engineering courses [1,17], the difficulties experienced by the students enrolled in a CRE course are very often attributed to the calculations rather than in understanding the phenomenon in question [22]. In this sense, depending on the ideal reactor model (i.e., continuous stirred-tank reactor, plug flow reactor [PFR], or bath reactor) and the operating conditions (i.e., isothermal, adiabatic, or with heat exchange), the reactor design would involve complex mathematical operations like the simultaneous solution of nonlinear algebraic equations, integral equations or ordinary differential equations (ODEs) [24]. All these contents are covered over the Bachelor's degree first course. Even so, students do not always attain meaningful learning because the unknowns are not associated with specific process variables, neither are the equations to specific engineering phenomena. Thus, contextualization into real engineering scenarios makes the students improve their understanding of the maths core principles and realize why they need them [22].
Over the last few years, the pedagogy of teaching engineering courses is continually advancing through the implementation of mathematical software packages, since they minimize the mathematical difficulties for students and allow teachers to devote more time addressing other engineering issues, making the course more student-friendly [4,8]. In this context, a significant number of computer-aided educational experiences by using mathematical software packages like Polymath [14,25,26], Mathcad [5,6,9,22], MATLAB [12,13,19], Maple [11,20], MATHEMATICA [2,10,15], or Excel [3,16,18,27] are found in the engineering field. Even though engineering competencies can be well acquired using any of this software, it is important to select one that does not involve an additional complexity given the limited time that can be dedicated to learning it in class [4]. Mathcad allows engineering students with limited or no formal training in programming to solve complex problems (i.e., simultaneous ODEs) thanks to the use of predefined functions [4,22]. Moreover, Mathcad offers students symbolic calculation, easy matrix operation, the use of units, creation of matrix, graphics, text, and tables into the same worksheet, as well as tutorials to assist them in their learning [4,7,22,23]. Therefore, in our case, given that most chemical engineering students were not trained in any operating mathematical software, Mathcad is a simple-to-use and intuitive software that help them to solve and understand the engineering problems by focusing exclusively on the phenomenon itself rather than on solving mathematical difficulties. Thus, it is intended that with the use of Mathcad the main obstacle to failure in the resolution of engineering problems is not their mathematical complexity.
The School of Engineering at the University of Huelva, thanks to various Teaching Innovation Projects, is currently implementing the use of Mathcad in some courses corresponding to a Bachelor's Degree in Industrial Chemical Engineering and a Master's Degree in Chemical Engineering [5,6,21,22]. As for the CRE area, we used Mathcad for the design of isothermal reactors involving multiple reactions (series and parallel reactions) [6]. Now, in this study, we consider the use of Mathcad to address the design of nonisothermal reactors, more precisely, nonisothermal PFR. Compared to isothermal reactors [6], the mathematical complexity increases since the molar and energy balances for nonisothermal PFR, as well as the energy balance for the heating/cooling fluid used for the heat exchange, are governed by ODEs, which should be simultaneously solved. Here, we present two case studies that were analyzed with the aid of Mathcad, and we provide guidelines on how Mathcad facilities the design of non-isothermal PFR involving simple and multiple reactions. In the first case study, we explore the effects that reactor operating conditions (adiabatic conditions or with a heat exchange) exert on the design of a PFR in which a simple liquid-phase reaction takes place. In the second one, the temperature and molar flow rate profiles are obtained for a PFR with multiple reactions (parallel reactions) and with a heat exchanger attached in a cocurrent mode. This paper not only deals with the fact of using Mathcad for solving problems related to reactor design but also in proposing a novel teaching methodology based on computational tools that favor students' problem-solving skills.

| Teaching methodology
The experience was carried out by 34 students (24 of them male and 10 were female) enrolled in the course "Chemical Reactor II" (reference: 606210217) corresponding to the 3rd year of the Bachelor's Degree in Industrial Chemical Engineering, at the School of Engineering, University of Huelva, Spain. This course is mandatory and was imparted over 60 h during the summer semester of the academic course 2020-2021. Its contents are clearly divided into two parts: the first one is a comprehensive study about the design of nonisothermal homogeneous chemical reactors, whereas the second one is dedicated to heterogeneous reactors. For this educational experience, which belongs to the first part of the course, the students' prerequisites, among others, include an adequate application of molar and energy balance on isothermal chemical reactors and the knowledge of chemical kinetics for homogeneous reactions.
The teaching methodology for the computer-based teaching methodology is illustrated in the flowchart of Figure 1. The teaching methodology was divided into two Mathcad workshops (5 h each) that took place in a computer room at the School of Engineering where a version of Mathcad Prime 4.0 is available. The first Mathcad workshop was mainly focused on solving ODEs by using the "Odesolve" Mathcad function, given that it is essential to address the design on nonisothermal PFR. It is important to note that this first workshop is not a training course on Mathcad, so students should already have assimilated the basic notions necessary to start using Mathcad (e.g., knowledge about the software-user interface, the definition of simple functions, or 2D graphical representation). This fact should not be an impediment for the chemical engineering students who participate in this educational experience, as they are familiar with Mathcad and used this software in courses such as "Fluid Flow" [5], "Chemical Reactor I" [6], or "Unit Operations of Chemical Engineering II" [22]. All of this means that the first Mathcad workshop is dedicated almost exclusively to learning how to apply the "Odesolve" function for solving ODEs, which is not an easy task and is completely new for students.
The second 5-h Mathcad workshop was devoted to solving two case studies about the design of nonisothermal PFR involving simple or multiple reactions. The case studies were selected to expand knowledge and understanding about this issue. During this workshop, the teacher proposed to the students the main guidelines to solve the case studies with Mathcad and was in permanent support to answer possible questions.
Finally, it is worth noting that this methodology was not intended to replace the traditional theory sessions, which result to be necessary to acquire the fundamentals behind the design of nonisothermal chemical reactors but to F I G U R E 1 Flowchart of the computerbased teaching methodology complement them reducing the mathematical complexity involved. Thus, the fundamentals related to the design of these reactors are taught in the theory sessions and leaving their practical applications for this educational experience. Finally, the student's opinion on the development of this educational experience and some authors' reflections are presented.

| Theoretical considerations
The case studies here proposed deal with the design of nonisothermal PFR. In this section, we will briefly summarize the main fundamentals that support their design. So, first of all, it is important to mention that basic design equations, rate laws, and stoichiometric relationships used for isothermal reactor design are still valid for nonisothermal reactors. Thus, the molar (1) and energy (2) balance for any type of ideal reactors can be written as [24]: The molar balance (1) states that the flow rate of a specie j that enters (F j0 ), leaves (F j ), and reacts (  r dV is equal to the rate of the heat flow (Q) that goes in the system from the surroundings, minus the rate of shaft work (Wṡ ) done by the system on the surroundings, and plus the difference between the inlet and outlet conditions of the summation of the molar flow rates (F j0 and F j , respectively) of each species multiplied by their corresponding molar enthalpies (H j0 and H j , respectively). For the case of a PFR, its ideal flow model (i.e., assuming that each reactor section is perfectly mixed but there is no mixing between sections in the axial direction) makes that the molar balance applied for a volume differential results in an ODE. Considering the compound A as a limiting reactant, we obtain the following ODE (3): Thus, the variation of the conversion (X A ) with the position in the reactor (i.e., the conversion profile) depends on its reaction rate ( r − A ) and inlet molar flow in Equation (3), it can be also written as Equation (4) allows determining the molar flow profile for each specie.
The molar balance is sufficient for the design of an isothermal PFR (without variation in temperature). However, a PFR can exchange heat with other heating/ cooling fluid, as can be illustrated in Figure 2A (in this case the heat exchange occurs at a cocurrent mode; this is, the reactant/product current flows through the reactor in the same direction as the heating/cooling fluid). Now, the energy balance (Equation 2) for the specie A is also converted into the following ODE: where U is the overall heat transfer coefficient, a is the heat exchange area per unit volume of reactor area, T a the temperature of the heating/cooling fluid, ∆H R the reaction enthalpy, I j is the ratio between the initial molar fraction of species j (y j0 ) and the initial molar fraction of specie A (y A0 ), and ∆C p the variation of specific heats of products and reactants. Recalling , where ν j is the ratio between the stoichiometric coefficient of specie j and specie A, it is possible to obtain Equation (6) from (5): If the temperature of the cooling/heating fluid (T a ) varies down the reactor, we must add the energy balance for this fluid, Equation (7): Graphical representation of (a) a PFR with heat exchange at a cocurrent mode, and (b) an adiabatic PFR. PFR, plug flow reactors dT dV where m f and C pf are the mass flow rate and specific heat of the heating/cooling fluid, respectively, and T a1 and T a2 are its initial and final temperature (see Figure 1A). In addition to that, a PFR can be thermally insulated from the surroundings and, thereby, operate at adiabatic conditions ( Figure 2B).
is zero in Equation (5), the energy balance for an adiabatic PFR is Finally, a nonisothermal PFR can involve more than one reaction (i.e., multiple reactions). Considering a PFR with the configuration shown in Figure 2A (heat exchange) and that involves multiple reactions, its corresponding energy balance will be written as Equation (9): where the subscript j refers to the species, the subscript i refers to the particular reaction, q is the number of independent reactions, and n is the number of species.

| Case study 1: Effects of operating conditions on a PFR involving a simple reaction
In this case study, a simple liquid-phase reaction is considered in a PFR working at different operating conditions: (i) adiabatically, (ii) with heat exchange and no variation in the temperature of heating/ cooling fluid, and (iii) with heat exchange and variation in the temperature of heating/cooling fluid. Thus, by analyzing both reactor temperature and conversion profiles, students can easily explore the effects that operating conditions exert on the reactor design and, consequently, make a decision to select the best option. The resolution of this case study with Mathcad entails the use of the "root" predefined function and the solution of a system of equations, which includes nonlinear algebraic and first-order ODEs, by using the "Odesolve" function. The problem statement is: The liquid-phase reaction A +B→C is carried out in a PFR with an internal diameter of 0.4 m. An equal molar feed in A and B enters at 300 K with a volumetric flow rate of 2 L/s. The initial concentration of each reactant is 0.1 mol/L. Most reacting systems involve more than one reaction and do not operate isothermally [24]. Consequently, it is considered to be one of the main issues that CRE has to address. More precisely, this case study considers a PFR with a heat exchanger attached, in a cocurrent mode, and occurring multiple reactions (parallel reactions). Its mathematical resolution implies the solution of a large system of equations including ODEs, which, without the use of a computational tool, would be very complex to solve. This is, therefore, a clear case in which to prove the enormous usefulness of the "Odesolve" Mathcad function for solving engineering problems. The problem statement is: The following liquid-phase reactions occur in a PFR (internal diameter of 0.5 m and length of 5.1 m): Pure A is feed at a rate of 150 mol/s, a temperature of 150°C, and a concentration of 0.1 mol/L. A heat exchanger at a cocurrent mode is added to the PFR. Determine the temperature and molar flow rate profiles down the reactor.
Data: To begin with, as for solving any other engineering problem with Mathcad, the first step is to include in the Mathcad worksheet the known parameters. Thus, the conditions of the feed current (volumetric flow rate, temperature, and reactants' concentrations), formation enthalpies at a reference temperature of 273 K, specific heats and kinetic parameters (activation energy, Ea, and the value of the reaction rate constant at 300 K, k_300K) are included in Figure 3 (framed with a dashed line). Item (a) asked students to calculate the reactor volume to achieve 90% conversion if the PFR operates at adiabatic conditions. Before applying the molar and energy balance (Equations (3) and (8), respectively) needed to solve this item, it is necessary to define other parameters that are included in them. Thus, the stoichiometric relationships (i.e., the reactants' concentration as a function of the conversion) and the dependence of reaction rate constant on temperature, k(T), should be implemented into the rate equation. The Arrhenius equation was used to determine the temperature dependence of the reaction rate constant. First, the pre-exponential factor (ko) is calculated from the value of the reaction rate constant at 300 K, (k_300K) and the activation energy (Ea). Once the preexponential factor has been calculated, the value of the reaction rate constant at any temperature, k(T), can be obtained. Finally, as can be observed in Figure 3, the reaction enthalpy is calculated. Note that, in this case, there is no dependence of this F I G U R E 4 A portion of the Mathcad worksheet for case study 1, item a) parameter on temperature, since the variation of specific heats of products and reactants (∆C p ) is zero. Now, after these preliminary calculations and the inclusion of the known parameters, it is the moment to apply the molar and energy balance in the worksheet (Figure 4). With the aid of Mathcad, the solution of the system of equations, which also includes ODEs, can be obtained straightforwardly with the application of the "Odesolve" function. So, after writing the command "Given," the ODEs corresponding to the molar and energy balance, as well as the equation describing the dependence of the reaction rate constant on temperature, are included in the Mathcad worksheet. To solve the system, the boundary conditions of the reactor position-dependent variables (conversion, temperature, and reaction rate constant) should be also included; these are, the initial conversion (zero), the initial temperature (300 K), and the initial value of the reaction rate constant at 300 K. Then, the solution variables are defined ("X_adiabatic," "T_adiabatic," and "k_adiabatic") and the system of equations is solved. Attending to the ease to insert graphs into the Mathcad worksheets, both conversion and temperature profiles are obtained (Figure 4). By analyzing them, it is possible to determine the volume required to achieve 90% conversion and the final temperature. However, the "root" Mathcad function, which allows a single nonlinear algebraic equation to be solved with no need for solve block, provides an accurate solution. To do so, the root of a new function ("f1(VVf)"), which is defined as the difference between the reactor conversion for each position ("X_adiabatic(V_reactor)") and the final desired conversion ("X_final"), gives the adiabatic PFR volume. Thus, the adiabatic reactor volume is 317.8 L (internal diameter of 0.4 m; length of 2.5 m).
F I G U R E 5 A portion of the Mathcad worksheet for case study 1, item (b) and (c) For the two next items, considering a reactor volume equal to the previously calculated for the adiabatic conditions (317.8 L), students are requested to determine both conversion and temperature profiles if a heat exchanger is added to the PFR, in a cocurrent mode, and the heating/cooling fluid is constant (item [b]) or varies (item [c]) down the reactor. Figure 5 shows the system of equations needed to solve both items by using "Odesolve" Mathcad function. Comparing them with the system of item (a), it is stated that: (i) the molar balance remains the same, (ii) the energy balance should include the term "Ua(Ta-T)", and (iii) the balance energy for the heating/cooling fluid should be also included as its temperature varies down the reactor (item [c]). Finally, since the volume reactor considered for the three items are the same, the effects of PFR operating conditions can be addressed by the conversion and temperature profiles displayed in Figure 6. As deduced, the best option would be to incorporate a heat exchanger varying the temperature of the coolant fluid; however, if its temperature is constant, there is a significant reduction in the conversion, which is even lower than that noticed for the adiabatic conditions. Thus, according to the energy balance, an increase in the conversion leads to a proportional increase in the reactor temperature if the exothermic reaction occurs in adiabatic conditions or keeping the temperature of the cooling/heating fluid constant; consequently, as seen in Figure 6, both conversion and reactor temperature profiles display the same shape. However, if the temperature of the heating/cooling fluid varies down the reactor, these profiles cease to be proportional and different zones can be observed: (a) in the first part of the reactor as the temperature of the heating/ cooling fluid (350 K) is higher than the feed current (300 K) there is heat exchange. It results in a slight decrease in the temperature of heating/cooling fluid and a more noticeable increase in the reactor temperature which is favored by the energy released by the exothermic reaction; (b) in the middle of the reactor, as conversion increases, the reactor temperature continues rising and the coolant temperature begins to increase; (c) finally, in the last part, the conversion tends to level off so that the effect of the energy released by the reaction on temperature profiles is smaller. Consequently, the temperature changes are attributed to the heat exchange F I G U R E 6 Conversion and temperature profiles for case study 1 between the heating/cooling fluid and the reactor, which results in a decrease in the reactor while the temperature of the heating/cooling fluid continues to rise.

| Computer-based solution of the case study 2
As in case study 1, the first step is to introduce the known parameters in the Mathcad worksheet ( Figure 7); so, the conditions of the feed current, enthalpies of reaction, specific heats, kinetic data, and information about the heat exchange are included. Afterward, it is time to add the equations needed to obtain both temperature and conversion profiles, as shown in Figure 8. First, the molar balances for all species are added. From the kinetic data, the reaction rate of specie A in reaction 1 and in reaction 2 (Equations 4′ and 5′, respectively) are used to express the reaction rate of each species. Thus, considering the following relationships, it is possible to write the reaction rate of each species as follows: So, molar balances for each species (Equations 1′, 2′, and 3′) are written in the Mathcad worksheet. Then, it is needed to know the concentration of each species as a function of the reactor position. Thus, considering the stoichiometry and the fact that feed current is solely composed of species A, the Equations (8′), (9′), and (10′) are added. In addition, the total molar flow is equal to the sum of the individual molar flow (Equation 11′). To F I G U R E 7 A portion of the Mathcad worksheet including all known variables for case study 2 finish with the system of equations, the energy balances for the reactor (Equation 12′) and for the coolant fluid (Equation 13′) are included. Thus, taking into consideration that there are multiple reactions, the energy balance for the reactor was written according to Equation (9) of the "Theoretical consideration" section. The boundary conditions for all the reactor position-dependent variables are also included: initial molar flow rates, initial values of reaction rates, initial values of reaction rate constants, initial concentrations, the initial value of the total molar flow rate, initial reactor temperature, and initial temperature of heating/cooling fluid.
After all this, the system of 13 equations including 5 ODEs can be easily solved by using the predefined "Odesolve" Mathcad function without the need to use any programming. Figure 9 shows the molar flow rate and temperature profiles for a volume reactor of 1 m 3 . As expected, there is a progressive decrease in F A while F B and F C are growing. When the specie A is totally consumed the reaction stops and the molar flows remain constant, which means that, beyond this point, the PFR behaves just as a heat exchanger. This fact is also reflected in the temperature profile. While the reaction occurs, the reactor temperature ("T_x") F I G U R E 8 A portion of the Mathcad worksheet including all equations for case study 2 undergoes an exponential increase which is accompanied by a slight decrease in the temperature of the coolant fluid; however, when the reaction has finished, the coolant fluid is able to cool the reactor.

| On the students' feedback and authors' reflections
The students' opinion on the development of this educational experience was obtained through a survey conducted at the end of the second workshop. The questions and their corresponding scores and standard deviations are gathered in Figure 10. The survey, which was answered by all the students (34 in total), had two purposes: on the one hand, to find out if the students have learned to use the "Odesolve" function (items Q1 and Q2) and, on the other hand, to know the impact that this Mathcad-based educational experience has on the teaching-learning process (items Q3-Q5). Finally, an overall rate of the educational experience was also requested. The main conclusions derived from the students' survey are: • Students are familiar with the use of Mathcad which allows them to easily learn how to use the "Odesolve" Mathcad function, as deduced from the answers to items Q1 and Q2. This is essential as the mathematical solution of the case studies is based on the usage of this predefined Mathcad function. • The mathematical complexity involved in the design of nonisothermal PFRs is minimized by using Mathcad (item Q3 displays the highest score of all), which clearly allows students to better understand the theoretical concepts (Q4). • It is quite satisfying to note that students would encourage other teachers to use this software (item Q5). This undoubtedly indicates the great usefulness of Mathcad. Therefore, it is the students themselves who believe it is necessary for teachers to be trained in Mathcad and use it in their classes.
F I G U R E 9 Molar flow rate and temperature profiles for case study 2 • The experience rating was 4.10, pointing out that Mathcad made the design of nonisothermal PFR more appealing, facilitated the understanding of the design process, and brought another dimension to the way the students perform complex calculations.
Based on the conclusions derived from the students' survey, as well as from our own perception, some authors' reflections are presented: • It would be highly interesting to integrate Mathcad software in the teaching-learning process of engineering degree courses. This would mean that students would minimize their mathematical difficulties, focusing their efforts on understanding the phenomena under study. • Compared to other software packages, Mathcad is a highly recommended option for those students who have no prior knowledge of programming (as in the case of chemical engineers). Thus, its natural maths notation and the use of predefined functions (e.g., "root" or "Odesolve") make it possible to solve complex mathematical operations. • The ease with which graphs can be inserted into the Mathcad worksheet facilitates a better analysis of the results obtained and, in the case of chemical reactor design, the selection of the best option.
• The School of Engineering at the University of Huelva is currently doing a lot of work in extending the use of Mathcad. The next step would be to make students aware that Mathcad is of huge benefit to their academic or working life. • As a further incentive to use Mathcad, it seems reasonable to think that in today's knowledge society those students who have acquired more skills for the management of mathematical software packages or programs have greater employability prospects.

| CONCLUSIONS
This paper shows how to integrate the Mathcad software into a CRE course to address the design of non-isothermal PFR. The educational experience, which was carried out by third-year chemical engineering students through 5-h Mathcad workshops, consisted of the resolution of two real case studies. In the first one, a simple liquid-phase reaction is considered in a PFR working at different operating conditions, whereas the second case evaluates a PFR taking place multiple reactions (parallel reactions) with a heat exchanger attached. The main findings derived from the students' survey revealed that the mathematical complexity involved in the design of F I G U R E 10 Student survey questions about this educational experience and their corresponding responses these reactors is minimized by using Mathcad, making it more appealing, and helping to understand the design process. In addition, students encourage other teachers to use Mathcad in their courses which highlights the usefulness of this software.
NOMENCLATURE t time (s) N j moles of specie j (mol) Ê sys energy within the system (J) F j0 molar flow rate in the feed current of specie j (mol s −1 ) F j molar flow rate in the output current of specie j (mol s −1 ) V reactor volume (m 3 ) r ij reaction rate of specie j in reaction i (mol m −3 s −1 ) k ij reaction rate constant of specie j in reaction i (variable) C j concentration of specie j (mol m −3 ) Q̇rate of the heat flow (W) Ẇs rate of shaft work (W) H j molar enthalpy of specie j (J mol −1 ) X j conversion of specie j (−) U overall heat transfer coefficient [W·m −3 ·K −1 ] a heat exchange area per unit volume of reactor area (m −1 ) T reactor temperature (K) T a1 initial temperature of heating/cooling fluid (K) T a2 final temperature of heating/cooling fluid (K) T 0 initial temperature of feed current (K) ∆H R reaction enthalpy (J mol −1 ) y j0 initial molar fraction of species j (−) y A0 initial molar fraction of species A [−] I j ratio between y j0 and y A0 [−] C pj specific heat of specie j (J mol −1 K −1 ) ∆C p variation of specific heats of products and reactants (J mol −1 K −1 ) ν j ratio between the stoichiometric coefficient of specie j and specie A (−) m f mass flow rate of heating/cooling fluid (g s −1 ) C pf specific heat of heating/cooling fluid (J mol −1 K −1 ) E a activation energy (J mol −1 ) k 0 pre-exponential factor (variable) ν 0 volumetric flow rate of feed current (m −3 s −1 ) R1 universal gas constant (J mol −1 K −1 )