Tolerance of ODE15s becomes too small

3 次查看（过去 30 天）

显示更早的评论

Mat Hunt 2024-2-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2079886-tolerance-of-ode15s-becomes-too-small

编辑： Torsten 2024-2-13

ode_sintering.m

在 MATLAB Online 中打开

I want to solve the equations:

with boundary conditions

, the initial conditions are

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