#*---------------------------------------------------------------------
* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
*
* This LIBRARY is free software; you can distribute it and/or modify
* it under the therms of the ALSOC FREE LICENSE as available at
* http://www.enq.ufrgs.br/alsoc.
*
* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
* from http://www.rps.eng.br Copyright (C) 2002-2004.
* All rights reserved.
*
* EMSO is distributed under the therms of the ALSOC LICENSE as
* available at http://www.enq.ufrgs.br/alsoc.
*
*----------------------------------------------------------------------
* 7. Diffusion with chemical reaction in a one dimensional slab
*----------------------------------------------------------------------
*
*   Description:
*      This problem is part of a collection of 10 representative
*	problems in Chemical Engineering for solution by numerical methods
*	developed for Cutlip (1998).
*
*   Subject:
*       * Transport Phenomena
*		* Reaction Engineering
*
*	Concepts utilized:
*		Methods for solving second order ODEs with 2 point boundary
*	values typically used in transport phenomena and reaction kinetics.
*
*	Numerical method:
*		* Simultaneous ODEs with split boundary conditions
*		* Resolved by finite difference method
*
*   Reference:
*		* CUTLIP et al. A collection of 10 numerical problems in 
*	chemical engineering solved by various mathematical software 
*	packages. Comp. Appl. in Eng. Education. v. 6, 169-180, 1998. 
*       * More informations and a detailed description of all problems
*	is available online in http://www.polymath-software.com/ASEE
* 
*----------------------------------------------------------------------
* Author: Rodolfo Rodrigues
* GIMSCOP/UFRGS - Group of Integration, Modeling, Simulation, 
*					Control, and Optimization of Processes
* $Id$
*--------------------------------------------------------------------*#
using "types";



Model problem
	PARAMETERS
outer N as Integer	(Brief="Number of discrete points", Lower=3);
	
	Co		as conc_mol		(Brief="Constant concentration at the surface");
	D		as diffusivity	(Brief="Binary diffusion coefficient");
	k		as Real			(Brief="Homogeneous reaction rate constant", Unit='1/s');
	L		as length		(Brief="Bottom surface");
	
	
	VARIABLES
	C(N+2)	as conc_mol		(Brief="Concentration of reactant");
	z(N+2) 	as length		(Brief="Distance", Default=1e-3);
	dz		as length_delta	(Brief="Distance increment");	


	EQUATIONS
	"Discrete interval"
	dz = (z(N+2) - z(1))/(N+1);
	
	for i in [2:(N+1)] do
	"Concentration of reactant"	
		(C(i+1) - 2*C(i)+ C(i-1))/(z(i) - z(i-1))^2 
			= (k/D)*C(i);
		
	"Discrete length"
		z(i) = z(i-1) + dz;
	end

	# Boundary conditions
	"Initial and boundary condition" # z = 0
	C(1) = Co;
	
	"Upper boundary" # z = L
	(C(N+2) - C(N+1))/(z(N+2) - z(N+1)) = 0*'kmol/m^4';
	
	
	SET
	Co= 0.2*'kmol/m^3';
	D = 1.2e-9*'m^2/s';
	k = 1e-3/'s';
	L = 1e-3*'m';
end



FlowSheet numerical_solution
	PARAMETERS
	N		as Integer;
	
	
	DEVICES
	reac	as problem;
	
	
	SET
	N = 10; # Number of discrete points
	
	
	SPECIFY
	reac.z(1) = 0*'m';

	reac.z(N+2) = reac.L;

	OPTIONS
	Dynamic = false;
end



FlowSheet comparative
	PARAMETERS
	N		as Integer;
	
	
	VARIABLES
	C_(N+2)	as conc_mol	(Brief="Concentration of reactant by analytical solution");
	
	r_		as Real		(Brief="Pearson product-moment correlation coefficient");
	
	Cm 		as conc_mol	(Brief="Arithmetic mean of calculated C");
	C_m		as conc_mol	(Brief="Arithmetic mean of analytical C");
	
	
	DEVICES
	reac	as problem;
	
	
	SET
	N = 10; # Number of discrete points
	
	
	EQUATIONS
	"Analytical solution"
	C_ = reac.Co*cosh(reac.L*sqrt(reac.k/reac.D)*(1 - reac.z/reac.L))
	/cosh(reac.L*sqrt(reac.k/reac.D));

	"Pearson correlation coefficient" # used by softwares like MS Excel
	r_ = (sum((reac.C - Cm)*(C_ - C_m)))/sqrt(sum((reac.C - Cm)^2)*sum((C_ - C_m)^2));
	
	"Arithmetic mean of C"
	Cm = sum(reac.C)/(N+1);
	
	"Arithmetic mean of C_"
	C_m = sum(C_)/(N+1);
	
	
	SPECIFY
	reac.z(1) = 0*'m';

	reac.z(N+2) = reac.L;

	OPTIONS
	Dynamic = false;
end


FlowSheet analytical_solution
	PARAMETERS
	Co		as conc_mol;
	L		as length;
	k		as Real(Unit='1/s');
	D		as diffusivity;
	
	
	VARIABLES
	C		as conc_mol		(Default=0.2);
	z		as length		(Default=1e-3);
	
	
	EQUATIONS
	"Change time in z"
	z = time*'m/s';
	
	"Analytical solution"
	C = Co*cosh(L*sqrt(k/D)*(1 - z/L))/cosh(L*sqrt(k/D));
	
	
	SET
	Co= 0.2*'kmol/m^3';
	D = 1.2e-9*'m^2/s';
	k = 1e-3/'s';
	L = 1e-3*'m';
	
	
	OPTIONS
	TimeStart = 0;
	TimeStep = 1e-6;
	TimeEnd = 1e-3;
end