Getting error while solving coupled PDEs

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

[Torsten]

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]

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]

There is no energy sink included in your equations.

Furthermore, velocity must be negative because the fluid flows from high to low pressure.

Answer 1

[Walter Roberson]

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.

[MANISH KUMAR GUPTA]

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]

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]

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.