Tolerance of ODE15s becomes too small
3 次查看(过去 30 天)
显示 更早的评论
I want to solve the equations:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1611696/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1611701/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1611706/image.png)
I decided to solve this using the finite volume method, so that's why I don't need a boundary condition for the specific volume. I am told that this is a "stiff system", hence I am using ode15s. Do I need to fiddle with the ode15s options to get this thing to work?
%this code solves the following system of equations:
%v_t=u_h
%u_t=(zeta/v*u_h)_h
%dX/dt=u
%h is defined by h_x=1/v, and X is the new position of the interfacial
%points
%The boundary conditions for u is u(t,0)=0, u_h(t,1)=-v(t,1)
%The BC's for v, can be derived from ones on u.
%here u is the velocity, and v is the specific volume defined as 1/rho, where rho is the density
%The system is in conservative form, and the finite volume method is the best method for these types of equations
clear;clc;
S=struct;
%---define physical constants
N=200; %This is the number of cells in the simulation
S.N=N;
S.alpha=0;
c_p=1; %heat capacity.
rho_half(1:N)=0.8; %This is the initial density
eta_0=80;
S.eta=eta_0;
S.nu_m_int=zeros(N+1,1);
%---Set up geometry
x=linspace(0,1,N+1); %Initial spacing of the cell interfaces and ends.
dx=x(2); %initial lengths of cells
L(1)=1; %Initial length
%---The solution of the system will be done via Newton-Raphson. The two
%variables v and u will be put into one vector y=[v u]'
nu(1,:)=1./rho_half; %Initial specific volume
%This interpolates the density from the cell centres to cell boundaries.
rho_0(2:N)=0.5*(rho_half(1,N-1)+rho_half(2:N));
rho_0(1)=1.5*rho_half(1)-0.5*rho_half(2);
rho_0(N+1)=1.5*rho_half(N)-0.5*rho_half(N-1);
h=cumtrapz(x,rho_0); %calculates h, which is used throughout the timesteps.
nu_m(1:N)=1; %This is the matal powder density, the maximum possible density.
S.nu_m=nu_m;
dh=diff(h); %This is used to approximate the derivative
S.dh=dh;
S.h=h;
IC=zeros(2*N,1);
IC(1:N)=nu(1,:)';
%tspan = [0 5];
tspan = [0 1.8];
[t,y] = ode15s(@(t,y) sintering_RHS(y,S),tspan,IC);
figure(1)
plot(t,y(:,180))
figure(2)
plot(t,y(:,380))
function y=flux(u,nu,S)
%This computes the fluxes required for the finite volume method
N=S.N;
dh=S.dh;
k=40;
a=1;
nu_face=0.5*(nu(1:N-1)+nu(2:N)); %This is specific volume on the cell interface
nu_L=1.5*nu(2)-0.5*nu(2);
nu_m=S.nu_m;
nu_m_int=S.nu_m_int;
nu_m_int(2:N)=0.5*(nu_m(1:N-1)+nu_m(2:N));
nu_m_int(N+1)=1.5*nu_m(N)-0.5*nu_m(N-1);
nu_m_int(1)=1.5*nu_m(1)-0.5*nu_m(2);
nu_end=1.5*nu(N)-0.5*nu(N-1);
y=zeros(2,N+1); %This vector will store the fluxes for the mass and momentum
y(1,2:N)=0.5*(u(2:N)+u(1:N-1)); %Fluxes for the mass
y(1,N+1)=u(N)-0.5*k*nu_end*dh(N);
y(2,2:N)=P_L(a,nu_face,nu_m_int(2:N))+(2*zeta(nu_face,nu_m_int(2:N),S)'./nu_face').*((u(2:N)'-u(1:N-1)')./(dh(2:N)+dh(1:N-1))); %The flux for the momentum equation
y(2,1)=P_L(a,nu_L,nu_m_int(1))+(zeta(nu_end,nu_m_int(1),S)/nu_end)*u(1); %the fluxes at the end
y(2,N+1)=-zeta(nu_end,nu_m_int(N+1),S);
end
function y=P_L(a,nu,nu_0) %This is the laplace pressure
x=1-nu_0./nu; %This is the porosity for an imcompressible metal powder
y=a*(1-x').^2;
end
function y=zeta(z,nu_0,S)
eta_0=S.eta;
x=1-nu_0./z;
y = 2*(1-x).^3*eta_0./(3*x);
end
function f=sintering_RHS(y,S)
%the vector f has the stencil for each cell for both the mass and momentum
N=floor(length(y)/2);
f=zeros(1,2*N);
flux_1=flux(y(N+1:2*N),y(1:N),S); %This is flux at timestep i
f_mass=flux_1(1,:);
f_mom=flux_1(2,:);
f(1:N)=f_mass(2:N+1)-f_mass(1:N); %These are the stencils coming from the system of equations
f(N+1:2*N)=f_mom(2:N+1)-f_mom(1:N);
f=f';
end
13 个评论
Torsten
2024-2-13
编辑:Torsten
2024-2-13
As far as I know, the book only treats scheme for hyperbolic systems of the form
du/dt + d(f(u))/dx = s(u,x,t)
I can't remember having seen schemes for 2nd order PDEs.
Your system is really special - a mixture of first and second order PDEs. I don't know either how to discretize it appropriately.
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!