Getting error while solving coupled PDEs

2021 12 24
1 个回答
更新时间：2021 12 25
5 次查看（30 天）
Ran in:

0 个投票

Ran in:
clc
clear all
m=50;
delx=1/(m-1);
h_int=0.0005e7;
final_time=1e7;
N_samples=(final_time/h_int)+1;
t0=0;
tf=final_time;
p_mat(1,:)=2*10^5;
T_mat(1,:)=287;
[t,p_mat] = ode15s(@GH,[t0 tf],p_mat(:,1));
Unable to perform assignment because the left and right sides have a different number of elements.

Error in solution>GH (line 45)
F_X(i)=[dp_dt(i) dT_dt(i)]';

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 152)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
[t,T_mat] = ode15s(@GH,[t0 tf],T_mat(:,1));
plot(t, p_mat);
% plot(t, T_mat);
xlabel('time'),ylabel('pressure');
%%
function F_X= GH(t,p,T,vg)
m=50;
delx=1/(m-1);
F_X = zeros(2,m);
dp_dt=zeros(1,m);
dT_dt=zeros(1,m);
k=0.98*10^-12;
mu=0.001;
phi=0.4;
K1=k/((1-phi)*mu);
Thcond=0.034;
Cp=2232;
K2=Thcond/(0.657*Cp);
pg=2*10^5;
ph=30*10^5;
p(1,:)=pg;
p(m,:)=ph;
T(1,:)=287;
T(m,:)=287;
vg(1,:)=0.01;
vg(m,:)=0;
for i=2:m-1
dp_dt(i)=K1*((p(i+1)-2*p(i)+p(i-1))/(delx)^2);% K1 is k/(1-phi)*mu
vg(i)=(k/mu)*(p(i+1)-p(i))/(delx);
dT_dt(i)=(K2*((T(i+1)-2*T(i)+T(i-1))/(delx)^2))-(vg(i)*((T(i+1)-T(i-1))/(2*delx)));% K2 is thermal diffusivity
F_X(i)=[dp_dt(i) dT_dt(i)]';
end
end

3 个评论
显示 1更早的评论 隐藏 1更早的评论

Torsten
Torsten 2021-12-24
Why do you include a heat balance although there is no source that affects temperature ? The temperature will remain at 287 K in the domain for all times.
MANISH KUMAR GUPTA
MANISH KUMAR GUPTA 2021-12-24
Temperature will decrease to 273 from 287 due to joule thomson cooling. I am unable to make function with couple odes with loop.
I am making mistake at this part, and unable to solve.
F_X(i)=[dp_dt(i) dT_dt(i)]'
Torsten
Torsten 2021-12-24
There is no energy sink included in your equations.
Furthermore, velocity must be negative because the fluid flows from high to low pressure.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Walter Roberson
Walter Roberson 2021-12-24

0 个投票

F_X = zeros(2,m);
That is a 2 x something array.
F_X(i)=[dp_dt(i) dT_dt(i)]';
You assign a pair of values to the same single location in the array, same as if you had assigned to F_X(i,1)
Perhaps you want
F_X(:,i)=[dp_dt(i) dT_dt(i)]';
if so then are you sure you would want the first and last columns to be 0 (because the loop does not assign to all columns) ?
Also remember that ode functions must return a column vector, never a 2D array or row vector.

3 个评论
显示 1更早的评论 隐藏 1更早的评论

MANISH KUMAR GUPTA
MANISH KUMAR GUPTA 2021-12-24
Hi
Thanks for the reply, I was trying to solve two coupled odes , dp_dt(i) and dT_dt(i), i am unable to assign function for these coupled odes.
dp_dt(i)=K1*((p(i+1)-2*p(i)+p(i-1))/(delx)^2);% K1 is k/(1-phi)*mu
vg(i)=(k/mu)*(p(i+1)-p(i))/(delx);
dT_dt(i)=(K2*((T(i+1)-2*T(i)+T(i-1))/(delx)^2))-(vg(i)*((T(i+1)-T(i-1))/(2*delx)));
Walter Roberson
Walter Roberson 2021-12-25
Ran in:
Ran in:
clc
clear all
m=50;
delx=1/(m-1);
h_int=0.0005e7;
final_time=1e7;
N_samples=(final_time/h_int)+1;
t0=0;
tf=final_time;
p_mat(1,:)=2*10^5;
T_mat(1,:)=287;
[t,p_mat] = ode15s(@GH,[t0 tf],p_mat(:,1));
Error using odearguments (line 95)
GH returns a vector of length 100, but the length of initial conditions vector is 1. The vector returned by GH and the initial conditions vector must have the same number of elements.

Error in ode15s (line 152)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
[t,T_mat] = ode15s(@GH,[t0 tf],T_mat(:,1));
plot(t, p_mat);
% plot(t, T_mat);
xlabel('time'),ylabel('pressure');
%%
function F_X= GH(t,p,T,vg)
m=50;
delx=1/(m-1);
F_X = zeros(2,m);
dp_dt=zeros(1,m);
dT_dt=zeros(1,m);
k=0.98*10^-12;
mu=0.001;
phi=0.4;
K1=k/((1-phi)*mu);
Thcond=0.034;
Cp=2232;
K2=Thcond/(0.657*Cp);
pg=2*10^5;
ph=30*10^5;
p(1,:)=pg;
p(m,:)=ph;
T(1,:)=287;
T(m,:)=287;
vg(1,:)=0.01;
vg(m,:)=0;
for i=2:m-1
dp_dt(i)=K1*((p(i+1)-2*p(i)+p(i-1))/(delx)^2);% K1 is k/(1-phi)*mu
vg(i)=(k/mu)*(p(i+1)-p(i))/(delx);
dT_dt(i)=(K2*((T(i+1)-2*T(i)+T(i-1))/(delx)^2))-(vg(i)*((T(i+1)-T(i-1))/(2*delx)));% K2 is thermal diffusivity
F_X(:, i)=[dp_dt(i) dT_dt(i)]';
end
F_X = F_X(:);
end
Walter Roberson
Walter Roberson 2021-12-25
Your function definition implies you are passing in T and vg, but you do not pass those in.
You are not working with two coupled differential equations: you are working with m*2 = 50*2 = 100 coupled differential equations. So you need 100 initial conditions. Those will be the initial p and T values. It looks to me as if you should not be passing T into the function, but should instead be using
function F_X= GH(t,boundaries,vg)
p = boundaries(1:2:end);
T = boundaries(2:2:end);
m = length(p);
delx=1/(m-1);
F_X = zeros(2,m);
dp_dt=zeros(1,m);
dT_dt=zeros(1,m);
k=0.98*10^-12;
mu=0.001;
phi=0.4;
K1=k/((1-phi)*mu);
Thcond=0.034;
Cp=2232;
K2=Thcond/(0.657*Cp);
pg=2*10^5;
ph=30*10^5;
p(1)=pg;
p(end)=ph;
T(1)=287;
T(end)=287;
vg(1,:)=0.01;
vg(end,:)=0;
for i=2:m-1
dp_dt(i)=K1*((p(i+1)-2*p(i)+p(i-1))/(delx)^2);% K1 is k/(1-phi)*mu
vg(i)=(k/mu)*(p(i+1)-p(i))/(delx);
dT_dt(i)=(K2*((T(i+1)-2*T(i)+T(i-1))/(delx)^2))-(vg(i)*((T(i+1)-T(i-1))/(2*delx)));% K2 is thermal diffusivity
F_X(:, i) = [dp_dt(i) dT_dt(i)];
end
F_X = F_X(:);
end
I am not sure why you are storing all of the values of vg. You only ever use the "current" vg, vg(i), and that only after storing into it. You could use an unindexed variable instead, and not pass in vg at all.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by