SPTP CHE 596F
Popular in Course
verified elite notetaker
Popular in Chemical Engineering
verified elite notetaker
This 203 page Class Notes was uploaded by Mr. Lilian Johns on Thursday October 15, 2015. The Class Notes belongs to CHE 596F at North Carolina State University taught by Staff in Fall. Since its upload, it has received 40 views. For similar materials see /class/223787/che-596f-north-carolina-state-university in Chemical Engineering at North Carolina State University.
Reviews for SPTP
Report this Material
What is Karma?
Karma is the currency of StudySoup.
You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!
Date Created: 10/15/15
CHE596M MultiScale Modeling of Matter Dr Erik E Santiso Lecture 26 Chemical Reactions and Activated Processes II 7 NC STATE UNIVERSITY Outlook Finding reaction mechanisms The ZeroTemperature String Method Calculating reaction rates and defining reaction coordinates in complex systems Introduction The BennettChandler method Blue Moon molecular dynamics Transition path sampling Aimless Shooting Likelihood Maximization Finite Temperature String Method Other methods 7 NC STATE UNIVERSITY 7 Outlook Finding reaction mechanisms The ZeroTemperature String Method Calculating reaction rates and defining reaction coordinates in complex systems Introduction The BennettChandler method Blue Moon molecular dynamics Transition path sampling Aimless Shooting Likelihood Maximization Finite Temperature String Method Other methods 7 NC STATE UNIVERSITY Finding reaction mechanisms The ZeroTemperature String Method 18 19 Main idea Minimize a reaction path string by evolving it according to the equation pl VVpl In order to keep resolution close to the transition state region after each timestep in the integration of the previous equation do a reparameterization step This is done by fitting a spline to the new points and redistributing them along the path 7 NC STATE UNIVERSITY 7 Finding reaction mechanisms The ZeroTemperature String Method 15 16 It can be shown that NEB is a particular case of the zero temperature string method However it is simpler and more efficient to parameterize the reaction path using arc length intrinsic reaction coordinate or energyweighted arc length Example httpwwwcims nvuedueve2strinqhtml j NC STATE UNIVERSITY Outlook Finding reaction mechanisms The ZeroTemperature String Method Calculating reaction rates and defining reaction coordinates in complex systems Introduction The BennettChandler method Blue Moon molecular dynamics Finite Temperature String Method Transition path sampling Aimless Shooting Likelihood Maximization Other methods 7 NC STATE UNIVERSITY Introduction Remember this slide Chemical Reaction A process that results in the interconversion of chemical species A chemical species is defined as An ensemble of chemically identical molecular entities that can explore the same set of molecular energy levels on the time scale of the experiment From IUPAC glossary ofterms used in physical organic chemistry available online at httpwwwchemqmuacukiupacgtpocl Key point Species are defined in terms of a time scale So which chemical species are you I don t know yet Give me some time 7 NC STATE UNIVERSITY 7 Introduction Separation of time scales By definition a chemical reaction implies a separation of time scales between molecular motion translation rotation vibration and the reactive process Usually the reaction time is many orders of magnitude longer Consider for example a simple reaction q H J L m itsa ck At 800 K the rate constant is 46x10394 5391 1 so the characteristic time for reaction is 35 min Fastest mode C H stretching is 3000 cm391 9x1013 5391 so the characteristic time for vibration is 11 fs If we use a timestep of 1 fs and a system with 1000 molecules we need to do 2x1014 MD steps to see a reaction event With current computers this would take 60000 years of CPU time NC STATE UNIVERSITY 7 J 9 6 Introduction Separation of time scales Most of the time in a normal MD simulation is spent exploring the nonreactive part of phase space This problem of separation of time scales is not only characteristic of chemical reactions Other processes such as diffusion in solids micelle formation protein folding etc have the same problem We need a way to sample the interesting part of the dynamics But first we need a statistical mechanical expression for the rate of a reaction J NC STATE UNIVERSITY 7 7 Reaction rates Reaction rates in chemistry are defined in terms of rate laws Rate laws are phenomenological equations and cannot be derived They arise as the result of experimental observations We can find statistical mechanical expressions for phenomenological coefficients by assuming that the rate law is correct and then relating it to microscopic variables HmmLooks like rAz kACA Concentratlon mol L gtMHH 0 05 1 15 2 Time min 7 NC STATE UNIVERSITY 7 Calculation of rate constants The rate law for a chemical reaction is a phenomenological equation We need to assume that it is correct and find an expression for the rate constant The derivation is in many ways similar to the derivation of the fluctuationdissipation theorem we saw for Langevin dynamics or Einstein s equation for the diffusion coefficient They relate the linear response of the system to a perturbation to the equilibrium fluctuations of microscopic quantities over time 2 6D tigg dd 02gt G 2kT We want to find a relation of this kind between the rate constant and the equilibrium fluctuations of the concentration Y 71 NC STATE UNIVERSITY 7 Calculation of rate constants Let us consider a simple reaction A lt gt B for which we know that the rate equations are amp k c kc fA rB g h 7fkch kcB Where kf and k are the fonlvard and backward rate constants and cA cB are the number densities of A and B We will assume that they are normalized ie CA CB 2 1 At equilibrium the concentrations are constant and we have eq k 0 eq I f 262 CA k NC STATE UNIVERSITY 1 Calculation of rate constants Now imagine that we disturb the system by adding a little bit of A If we define AcAt cAt cff we have kch 16ch 2 kfcf ACAkr1 C1q ACA gt dt dt gt kf krAcA k kfcjf krcjfl gt gt kf krAcA k krcgq krcjfl kf krAcA 263 We can solve this equation to get ACA t ACA 06Xpl kf lg ACA 0exp lTR 264 NC STATE UNIVERSITY Calculation of rate constants Equation 264 will be our connection with the microscopic variables The characteristic reaction time ER is defined as 1 1 K 1 CalCa 66quot TR E Z e61 B eq A W Z i 265 kfkr kf1Keq 191ch cAl kf The next thing we need is a microscopic definition of cA Since it is normalized we can interpret it as the probability that a given molecule is a molecule of species A given that the system is large enough Now let us consider the variable 2 1 33 eA e i 0 b 266 NC STATE UNIVERSITY 7 1 Calculation of rate constants The reaction coordinate a is a function of the configuration The value in corresponds to the transition state between A and B and defines a dividing surface in configurational space that separates reactants from products If we now average 6A over the canonical ensemble we obtain the probability that the configuration is a reactant configuration which is equal to cA lt9 gt ldF9AzeXpBH A Q Here H is the Hamiltonian of the system Although we are deriving the expression for the canonical ensemble the final result is actually the same for other ensembles as well For derivations see 23 2 CA 267 7 NC STATE UNIVERSITY 7 Calculation of rate constants Now we would like to find an expression for ACA t the relaxation ofthe concentration upon an initial perturbation so that we can relate it to equation 304 If we perturb the Hamiltonian so that the states with Q lt E have a lower energy the equilibrium concentration of A will be higher We can do this by lowering the energy of these states HHO xeA Unperturbed Hamiltonian Small perturbation 268 If the perturbation is small we can expand ACA on a Taylor series and keep only the leading term 0 mo NC STATE UNIVERSITY 7 7 Calculation of rate constants Thus we have all 9 eBHO739eA ax H ax H ax are 0 A H Jdmez e mHo xeA B awe e MHo xeAlf 7L J df eiBH07 eA UdFeZHOAGlqz A20 BAkeilgt0 lt9Agt 2610 Or since 92 9A ACA Lukego 933 BAlt9Agt01 lt9Agt0 mech AcA 0 2611 NC STATE UNIVERSITY Calculation of rate constants If we turn off the perturbation now at t O the system will relax back to the equilibrium concentration As before we can expand on a Taylor series in 7 and keep the linear term This way we obtain magic llgloreAltrgteMH0W 2 57L 57L J dFe BlHo MA 0 J dFBeAOpAteBHO7 eA0 7L J dFeBlHoWAOH J39dl GA teBHO7quoteA 0HdFBGA OeBHO739eA 01 UdFe BH0K9AZ 70 BikeAwait 0 lt9A0gt0lt9Atgt0 2612 1 NC STATE UNIVERSITY 7 Calculation of rate constants Be patient The ugly math is almost 211 the manipulations we got qg g A6140 2 Bl lt9A 0914 tgt0 lt9A 0 GA 00 The first term is a time autocorrelation function The second term is the product of the equilibrium concentration of A at t 0 times the same concentration at time I as A gt 0 In this limit both of these terms are equal to the equilibrium concentration of A in the unperturbed system and hence ACA 1 mk9A0911 0 02qu 2613 5 Calculation of rate constants Now from equation 2611 we had AcA0 Bkcif c q 2 BA A39 CA CB Substituting this into 2613 we get A620 AcA0 lt9A09Ar0 05f eq eq CA CB 2614 From the phenomenological rate law we had mt Aclltogtexpi rm TR 3 f NC STATE UNIVERSITY Calculation of rate constants Comparing the last two equations we get exp tIR lt9AltOgt92SZgt 2 Giff 2615 This equation is actually valid only for long times As is the case with Einstein s equation for the diffusion coefficient and other autocorrelation function expressions the shorttime behavior of the function can be quite different Differentiating 2615 with respect to time we get 1 9 09 t T eXP tTR lt Acme gt0 2616 R A B 7 NC STATE UNIVERSITY 7 Calculation of rate constants If we now assume that the reaction is truly a rare event the time sampled by our simulation will be much smaller than TR Then we can assume that eXp tIR1 and we get the very important result 1 lt9A 061 tgt0 lt6A 0914 tgt0 2 2617 T R 01 cf 01 Off The last equality can be obtained by changing time tto t 0 and reversing time This only changes the sign of the autocorrelation function If we now recall that ER eggcf we get the Bennett Chander equation lt6AltogteAltrgtgt0 for T f eq m0 CA ltlttltlt 1R 2618 1 NC STATE UNIVERSITY Calculation of rate constants Recalling that 9At 9k it we can rewrite the Bennett Chandler equation in the more common form 23 3539 3 kf lt5 0 a a gt for 1mg ltlttltlt 1R 2619 There are two common ways of factoring this equation and obtaining rate constants from simulations We can write again that the equilibrium concentration is 02quot 6A and factor it into an Arrheniuslike equation 8a0 ail 0elar ail 8a0 ail k 2620 f lt5l 0 933w lt9Agt Relative probability of an activated Conditional average dynamic configuration equilibrium 7 NC STATE UNIVERSITY 7 Calculation of rate constants In order to obtain the other factorization we need to look at the limit when t gt 0 2 2 E3 lt5la0 claim ail lt5lamp0 ampl 0el 0lgt 262 This is a rate constant for which we count as reactive all trajectories that pass the dividing surface in the direction of the products It can be shown 345 by changing to a coordinate system where the one of the coordinates is the reaction coordinate that the righthand side is equivalent to the transition state theory TST rate constant kT kT A160 7K 7e W 2622 Here K is an equilibrium constant between the reactant and the transition state k kTST 7 NC STATE UNIVERSITY 7 Calculation of rate constants Using this result we can factor the BennettChandler equation as ltsialtogt ail 0elar ail ltsialtogt ail 09amp0lgt kf KkTST 2623 lt6la0 ail oelaoll of Here K is the transmission coefficient This is analogous to equation 2620 in that one factor is an equilibrium quantity and the other is a dynamic quantity The TST constant assumes that everything that reaches the top of the barrier ultimately reacts The transmission coefficient accounts for the possibility of recrossings of the barrier and is always less than or equal to 1 You may have seen equation 2623 with an extra factor which accounts for quantum tunneling through the barrier This is a non classical term and cannot be obtained from this theory 2 NC STATE UNIVERSITY 7 The BennettChandler Method There are two ways we can use the BennettChandler equation to calculate rate constants If we use equation 2620 we can obtain the second equilibrium term by realizing that lt5i zb Probability density that the system is at the top of the barrier GA Probability that the system is on the quotreactantquot side of the barrier The calculation of probabilities is equivalent to the calculation of free energies remember the coarsegrained Hamiltonian There are standard techniques to deal with this problem such as umbrella sampling or thermodynamic integration We will see shortly an example of how to do this 7 NC STATE UNIVERSITY 7 The BennettChandler Method The first term in equation 2620 is a conditional average lt5E O gik09E t ab Count only trajectories that lt5EO ED start at the top of the barrier Qt We can obtain this quantity by using molecular dynamics We sample only dynamic trajectories that start at the top of the energy barrier This condition however introduces a bias in the initial distribution function and we will need to correct for this NC STATE UNIVERSITY The BennettChandler Method Another option is to use equation 263 We can obtain the equilibrium constant between the reactants and the transition state using for example the Reactive Monte Carlo method From this we can obtain the TST rate constant If necessary the TST rate constant can then be corrected by evaluating the transmission coefficient lt5a0 ai 09ar aigt lt6la0 ail oel olgt KE for 1mg ltlt t ltlt 1R These averages can also be obtained from a molecular dynamics simulation sampling from trajectories that start at the top of the barrier 7 NC STATE UNIVERSITY 7 TST RXMC Example Decomposition of HI in the bulk in solution and in carbon nanopores 6 0 25 Solvent 020 ccta I C02 NH3 p o 0 UI Relative Reaction Rate Iogmklka 2 U1 0 50 1 00 1 50 200 250 Pressurebar FIG 6 TSRXIVIC simulations of hydrogen iodide decomposition in differ ent inert solvents The temperature is 57315 K te reaction is retarded by N113 N1 and CH4 while C03 and 14 solvents are predicted to increase the rate somewhat Relative Reaction Rate logmkk0 o 40 5 Initial HI Concentration moiL P400 bar P200 bat Relative Reaction Rate Iogmkko 0 I 10 20 30 40 50 00 70 80 Initial HI Concentration moiL FIG 2 TSRXMC simulations lot and experimental measurements Ref 30 39 iiie 39 quot quot i39wctinn The temperature is int 57315 K lb 59455 K The simulated pressure is shown corresponding to different concentrations NC STATE UNIVERSITY TST RXMC Example Decomposition of HI in the bulk in solution and in carbon nanopores 6 50 m b W CARBON NANOTUEES 9 9 a m 0 0 5 Relative Reaction Rate logm tnbsKV F39 9 Relative Reaction Rate iogmkk Pressureoar PressureBat FIG 7 TSRxMC simulations ofthe hydrogen iodide decomposition rate at pp 9 qumMc sunnlauons cfme hydrogen iodide lecolnposi on rate at g t y 7 t v J I in 57315 K m the bulk phase Oi and in carbon nanombes of different ham 100 um 1D 125 Inn A 150 1111110 and l 75 11111 9 eter 188 nauombe Di and 1010 nanotube O i NC STATE UNIVERSITY Outlook Finding reaction mechanisms The ZeroTemperature String Method Calculating reaction rates and defining reaction coordinates in complex systems Introduction The BennettChandler method Blue Moon molecular dynamics Transition path sampling Aimless Shooting Likelihood Maximization Finite Temperature String Method Other methods NC STATE UNIVERSITY Blue Moon Molecular Dynamics From the point of view of the molecules reactions happen only once in a blue moon actually blue moons are comparatively much more common they happen every two or three years not 60000 Last one was May 30 2007 the next one will be December 2009 The Blue Moon MD method is a very elegant way to obtain both terms in equation 2620 from molecular dynamics The main idea is to 1 Use constrained dynamics to generate configurations at the top of the barrier thermodynamic integration to obtain the probability ratio SE gibGA then 2 Use these configurations as starting points to run unconstrained dynamics and estimate the dynamical term lt5 0 gik0e t giDlt6 0 i Picture from Michael Myers blue moon web page at httpwwwnetaxscommhmyersbluebluehtml 72 NC STATE UNIVERSITY Blue Moon Molecular Dynamics Schematically a Blue Moon MD simulation goes like this 7 Reactants ark Products The thick line represents a constrained MD simulation in which the system is kept at the top ofthe barrier and hence ii 0 This can be done by a standard constrained MD algorithm such as SHAKE or RATTLE 89 The thin lines are unconstrained simulations with ii 0 that start at the top ofthe barrier and evolve to become reactants or products 2 NC STATE UNIVERSITY Blue Moon Molecular Dynamics In order to find the dynamic part of the reaction rate constant we will be sampling over trajectories that start at 0 if It is possible to relate the canonical probability distribution of the constrained system and the real system see eg Frenkel and Smit chapter 15 and find the effect of the constraint on the averages This way we obtain the correction for the bias lt5l lt0gt ail ltoelar ail ltZ 2 0elar ailc Mamfil lizrfgtc The C indicates that the average is done on the ensemble where initial conditions are constrained and the correction factor is Lg lZl lzmi 5L 2624 2 2625 NC STATE UNIVERSITY 7 71 Blue Moon Molecular Dynamics The equilibrium part of the rate constant can be obtained using thermodynamic integration First we rewrite it as sli ail z lewd 6er W P w1th P W We can differentiate 1nP with respect to Sam change variables to a new coordinate system in which i is one of the coordinates eg the eigenbasis of the Hessian This way we obtain 2626 61nPi lt6 iWgtltIZIIZWgtC 26 27 a6 662a IZI lZlc 39 Here J is the Jacobian of the coordinate transformation This equation can be integrated numerically to get PE 7 NC STATE UNIVERSITY 7 Blue Moon Molecular Dynamics In summary a Blue Moon MD simulation goes as follows 1 Run constrained dynamics with values of Z ranging from EA the reactant side to if the transition state Use these simulations to evaluate the averages on the righthand side of eq 2627 Integrate this equation numerically to find the equilibrium part of the rate constant Using as initial conditions the configurations from the constrained MD with i ii run normal MD simulations to calculate the averages on the righthand side of eq 2624 This gives the dynamic part of the rate constant 3 Multiply the two terms and voila you have the rate constant 2 71 NC STATE UNIVERSITY 7 Blue Moon MD Example Simulation of the 8N2 reaction F CH3C1 gt C1 CH3F using ab initio molecular dynamics in the Blue Moon ensemble 11 250 00 7 g 750 g E g 2D0 3 a a E u LIJ I5 4250 3 I 400 4750 600 39 15 05 05 15 25 10 00 20 30 ch on coordinate A Reaction Coordinate A FIG 5 Energy pro le along the reaction path at 0 and 6300 K computed FIG 9 Free energy as a function of the reaction coordinate in the Blue Moon ensemble Energies relative to the isolated reactants at 0 139C7Cl7rC7F K NC STATE UNIVERSITY Outlook Finding reaction mechanisms The ZeroTemperature String Method Calculating reaction rates and defining reaction coordinates in complex systems Introduction The BennettChandler method Blue Moon molecular dynamics Transition path sampling Aimless Shooting Likelihood Maximization Finite Temperature String Method Other methods NC STATE UNIVERSITY Transition Path Sampling The problem with schemes based on the BennettChandler equation is that they require defining a reaction coordinate There are many reactions for which it is not evident what the reaction coordinate is and we need to guess As we saw in the previous class we can check to see if it is correct by looking for hysteresis in the free energy diagram We can find the free energy diagram for example in Blue Moon MD by doing the thermodynamic integration starting from the reactants and the products and then checking for hysteresis 10 However this is an awful lot of work and if we guessed wrong we have to start all over again I 39 lt 60 BACK TO SQUARE ONE 7 Transition Path Sampling In many cases we simply cannot define a reaction coordinate because the process is too complex eg nucleation protein folding We need a method to evaluate rate constants that does not require introducing a reaction coordinate Figure from ref 12 From the BennettChandler method we sawthat the dynamic part of the rate constant can be obtained by sampling only over reactive trajectories If we know what the reactants and products are we only need trajectories connecting the two 4 NC STATE UNIVERSITY Transition Path Sampling We need to write the BennettChandler equation in a way that does not require defining a reaction coordinate One possibility is kf 610 with Cg W The correlation function Ct counts how many trajectories that start at the reactants at t 0 end up being products at time t This still requires us to define what the reactants and products are but not a reaction coordinate The correlation function can be rewritten as a ratio of path integrals 2628 co 0 2629 NC STATE UNIVERSITY 7 71 Transition Path Sampling The integration in 2629 is over the space of all dynamic trajectories This is a bit different from the sampling over configurations that we are used to The Transition Path Sampling TPS method attempts to evaluate Ct and hence the rate constant by doing Monte Carlo sampling on the space of trajectories that connect reactants with products A trajectory is in principle an infinite collection of configurations In practice it is approximated b a discrete set of points in configuration space to 01v 611 zM 215 The probability of a given trajectory 7 N t is calculated differently depending on the dynamic method use to generate it eg Langevin dynamics standard MD kinetic Monte Carlo It is usually written as a product of transition probabilities from one point in the trajectory to the next times the probability density of the initial point NC STATE UNIVERSITY 7 7 Transition Path Sampling In order to carry out a TPS simulation we need several things An initial trajectory that connects the reactants and the products There are several ways to obtain this eg run a simulation at a much higher temperature and then anneal the trajectory obtained to the desired temperature do a minimization of the classical action functional ref 13 use the Finite Temperature String Method etc Choose the length of the trajectories the total time spanned This can be done by trialanderror One way to estimate it is finding how long it takes for two trajectories to become uncorrelated Of course defining the reactants and products regions in configurational space This sounds easy but in complex systems it can be very challenging 7 NC STATE UNIVERSITY 7 Transition Path Sampling How do we sample traiectorv space There are several possible Monte Carlo moves we can use to generate new trajectories The two most common are H NC STATE UNIVERSITY Transition Path Sampling How do we sample trajectory space There are several possible Monte Carlo moves we can use to generate new trajectories The two most common are Shooting moves Change the momenta a little bit at an intermediate point and integrate fonNard and backwards in time From 12 E NC STATE UNIVERSITY Transition Path Sampling How do we sample trajectory space There are several possible Monte Carlo moves we can use to generate new trajectories The two most common are Shifting moves Delete a part of the trajectory at the beginning or end and integrate fonNard backwards in time from the end beginning A B From 12 2 NC STATE UNIVERSITY Transition Path Sampling There are other possible trial moves eg configurational biaslike moves to generate new trajectories So far the combination of shooting and shifting seems to be the most efficient Note that we did not get away with calculating the reaction rate without finding free energies In the equation for the reaction rate in the transition path ensemble both the numerator and denominator are probabilities This can be related to the average reversible work needed to take a trajectory that ends anywhere and make it end in the products region This is in practice obtained using thermodynamic integration in trajectory space 7 NC STATE UNIVERSITY 7 Transition Path Sampling Example Autoionization of liquid water using TPS ab initio MD 14 Fig 1 A through F Snapshots from an au toionization trajectory for 32 water moleculs periodically replicated in space The illustrat ed configurations are separated by 30 fs Dy 39u 1 temperature and with a force field deter mined by the BLYP lectronic density functional theory The trajectory was gener ated by the CPMD al gorithm and was har vested by transition path sampling R d spheres represent oxy gen atoms and white l a d blue gen atoms highlight the hyd and droxid Ions re p c show hydrogen bond wires along which ionic species move relatively quickly 3 NC STATE UNIVERSITY Outlook Finding reaction mechanisms The ZeroTemperature String Method Calculating reaction rates and defining reaction coordinates in complex systems Introduction The BennettChandler method Blue Moon molecular dynamics Transition path sampling Aimless Shooting Likelihood Maximization Finite Temperature String Method Other methods 7 NC STATE UNIVERSITY Aimless Shooting and Likelihood Maximization As we discussed earlier reaction coordinates are often found through trialanderror The problem with the trialanderror approach is that it requires either looking for hysteresis in a free energy profile or generating a committor probability distribution for the trial reaction coordinate both of which are very computationally expensive Alternatively one can generate an unbiased set of reactive trajectories and analyze the data afterwards to determine suitable reaction coordinates This is the idea of the Aimless Shooting Likelihood Maximization LMaX approach 20 21 7 NC STATE UNIVERSITY 7 Aimless Shooting and Likelihood Maximization Aimless shooting generates an unbiased set of reactive trajectories by sampling the transition state ensemble The procedure is as follows 1 Find one point x0 close to the transition state hypersurface and build a reactive trajectory i e one that connects the basins of reactants A and products B by shooting fonrvard and backwards in time from that point up to time T2 and T2 where Tis the total trajectory time 2 Select with equal probability one of the points X0XiA2 where At is the timestep Call this point the new initial t 0 point and shoot fonvard and backwards in time again with a new set of velocities sampled from a MaxwellBoltzmann distribution 3 If the new trajectory connects A and B accept it othenvise reject it and start from the old trajectory Go back to step 2 Store both the accepted and the rejected trajectories 7 NC STATE UNIVERSITY 7 Aimless Shooting and Likelihood Maximization FIG 2 Aimless shooting creates a sequence of interconnected trajectories The three points on each trajectory are separated by a short time Ar 7 NC STATE UNIVERSITY Aimless Shooting and Likelihood Maximization After collecting a large number of trajectories we can test different models for the reaction coordinate by using Likelihood Maximization LMax The procedure is as follows Identify a set of relevant parameters qk that may be important to describe the reaction Postulate a functional dependence of the reaction coordinate on the parameters For example for a linear model quotQ 205qu 050 k where the 06 are free parameters Define the likelihood of the trajectories harvested as B B La tanh tanh r k1 accepted rejected this is only one possible model for the likelihood see 20 7 NC STATE UNIVERSITY 7 Aimless Shooting and Likelihood Maximization The best set of parameters at for the current model can be obtained by finding the values or that maximize the log likelihood a1nLa Let s say we have M parameters qk How do we know how many and which of them we need to get a good reaction coordinate The answer comes from the Bayesian information criterion BIC According to this if Emis the maximum log likelihood of all possible combinations of m parameters then a model with m1 parameters will be significantly better if 6 m1 m gt 1nNR where NR is the number of trajectories harvested 7 NC STATE UNIVERSITY 7 Aimless Shooting and Likelihood Maximization In summary the algorithm to find a reaction coordinate using Aimless Shooting LMax is 1 Harvest NR independent trajectories using aimless shooting 2 Propose a set of M relevant parameters start using models for the reaction coordinate with one parameter m 1 3 For each of the possible combinations of m parameters maximize the log likelihood E Let Em be the maximum of all the values of E obtained 1 4 If m1Em lt 1IINR ofm M stop If m1Em gt 1IINR orm 1 let m m1 and go back to step 3 This process finds the best reaction coordinate for the given set of important parameters 7 NC STATE UNIVERSITY 7 Aimless Shooting and Likelihood Maximization An example polymorph transformation in terephtalic acid 22 WW blur 120 A j 39 gt hemermmmm m Wilma blc100 quot Figure 9 Snapshots ow the hvaxis and minds news along the initial trajectory for the 343 molecule system obtained from MD umbrella sampling These configurations correspond to the PM we shown in Figure 8 From top to bottom the con gurations are taken at bc 160 preuansformatxon bi 140 at the onset of nucleation bc 140 as nucleation ls occunmg at the peak of the PM curve hc 120 as growth is occurring and at We 100 in the form II basin None that the nucleation event takes place at the top le comer as seen from the axis View and propagates to the bottom right comer again on the caxis View a NC STATE UNIVERSITY Outlook Finding reaction mechanisms The ZeroTemperature String Method Calculating reaction rates and defining reaction coordinates in complex systems Introduction The BennettChandler method Blue Moon molecular dynamics Transition path sampling Aimless Shooting Likelihood Maximization Finite Temperature String Method Other methods NC STATE UNIVERSITY Finite Temperature String Method The FSTM can be used to obtain a minimum free energy path between two metastable basins at nite temperature 23 24 The main idea is to construct a chain of states string that evolve in time guided by the negative gradient of the ee energy with respect to some collective variables the mean force in analogy with chainof states methods to obtain minimum energy paths where images are guided by the forces In order to prevent the images from clustering at the stable basins the string is reparameterized after each time step like in the ZTSM The reparameterization involves two steps Fitting a continuous curve to the images after the time step e g by using cubic splines Use the curve to interpolate between images based on some criterion e g equal separation in arc length or energyweighted arc length The mean force at each image can be obtained for example by using a constrained dynamics approach e g Blue Moon MD or restrained Langevin dynamics NC STATE UNIVERSITY 7 71 Outlook Finding reaction mechanisms The ZeroTemperature String Method Calculating reaction rates and defining reaction coordinates in complex systems Introduction The BennettChandler method Blue Moon molecular dynamics Transition path sampling Aimless Shooting Likelihood Maximization Finite Temperature String Method Other methods NC STATE UNIVERSITY Other methods There are many other techniques beyond those discussed here to find rates or mechanisms of rare events as this is currently a very active field of research Some other methods are Parallel replica dynamics This is useful mostly for classical MD simulations of complex systems It consists on running a set of MD calculations in parallel until a rare event occurs Then the simulation clock is advanced by the time taken by all simulations 15 Hyperdynamics Similar in spirit to metadynamics and umbrella sampling this method increases the energy of the reactants to reduce the activation barriers 16 Temperatureaccelerated dynamics TAD As its name indicates this involves running at higher temperature where the rare event is not so rare By assuming a form of the temperature dependence of the rate constant this can be used to obtain the rate at lower temperatures 17 i NC STATE UNIVERSITY References 1 D Masson C Richard R Martin Int J Chem Kinet 8 37 1976 2 D Chandler J Chem Phys 68 2959 1978 3 D Chandler Introduction to Modern Statistical Mechanics Oxford University Press Oxford 1987 4 BH Mahan J Chem Educ 51 709 1974 5 WH Miller Acc Chem Res 9 306 1976 6 CH Turner JK Brennan JK Johnson KE Gubbins J Chem Phys 116 2138 2002 7 G Ciccotti M Ferrario J Molec Liquids 89 1 2000 8 JP Ryckaert G Ciccotti HJC Berendsen J Comput Phys 23 237 1977 9 HO Andersen J Comput Phys 52 24 1983 10 See for example EJ Meijer M Sprik J Am Chem Soc 120 6345 1998 or B Ensing EJ Meijer BE Bl60hl EJ Baerends J Phys Chem A 105 3300 2001 11 M Mugnai G Cardini V Schettino JChem Phys 118 2767 2003 12 C Dellago PG Bolhuis PL Geissler Adv Chem Phys 123 1 2002 13 D Passerone M Parrinello Phys Rev Lett 87 108302 2001 14 PL Geissler C Dellago D Chandler J Hutter M Parrinello Science 291 2121 2001 15 AF Voter Phys Rev B 57 R13985 1998 16 AF Voter J Chem Phys 106 4665 1997 AF Voter Phys Rev Lett 78 3908 1997 17 MR Sorensen AF Voter J Chem Phys 112 9599 2000 18 WN E WQ Ren E VandenEijnden Phys Rev B 66 052301 2002 19 WN E WQ Ren E VandenEijnden J Chem Phys 126 164103 2007 20 B Peters and BL Trout J Chem Phys 125 054108 2006 21 B Peters GT Beckham BL Trout J Chem Phys 127 034109 2007 22 GT Beckham B Peters C Starbuck N Variankaval BL Trout J Am Chem Soc 129 4714 2007 23 WE W Ren J Phys Chem B 109 6688 2005 24 L Maragliano A Fischer E VandenEijnden G Ciccotti J Chem Phys 125 024106 2006 7 NC STATE UNIVERSITY 7 fCHE596M Multi Scale Modeling of Matter Instructor Keith E Gubbins Lecture 16 Phase equilibria bulk phases and confined systems 7 NC STATE UNIVERSITY Outlook Thermodynamic integration solidfluid equilibria Introduction to adsorption and confined systems Motivation Porous materials Gasliquid transitions capillary condensation Layering transitions Freezingmelting transitions 71 NC STATE UNIVERSITY 7 Thermodynamic Integration References Frenkel amp Smit pp 168172 241267 and references therein Allen amp Tildesley pp 4950 and references therein Experiments always determine a derivative of the free energy eg j 2 p 11 6V NT aAT g M W T2 12 Therefore to compute the free energy of a system at given temperature and density we should find a reversible path in the VT plane linking the state under consideration to a state of known free energy Then the change in A along that path can be evaluated by W W of 11 and 12 There are only a few thermodynamic states for which the free energy of a substance is known One of them is the ideal gas phase another is the lowtemperature harmonic crystal NC STATE UNIVERSITY 7 Thermodynamic Integration Let us rewrite 11 in terms of the excess Helmholtz free energy Aexp Ap Aidp where Aid is the Helmholtz free energy of an ideal gas reference state at the same density and temperature as the system under study Then ex id r e r 8V NJ 8V NJ 8V NJ aAex P ka 8V NT 13 14 where ka is the pressure of an ideal gas at the same density as the system under study NC STATE UNIVERSITY 7 7 Thermodynamic Integration Therefore to determine the Helmholtz free energy of a liquid phase with density p we can integrate 14 from an ideal gas phase for which A 0 to the state of interest keeping the temperature and N constant AND V 1 p 1 619 0 P kT dV P kT NkT iNkT p lax p pz Aex A Aid 1 P P P P pJ39 Edp 15 NkT NkT NkT JO 9 An important condition is that the integration path in 15 should be reversible If the integration path crosses a firstorder transition hysteresis may occur and 15 can no longer be used NC STATE UNIVERSITY 7 Thermodynamic Integration For a liquid phase at pand T this problem can be avoided by performing the integration in two steps First at a temperature T well above the critical use 15 to integrate along this isotherm from the ideal gas phase to the supercritical fluid state at p and TH id P I ApT NkT A p 1 PTPfT jab 16 0 NkT H 0 Then integrate 12 along a line of constant density cooling the system from T to the temperature of interest T ApT AlpT l T U T T11 IFC T 7 TH 7 NC STATE UNIVERSITY 7 Thermodynamic Integration It is not really necessary to go all the way to the ideal gas One can start from a state p1 T1 for which we can compute the free energy In practice p1 should be a sufficient low density to ensure that the free energy can be computed accurately eg by using the Widom method to determine the excess chemical potential Then one can determine the Helmholtz free energy of a liquid phase at density p2 and temperature T2 by integrating 11 A A a wjdV VyJ NkT 2 Nle V1 NkT V kT p2 p1 18 Again the integration should be done along a reversible path linking the States P1 T1 and P2 T2 NC STATE UNIVERSITY 7 SolidLiquid Equilibria and Thermo Integration The solidliguid phase equilibria problem requires different approaches For a review see eg P A Monson and D A Kofke in Advances in Chemical Physics ed Prigogine and S A Rice vol 115 pp113179 Wiley New York 2000 Also see Frenkel amp Smit ch 10 The main problem is that the solidliquid coexistence curve itself does not end in a critical point and hence there exists no natural reversible path from the solid to the ideal gas that does not cross a firstorder phase transition It is usually possible however to use thermodynamic integration and construct a reversible path to other states of known free energy The basic idea of one of the methods is to transform the solid under consideration reversibly into some reference system for which the free energy can be evaluated analytically For a brief description of the reference systems used for solids see the review paper by Monson amp Kofke i NC STATE UNIVERSITY SolidLiquid Equilibria and Thermo Integration One of the most widely used reference systems is the noninteracting Einstein crystal 0X yo gt40 Om x marks positions of the noninteracting X particles which are attached to sites in a O 0 O 039 lattice indicated by the small circles via X Hookean springs indicated by the dotted XE xiv lines A simple square lattice is shown for gt9 9 q illustration X X Figure taken from P A Monson and D A Kofke D X9 gt3 X9 t l in Advances in Chemical Physics ed Prigogine X and S A Rice vol 115 pp113179 Wiley New 3 x 3lt o X York 2000 Suppose we have an atomic solid system where the atoms interact with a continuous potential MEN We shall use thermodynamic integration to relate the free energy of this system to that of a noninteracting Einstein crystal J Q Broughton and G H Gilmer J Chem Phys 79 5095 1983 D Frenkel and A J C Ladd J Chem Phys 81 3188 1984 SolidLiquid Equilibria and Thermo Integration During the thermodynamic integration we switch on the spring constants and switch off the intermolecular interactions To this end we consider a potential energy function Morale lt1 2gtTUKNgtur v Mice 21 ill2 am where 110 is the lattice position of atom z and LIQON is the static contribution to the potential energy ie the potential energy of a crystal with all atoms at their lattice positions A is the switching parameter and a is the Einsteincrystal spring constant coupling atom ito its lattice site Note that for A O we recover the original interactions For A 1 we have switched off the intermolecular interactions completely except for the constant static term and the system behaves like an ideal non interacting Einstein crystal 7 NC STATE UNIVERSITY SolidLiquid Equilibria and Thermo Integration The partition function for a system with a potential energy function that corresponds to a value of 2 between 0 and 1 is LNszN exp Zll 20 QNVTl Am The derivative of the Helmholtz free energy AOL with respect to A can be written as an ensemble average Mg JDLM 1nQNVTA N aupl 1 aguvym laxityfxpl w QNVT1 61 IdzNeXp 8Ll 6A A 61 A NVT ltgtl 7 NC STATE UNIVERSITY 7 SolidLiquid Equilibria and Thermo Integration where denotes an ensemble average for a system with a potential energy function 011 such as 19 Therefore the free energy difference between the systems with A 0 original system and A 1 noninteracting Einstein crystal can be obtained by integrating 21 At 0 2 Am mfg0X am ldrliaxnweir turw uri The configurational free energy of the noninteracting Einstein crystal is given by d N where Am 2 HQ Zln L d dimensionality z 391 ai 22 23 7 NC STATE UNIVERSITY SolidLiquid Equilibria and Thermo Integration Important implementation details of this technique which is called the FrenkelLadd technique vary with the nature of the system being studied including whether the intermolecular potential is continuous or discontinuous and whether the molecules are monatomic or multiatomic See Frenkel amp Smit ch 10 for a thorough discussion and implementation details for each one of these cases 7 NC STATE UNIVERSITY Phase Transitions in Confined Systems Introduction Review papers on Phase Transitions in Confined Systems L D Gelb K E Gubbins R Radhakrishnan and M SliwinskaBartkowiak Rep Prog Phys 62 15731659 1999 C AlbaSimionesco B Coasne G Dosseh G Dudziak KE Gubbins R Radhakrishnan and M SliwinskaBartkowiak Effects of Confinement on Freezing and Melting Journal of Physics Condensed Matter 18 R15R68 2006 All figures and tables presented in this lecture came from these references unless othenNise noted Motivation Molecules confined within narrow pores can exhibit a wide range of physical behavior Wall forces and competition between fluidwall and fluidfluid forces can cause Phase transitions not found in the bulk eg Layering Wetting Commensurateincommensurate transitions 7 NC STATE UNIVERSITY 7 Introduction Motivation Shifts in transitions that also occur in bulk systems Liquidsolid Freezing Gasliquid Liquidliquid Understanding of these phenomena is necessary for many industrial and geophysical operations Pollution control mixture separation catalysis Lubrication adhesion Tertiary oil recovery gas field technology removal of pollutants from ground water and soils frost heaving Fabrication of nanomaterials 7 NC STATE UNIVERSITY 7 Introduction Porous Materials able 1 Nmmpmous mmcrials Dun lakcn in pun from m7 Pore uhh Mamm Sui1m Pom gt1me mm A Cry Mina mgulzu Aluminusilicmcs 0 AL Si Cyhmch cage 0340 Ahimilmphosphmc 0 AL P Cylnulcl 0843 MCMAI o ALS Cylinder 1 11 Carbon uanombc C Cylinder 271 U B Anmrphaus Porous gm 0 5 ylindcx Sil39 Dudes 0 3 cit Cylimlcr 0 5 oid 0mm SmpHurs B N H 5m C 5m C SH mid Figure 2 Scanning eleclrnn micrograph 0f 3 CFO of 300 nm mean pore diameter From 128 Figure 1 Transmission electron micro 39lph TEM01 MCM 41 showing uniform pores of hexagonal cross section From I49 39 rug17E w WW 3 Introduction Porous Materials Figure 3 TEM of a carbon aerogel magni cation Figure 4 TEM of a pitchbased activated carbon bre 2200 000 From 132 White areas are poreg grey areas are pore walls dark areas are thicker pore walls 271 NC STATE UNIVERSITY MC Flag E I LLMWER 1 Y From 2 increasing numbers ol concenl ric tubes one to ve layers in Figure 6 Highresolution TEM images of carbon nanotubes with uHc anthracene soul sample From 277 1 image analysis of a graphitized Figure 5 Hi2 l39l 1 63 olution TEM imn 2 e before a and after 10 nm 5WW j mgmg UyWWE pmnnk gwwwgwme mg 4Wuy f r j I My 4 x x fill 1H amn ymw9 Introduction Simulation and Theoretical Methods Both experimental and theoreticalsimulation work in this field face significant difficulties Experiments Accurate characterization of porous materials is difficult Molecular structure Morphology Distribution of pore sizes and shapes Pore connectivity surface chemistry Metastability usually can tdon t get free energies Changes in surfaces and pore structure due to temperature and pressure changes Trace amounts of impurities in the adsorbate may preferentially adsorb on pore walls NC STATE UNIVERSITY 7 1 Introduction Simulation and Theoretical Methods Theoretical and simulation methods For amorphous porous materials don t know exact morphology of the material Mimic fabrication process of real material in simulations Build a model matching experimental structural data of real material Metastability Determine free energies Uncertainties in intermolecular potentials Walls are often assumed to have a negligible influence on fluidfluid potentials Importance of electrostatic induction three and many body interactions is often unknown Pore swelling can occur in some cases on adsorption NC STATE UNIVERSITY 39 7 Introduction Simulation and Theoretical Methods Simulation and theoretical methods commonly used Grand Canonical Monte Carlo GCMC Gibbs Ensemble Monte Carlo GEMC Semigrand Monte Carlo SGMC Quench molecular dynamics Histogram reweighting Landau free energy methods using bond orientation parameters Lattice models Density Functional Theory DFT NC STATE UNIVERSITY GasLiquid Transition in Pores Capillary Condensation Extensively studied by experiment DFT and molecular simulation Transition occurs when T is below the pore critical temperature and pore width is greater than a few times the adsorbate molecular diameter Transition appears in adsorption isotherm as a sudden and large jump in the amount adsorbed at the capillary condensation pressure All pores of same size and shape homoqeneous walls sharp transition sudden vertical jump in amount adsorbed Real system steeply rising adsorption continuous filling due to distribution of pore sizes and shapes pore connectivity and wall heterogeneity Usually accompanied by hysteresis desorption starting from pores filled with a dense fluid phase occurs via a different path than adsorption 7 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation Experimental studies DH Everett eta1I Reversible region at low pressures Sharp increase in adsorption for a narrow range of pressures on E E E n accompanied by an hysteresis loop Capillary condensation pressure increases with T Hysteresis loop becomes narrower as T increases disappears at some hysteresis critical temperature Tch below the bulk Tc anmol 9 1 Hysteresis arises from metastability Difference in hysteresis loop shapes due to pore shape lquot i g 4 m um um 202K Zl UK 22 v 233K 253K 242K 255K 273K 3 0 1 0 log tatrn a I 1 wax 213K 233K 242K 25 u 39 7 Adsorption isotherms For xenon plotted us mole Idwibcd Pk 39 of ELM o 155 20 2n fatm b adsorbent crsm lo IClI ol xenon in thc bulL gm plmxc L1 in Vycot39 glass II in urine carbon From 4h GasLiquid Transition in Pores Capillary Condensation Experimental studies 2 WW 250 250 HystereSIs phase 240 diagram gt plot amount 230 g o T p adsorbed at the extremes 200 of hysteresis loop vs T 20 quot I I 190 l70 W I HystereSIs vanishes 30 40 5 130 m 150 150 70190 above the hysteresis n lmmol 91 nvmmolgquot critical temperature Tch Hrai L IV f m f H w xiciwupmc m m in gun cmi vi Different shape due to Lquotl3l ill 1l39sElli irimtutti different pore morphology W m When plotted in reduced units 39 o 6 hysteresis phase diagrams for is XeNycor and COZNycor coincide E 0839 Within experimental error figure 9 I igurc i L ompzu39ixuu ul lnxlcrceix phuxc diugmmx for 035 39 39 39 13 10 o mm and X9 new in reduced mm mm 40 nlnclH GasLiquid Transition in Pores Capillary Condensation Exgerimental studies GH Findenegg 32 bulk coex Comparison of hysteresis With m bulk phase diagram W coexistence curve in pore is 36 3 narrowed hysteresis critical point is depressed relative to 3394 39 24am bulk critical point 32 Density of liquid phases is similar in pore and in bulk m M napWWW 390 392 higher denSity for the gas Figure 1 Elpllll39lk1llltl cucxrstcncc CLIHC for bulk phase in the pore due to xulfur hexu uoridc and ll3939tC39Cl39 phase diagrams for SE in CFO of mean pore adsorptlon of dense layers 0 diameter 3 Ind 24 nm From 354 the pore walls gt higher average values forp Shift in hysteresis critical point with pore width H is consistent with ATC bH y y N 25 gt ch changes with pore size GasLiquid Transition in Pores Capillary Condensation Simple pore geometries means Regular shape slit cylindrical Homogeneous surfaces Pores not interconnected Theoretical macroscopic treatment the Kelvin equation Valid for sufficiently large pores For a pore of general shape the exact differential of the grand free energy Q is given by d z SdT Pde ltNgtdyydA 1 where S is the entropy Pb is the bulk phase pressure y is the solidfluid interfacial tension and A is the surface area 1 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation Theoretical macroscopic treatment the Kelvin equation For slitshaped pores 1 is often written as do SdT 19de ltNgtdy ZydA Ade 2 where A is now the surface area of one of the walls H is the pore width andfis the solvation pressure PH Pb PH is the pressure exerted on the pore walls by the adsorbate For sufficiently large pores oscillations in the density profile due to fluidwall forces can be ignored By integrating 2 using Euler s theorem we get 2 PbV 232A 3 1 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation Theoretical macroscopic treatment the Kelvin equation By writing equations for Q of the confined gaslike and liquidlike phases assuming that the liquid wets the walls completely contact angle is zero and equating the grand free energies it is possible to get an equation for the pressure at which capillary condensation occurs the Kelvin eguation 2 1350 RTPIH where Pl0 is the saturated vapor pressure of the bqu fluid 7g is the gasliquid surface tension of the bqu fluid R is the gas constant and 0 is the density of the bulk liquid Kelvin equation is valid for large pores and for temperatures well below the pore critical temperature 7 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation Macroscopic treatment the Kelvin equation Assumptions made in the derivation of the Kelvin equation Density oscillations due to pore walls can be neglected System is large enough that a surface tension can be defined and its value is that for the bulk fluid at this temperature Gas phase can be treated as an ideal gas Liquid phase is incompressible Therefore the Kelvin equation is valid only for large pores and temperatures well below the pore critical temperature Modified Kelvin equations have been proposed to take into account strongly adsorbed layers on pore walls nonideality of the gas phase liquid compressibility and nonzero contact angle 7 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation Theoretical macroscopictreatment too the Kelvin eguation Inaccurate for pore sizes below 10 7510 nm for small molecules such as N2 when liquid wets the 0 surfaces Kelvin equation overestimates capillary condensation pressure for a given pore size errors become larger for smaller pore sizes N o amp l1 6 w tttmut tttttttt lllllll tittttii tttttttt llllllnl tttttttt lllllllu tttttttt tttttt 0 1039 Hui iigtttcis 139 t i quotquot GC Homih sznzoc quiulion Hlxl DH tutti Inolccllldl qmuiunott liuliihl int i 31mph mOtlL ni ntitogcn adsorption in L lti cattlmit pate u 77 K Hctc II is dk llCLI l the HAIRY beittccn 1 mt Lict39 at cutit iii the tittncx Mitt lIu Ullgh the mum untit ti imit mutt imittttt 39i i wolf n tilt Dmimimh it FT and Simulation Itsulh Sim ntptiiitt minimisingquot rat the l 7 i c In llh mini Il NVIHL 1mlc math The huip tutmm in llh clll w int 1 tttoiccttittt IL ing aim rim i 107i GasLiquid Transition in Pores Capillary Condensation Slit pores theory and molecular simulation pure fluids Methods DFT Lattice gas model calculations Integral equation theory Molecular simulation Adsorption isotherms from theorysimulation are qualitatively very similar to those obtained in experiments For mesopores H2 5O nm and subcritical temperatures As the pressure increases layers of adsorbate first form on the pore walls rapid increase in adsorption At higher pressures there is a sudden discontinuous jump in adsorption capillary condensation transition After transition pores are filled with an incompressible liquid quotke phase a GasLiquid Transition in Pores Capillary Condensation Slit pores theory and molecular simulation pure fluids For mesopores and subcritical temperatures Since pores are monodisperse the jumps in adsorption are sharp discontinuities Starting from filled pores and decreasing the pressure capillary evaporation is delayed and occurs at a lower pressure than that observed for condensation Hysteresis effects are due to metastability As we increase the temperature capillary condensation moves to higher pressures and hysteresis loops decrease in size until they eventually disappear Can determine the true thermodynamic transition by thermodynamic integration see eg B K Peterson and K E Gubbins Mol Phys 62 215 1987 L D Ger eta Rep Prog Phys 62 1573 1999 and references therein 7 GasLiquid Transition in Pores Capillary Condensation Slit pores theom and molecular simulation pure uids At subcritical temperatures the capillary condensation pressure decreases as H decreases Pore critical temperature Top also decreases as H decreases For small pores there is no longer a capillary condensation pressure it is replaced by a steep but continuous filling 08 06 Q0 j l M 10quot 10394 103 102 10 10 0 PP Fing m Adwlplion modicum m L mlch cfninngcn m mmlmped gluphmr Calb n pom th llll 3 nulls m H K 1mm m L unu l 39m of Dl39 I Raiding lmm ion In l39lglt the pore ldllh luc ML l L x ll l4 and 4th lhcrc n is the LJ lllllmctcl nftlic Nmmlccul 39 a l 39l l From 1 r NC STATE UNIVERSITY GasLiquid Transition in Pores Capillary Condensation Slit pores theory and molecular simulation pure fluids Capillary condensation is very sensitive to the pore width This is the basis of commonly used methods for measurement of the pore size distribution PSD of porous materials Assume pore structure is represented by a collection of non interconnected pores All pores have same geometry slit cylinder and surface structure but they can have different pore sizes Total adsorption F can be represented by F jpPHfHdH 5 where pPH is the mean density of fluid in the pore of width H and pressure P H is the PSD and the integration is over all possible values of H 7 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation Effect of system variables on adsorption behavior can be studied in a systematic manner using DFT or simulations l igurc 19 A and rcl ur lo quiililuinc clmnu 39 llllc l1lt f a leaol39l fl ll re Clasxi cmion of m m l tloll isotherms for xiii porclt for mfmj 5 w 39I39 ix 139 in the iiilxni pnoii belim ioiir Si Fri 1 is bola Ihc wetting tompcratui cl Slipchcript I I B MlCROPORES NON POROUS III IE 390 WEAK E SUBSTRATE 396 w a c mssuponss 3 CAPILLARY con o DENSATION E E Y If 3 lt1 WEAK SUBSTRATE LAYERING Relative pressure PP Figure 1X The six muss ol39 stui39plmn imihcrm according to 11 1985 lUPAC Cldssl cillmll From 14 li om DFT The lint in thc iilc L st 39supcrcriiicul 1739 is IbOlC tlic pore n iiicul lcmpciiiliircl From 1 l 10402 ml 1 C zii39c ippruxiinmc l l5 3 IE 5 n E B 2H m I m is 2 5 i S Y 4 T OB IZ Vlf v H Vii T 6 iv W E 0 ol PP b n l l 02 o 06 sf lf 2 J r 50 u H p 6 I l C 0 1 mcdm ISC Ic 0392 sf ff 0394 GasLiquid Transition in Pores Capillary Condensation Cylindrical pores theory and molecular simulation pure fluids Confinement effects are greater than in slit pores since fluid is now confined in two dimensions rather than one Capillary condensation occurs at lower pressures and TCP is also lower than for slit pores of same width For narrow cylindrical pores the system approaches one dimensional behavior and so cannot exhibit true phase transitions at temperatures above 0 K Any firstorder phase transition in a cylinder must be rounded by finitesize effects In both simulations and DFT for pores as small as R Rc7 25 capillary condensation transition appears sharp Single domains of gas and liquid phases are not observed inside a cylindrical pore but rather a series of alternating domains of each phase average size um 7 GasLiquid Transition in Pores Capillary Condensation Mixtures in pores of simple qeometry theory and molecular simulation New feature is the change in composition due to confinement pore walls can attract one of the components more strongly than the others Change in composition of confined phase is often expressed as a selectivity bbulk ppore S prpr 6 l xbixbj NC STATE UNIVERSITY GasLiquid Transition in Pores Capillary Condensation Mixtures in gores of stmgle geometm theom and molecular simulation mun zz Sclccliiny or cthunc m u were mm 11 uml mcthuncl 1 uinmcm mixtth in a slit mam pole at Iquot lrm 15 w Rmm arc shoun For Um promu x l w t I pm 395 71H lsolitl cum ml 02 Lithhcll curvcl Vcitmaljumps m ducmcaplllary cm Llcnmuon From 340 I v I 0 00 002 00a o 05 006 0 o P Figurr ZJ Saleem ll Isollmms llu39 llic mulk lcm J m gurr or z x and l l i urn AH C alumnamum l mm g NC STATE UNIVERSITY Mixtures in pores of simple qeometrv GasLiquid Transition in Pores Capillary Condensation Buildup of 11 classification of selectiVity isotherms multiple layers Class I found at temperatures well above Top 0 Class II corresponds to temperatures slightly above Top Capillary n Layering 1v condens Class III occurs at temperatures somewhat below Top and exhibits N p capillary condensation Class IV found at temperatures well below Top where layering transitions filj ft f 2 f m mm can occur in addition to capillary condensation NC STATE UNIVERSITY GasLiquid Transition in Pores Capillary Condensation Complex pore geometries theory and molecular simulation Simple pore geometries omit several factors which make it difficult to make direct comparison with experiments Porepore correlation effects Surface chemical heterogeneity Geometric heterogeneity variations in pore geometry networking of pores etc Porepore correlation effects Systems confined within very narrow cylindrical pores one dimensional systems cannot exhibit true phase transitions at T gt OK However if several narrow cylindrical pores are close together a gasliquid transition can be observed at low temperatures See eg L D Gelb eta Rep Prog Phys 62 1573 1999 and references therein 7 7 GasLiquid Transition in Pores Capillary Condensation Complex pore qeometries theory and molecular simulation Chemically heteroqeneous surfaces Most studies in slit pores two main types of chemical heterogeneities Pore walls with stripes of strongly and weakly attractive regions gt behavior depends on dH d stripe width Place individual chemical groups on the pore walls gt behavior depends on nH n site density on surface See eg L D Gelb eta Rep Prog Phys 62 1573 1999 and references therein 7 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation Complex pore qeometries theory and molecular simulation Geometric heteroqeneitv Most studies have used independent pores of regular shape Most amorphous porous materials have interconnected pores that may be tortuous and of varying size Realistic models of oxides and glasses have been developed eg Langmuir14 2097 1998 Simulate manufacturing process eg spinodal decomposition templating solgel process Use quench molecular dynamics lattice MC simulation 7 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation Complex pore 39 theory and molecular Geometric heterogeneity quench time D T 225 I 390 r 450 I Pore size is controlled by varying length of quench period Porosity is controlled by varying initial composition D nclt molecular a 1 28 Preparation ofmodcl porous r rm Quc r mixture ltopl ptoducc a networked attucmtc Snupxltots tttltcn at in 33 nm average pore size in simulations tor a 39 l m t to plmc tr 1 Simulation 1 m plum 39 39 is cqtrimlcunoo781 psi From 1 l0 GasLiquid Transition in Pores Capillary Condensation Complex pore qeometries theory and molecular simulation Geometric heteroqeneitv Those porous glass models can be used to test the accuracy of commonly used characterization methods BET Brunauer Emmett Teller gt surface area Accurate for pore diameters larger than 5 nm overestimates area for smaller pores 2025 for pore diameters of 2 nm BJH BarrettJoyner Halenda gt PSD Gives PSDs that are narrower maximum pore diameters are 1 nm smaller than reality See results on next slide 7 NC STATE UNIVERSITY 7 GasLiquid Transition in Pores Capillary Condensation 60 008 000 A quotAll 391 MB E r g 004 0 m n 0 o to Excess adsorption mmolg 0 000 39 0 1 2 3 4 5 6 7 7 008 00 A 006 939 3 004 00 D 00 02 04 06 08 10 g PF relative pressure 002 Figure 2 Adsorption ilt0rhcrms for 1 model 000 ofnnrogcn in the moch glilrmLK depiclcd in 0 1 2 3 4 5 6 7 0 1 2 a 4 5 6 7 gure H017 K Jgturpnon dnlu are 5110 n diameter nm diameter nm 15 Cll39Clk a39 dcsm39plion dam shown 0111 For rho A and D 3 mm as 11 Llul39Cs The isothcrrm arc c m Figure 30 PSDs from geometric analysis wing Spherical probes gt0le cum 23 together W1th PSD Chum and lines are a estimzncs from BJH analysis Llushcd CLH39VCS of the nitrogen adsorption isotherms of gure 20 Uncci unnlm 11 the Ll mullcr than the From 1 10 symbols Glilabcgt A B C and D had mcmgc pore 31203 of approximately 33 3 45 uml 50 nm all glmsm haul u poms CI time to 30quot From 110 5 NC STATE UNIVERSITY quot dlx39plncerl For I Arbitrary Units Layering Transitions At relatively low temperatures adsorption isotherms on smooth surfaces can show stepwise behavior Adsorbing gas prefers to complete each successive monolayer before beginning the next one Steps in isotherms are often quite sharp gt these systems undergo a series of firstorder phase transitions as layers are added Typical experimental adsorption isotherm on planar surfaces Ar on graphite Figurr 35 Ellllvwmctric cmcrugc Isolltcrmx for argon on graphite at our Imiw unlrc The abscissa l 1 function al the ICLIHCCLI 7lcgtll39t gncn w M In quot nipmlmll Mth ilic llwqml lcr r l V r u n i ll mu 39 chimluuml from 3 14 lth 12 ILL UOIquot 3 KlVa NC STATE UNIVERSITY Layering Transitions For a stronglyadsorbing surface surfacefluid interaction is large compared to fluidfluid Adsorption at low temperatures occurs via a series of layer wise firstorder transitions Film thickness becomes infinitely large at liquidvapor coexistence Each of these transitions has an associated critical layer temperature at which the layer transition becomes smooth These layer critical points approach the surface rouqheninq temperature at liquidvapor coexistence Above such temperature there cannot be layerwise phase transitions See Pandit et al Phys Rev B 26 5112 1982 for more details 2 NC STATE UNIVERSITY 7 Layering Transitions Another surface phase transition is the prewettinq transition which corresponds to a discontinuous but finite jump in the adsorbed layer thickness typically of several monolayers size Such transition is first order with an associated critical point and occurs for some intermediatestrength substrates In porous materials we can find adsorption isotherms that consist of Layerwise regimes at low pressures Smooth multilayer adsorption at intermediate pressures Capillary condensation at higher pressures However nontrivial geometry and surface heterogeneity of most porous materials tends to smooth out layerwise behavior 7 NC STATE UNIVERSITY 7 Layering Transitions Layering transitions have been studied using molecular simulation in a number of geometries see L D Gelb et al Rep Prog Phys 62 1573 1999 and references therein DB I 08L Several layertransitions i 39 occur at high pressures l j i 3 within the capillary Nquot a 39 quot396 l i condensation hysteresis a i i a i 2 loop 0 0 f Number of layering FT 3 transttions depends on pore c2 1W 02 FT 352 7 widthH m M Competition between 0 4 g A A 0 4 7 L layering caused by surface V and capillary condensation caused by confinement licon Figurv 36 GCMC results m i lsmption isolhcrim or mcthzmc while slit pom in 39r 05 I e i when 39r ATe Jud 11 7 o iimrpiion ix lcnorctl ind Human transitions 4x dashed line FT is the livczing by cirrlcx Llcmlptlon b l transition chimluccd mm NC STATE UNIVERSITY FreezingMelting Transitions Macroscopic treatment the GibbsThomson equation A classical thermodynamic treatment can produce the analogue of the Kelvin equation for the solidfluid transition See D Nicholson and N G Parsonage Computer Simulation and the Statistical Mechanics of Adsorption London Academic 1982 o R Evans and U Marini Bettolo Marconi J Chem Phys 86 7138 1987 The GibbsThomson equation is then obtained ATfZTfT J 2sz szv 7 T T Hl where T9 T are the pore and bulk freezing temperature yws and ywl are the corresponding wall solid and wall fluid surface tensions v is the molar volume of the liquid 1 is the latent heat of melting in the bulk and His the pore width 7 i FreezingMelting Transitions Many of the early experimental studies of freezing in porous media employed silicabased materials and showed A depression of Tf the pore freezing temperature relative to the bulk fluid As H decreased Tfalso 200 A decreased Different structures forthe 00 confined phase near pore Tff fgeIgass walls and center of pore g m A However increases in T have E K 043555 been observed experimentally 4007 66 OCTMS on mica surface ng force apparatus gt603 y CCI4 benzene on activated 03900 m Hquot In W 0 carbon bers Figurc 45 Shift in l rcczmg tcmpuidtllrc from experiment 0 in solegc CPO 304 q clolmwnc bot ccn mica plums in u u r lcc force apparatus I N The thlxhcd lmcs Lll39C 1 gulth to 1110 c c t 370 indium in Ch 1n 39I G 35 5 FreezingMelting Transitions There is experimental evidence that the GibbsThomson equation breaks down for small pores where the confined phases are highly inhomogeneous 1500 7 I 1000 39 II I M oQE I l I a 500 033 I 39 Q 5 a use snica 40quot r o I DSC Activated Carbons 9 c Simulation Activated Carbon 1ooo 39 7 000 050 1 you 1 o H nm FiguIc 46 51m m rm n Human 11 L39 l4 mm WWW uml Ilnululiun Fur lhc unh mlsoxbcd unl m It the man at the hm dcpenlls on In nmm oi the uid all mcmuon m l nbxcncd ll the me nl mmcmc gldph t mm and A ncgauxc hin nbscncll m tlcncc um 6 union mm Immiil m AT mun x 1 p0 c lhc um ui pomm NlIICtL The plni in mom I II in me Inlciopm oih icgimu no 13053 NC STATE UNIVERSITY FreezingMelting Transitions 39 39 39 39 39 39 39 39 39 Phase diagram for CO2 in 1 Vycor glass with a mean E pore diameter of 7 nm obtained by positronium g s annihilation spectroscopy 3 Liquidsolid transition shifted 3 4 to lowertemperatures and gasliquid transition shifted 2 to higher temperatures than 00 23905 I 20 25 39 220 I 2 25 I 20 I 25135 I 24 in the TemperatureK mrc 47 mm Imgl39um Dru3 m um39 gm the mm m the ma cum a rctbi m the imik I r in 11w con rth pow icm Dmhct Iquoti m 11 ml nu ymtmk to i and the Mud mum 4bctlctP IpotcmndcnsatmnlnndPqurchcczmu n mm c Bt dcnnm cm w the hull H39iplc pnim uluic PT milm in 111 39quaar39 mpic mun oi the pmc Hind Fun 7 NC STATE UNIVERSITY FreezingMelting Transitions Studies of freezing of LJ methane in slit pores via GCMC simulations Miyahara and Gubbins KE J Chem Phys 106 2865 1997 showed Magnitude and sign of ATf Tfi T mlk depended on the strength of the 10 adsorbatewall attractive interaction 33 Aquot 39 quot 7395 For weakly attractive walls eg 0 9 mm silica Awaas negative as found in g 0 most experiments Ls However for strongly attractive walls 5 M11 eg carbons Awaas positive SE 07 N For walls in which the fluidwall interaction was similar to the uid 06 39 fluid interaction and the density of wa atoms was similar to the 0560 8390 100 go 140 adsorbate density little or no change TemperatuxeT m in waas observed Flo 1x L39l fcu or pmc wall on li cczing bclmuoi mutual m mcrull 75 llCllxllV clmngc against lcmpcraturu llur mic itlllix ur ll FreezingMelting Transitions Thermodynamic integration breaks down for attractive pores To overcome this difficulty we can use a method based on an order parameter formulation to calculate the free energy from GCMC simulations in a slit pore Radhakrishnan and Gubbins Mol Phys 96 1249 1999 Method relies on calculation of the Landau free energy as a function of an effective bond orientational order parameter CD The Landau free energy can be defined as ACD kT1nPCD constant 8 where the probability distribution function P is calculated in a GCMC simulation with the help of umbrella sampling The grand freeenergy Q is then related to the Landau free energy by exp 139 I dlt1gtexp BAltD 9 7 NC STATE UNIVERSITY 7 Order Parameters 2 d D lNb 2k eXpi69k Mermin 1968 k1 Nb sum over nearest neighbor bonds at position r y 00 9 Landau free energy 00 OO ADkT1nPD Orientational pair correlation G6r G6ltrgt lt ltIgtlt0gtltIgtltrgt gt 7 NC STATE UNIVERSITY FreezingMelting Transitions Slit pore with two molecular layers of adsorbate a each layer behaves as a twodimensional system Three phases can be observed 9 04 Gel Liquid positional and 30 03 orientational disorder h Hexatic positional disorder g 0 and quasilongrange T34o K 04 orientational order 10 02 WNW 00 00 39 Crystal quasilongrange 50 0 positional order and long range 3 W orientational order iig y mm 0 2 00 00 Meltmg of a 2D crystal occurs 00 100 200 300 00 1026200 300 n H l39 6 VIa the KTHNY Kosterlltz ThoulessFHalperinNelson Fro l mm and calyx123Lilian l 39 Slxll k lb Cllllc l ll l plmc it 1140 K c anmllinc l quotl phaxc Young mechanism D R m 729 K Nelson and B Halperin R Radhakrishnan K E Gubbins and M SIiwinska Bartkowiak J Chem ys ev 3 Phys 116 1147 2002 mn ncd in a 150 i T 5 E o 9 FreezingMelting Transitions CCI4 in ACF Hexatic phase is the thermodynamically stable phase at 335 K At 290 K the hexatic and crystalline phases coexist o 02 39o4 390e o 3902 3904 3906 Liquidhexatic transition 761 V39 vl occurs at 347 The lm onlci Aimmutton inmons or the Landau rm cncl g rot Hexatic39cryStal transmon ill 1n5 K m Tl lll K for LJ cl in a gl39aplmc pni c llquot occurs at 290 K Both transitions are first order at least forthis 15000 L particular system size 5 H a 5000 5000 R Radhakrishnan K E Gubbins and 25 390 2900 3quot 3700 M SliwinskaBartkowiak J Chem Phys 116 1147 2002 c tum cc on r function of temperature 1 liquid 1 3 FIG 3 The l A ltchlir Ami cigsmllinc pit fol LJ 1 in a graphite mic II 7 NC STATE UNIVERSITY Tm I39l39mk FreezingMelting Transitions 17 15 L 12 l2 3 quotJme 09 09 C 06 06 1 I7 n a I I I I I I 0 l 2 3 0 1 2 3 l FIG H II Global phase IIILIgIAm nfzI uid In Ll Iit pore ofw idth I 3 UN from imIIlation Three different 39IIC 01153 i I I CVHIF IIIHiIIc I The d DUI I of 1 mm in u 31 pm of mm 117 aquot ram I mm The cmumcms are for aI39IoIIs IIImrbmcs c n ned 1min acIIVutcd when bc LAC mean pom widm IA mm I H10 Rd I3 2 5H5No2 IR III 3 c5 5 H2 c I5 Ref13y5 l4 Rm H lb 6 5He Ru 1 o R Radhakrishnan K E Gubbins and M Sliwinsk Phys 1 aBartkowiak J Chem 16 1147 2002 NC STATE UNIVERSITY Concluding Remarks Capillary condensation and layering transitions have been extensively investigated Much work being done on freezing in such systems a clearer picture is beginning to emerge In experiments we need To synthesize wellcharacterized porous materials with desired and known features pore size and shape surface chemistry connected or unconnected pores etc Much work now in this area mesoporous carbons highly ordered silicas organosilicas MOFs Improved and more accurate methods for characterization of porous materials In modeling we need To develop better models for porous materials and to use these to gain a better understanding of phase equilibria in real materials Improved intermolecular potentials 7 NC STATE UNIVERSITY 7 Concluding Remarks Difficulties found in experiments and simulations are of complementary nature However it is difficult to bring experiment and theorysimulation together in a fruitful way Difficulty of characterization of materials on the experimental side Difficulty of developing adequate models for use in theory and simulation Situation will improve as increased computational power makes more sophisticated models a reality and as new materials become available Only covered simple fluids small neutral molecules adsorbed on wellcharacterized porous materials Complex fluids were not covered ionic fluids polymers colloidal suspensions amphiphiles etc Not covered liquidliquid transitions wetting solidsolid transitions 7 7 3 NC STATE UNIVERSITY I N From 240 GasLiquid Transition in Pores Capillary Condensation Experimental studies Hysteresis critical temperature Tm temperature at which hysteresis disappears Pore critical temperature Top temperature where the sharp jump in adsorption due to capillary condensation disappears There is evidence that these two temperatures are different with TCPgt Tm see figure 13 7 Top also changes with pore size 0 ml the 1mm Iopc who in Kiliccmis MC Mull Hrlt39 slnpcathighci pmsummlimumlonmcmm mh I Inverse slope arb units ca 0 I m o Arrp1 2nml Tcp98K BlocPavr to alogl3V1 A if u of LGm The ogllEII39 randlfr 90 the pmm t TK NC STATE UNIVERSITY log mum GasLiquid Transition in Pores Capillary Condensation bx er ental stud es rpmoleculardiameter Plots of T07 TCPTC and 5 i 3 TC 7 T011 C vs drp should be linear according to simple theories of adsorption This is observed experimentally with deviations for small pores see figure 10 14 Shifts in the msmm and pure 1 a Cmpci uli uluccd invcisc pol39c l ttzlllls rm m l inn mm o w x it x mt39 n rcxulrx m rd rm ncmorkctl mmpmuus mlxoibcntx V360 40 and 11 m the eye Hum 240 l t us mm mm Open mu all at pomm s It and 7 our open 1mm at smaller alucs Drul w inch lie well bclon m lino m the lcmpcrmm39cofMCh Jl r Silica gel 2M Lincx a NC STATE UNIVERSITY GasLiquid Transition in Pores Capillary Condensation LJ uid adsorbed on a random porous solid composed ofa matrix of nonoverlapping spheres see inset Results were obtained by thermodynamic integration of GCMC simulation results by Monson and coworkers rigmc 33 1cmmxmmlimbm phase dldgl nm rm 3 LJ uid m 4 11d made up ofa mnu n or wheres Tlic nnixtcncc rune rm the lull uid n slum m Ihc mini un39c Fillnl circles urc he mum m lhc random and urnunrm upping wiltm my u mime Kim nrrcspulhlillg m u 5mm ciogcl um i uid modelled on mctlizinc own are iltl r am l as IOI tin nixlciml Ict mulrn lore p 39 u ICI39C I an IT two in ll quot lllt Ll parameters 1 um thc permit u illi ll being the fraction Ilrlllk nlumc UCClIl lcd I ll mlnl plume Frnm 313 NC STATE UNIVERSITY Layering Transitions Complex adsorbed phase diagram can be determined from experiments of adsorption on surfaces can nd as many as seven layering transitions at low temperatures eg for Ar on graphite 0 VquotUzi1mmWWAL mm 34 Pliztsctimgrnnt for m mum n5 y c 0 woxdcr phase mmmon m7 uwwwxwmkaakc i mhcmt mm wluic the min 35 tum Licnotc poxqbic antinuou l39JnMIiDn 30 39 iinwikinksinthcisothcrim Thu 1 AWWquotWIWWW cm F rcp sun sand gwnd nu E 356 mpmivc The uid or imdcr d i a mm mm cxm in a legion Mm them i 254 L 50 quot dc n J L tics m it rosewood r z 7 K 16 39 has i 05 k 391 m 072 i 02 ii and the 236 mm triple point 41 m2 7 ms i 05 K um m3 672 i 05 Ii Reproduced 90 i I I 39mni li NC STATE UNIVERSITY Layering Transitions For intermediatestrength surfaces surfacefluid interaction is comparable to fluidfluid Picture is complicated by the additional presence of a non zero wettinq temperature Below the wetting temperature adsorbed layer is of finite thickness even at bulk liquidvapor coexistence Above the wetting temperature the thickness of the adsorbed layer diverges Wetting temperature is lower than roughening temperature and the layerwise transition lines intersect the bulk liquid vapor curve Get adsorption isotherms at intermediate temperatures that show several steps at low pressures followed by a smooth divergence at high pressures See Pandit et al Phys Rev B 26 5112 1982 for more details 7 FreezingMelting Transitions In addition the contact layers those adjacent to the two pore walls behaved differently than the rest of the fluid Contact layers appeared to freeze at a different temperature than the rest of the fluid in the pore Attractive walls higher freezing temperatures that the rest of fluid in pore Repulsive walls lower freezing temperatures than the rest of fluid in pore Maddox and Gubbins J Chem Phys 107 9659 1997 studied freezing of LJ methane inside cylindrical pores using GCMC and MD simulations Similar conclusions about the nature of the freezingpoint shifts as in slit pores Found important differences because of the increased confinement in a cylindrical geometry smaller values in the pore freezing temperatures when compared with those obtained within a slit geometry of the same pore size 7 NC STATE UNIVERSITY 7 FreezingMelting Transitions Freezing and melting in the previous simulation studies displayed larqe hysteresis loops 1030 K Cannot precisely locate the true equilibrium freezingmelting point Dominguez et a Mol Phys 96 209 1998 calculated free energies of the solid and fluid phases in a slit pore used thermodynamic integration Method needs a reversible path of integration ie path must not cross any firstorder transition gt method fails for strongly attractive pores Gibbs free energy was calculated for purely repulsive and weakly attractively pores Thermodynamic transition temperature was located Freezing transition was first order Authors calculated the condensation line and mapped out the complete phase diagram for weakly attractive pores figure 50 7 NC STATE UNIVERSITY 7 FreezingMelting Transitions 002m I W I 00151 r ll I l l l I 393 0010 I 9 Liquid I I I I I I x o Solid 1 I I X r I x 00051 y c x I 5 39 I If gt IL r 9 IA XX 3 0000 i M 4 06 as T bulk LJ uid III III pom mm 11 10w re re whereas for bulk k is the bqu p IIIc Fith n 50 Coexistence cuncs For Con ned and NoIthIII oI he pom P Ammo Ihc normal pres Ion I cplc m D hrlulchallLI Thmym 5qu o Uquhlrgns chIIIquII oIIquIigaI AIIquIamIIII ouu chIIlschwaH inquidr 1iquidr 0H N 300 wma IIIIIIICIIIMIII x HqIIiI gm liquidaolid Nsm n ofmuicculcs NI in he sImIIIqun From 09 39 Is hC number NC STATE UNIVERSITY CHE596M MultiScale Modeling of Matter Dr Erik E Santiso Lecture 27 Continuum modeling Conservation principles and phenomenological equations 7 NC STATE UNIVERSITY OUTLINE Why continuum modeling The continuum hypothesis What s in a continuum model Generalized transport equation Numerical methods Discretizations The Finite Differences Method The Finite Control Volumes Method 2 NC STATE UNIVERSITY 7 Why continuum modeling We can often obtain information about large systems without looking at the atomic level of detail just by applying macroscopic conservation principles For many macroscopic systems this is all that is necessary or possible Most problems in civil mechanical materials and chemical engineering are typically solved using continuum models 7 NC STATE UNIVERSITY 7 The Continuum Hypothesis The fundamental assumption underlying continuum models is the continuum hypothesis which may be stated as quotMatter can be subdivided inde nitely into smaller portions An equivalent statement is to say that matter is a continuum We know since the times of Dalton that this is not correct However it is a very good approximation at length scales that are much greater than the mean free path of the molecules in the system Not to be confused with Cantor s continuum hypothesis 7 NC STATE UNIVERSITY The Continuum Hypothesis A A A consequence of the continuum hypothesis is that all physical properties of a system can be treated as elds ie they can be assigned a value at any point within the system How can we define these fields if real systems are not continuous Answer Use a limiting process W l7 Tempa39zmre prn le nn 2 lasarhmted surface nhmined with the nite demerits methnd Frmn sM Ra39zdh ksha p Michaleris Int J Numer Meth Eng 47 mm min NC STATE UNIVERSITY The Continuum Hypothesis A property at a point in a continuum model is defined by a limiting process We first define 5V0 as the smallest volume compared to the size of the system that is still large compared to the molecular dimensions The point value of a property in a continuum model is taken as the average value within this volume Region R Volume W1 y 53 Example Definition of density lim 5 5Vgt5V0 y Figure from the Aerodynamics for students web book by DJ Auld and K Srinivas the University of Sydney URL httpwwwaeromechusydeduauaero 7 NC STATE UNIVERSITY 7 The Continuum Hypothesis There are two main consequences of working within the continuum hypothesis Any properties or relations that result from the nature of the molecules that comprise the system cannot be derived from a continuum theory Such relations are called phenomenological equations and the corresponding properties are called phenomenological coef cients On the other hand any properties or relations that result from conservation principles or the second law of thermodynamics can be obtained from a continuum theory Thus continuum models often involve applying one or more conservation principles and the use of external information in the form of phenomenological equations and coefficients is NC STATE UNIVERSITY Continuum Modeling Examples of phenomenological equations and coefficients Relations between fluxes and driving forces eg Newton s law i 6 7 Fourier s law cf 2 k T Fick s law 7 Dvc Equations of state eg the ideal gas law van der Waals equation the PengRobinson equation etc These equations complement the conservation principles so that the resulting equations can be solved The phenomenological coefficients eg u k D must be obtained from experiments or molecular models 7 7 NC STATE UNIVERSITY A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet 4nr2Nr r Ar 475r2rA Consumed by reaction NV Flux in the radial direction 116 moles per unit area and time that cross a spherical control surface rA Volumetric reaction rate 1 e moles per unit volume and time consumed by reaction NC STATE UNIVERSITY 7 71 A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet Conservation of mass requires that 2 Moles accumulated 47W NV Ar within the control 2 volume Moles Moles Moles consumed commg 1n g01ng out by react10n 4TH 2rA Consumed by If the system 1s at steady state the reaction accumulation term is zero Thus we get 0 4nr2Nrr 4nr2Nr rAr 4nr2ArrA 7 NC STATE UNIVERSITY 7 A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet From conservation of mass 4nr2NrrAr O 4nr2Nrr 4nr2NrrAr 4nr2ArrA We can divide by 475m and rearrange to 2 errr r2NrrAr 475r2rA rAr Ar Consumed by Taking the limit as Ar gt O and using reaction the de nition of a derivative we get a differential equation r2rA dir2Nr r 7 NC STATE UNIVERSITY 7 A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet 4nr2Nr r Ar 475r2rA Consumed by reaction From conservation of mass we got the differential equation r2rA r2Nr This is as far as we can go by using the principle of conservation of mass alone However this is not enough We have one differential equation but two unknown functions rA and NV We need a relation between these functions to solve the problem Here s where phenomenological equations come 1n NC STATE UNIVERSITY 7 7 A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet From conservation of mass we got the differential equation 4751 2Nr d rAr 2 2 r r r N A dr r If we know the rate law we can relate the reaction rate to the local concentration Let s assume that this is a firstorder reaction 475r2rA rA akC Consumed by reaction If we know the diffusion regime we can relate the ux to the concentration gradient Let s assume that it is Fickian diffusion NV 2 D dr 7 NC STATE UNIVERSITY 7 A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet 4nr2Nr r Ar 475r2rA Consumed by reaction Thus the conservation of mass equation d r2r rZN A dr r Becomes assuming constant D rzakC Dip2 dr dr The phenomenological coef cients D a and k the effective diffusion coefficient speci c surface area of the catalyst and the rate constant must be supplied to this model They cannot be obtained from a continuum theory only om experiments or a molecularlevel theory NC STATE UNIVERSITY 7 7 A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet 4nr2Nr r Ar 475r2rA Consumed by reaction Finally we need two boundary conditions If we assume that the external liquid or gas is well mixed we have a xed concentration on the surface of the catalyst And because of the spherical symmetry of the system there cannot be a ux at the center dC Z 0 dr r0 NC STATE UNIVERSITY 7 7 A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet 4nr2Nr r Ar 475r2rA Consumed by reaction Thus our continuum model for this system is rzakC Dip2 dc dr E C RCs 0 dr r0 In order to obtain a numerical solution to this model we need values of the phenomenological coef cients D and ak These can be obtained from experiments or molecularlevel theories simulations NC STATE UNIVERSITY 7 7 4nr2Nr r Ar 475r2rA Consumed by reaction A Continuum Model Example Steady state diffusion and reaction in a spherical catalyst pellet In this particular case the model can be solved exactly The analytical solution is Cr CS Esmh hrR r smhd Where I E RwlakD is the Thiele modulus a parameter that compares the reaction rate to the diffusion rate From this we can obtain useful information such as the effective reaction rate CS I 1 tanh I NC STATE UNIVERSITY 7 47ch dC r 3D e 47cR33 dr 2 rR R 7 Continuum Modeling The previous example shows the main features common to all continuum models Apply one or more conservation principles to a control volume Take the limit of the conservation equations as the volume goes to zero to get a differential equation Use of one or more phenomenological equations to reduce the number of unknowns Choose appropriate boundary conditions Solve the differential equations subject to the boundary conditions It is possible to formulate conservation principles in integral form as well the two formulations are related by Gauss theorem fly de dl V 6V We will deal mostly with the differential formulation in this lecture 7 NC STATE UNIVERSITY Continuum Modeling So why do we need computer simulations in continuum modeling 7 NCSTATEUNIVERSITY 7 Continuum Modeling So why do we need computer simulations in continuum modeling Unfortunately most continuum models cannot be solved exactly Thus in most cases continuum modeling requires finding approximate solutions to systems of partial differential equations We will go over the most common methods to solve this kind of problem later 7 NC STATE UNIVERSITY 7 Continuum Modeling An example of a more realistic therefore complicated continuum model an evaporating water droplet containing latex particles in a fluorinated oil air interface ST Chang and OD Velev Langmuir 22 1459 2006 V V 0 Conservation of mass 195K 6V VP V VV VVZV Conservation of momentum 61 p 6T 5 V VT ocva Conservation of energy V 0 Airwater interface kwatzrVT 3911 hT Too hl39JE 29 1011 T WSWa t Airoil interface k VTn kaT n L711 V0 TToo 17ml Twatzr V6 t allm Wateroil interface k VTn k VTn watzr L711 NC STATE UNIVERSITY Container edges Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as M a Accumulation Convection Diffusion VpovVpDVqgt5 Source generation In this equation d A specific per unit mass or mole property Examples velocity momentum per unit mass specific energy mole fraction 9 Density or molar density v Veocity D Diffusivity tensor 6 V del operator Vi jk ax ay 52 NC STATE UNIVERSITY 7 7 Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as Vp vV39pDVCDS Some examples NC STATE UNIVERSITY 7 Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as Vp vV39pDVCDS Some examples q 1 gt Z V pv 0 Continuity equation I Total mass per unit mass 1 NC STATE UNIVERSITY 7 Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as Vp vV39pDVCDS Some examples a 21 gt V pv 0 Continuity equation Total mass per unit mass This equation is often combined with the generalized transport equation to give ppvV VpDVCDS NC STATE UNIVERSITY 7 71 Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as Vp vV39pDVCDS Some examples q 1 gt 2 SV pv 0 Continuity equation q 2 V V pVV V TV VP pf Momentum equation Velocity momentum per unit mass NC STATE UNIVERSITY 7 1 Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as Vp vV39pDVCDS Some examples 69 21 gt a V39pvO Continuity equation I a v qD v 39 VDVVV39TV VPpf Momentum equation Velocity momentum per For a Newtonian uid of constant density and viscosity unit mass this can be rewritten as the NavierStokes equation ppvVvuV2v VPpf NC STATE UNIVERSITY 7 7 Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as Vp vV39pDVCDS Some examples 69 q 2 EV pv 0 Continuity equation d v gt V pvv V TV VP pf Momentum equation q e gt V Dev V pa Ve Energy equation Specific energy hTw T 5T T4 NC STATE UNIVERSITY 7 71 Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as Vp vV39pDVCDS Some examples 69 q 1 gt a V pv 0 Continuity equation I 82 TV 39 9W V TV VP pf Momentum equation q 2 e V pev V pa Ve Energy equation hTw T 5T T4 qDZXA gt VVCAVpDVxArA Mass equation Mole fraction of component A 7 NC STATE UNIVERSITY 7 Generalized transport equation Most continuum methods are designed to solve some form of generalized transport equations Such an equation can be written as Vp vV39pDVCDS For simplicity we will restrict our discussion to the case where the medium is isotropic ie properties do not depend on direction and has constant physical properties Under such conditions the generalized transport equation becomes VV DV2 S 2 2 2 Where V2 EVV is the Laplacian operator V2 a a a and s Sp 6x2 aJZ a 2 Z 7 NC STATE UNIVERSITY 7 Numerical methods So now we are faced with the question how do we find an approximate solution to one or more generalized transport equations VV DV2 S The most widely used methods are The method of finite differences FDM The method of finite control volumes FVM The method of weighted residuals and in particular The finite elements method FEM These methods deal with the spatial components of the transport equation The time dependence is usually treated separately We will discuss how to deal with timedependent problems at the end NC STATE UNIVERSITY 7 Numerical methods For now let us consider a steadystate time independent problem vV DV2 S There are two fundamental ways to approach a problem in which the unknown is a function Discretize the function s domain and deal only with the values of the function at certain points Expand the unknown function on a series of independent functions and find the expansion coefficients The former is actually a particular case of the latter but we will consider it separately NC STATE UNIVERSITY 7 7 Discretizations By discretizing the domain we replace the unknown function by the set of values of the function at a finite number of points This turns the differential problem into an algebraic problem that can be solved on a computer For simple geometries and relatively small gradients a simple grid is usually good enough For complex geometries or problems where large gradients are concentrated in small regions it is necessary to use more complex andor adaptive meshes 7 NC STATE UNIVERSITY 7 Discretizations Simple discretizations Canes ian Polar Relatively easy to implement Not efficient if gradients vary widely within the domain NC STATE UNIVERSITY Discretizations More elaborate discretizations v Left femur smoothed decimated and smoothed again Number of Nodes 2482 Surface Triangles 4935 More realistic no spherical cow More dif cult to implement Figure from The Visible Human Project s web page httpWWW nlm nih quot quotnnAUTHORSSULLIVANFIGURE9HTME Discretizations Once we have a discretized space how do we represent The system s properties The differential operators that appear in the transport equations Two possibilities Use the values of the function at the grid points and define derivatives at the same points Finite differences Associate a single function value to each grid cell and define fluxes at the cell boundaries Finite control volumes Derivatives are usually obtained using finite difference formulas 7 NC STATE UNIVERSITY 7 Finite Difference Formulas Finite difference formulas can be derived from Taylor expansions for example fxh fxhoh2 3 fxh2fx0h More precise formulas can be obtained by using more expansions and more terms in each expansion d h2 d2 h3 d3 4 fxhfxhd J7 dx dx0h 561215de 4 fx h fx hdx 2 dxz 6 dx3 Oh xhfxh0h2 dx 2h d2ffxh2fxfx h dx2 h2 gt 0012 7 NC STATE UNIVERSITY 7 The Finite Differences Method In the Finite Differences Method FDM both the properties and the fluxes are defined at the grid points The derivatives are replaced by finite difference formulas This turns the differential transport equation and the boundary conditions into a set of algebraic equations i1j 0 These equations can be solved to obtain the properties at the grid points More precise solutions can be obtained by using finer grids NC STATE UNIVERSITY 7 The Finite Differences Method For example if we work on a 2D uniform Cartesian grid we have vV DV2 S Ax q gtqij s gtsl sCDZj I 1144 V Ay 8x 8y 11414 11214 1313141 11314 c1gt c1gt cp a M 2m 1 My 121 z39j i1j 2 2 qu3 a 33 5 gt 6x 8y Did1 DHM 21 1317111 1311141 21 131111 NY NY 1 NC STATE UNIVERSITY 7 The Finite Differences Method For example if we work on a 2D uniform Cartesian grid we have I i1j vV DV2 s gt gt i1j quelj39 DEW 2c1gtw c1gt AxY Set of algebraic equations for the Did 3 Lqij1 2121 I J 2Ay iilj CD 2 QD 7 11 1 111 S 1 NC STATE UNIVERSITY 7 M The Finite Differences Method As an example let us consider the heat transfer problem vV DV2CDS gt 6T 6T aZT aZT Conduction convection gt u v 2 2 qs ax ay ax 5y surface heat generation In 2D T O x OVy T20 xZI Vy We ll assume constant velocities heat 6T a 0 y 0Vx diffuswlty and generation T O y 1Vx NC STATE UNIVERSITY 7 The Finite Differences Method Discretizing the differential equation inner nodes Ax u va T 0L aZT aZT gt ax ay 6x2 832 qs Til1 Ay gt u ZHJ 771j v 72141 7j71 2Ax 2Ay a71j 77139 72141 713971 qs 7271 21 241 Ax Ay Tm 1 2Nx 1 j 2Ny 1 1 NC STATE UNIVERSITY 7 The Finite Differences Method Discretizing the differential equation inner nodes Ax u va T or aZT aZT gt ax ay 6x2 832 qs Til1 Ay gt u ZHJ 771j v 72141 7j71 2Ax 2Ay a71j 77139 72141 713971 qs 7271 21 241 Ax Ay Tm i2Nx 1 jzzjjNy1 N x 2Ny 2 equations 1 NC STATE UNIVERSITY 7 The Finite Differences Method Discretizing the boundary conditions Ax T20 xOVy i1j1uNy T20 XZLVy gt ZJZO iZNxj17HaNy Zj1 T20 y1Vx JZNyi27 39aNx1 Ay 2717 Zd Ti1j 2Ny Nx 2 equations 7 NCSTATEUNIVERSITY 7 The Finite Differences Method Discretizing the boundary conditions Ax T20 xOVy i1j1uNy T20 XZLVy gt ZJZO iZNxj17HaNy Zj1 T20 yZLVx jZNyl27HaNx 1 a TO yOVxgt 5y 741 721 gt 20 for 1l2Nx 1 2Ay NC STATE UNIVERSITY 7 The Finite Differences Method Discretizing the boundary conditions Ax T20 xOVy i1j1uNy T20 XZLVy gt ZJZO iZNxj17HaNy Zj1 T20 yZLVx jZNyl27HaNx 1 a TO yOVxgt 5y 7 41 7 21 gt 1 0 for 1l2Nx 1 2Ay I Til1 7 Noj0 NC STATE UNIVERSITY 7 The Finite Differences Method Discretizing the boundary conditions Ax T20 xOVy i1j1uNy T20 XZLVy gt ZJZO iZNxj17HaNy Zj1 T20 yZLVx jZNyl27HaNx 1 a TO yOVxgt 5y 7 41 7 21 gt 1 0 for 1l2Nx 1 2Ay I Til1 7 Noj0 gt Use mirror method or a fonvard difference 1 NC STATE UNIVERSITY 7 The Finite Differences Method T T20 x1Vy 3 Tw0 T Discretizing the boundary conditions 0 xOVy O y1Vx a TO yOVx gt 5y Ax i1j1Ny iNxj1Ny Tim ijyi2Nx 1 Ay 2111 21 Ti1j T 7 Z141 7 j 1 1 qs NC STATE UNIVERSITY 7 1 The Finite Differences Method T T20 x1Vy 3 Tw0 T Discretizing the boundary conditions 0 xOVy O y1Vx a TO yOVx gt 5y Ax i1j1Ny iNxj1Ny Tim ijyi2Nx 1 Ay 2111 21 Ti1j T 7 71 7j 1 11 qs gt Nx 2 equations NC STATE UNIVERSITY 7 1 The Finite Differences Method In total we have Nx 2Ny 2 equations inner nodes 2Ny Nx 2 equations left right and top nodes Nx 2 equations bottom nodes 72 NC STATE UNIVERSITY The Finite Differences Method In total we have Nx 2Ny 2 equations inner nodes 2Ny Nx 2 equations left right and top nodes Nx 2 equations bottom nodes NxNy ZNx ZNy42NyNx 2Nx 2 NxNy equations for NxNy unknown temperatures i NC STATE UNIVERSITY 7 The Finite Differences Method For example here s the solution for ca 1Nx Ny 31u v 1 qs 10 T m m If T 5 NC STATE UNIVERSITY The Finite Differences Method A trickier problem on 1Nx Ny 31u 50v1 61S 10 mmmmmmmmmmmm i i 0 T iTEEEE NCSTATEUNNERSWY The Finite Differences Method If we reduce the number of nodes 0L 1Nx Ny 6u 50v1qs 10 T 0 02 04 06 08 7 773 The Finite Differences Method If we reduce the number of nodes or 1Nx Ny 6u 50v1qs 10 Unphysical This is the main defect of FDM for coarse grids it does not guarantee conservation in this case of energy T 5 NC STATE UNIVERSITY The Finite Differences Method Advantages of FDM Simple Easy to implement for simple geometries Disadvantages Does not guarantee conservation for coarse grids Difficult to implement for complex geometries These problems are solved to some extent by FCV and FEM NC STATE UNIVERSITY 7 3 NC STATE UNIVERSITY The Finite Control Volumes Method A big problem with FDM is that it may give unphysical solutions for coarse grids How can we guarantee that the solution does not violate any conservation principles NC STATE UNIVERSITY 7 The Finite Control Volumes Method A big problem with FDM is that it may give unphysical solutions for coarse grids How can we guarantee that the solution does not violate any conservation principles One solution apply the conservation principles to the grid elements NC STATE UNIVERSITY 7 1 The Finite Control Volumes Method For example on a 2D uniform Cartesian grid we have vV DV2 S Ax lt gt CIN Ay I gt I gt I T Tv T CD CD C2Wgt P E Ll T W Tv e T CDS l gt I gt I NCSTATEUNIVERSITY 7 The Finite Control Volumes Method For example on a 2D uniform Cartesian grid we have vVCDDV2CDS 3 VvoDV2qgts Ax lt gt c1 From continuity equation Vv 0 Ay l gt IN gt I T TV T lt2 lt22 2 r W Tv 6 r CDS l gt I gt I NC STATE UNIVERSITY 7 The Finite Control Volumes Method For example on a 2D uniform Cartesian grid we have VV DV HwgtVWMDV Hw A CDN From continuity equation Vv 0 Ay l gt I gt I l TV l I I n I Integrate over grid cell to find lt2 lt22 2 T w TV e T ueqe uwCDwAyanDn VSCDSAx I 23 gt I D 62 Ay 62W 162 Ax sAxAy 6x 9 6x w Kay tam 7 NC STATE UNIVERSITY 7 The Finite Control Volumes Method For example on a 2D uniform Cartesian grid we have VV DV HsgtVWMDV Hw A CDN From continuity equation Vv 0 Ay l gt I gt I l TV I I I n I Integrate over grid cell to find lt2 lt22 2 T w TV e T ueqe uwCDwAyanDn VSCDSAx I 23 gt I D 62 Ay 62W 162 Ax sAxAy 6x 9 6x w Kay tam Define n E VQD DVQD then 111 tn1 SA Same as conservation W m law applied to grid cell 7 NC STATE UNIVERSITY 7 The Finite Control Volumes Method For 3D case the idea is similar vVCDDV2CDS 3 VvoDV2qgts From continuity equation V v O 2 NC STATE UNIVERSITY The Finite Control Volumes Method For 3D case the idea is similar vVCDDV2CDS 3 VvoDV2qgts From continuity equation V v 0 Integrate over volume of grid cell to find From Gauss theorem mVvqgtdVDJIIVVqgtdimst 3 IJIV39Fde F39quotSdl vqgt nSdA g wop nSdA Vmst W W V Define n E vlt1gt DWI then dA 07 Same as conservation n us His 6V law applied to grid cell 7 NC STATE UNIVERSITY 7 The Finite Control Volumes Method Back to the 2D uniform Cartesian grid Ax lt gt CIN Ay I gt I gt I T Tv T CD CD C2Wgt P E Ll T W Tv e T CDS l gt I gt I i7 NC STATE UNIVERSITY 7 The Finite Control Volumes Method Back to the 2D uniform Cartesian grid We had a uEQE N D 62 62 lining 163le MW axe 8x w Kay tam Let us first consider a purely diffusive problem ie u v 0 Everywhere Then Diirigl ti lwlAwlt35lngllmr sAxAy O Ay R I l n m I gt 9 gt I Q l gt I ve gt 9 gt 9 gt I 7 NC STATE UNIVERSITY 7 The Finite Control Volumes Method Back to the 2D uniform Cartesian grid Let us first consider a purely diffusive problem Ax ie u v 0 Everywhere Then c1 Ay gt IN I D Ay 62 62 Ax 5x 9 5x w 6y 8y l ivquot r quot S D D D sAxAy O Tl 7 We can use a finitedifference approximation for l Vs I the derivatives at the edges to get CDS l gt I gt I DttJ tlAy Dir u jlmmyzo KAyJKAy 7 NC STATE UNIVERSITY 7 The Finite Control Volumes Method Back to the 2D uniform Cartesian grid Can simplify to Ax lt gt CIDE ZCDPCDW CDN ZCDPCDS D 2 2 s0 Ay 9 s 9 Ax M T iv T lt2W lt22 2 T W T CDS l gt I gt I NC STATE UNIVERSITY 7 The Finite Control Volumes Method Back to the 2D uniform Cartesian grid Can simplify to Ax lt gt DCIgtE ZCIgtPCDWCDN 2CDPCIgtSHSO 2 2 Ay 9 s 9 Ax M l lv l Eu CE 2 l W Tv 6 l CD Set of algebraic equations for the s l gt IS gt I NOTE In the absence of convection and if we use the central difference approximation we get the same equations as FDM 7 NC STATE UNIVERSITY 7 The Finite Control Volumes Method Back to the 2D uniform Cartesian grid Can simplify to Ax lt gt D E 2CIPCDW N 2 P 5HSO 2 2 Ay 9 s 9 Ax M l lv l Bus C22 2 l W Tv 6 l CD Set of algebraic equations for the s l gt IS gt I NOTE In the absence of convection and if we use the central difference approximation we get the same equations as FDM What if we do have convection next class 7 NC STATE UNIVERSITY 7
Are you sure you want to buy this material for
You're already Subscribed!
Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'