33节点DG最优分布_粒子群选址_粒子群_粒子群选址_粒子群算法_33节点

%i;exec('matlab2scilab.m',-1) function [resode,respde,resanalytic]=sensc(methode,F,ploton,dispon) % // SENSC simule un transfert RD+RK+RH lors d'un contact solide fluide avec des isothermes arbitraires (BET, GP...) % // Syntaxe: F = sensc() % // Syntaxe: [RESODE,RESPDE,RESANALYTIC] = SENSC(MODE,F[,PLOTON,DISPON]) % % Values for MODE : % 'init' : constructor of F (default value of MODE). % 'auto' : tout. % 'ode' : Algebraic solution based on Vitrac Hayert 2006. Result in resode. Attention, ne marche pas ! % 'pde' : Finite element solution. Result in respde. % 'analytic' : Analytic solution based on Vitrac Goujot 2008 (in preparation). Result in resanalytic. If first/second letter is uppercase, use only first/second matrix product at corners. If third letter is uppercase of if ploton is true, combine F.t and F.x in resanalytic.tC % voir adsroptionv70.doc % 'odepde' : ode and pde. % 'all' : ode and pde and analytic. % % Values for F : % sol : le solide charge initialement avec la concentration C0sol mol/m^3 initiale. % flu : le fluide charge initialement avec la concentration C0flu mol/m^3 initiale. TODO BUG: Tres mal gere si non nul. % % Values for ploton : % false : no figure (default). % true : various plots for visual check. % 2 : various plots for visual check, do not close or clear figures. % % Values for dispon : % false : no comment (default). % true : various output showing progress. % % Debug is triggered if some letters of 'analytic' are upper case, with following meaning: % 'Analytic' : include alternate matrix product. % 'ANalytic' : exclude stable matrix product. % 'anAlytic' : include profiles for all user timesteps F.t and some additional times around corners (automatic if ploton is set). % 'anaLytic' : use Sagiv product with inverse matrix. % 'anaLYtic' : compare Sagiv with stable matrix product. % 'anaLyTic' : Use Sagiv enhanced with linear inversion. % 'analytIc' : resanalytic.err gives clues about precision near t=t0, but no more around corners in profiles. % 'anaLytiC' : use Sagiv product, with inv instead of \ ; also give plot debug about eigenvalue computatoin; % SENS 1.0 - 22/02/07 - Olivier Vitrac - rev. % SENSC 1.1 - 09/05/07 - Daniel Goujot - rev. % SENSC 1.2 - 11/10/07 - Daniel Goujot - rev. % % revision history % SENSC 1.1 - 09/05/07 - isotherme a coin. % SENSC 1.2 - 11/10/07 - solution analytique. %i;// a utiliser avec .PvsatsurRT %i;// appel sous matlab : cd K:/recherCHE/matlab/senscrep,F=sensc,betiso=feval(F.iso.helperfunctions.baldev02isotherms,'BET',feval(F.iso.helperfunctions.baldev02STARCHandPE,'BET',.4));betiso.awminawmax(3)=-2-1/10;F.iso.S=feval(F.iso.helperfunctions.approximate_isotherm,betiso);F.C0sol=5;F.L=.1;[resode,respde,resanalytic]=sensc('anaLYTic',F,true,true) ; OK, marche 20080611 %i;// 'Pvsat(273.15+27) et 1atm en pascal :';pvsat27=3564.9138;atm=101338.52; %i;// density of LDPE http://www.freepatentsonline.com/EP0152865A1.html 922 kg/m^3. %i;// wikipedia density air: density of air at 1atm at 20�C under 1atm: 1.2041 kg/m^3 %i;// wikipedia density air: \rho~_{_{humid~air}} = \frac{(p-aw*pvsat)}{RsurMair \cdot T} + \frac{aw*pvsat}{RsurMeau \cdot T} en kg/m^3. avec RsurMair=287.05 et RsurMeau=461.495. %i;// /usr/bin/bc : 101338.52/287.05/(273.15+20) donne 1.20427898432070391574 %i;// /usr/bin/bc : 101338.52/287.05/(273.15+27) donne 1.17619318425325454905 %i;// /usr/bin/bc : (101338.52-.2*3564.9)/287.05/(273.15+27)+.2*3564.9/461.495/(273.15+27) donne 1.17306513937994601324 %i;// %i;// http://www.gpa.uq.edu.au/UQweather/airdensity.htm : Air Density (in Kg/m3) = 1.2929 X 273.13 X ( AP - ( SVP x RH )) ( T + 273.13) 760 %i;// %i;// appel sous matlab : cd K:/recherCHE/matlab/senscrep,F=sensc,betiso=feval(F.iso.helperfunctions.baldev02isotherms,'BET',feval(F.iso.helperfunctions.baldev02STARCHandPE,'BET',.2));betiso.awminawmax(1)=0;betiso.awminawmax(3)=-1/10;F.iso.S=feval(F.iso.helperfunctions.approximate_isotherm,betiso);F.C0sol=1;F.L=.01;F.Bi=10;F.t=[1 10]';[resode,respde,resanalytic]=sensc('anaLYTic',F,true,true) %i;// figure(1);clf;hold on;F=sensc;for i=.1:.1:.5;betiso=feval(F.iso.helperfunctions.baldev02isotherms,'BET',feval(F.iso.helperfunctions.baldev02STARCHandPE,'BET',i));betiso.awminawmax(3)=-1/10;F.iso.S=feval(F.iso.helperfunctions.approximate_isotherm,betiso);plot([-.1,F.iso.S.aws],[i,F.iso.S.Cs/F.iso.S.Cs(end),i]);end; %i;// figure(2);clf;hold on;densiteLDPE=922;pvsat27=3564.9138;RsurMeau=461.495;F=sensc;for i=.1:.1:.5;betiso=feval(F.iso.helperfunctions.baldev02isotherms,'BET',feval(F.iso.helperfunctions.baldev02STARCHandPE,'BET',i));betiso.awminawmax(3)=-1/10;F.iso.S=feval(F.iso.helperfunctions.approximate_isotherm,betiso);plot(F.iso.S.aws,[pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K/densiteLDPE,i/100]);end; %i;// clear;figure(10);clf;liste={'bet' 'smith' 'halsey' 'caurie' 'oswin'};F=sensc;for noliste=1:length(liste);hold on;t=liste{noliste};for i=.1:.1:.5;subplot(221);title(t);hold on;betiso=feval(F.iso.helperfunctions.baldev02isotherms,t,feval(F.iso.helperfunctions.baldev02STARCHandPE,t,i));betiso.awminawmax(1:2)=[.01 .9];betiso.awminawmax(3)=-1/30;F.iso.S=feval(F.iso.helperfunctions.approximate_isotherm,betiso);plot([-.1,F.iso.S.aws],[i,F.iso.S.Cs/F.iso.S.Cs(end),i]);subplot(222);hold on;plot([-.1,F.iso.S.aws],[i/10,F.iso.S.Cs/100,i/10]);subplot(223);hold on;densiteLDPE=922;pvsat27=3564.9138;RsurMeau=461.495;plot(F.iso.S.aws,[i/100,pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE,i/100]);subplot(224);hold on;plot(F.iso.S.aws,[i,[pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE]./max([pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE]),i]);end;end; %i;// clear;close all;liste={'bet' 'smith' 'halsey' 'caurie' 'oswin'};F=sensc;for noliste=1:length(liste);figure(noliste);hold on;t=liste{noliste};for i=.1:.1:.5;subplot(221);title(t);hold on;betiso=feval(F.iso.helperfunctions.baldev02isotherms,t,feval(F.iso.helperfunctions.baldev02STARCHandPE,t,i));betiso_awminawmax(1:2)=[.01 .9];betiso.awminawmax(3)=-1/30;F.iso.S=feval(F.iso.helperfunctions.approximate_isotherm,betiso);plot([-.1,F.iso.S.aws],[i,F.iso.S.Cs/F.iso.S.Cs(end),i]);subplot(222);hold on;plot([-.1,F.iso.S.aws],[i/10,F.iso.S.Cs/100,i/10]);subplot(223);hold on;densiteLDPE=922;pvsat27=3564.9138;RsurMeau=461.495;plot(F.iso.S.aws,[i/100,pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE,i/100]);subplot(224);hold on;plot(F.iso.S.aws,[i,[pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE]./max([pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE]),i]);end;end; %i;// clear;close all;liste={'bet' 'smith' 'halsey' 'caurie' 'oswin'};F=sensc;for noliste=1:length(liste);figure(noliste);hold on;t=liste{noliste};for i=.1:.1:.5;subplot(221);title(t);hold on;betiso=feval(F.iso.helperfunctions.baldev02isotherms,t,feval(F.iso.helperfunctions.baldev02STARCHandPE,t,i));betiso_awminawmax(1)=betiso.awminawmax(1)-.05;betiso.awminawmax(3)=-1/30;F.iso.S=feval(F.iso.helperfunctions.approximate_isotherm,betiso);plot([-.1,F.iso.S.aws],[i,F.iso.S.Cs/F.iso.S.Cs(end),i]);subplot(222);hold on;plot([-.1,F.iso.S.aws],[i/10,F.iso.S.Cs/100,i/10]);subplot(223);hold on;densiteLDPE=922;pvsat27=3564.9138;RsurMeau=461.495;plot(F.iso.S.aws,[i/1000,pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE,i/1000]);subplot(224);hold on;plot(F.iso.S.aws,[i,[pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE]./max([pvsat27/RsurMeau/(273.15+27)*100./F.iso.S.K(2:end)/densiteLDPE]),i]);end;end; %i;// appel sous matlab : cd K:/recherCHE/matlab/senscrep,F=sensc,betiso=feval(F.iso.helperfunctions.baldev02isotherms,'smith',feval(F.iso.helperfunctions.baldev02STARCHandPE,'smith',.4));betiso.awminawmax(3)=-2-1/5;F.iso.S=feval(F.iso.helperfunctions.approximate_isotherm,betiso);F.C0sol=5;F.L=.1;