How to Replace Derivative Terms in Equations
2 次查看(过去 30 天)
显示 更早的评论
clc,clear
%已知氢气在空气中浓度达到4.0%~75.6%就会爆炸
%假设初始条件:t = 0; H2 : O2 : N2 = 0.1 : 0.198 : 0.702
c = sym('c',[3,1]);%the units of concentration should be volume fraction
syms t z
c(1) = str2sym('c1(t,z)');%h2
c(2) = str2sym('c2(t,z)');%o2
c(3) = str2sym('c3(t,z)');%n2
c0 = [0.1;0.198;0.702];
% c = [c1;c2;c3];
%%
%parameters
%units of mass:kg; units of distance:km;
%moleculer weight
m_H = 1.67E-27;%kg
m = [2 * m_H;32 * m_H;28 * m_H];
%reduced mass
for i = 1:3
for j = 1:3
if i~=j
mu(i,j) = m(i) * m(j) / (m(i) + m(j));
else
mu(i,j) = 0;
end
end
end
%%
%the sum of kinetic raddi
k_r = [289E-15;364E-15;346E-15];%km
for i = 1:3
for j = 1:3
if i~=j
sigma(i,j) = k_r(i) + k_r(j);
else
sigma(i,j) = 0;
end
end
end
%%
%Boltzmann constant
k = 1.381E-23;%J/K
%T-z
dTdz = -6.5;%K/km
T0 = 293.15;%K
T = T0 + dTdz * z;
%scale height
g = 9.8;%m/s^2
R = 8.314;%J/mol/K
M = [2;32;28];
sum_c = sum(c);
sum_c0 = sum(c0);
%%
for i = 1:3
%thermal diffusion factor
alpha(i) = -log((c(i) ./ c0(i)) .* ((sum_c0 - c0(i)) ./ (sum_c - c(i)))) ./ log(T ./ T0);
%individual height
Hi(i) = (R * T) / (M(i) * g);
for j = 1:3
if j ~= i
D_d(i,j) = c(j) .* (16/3) .* (mu(i,j)./m(i)) .* sigma(i,j).^2 .* (pi.*k.*T./(2.*mu(i,j))).^0.5
end
end
end
for i = 1:3
% moleculer diffusion coefficent
sumD_d(i) = sum(D_d(i,:))
D(i) = Hi(i) .* g ./ sumD_d(i)
%velocity-moleculer diffusion
v_m(i) = D(i) .* ((1 ./ c(i))*diff(c(i),z) + (1 + alpha(i)) ./ T .* dTdz + 1 ./ Hi(i))
end
for i = 1:3
%eddy diffsion coefficient--< Transport phenomena and the chemical composition in the mesosphere and lower thermosphere>
K_M = 3.8E6;%cm2s-1
K_m = 6E5;
S2 = 0.18;%km-2
S3 = 0.098;%km-1
z0 = 102;%km
K = (K_M - K_m) * exp(-S2 * (z - z0).^2) + K_m * exp(-S3 * (z - z0));
M_air = 29;
H = (R * T) / (M_air * g);
%velocity-eddy diffusion
v_e(i) = K .* ((1 ./ c(i)) .* diff(c(i),z) + dTdz ./ T + 1 ./ H)
end
v = -v_m - v_e;
%%
for i = 1:3
eqn(i) = diff(c(i),t) == -diff(c(i) .* v(i),z)
end
%%
syms delta_z
C = [];
% for i = 1:3
eqn1 = subs(eqn(1),diff(c1(t, z), z),(C(i+1,n) - C(i,n)) ./ delta_z)
Matlab discretes the differential equation and replaces the differential term with the differential term. How can we do it ?The subs is not useful.
That is in an equation, replace the differential term on the left with the difference term on the right.
3 个评论
Torsten
2022-4-17
What do you intend with your symbolic setup ?
If in the end, you replace the differentials with finite-difference approximations and calculate the concentrations numerically, I'd skip all the symbolic calculations and set up the problem numerically right from the beginning.
回答(0 个)
另请参阅
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!