GreenPolyProducts/cantera
5
0
Fork 0
mirror of https://github.com/Cantera/cantera.git synced 2025-02-25 18:55:29 -06:00
Code Issues Packages Projects Releases Wiki Activity

[Matlab] Fixed issue causing crash when getting diffusion coefficients

Browse Source 
Other miscellaneous updates
This commit is contained in:
ssun30 2022-12-09 19:25:31 -05:00 committed by Ray Speth
parent 5adbb0ccf9
commit 4e122e443d
24 changed files with 165 additions and 70 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
interfaces/matlab_experimental/1D/CounterFlowDiffusionFlame.m
View File

@ -167,7 +167,7 @@ function flame = CounterFlowDiffusionFlame(left, flow, right, tp_f, tp_o, oxidiz
% Set the profile of the flame with the estimated axial velocities,
% radial velocities, temperature, and mass fractions calculated above.
flame = Stack([left flow right]);
flame.setProfile(2, {'u', 'V'}, [zrel; u; v]);
flame.setProfile(2, {'velocity', 'spread_rate'}, [zrel; u; v]);
flame.setProfile(2, 'T', [zrel; t] );
for n = 1:nsp
nm = tp_f.speciesName(n);

12
interfaces/matlab_experimental/1D/Domain1D.m
View File

@ -205,7 +205,7 @@ classdef Domain1D < handle
end
end
function n = componentName(d, index)
function s = componentName(d, index)
% Get the name of a component given its index.
%
% n = d.componentName(index)
@ -219,7 +219,7 @@ classdef Domain1D < handle
% Cell array of component names.
%
n = length(index);
s = cell(m);
s = cell(1, n);
for i = 1:n
id = index(i)-1;
output = callct2('domain_componentName', d.domainID, id);
@ -650,7 +650,7 @@ classdef Domain1D < handle
nc = d.nComponents;
for ii = 1:nc
callct('domain_setSteadyTolerances', ...
d.domainID, ii, rtol, atol);
d.domainID, ii-1, rtol, atol);
end
elseif iscell(component)
nc = length(component);
@ -662,7 +662,7 @@ classdef Domain1D < handle
else
n = d.componentIndex(component);
callct('domain_setSteadyTolerances', ...
d.domainID, ii, rtol, atol);
d.domainID, n, rtol, atol);
end
end
@ -686,7 +686,7 @@ classdef Domain1D < handle
nc = d.nComponents;
for ii = 1:nc
callct('domain_setTransientTolerances', ...
d.domainID, ii, rtol, atol);
d.domainID, ii-1, rtol, atol);
end
elseif iscell(component)
nc = length(component);
@ -698,7 +698,7 @@ classdef Domain1D < handle
else
n = d.componentIndex(component);
callct('domain_setTransientTolerances', ...
d.domainID, ii, rtol, atol);
d.domainID, n, rtol, atol);
end
end

4
interfaces/matlab_experimental/1D/Sim1D.m
View File

@ -355,12 +355,12 @@ classdef Stack < handle
for j = 1:np
ic = d.componentIndex(c{j});
callct('sim1D_setProfile', s.stID, ...
n - 1, ic - 1, sz(1), p(1, :), sz(1), p(j+1, :));
n - 1, ic - 1, sz(2), p(1, :), sz(2), p(j+1, :));
end
elseif sz(2) == np + 1;
ic = d.componentIndex(c{j});
callct('sim1D_setProfile', s.stID, ...
n - 1, ic - 1, sz(2), p(:, 1), sz(2), p(:, j+1));
n - 1, ic - 1, sz(1), p(:, 1), sz(1), p(:, j+1));
else
error('Wrong profile shape.');
end

18
interfaces/matlab_experimental/Base/Solution.m
View File

@ -24,12 +24,11 @@ classdef Solution < handle & ThermoPhase & Kinetics & Transport
% >> s = Solution('input.yaml'[, phase_name[, transport_model]])
%
% constructs a :mat:func:`Solution` object from a specification contained in
% file ``input.yaml``. Optionally, the name of the phase to be imported
% can be specified with ``phase_name``. If a :mat:func:`Transport` model is
% included in ``input.yaml``, it will be included in the :mat:func:`Solution`
% instance with the default transport modeling as set
% in the input file. To specify the transport modeling, set the input
% argument ``trans`` to one of ``'default'``, ``'None'``, ``'Mix'``, or ``'Multi'``.
% file ``input.yaml`` with the name of the phase to be imported specified with
% ``phase_name``. If a :mat:func:`Transport` model is included in ``input.yaml``,
% it will be included in the :mat:func:`Solution` instance with the default transport modeling as set
% in the input file. To specify the transport modeling, set the input argument
% ``trans`` to one of ``'default'``, ``'None'``, ``'Mix'``, or ``'Multi'``.
% In this case, the phase name must be specified as well. Alternatively,
% change the ``transport`` node in the YAML file, or ``transport``
% property in the CTI file before loading the phase. The transport
@ -47,16 +46,15 @@ classdef Solution < handle & ThermoPhase & Kinetics & Transport
% :param src:
% Input string of YAML file name.
% :param id:
% Optional unless ``trans`` is specified. ID of the phase to
% import as specified in the YAML file.
% ID of the phase to import as specified in the YAML file.
% :param trans:
% String, transport modeling. Possible values are ``'default'``, ``'None'``,
% ``'Mix'``, or ``'Multi'``. If not specified, ``'default'`` is used.
% :return:
% Instance of class :mat:func:`Solution`
if nargin == 1
id = '-';
if nargin < 2 || nargin > 3
error('Solution class constructor expects 2 or 3 input arguments.');
end
tp = ThermoPhase(src, id);
s@ThermoPhase(src, id);

16
interfaces/matlab_experimental/Base/ThermoPhase.m
View File

@ -68,18 +68,15 @@ classdef ThermoPhase < handle
% t = ThermoPhase(src, id)
%
% :param src:
% Input string of YAML, CTI, or XML file name.
% Input string of YAML file name.
% :param id:
% ID of the phase to import as specified in the input file. (optional)
% ID of the phase to import as specified in the input file.
% :return:
% Instance of class :mat:func:`ThermoPhase`
%
checklib;
if nargin > 2
error('ThermoPhase expects 1 or 2 input arguments.');
end
if nargin == 1
id = '-';
if nargin ~= 2
error('ThermoPhase expects 2 input arguments.');
end
tp.tpID = callct('thermo_newFromFile', src, id);
tp.basis = 'molar';
@ -87,12 +84,9 @@ classdef ThermoPhase < handle
%% Utility methods
function display(tp, threshold)
function display(tp)
% Display thermo properties
if nargin < 2 || ~isnumeric(threshold)
threshold = 1e-14;
end
buflen = 0 - calllib(ct, 'thermo_report', tp.tpID, 0, '', 1);
aa = char(ones(1, buflen));
ptr = libpointer('cstring', aa);

4
interfaces/matlab_experimental/Base/Transport.m
View File

@ -123,7 +123,7 @@ classdef Transport < handle
% symmetric: d(i, j) = d(j, i). Unit: m^2/s.
nsp = tr.th.nSpecies;
xx = zeros(1, nsp);
xx = zeros(nsp, nsp);
pt = libpointer('doublePtr', xx);
callct('trans_getBinDiffCoeffs', tr.trID, nsp, pt);
v = pt.Value;
@ -137,7 +137,7 @@ classdef Transport < handle
% diffusion coefficients. Unit: m^2/s.
nsp = tr.th.nSpecies;
xx = zeros(1, nsp);
xx = zeros(nsp, nsp);
pt = libpointer('doublePtr', xx);
callct('trans_getMultiDiffCoeffs', tr.trID, nsp, pt);
v = pt.Value;

12
samples/matlab_experimental/catcomb.m
View File

@ -132,10 +132,10 @@ surf.T = tsurf;
sim1D = Stack([inlt, flow, surf]);
% set the initial profiles.
sim1D.setProfile(2, {'u', 'V', 'T'}, [0.0, 1.0 % z/zmax
0.06, 0.0 % u
0.0, 0.0 % V
tinlet, tsurf]); % T
sim1D.setProfile(2, {'velocity', 'spread_rate', 'T'}, [0.0, 1.0 % z/zmax
0.06, 0.0 % u
0.0, 0.0 % V
tinlet, tsurf]); % T
names = gas.speciesNames;
for k = 1:gas.nSpecies
@ -207,11 +207,11 @@ sim1D.plotSolution('flow', 'T');
title('Temperature [K]');
subplot(3, 3, 2);
sim1D.plotSolution('flow', 'u');
sim1D.plotSolution('flow', 'velocity');
title('Axial Velocity [m/s]');
subplot(3, 3, 3);
sim1D.plotSolution('flow', 'V');
sim1D.plotSolution('flow', 'spread_rate');
title('Radial Velocity / Radius [1/s]');
subplot(3, 3, 4);

7
samples/matlab_experimental/diff_flame.m
View File

@ -5,13 +5,14 @@
%% Initialization
help diff_flame
clear all
close all
cleanup
clc
tic % total running time of the script
help diff_flame
runtime = cputime; % Record the starting time
%% Parameter values of inlet streams
@ -139,3 +140,5 @@ yyaxis right
plot(z, y(j, :), 'r', z, y(k, :), 'g', z, y(l, :), 'm', z, y(m, :), 'b');
ylabel('Mass Fraction');
legend('T', 'O2', 'H2O', 'C2H6', 'CO2');
toc

10
samples/matlab_experimental/equil.m
View File

@ -3,12 +3,19 @@ function equil(g)
%
% This example computes the adiabatic flame temperature and equilibrium
% composition for a methane/air mixture as a function of equivalence ratio.
clear all
close all
cleanup
clc
tic
help equil
if nargin == 1
gas = g;
else
gas = Solution('gri30.yaml');
gas = GRI30;
end
nsp = gas.nSpecies;
@ -60,4 +67,5 @@ function equil(g)
ylabel('Mole Fraction');
title('Equilibrium Composition');
toc
end

6
samples/matlab_experimental/flame.m
View File

@ -45,9 +45,9 @@ function f = flame(gas, left, flow, right)
else
t1 = right.T;
end
f.setProfile(2, {'u', 'V'}, [0.0 1.0
mdot0/rho0 -mdot1/rho0
0.0 0.0]);
f.setProfile(2, {'velocity', 'spread_rate'}, [0.0 1.0
mdot0/rho0 -mdot1/rho0
0.0 0.0]);
f.setProfile(2, 'T', [0.0, 1.0
t0, t1]);

9
samples/matlab_experimental/flame1.m
View File

@ -5,13 +5,14 @@
%% Initialization
help flame1
clear all
close all
cleanup
clc
tic
help flame1
t0 = cputime; % record the starting time
%% Set parameter values
@ -112,7 +113,7 @@ subplot(2, 2, 1);
plotSolution(fl, 'flow', 'T');
title('Temperature [K]');
subplot(2, 2, 2);
plotSolution(fl, 'flow', 'u');
plotSolution(fl, 'flow', 'velocity');
title('Axial Velocity [m/s]');
subplot(2, 2, 3);
plotSolution(fl, 'flow', 'H2O');
@ -120,3 +121,5 @@ title('H2O Mass Fraction');
subplot(2, 2, 4);
plotSolution(fl, 'flow', 'O2');
title('O2 Mass Fraction');
toc

11
samples/matlab_experimental/flame2.m
View File

@ -4,13 +4,14 @@
%% Initialization
help flame2
clear all
close all
cleanup
clc
tic
help flame2
t0 = cputime; % record the starting time
%% Set parameter values
@ -119,8 +120,10 @@ subplot(2, 3, 4);
plotSolution(fl, 'flow', 'CH');
title('CH Mass Fraction');
subplot(2, 3, 5);
plotSolution(fl, 'flow', 'V');
plotSolution(fl, 'flow', 'spread_rate');
title('Radial Velocity / Radius [s^-1]');
subplot(2, 3, 6);
plotSolution(fl, 'flow', 'u');
plotSolution(fl, 'flow', 'velocity');
title('Axial Velocity [m/s]');
toc

7
samples/matlab_experimental/ignite_hp.m
View File

@ -2,6 +2,12 @@ function ignite_hp(gas)
% IGNITE_HP Solves the same ignition problem as 'ignite', but uses
% function conhp instead of reactor.
clear all
close all
cleanup
clc
tic
help ignite_hp
if nargin == 0
@ -35,4 +41,5 @@ function ignite_hp(gas)
title('OH Mass Fraction');
end
toc
end

7
samples/matlab_experimental/ignite_uv.m
View File

@ -2,6 +2,12 @@ function ignite_uv(gas)
% IGNITE_UV Solves the same ignition problem as 'ignite2', except
% function conuv is used instead of reactor.
%
clear all
close all
cleanup
clc
tic
help ignite_uv
if nargin == 0
@ -35,4 +41,5 @@ function ignite_uv(gas)
title('OH Mass Fraction');
end
toc
end

8
samples/matlab_experimental/isentropic.m
View File

@ -4,6 +4,12 @@ function isentropic(g)
% In this example, the area ratio vs. Mach number curve is
% computed for a hydrogen/nitrogen gas mixture.
%
clear all
close all
cleanup
clc
tic
help isentropic
if nargin == 1
@ -56,4 +62,6 @@ function isentropic(g)
ylabel('Area Ratio');
xlabel('Mach Number');
title('Isentropic Flow: Area Ratio vs. Mach Number');
toc
end

24
samples/matlab_experimental/lithium_ion_battery.m
View File

@ -1,3 +1,5 @@
% LITHIUM_ION_BATTERY
%
% This example file calculates the cell voltage of a lithium-ion battery
% at given temperature, pressure, current, and range of state of charge (SOC).
%
@ -17,22 +19,25 @@
% open-source software Cantera,” Electrochim. Acta 323, 134797 (2019),
% https://doi.org/10.1016/j.electacta.2019.134797
% -----------------------------------------------------------------------------
% Input
% -----------------------------------------------------------------------------
%% Initialization
clear all
close all
cleanup
clc
% Operation parameters
tic
help lithium_ion_battery
%% Operation parameters
SOC = 0:0.02:1; % [-] Input state of charge (0...1) (can be a vector)
I_app = -1; % [A] Externally-applied current, negative for discharge
T = 293; % [K] Temperature
P = oneatm; % [Pa] Pressure
% Cell properties
%% Cell properties
inputFile = 'lithium_ion_battery.yaml'; % Cantera input file name
R_elyt = 0.0384; % [Ohm] Electrolyte resistance
S_ca = 1.1167; % [m^2] Cathode total active material surface area
@ -46,9 +51,7 @@ X_Li_an_1 = 0.75; % [-] anode Li mole fraction at SOC = 100 %
X_Li_ca_0 = 0.99; % [-] cathode Li mole fraction at SOC = 0 %
X_Li_ca_1 = 0.49; % [-] cathode Li mole fraction at SOC = 100 %
% -----------------------------------------------------------------------------
% Calculations
% -----------------------------------------------------------------------------
%% Calculations
% Calculate mole fractions from SOC
X_Li_an = (X_Li_an_1-X_Li_an_0)*SOC+X_Li_an_0; % anode balancing
@ -102,10 +105,9 @@ xlabel('State of charge / %')
ylabel('Cell voltage / V')
set(gca, 'fontsize', 14)
toc
%--------------------------------------------------------------------------
% Helper functions
% -----------------------------------------------------------------------------
%% Helper functions
% This function returns the Cantera calculated anode current (in A)
function anCurr = anode_curr(phi_s, phi_l, X_Li_an, anode, elde, elyt, anode_interface, S_an)

11
samples/matlab_experimental/periodic_cstr.m
View File

@ -1,6 +1,6 @@
function periodic_cstr
%
% A CSTR with steady inputs but periodic interior state.
% PERIODIC_CSTR - A CSTR with steady inputs but periodic interior state.
%
% A stoichiometric hydrogen/oxygen mixture is introduced and reacts to produce
% water. But since water has a large efficiency as a third body in the chain
@ -18,8 +18,13 @@ function periodic_cstr
% conditions needed to see the oscillations came from Bob Kee,
% Colorado School of Mines
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
close all
cleanup
clc
tic
help periodic_cstr
% create the gas mixture
@ -102,4 +107,6 @@ function periodic_cstr
plot(tm,y)
legend('H2','O2','H2O')
title('Mass Fractions')
toc
end

7
samples/matlab_experimental/PFR.m → samples/matlab_experimental/plug_flow_reactor.m
View File

@ -21,13 +21,17 @@
% Virginia Tech
%% Clear all variables, close all figures, clear the command line:
clear all
close all
cleanup
clc
tic
help PFR
%% Operation Parameters
% Temperature of gas, in K
T0 = 1473;
% Pressure of gas, in Pa
@ -54,6 +58,7 @@ gas_calc.TPX = {T0, P0, x};
gas_calc.equilibrate('HP');
%% Calculation of properties along the reactor length
% The Dimensions and conditions of the reactor are given below
% Inlet Area, in m^2
@ -166,3 +171,5 @@ plot(x_calc, P_calc)
xlabel('X-Position (m)')
ylabel('Pressure (Pa)')
title('Pressure Variation')
toc

8
samples/matlab_experimental/prandtl1.m
View File

@ -5,13 +5,14 @@ function prandtl1(g)
% hydrogen / oxygen mixture in chemical equilibrium for P = 1
% atm and a range of temperatures and elemental O/(O+H) ratios.
help prandtl1
clear all
close all
cleanup
clc
tic
help prandtl1
if nargin == 1
gas = g;
else
@ -78,4 +79,7 @@ function prandtl1(g)
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Thermal Conductivity');
toc
end

6
samples/matlab_experimental/prandtl2.m
View File

@ -4,13 +4,15 @@ function prandtl2(g)
% This example does the same thing as prandtl1, but using
% the multicomponent expression for the thermal conductivity.
%
help prandtl2
clear all
close all
cleanup
clc
tic
help prandtl2
if nargin == 1
gas = g;
else
@ -76,4 +78,6 @@ function prandtl2(g)
xlabel('Elemental O/(O+H)');
ylabel('Temperature (K)');
zlabel('Thermal Conductivity');
toc
end

8
samples/matlab_experimental/reactor1.m
View File

@ -4,6 +4,13 @@ function reactor1(g)
% This example illustrates how to use class 'Reactor' for zero-dimensional
% kinetics simulations. Here the parameters are set so that the reactor is
% adiabatic and very close to constant pressure.
clear all
close all
cleanup
clc
tic
help reactor1
if nargin == 1
@ -82,4 +89,5 @@ function reactor1(g)
clear all
cleanup
toc
end

8
samples/matlab_experimental/reactor2.m
View File

@ -5,6 +5,12 @@ function reactor2(g)
% zero-dimensional kinetics simulations. Here the parameters are
% set so that the reactor is adiabatic and constant volume.
clear all
close all
cleanup
clc
tic
help reactor2
if nargin == 1
@ -57,4 +63,6 @@ function reactor2(g)
ylabel('H2 Mole Fraction (K)');
clear all
cleanup
toc
end

12
samples/matlab_experimental/surfreactor.m
View File

@ -4,17 +4,21 @@
% zero-dimensional simulations including both homogeneous and
% heterogeneous chemistry.
help surfreactor
%% Initialization
clear all
close all
cleanup
clc
tic
help surfreactor
%% Set the initial conditions
t = 870.0;
gas = Solution('ptcombust.yaml','gas');
% set the initial conditions
gas.TPX = {t, oneatm, 'CH4:0.01, O2:0.21, N2:0.78'};
% The surface reaction mechanism describes catalytic combustion of
@ -78,6 +82,8 @@ for n = 1:nSteps
end
disp(['CPU time = ' num2str(cputime - t0)]);
%% Plotting
clf;
subplot(2,2,1);
plot(tim,temp);
@ -100,3 +106,5 @@ ylabel('Mole Fractions');
legend(names);
clear all
cleanup
toc

18
samples/matlab_experimental/test_examples.m
View File

@ -1,12 +1,28 @@
% runs selected examples without pausing
LoadCantera
clear all
close all
cleanup
equil();
isentropic();
reactor1();
reactor2();
surfreactor;
periodic_cstr;
plug_flow_reactor;
lithium_ion_battery
rankine(300.0, 2.0*oneatm, 0.8, 0.7);
prandtl1();
prandtl2();
flame1;
flame2;
catcomb;
exit;
diff_flame;
ignite_hp;
ignite_uv;
clear all
close all
cleanup
UnloadCantera