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
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.
Hou X.Y
Hou X.Y 2022-4-18
Yeah,I want to use finite-difference approximations to calculate the concentrations numerically,and the initial equation is complicated,so I use the symbolic calculations to deduce the final term of the diffrential function which can be used for calclations numerically.

请先登录,再进行评论。

回答(0 个)

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by