Code has numerical issues that I don't understand
1 visualizzazione (ultimi 30 giorni)
Mostra commenti meno recenti
I managed to get a system of PDEs down to a large system of ODEs via the method of lines. I managed to code it up using ode45 but I still have numerical errors when the code is run, and I don't quite understand why. The basic physics of the model is that of a porous 1D bar where the porosity and length decreases, whereas the velocity and temperature increases. I think I have reasonable initial conditions and I've checked my discritisation and it looks okay. I get te error:
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.
> In odemassexplicit>ExplicitSolverHandleMass34 (line 67)
In odefcncleanup>@(t,y)oldFcn(t,y,inputArgs{:}) (line 11)
In ode45 (line 315)
In sintering (line 12)
clear;clc;
N=10;
c=1;
h=1/N;
tspan=[0 20];
y_0=zeros(1,3*N+1);
y_0(1:N)=1;
y_0(N+1:2*N)=0;
y_0(2*N+1:3*N)=1;
y_0(3*N+1)=1;
opts = odeset('Mass',@(t,y) mass(y,h,c));
[t,y] = ode45(@(t,y) sint(N,y),tspan,y_0,opts);
function dydt = sint(N,y)
%This is the RHS of the sintering equations
g=9.81;
alpha=0.1;
h=1/N;
kappa=1;
x=linspace(0,1,N);
dydt=zeros(3*N+1,1);
dydt(3*N+1)=y(2*N);
%These are the bulk equations
for i=2:N-1
dydt(i)=y(2*N)*(y(i+1)-y(i-1))/(2*h*y(3*N+1))-(y(3*N+1)/(2*h))*(y(i+1)*y(N+1+1)-y(i-1)*y(N+i-1));
A=(y(2*N)*y(i)/(2*h*y(3*N+1)))*(y(i+1)-y(i-1))-(y(i)*y(N+i)/(2*h*y(3*N+1)))*(y(N+i+1)-y(N+i-1))+1/(2*h*y(3*N+1))*(P_L(y(i+1))-P_L(y(i-1)));
C=(1/(h^2*y(3*N+1)^2))*0.5*(zeta(y(i+1))+zeta(y(i)))*y(N+i+1)-0.5*(zeta(y(i+1))+zeta(y(i))+zeta(y(i))+zeta(y(i-1)))*y(N+i)+0.5*(zeta(y(i+1))+zeta(y(i-1)))*y(N+i-1);
D=-g+alpha/(2*h*y(3*N+1))*(y(2*N+i+1)-y(2*N+i-1));
dydt(N+i)=A+C+D;
E=-((y(N+i)/y(3*N+1))-y(2*N)/y(3*N+1).^2)*((y(N+i+1)-y(N+i-1))/(2*h));
F=(kappa*y(i)/(h^2*y(3*N+1)))*(y(2*N+i+1)-2*y(2*N+i)+y(2*N+i-1))-(P_L(y(i))/y(3*N+1))*((y(N+i+1)-y(N+i-1))/(2*h));
G=(alpha*y(i)*y(2*N+i)/h^2)*(y(N+i)/y(3*N+1)^2-y(2*N)/y(3*N+1))*(y(N+i+1)-2*y(N+i)+y(N+i-1));
dydt(2*N+i)=E+F+G;
end
%This is the boundary condition for the density;
dydt(1)=0;
dydt(N)=y(2*N)*(y(N)-y(N-1))/(h*y(3*N+1));
%This is the boundary condition for the velocity
dydt(N+1)=(P_L(y(2))-P_L(y(1)))/(h*y(3*N+1))+(y(N+2)*(zeta(y(2))-zeta(y(1))))/(h^2*y(3*N+1)^2)-g+alpha*(y(2*N+2)-y(2*N+1))/(h*y(3*N+1));
A=y(2*N)*y(N)*(y(N)-y(N-1))/(h*y(3*N+1))+2*y(N)*y(2*N)*(P_L(y(N))+alpha*T_a)/(y(3*N+1)*zeta(y(N)))+(P_L(y(N))-P_L(y(N-1)))/(h*y(3*N+1));
B=(0.5*(3*zeta(y(N))-zeta(y(N-1)))*(y(2*N-1)-2*h*(P_L(y(N)))+alpha*T_a)/zeta(y(N))+0.25*y(2*N)*(7*zeta(y(N))-zeta(y(N-1)))+0.5*y(2*N-1)*(zeta(y(N))+zeta(y(N-1))))/(h^2*y(3*N+1)^2);
C=-g+alpha*(y(3*N)-y(3*N-1))/(h*y(3*N+1));
dydt(2*N)=A+B+C;
%This is the boundary condition for Temperature
dydt(2*N+1)=y(1)*(kappa*y(3*N+1)-2*kappa*T_a+kappa*(2*T_a-y(3*N+1)))/(h^(2)*y(3*N+1)^2)-P_L(y(1))*y(N+2)/h;
dydt(3*N)=y(N)*(kappa*(2*T_a)-y(3*N-1)-2*kappa*T_a+kappa*y(3*N-1))-P_L(y(N))*(y(2*N)-y(2*N-1))/(h*y(3*N+1));
end
function M=mass(y,h,c)
N=length(y);
n=floor((N-1)/3);
l=ones(1,n);
M=zeros(N,N);
%Insert the conservation of mass terms
M(1:n,1:n)=diag(l,0);
%Insert the coefficients for the conservation of momentum
M(n+1:2*n,n+1:2*n)=diag(y(1:n),0);
%Insert the coefficients for the conservation of energy
M(2*n+1:3*n,2*n+1:3*n)=diag(c*y(1:n),0); %Diagonal elements
%Compute the off diagonal elements
l_sub=ones(1,n-1)/(2*h);
M_sub=diag(l_sub,-1)-diag(l_sub,1);
M_sub(1,1)=1/h;
M_sub(1,2)=-1/h;
M_sub(n,n-1)=1/h;
M_sub(n,n)=-1/h;
M((2*n+1):3*n,(n+1):2*n)=M_sub;
M(N,N)=1;
end
function y=P_L(x)
gamma=1;
r_0=1;
y=3*gamma*(1-x).^2/r_0;
end
function y=zeta(x)
eta_0=1;
y=2*(1-x).^3*eta_0/(3*x);
end
function y=T_a(t)
y=1;
end
Can you suggest might be wrong?
11 Commenti
Torsten
il 24 Mar 2023
Modificato: Torsten
il 24 Mar 2023
I know how to discretise these equations, they should work.
Believe me: these equations are very difficult to solve and adequate discretization schemes fill books in the literature.
One is tempted to use standard discretization schemes, but they won't work.
As long as velocity and pressure are given, standard discretization schemes work (e.g. for species or temperature). As soon as you try to add the continuity and momentum equations with a naive approach, the results become unphysical.
But the best way to learn is to fail - and finally to know why one failed.
Risposte (0)
Vedere anche
Categorie
Scopri di più su Boundary Conditions in Help Center e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!