How to Replace Derivative Terms in Equations
1 次查看(过去 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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/967280/image.jpeg)
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 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Boundary Conditions 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!