Hey everyone, long time no post :). It’s also been some time since I have been working in Matlab, but today I found enough time to waste. I am working on a UNIQUAC model for predicting binary solution activities, and also have been trying to generate a X-Y plot for distillation. I ran into a problem however, I have a function:

%Uniquac.m function gamma = Uniquac(r_matrix, q_matrix, a_matrix, T, x1) ... (More code here)

The extra parameters in the function above are really nice, I might have to call several different binary components with the same function. While I could hard-code the constants r_matrix, q_matrix, and a_matrix… I really feel cheap doing so. Why does it matter so much? I don’t want to corrode my perfectly general code with very specific constants. Consider the very real possibility that you might want to use this function in fsolve, a numeric solver in Matlab. Fsolve requires:

1. That the first input be a function.

2. The second input must be a ‘guess’ for a possible solution.

The problem is that when you use a function like Uniquac.m, fsolve(@Uniquac, x0) says to Matlab: “Use the matrix x0 as the starting point for meeting my constraints.” But what is x0?

Well, it’s the inputs of Uniquac.m… which are r_matrix, q_matrix, a_matrix, T and x1. But hold on… maybe I only want to allow T to change and if r_matrix, q_matrix, a_matrix, and x1 are supposed to be constants, then it would be more than desirable if Matlab didn’t change them. But in this configuration Matlab cannot distinguish between the ‘constants’ and the variables that you want to change.

So what then? How to fix this problem? Well, you *could* go back and ruin all that beautiful code you just wrote by hard-coding constants and making global variables… but you could also go back to dragging your fists and beating things with a club.

OR we could do something VERY hip and gets rid of the bad vibes. For this we pull out a seemingly useless friend: the anonymous function. We covered these previously but here’s the angle:

%Some other script than Uniquac.m shortUniquac = @(T) Uniquac(r_matrix, q_matrix, a_matrix, T, x1)

Yea! Now you’ve got a good as new function shortUniquac(T) that depends only on T… without a bunch of band-aids. Now for fsolve:

%fsolver answer = fsolve(shortUniquac, T0)

If you would like the files for the Uniquac stuff, leave a comment with your email (do not put “example@gmail.com”, rather “example \ gmail” to avoid spam) and what you intend to use it for and I’ll send you a link.

——-Updated 04/23/2012———

User List

- Changsoo Lee — Yonsei Univ. Korea
- Bobby J — The University of Edinburgh
- Santiago — University of Los Andes; Bogotá, Colombia

October 24th, 2011 at 6:07 am

I would like mfile for unifac model

October 24th, 2011 at 12:58 pm

Aidin,

Because I release the Unifac model as a non-profit, non-commercial use application I require that you send me an email that is associated with your learning institution (preferably a .edu) and what your intended use for the program is. Provided with these I can release the code to you, please keep in mind that it comes with ABSOLUTELY NO warranty.

Take care,

pavlovsr

November 29th, 2011 at 12:32 pm

I would like m file 🙂

I’m learning chemical engineering in Yonsei Univ. Korea.

Plz help me T T

ch1ch2ch3@naver.com

November 29th, 2011 at 12:51 pm

oh, sorry, ch1ch2ch3/naver.com plz T T

December 1st, 2011 at 12:28 pm

Dear Mr. Lee,

Is it possible to set the a_matrix as variable, so we ask fsolve to adjust the interaction parameters having a set of experimental data?

Best regards,

Elias.

—

emonteirofilho\ufsj.edu.br

December 2nd, 2011 at 4:15 am

@Elias,

You could!

%knowing r_matrix, q_matrix, T, x1

shortUnifac = @(a_matrix) Unifac(r_matrix, q_matrix, a_matrix, T, x1) - gamma_exp_matrix;

`fsolve(shortUnifac, a_matrix);`

Although I wouldn’t use this approach if you don’t have values to compare to… simply because the equations are very non-linear and the value for which you are trying to find is embedded in two exponentials. If you try it let us know!

February 18th, 2012 at 6:03 pm

Hi there.

I’ve been trying for a few weeks now to get UNIFAC to model a system for me through matlab. So far the best I’ve managed is to get numbers that are in the right ballpark (ie. not out by orders of magnitude from the actual values) but I am having real trouble getting past that.

It’s for the purpose of a university design project. Please let me know if you can help asap!

Bobby

s0828614 / sms.ed.ac.uk

February 18th, 2012 at 8:50 pm

I should also add that what I’m trying to do is model a 7-component system, so if your program is capable of handling that, I would greatly appreciate having a chance to try it out.

Bobby

February 18th, 2012 at 10:48 pm

So there could be a couple of issues:

1. If you’re using fsolve it could be that you need to set your tolerance option when calling it (ie you need to set parameters ‘TolX’ or ‘TolFun’ to > 1E-6… the default).

2. Make sure you are not using VLE binary interaction parameters for LLE situations and vis versa! They can show similar behavior but have, say an order of magnitude accuracy.

3. Your system may not be well estimated by UNIFAC. If I remember correctly UNIFAC does not work well for systems with lots of strong bonding forces (ex: Hydrogen bonding) and works best for non-polar systems. Some items you might find interesting are here:

http://www.aiche.org/uploadedfiles/students/departmentuploads/thermodynamics.pdf

http://people.clarkson.edu/~wilcox/Design/BYUprop.doc

http://people.clarkson.edu/~wwilcox/Design/thermodl.htm

Those can help you decide if UNIFAC is right. In any case I will send my code and add you to the user’s list.

Good luck,

Ryan

February 18th, 2012 at 10:52 pm

A seven component system is possible with this code but it’s not built in… codes like HYSYS and ASPEN I believe just do several binary system estimations for things like this (ie seven components = 49-7 interaction coefficients since there are 7 components and because an interaction coefficient b(i,j) does not necessarily equal b(j,i) ).

February 19th, 2012 at 12:48 pm

Do you mean they take into account that many binary interactions? [ b(i,j)s and b(j,i)s ] or that they only calculate activity coefficients for binary systems and then estimate activity within the mixture? If it’s the latter case that would make modelling this system considerably easier, but I am not sure how I would go about estimating activity in the mixture from activity in binary systems.

If this is what you meant, could you go into some more detail about how activity coefficients for a system can be built up from those calculated for binary systems?

February 20th, 2012 at 9:43 pm

So I have never done this myself but something that illustrates how you might go after this is listed in the literature:

http://onlinelibrary.wiley.com/doi/10.1002/aic.690180208/pdf

The section you may be most interested in starts at THREE ION SYSTEM part. This may not be what you want to do since the mixing rule that this paper uses to combine the activity coefficients may not be the best description of the system… this is something that you might do at a PhD or Master’s level to get the details right and it would probably be best to consult your thermodynamics professor about this topic. If you do get somewhere I would be interested so let us know!

February 28th, 2012 at 11:08 pm

Thanks for the info, but with deadlines coming up we’ve abandoned trying to get UNIFAC to work through matlab for this system (and found ways around needing activity coefficients – or use UNISIM to simulate non-idealities).

April 23rd, 2012 at 1:21 am

I would like the m file for the UNIFAC method for matlab. Thanks in advance!. santiagorodriguez23 / live.com.ar

April 23rd, 2012 at 3:10 am

Hey Santiago,

I can give you the code however the code comes with absolutely no warranty and I need to ensure that it is only used for educational purposes. If you can supply me with the name of your learning institution/how you intend to use the code, I can provide you with a copy.

April 23rd, 2012 at 10:53 am

Yeah, I study at the University of Los Andes, in Bogotá, Colombia. Right now I’m doing a class (homework) project where we develop the industrial process of Isoamyl acetate manufacturing and scaling, so we first need to find for the quaternary system (acetate+acetid acid+alcohol+water) all the activity coefficients by UNIFAC based in papers and articles that use it. It’ll be great if you could provide me with a copy. as you see is intended to be used for EDUCATIONAL purposes ONLY, cause is just a class project. Thanks for the fast reponse!!

April 23rd, 2012 at 4:57 am

Yeah, I study at the University of Los Andes, in Bogotá, Colombia. Right now I’m doing a class (homework) project where we develop the industrial process of Isoamyl acetate manufacturing and scaling, so we first need to find for the quaternary system (acetate+acetid acid+alcohol+water) all the activity coefficients by UNIFAC based in papers and articles that use it. It’ll be great if you could provide me with a copy. as you see is intended to be used for EDUCATIONAL purposes ONLY, cause is just a class project. Thanks for the fast reponse!!