ODE15s doesn't seem to be working, one of the equations seem to be unchanged

Question

Mat 2024-11-12，16:38

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2165674-ode15s-doesn-t-seem-to-be-working-one-of-the-equations-seem-to-be-unchanged

评论： Mat 2024-11-14，10:11

sintering_test.m

在 MATLAB Online 中打开

I want to solve the following set of equations:

via the finite volume method, so I integrate the above equations over an inteval [h_{i},h_{i+1}] to get a set of ODEs which I then proceed to solve using ODE15s, as this appears to be a stiff set of equations. The attached code seems to be getting the change in u, which is good, BUT, gradients in u should give rise to a RHS for the first equation and therefore evolves nu, but it doesn't. I'm at a loss why. I've checked the fluxes, which should be the source of the issue but they seem fine.

Can anyone point out any trivial errors?

%This is a code to solve the isothermal sintering equations.

S=struct;

%%---Physical parameters---

S.N = 251; %This is the number of cells-1

S.L = 1; %This is the initial length of the rod

S.g = 9.81; %Acceleration due to gravity.

S.K = 1e-2; %This is the Laplace constant from the sintering potential.

S.eta_0 = 0.05; %The shear stress of the skeleton

S.T = 8; %This is the time of sintering.

S.rho_m = 1; %This is the density of the individual metallic particles.

%%---Initial conditions

%Length

S.X=linspace(0,S.L,S.N); %This is the partition of the initial length

X = S.X;

%Density

rho_0_int=0.5-0.1*exp(-500*(X-0.5).^2);

rho_0=0.5*(rho_0_int(1:S.N-1)+rho_0_int(2:S.N));

%Amount of mass as a function of position

h_interface=cumtrapz(X,rho_0_int);

S.dh=diff(h_interface);

%Mass at the midpoint;

S.h=0.5*(h_interface(1:S.N-1)+h_interface(2:S.N));

%Total mass

S.M = h_interface(end);

U=S.M/(S.rho_m*S.T); %Scaling for velocity

a_1=S.T/(U*S.M);

a_2=S.eta_0*S.T*S.rho_m/S.M^2;

a_3=S.g*S.T/U;

%Initial velocity

u0=zeros(S.N-1,1);

% Define the system of equations as a function

y0 = [1./rho_0'; u0];

% Finite Volume Method solution loop

tspan = [0 1]; % You can adjust this based on the time range you're interested in

[t, y] = ode15s(@(t,y)(ode_system(t,y,S)), tspan, y0);

rho=1./y(:,1:S.N-1);

u=y(:,S.N:2*S.N-2);

Tspan=length(y(:,1));

disp('PDEs solved, now computing new grid');

PDEs solved, now computing new grid

figure()

hold on

for i=1:Tspan

plot(u(i,:))

end

hold off

function dydt = ode_system(t, y, S)

U=S.M/(S.rho_m*S.T); %Scaling for velocity

a_1=S.T/(U*S.M);

a_2=S.T*S.rho_m/S.M^2;

a_3=S.g*S.T/U;

nu = y(1:S.N-1)'; % specific volume

u = y(S.N:2*(S.N-1))'; % velocity

%Compute the left and right specific volumes:

nu_left=nu(1); %15*nu(1)/8-5*nu(2)/4+3*nu(3)/8;

nu_right=nu(S.N-1); %15*nu(S.N-3)/8-5*nu(S.N-2)/5+3*nu(S.N-1)/8;

%Compute the fluxes on the interfaces

nu_half=NaN(1,S.N);

nu_half(2:S.N-1)=0.5*(nu(2:S.N-1)+nu(1:S.N-2));

nu_half(1)=nu_left;

nu_half(S.N)=nu_right;

%Fluxes for nu equation

nu_flux=zeros(S.N,1);

nu_flux(2:S.N-1)=0.5*(u(1:S.N-2)+u(2:S.N-1));

du_dh_left=-P_L(S.K,nu_left)*nu_left/zeta(nu_left);

du_dh_right=-P_L(S.K,nu_right)*nu_right/zeta(nu_right);

a=P_L(S.K,nu_right);

b=zeta(nu_right);

nu_flux(1)=-3*du_dh_left*S.h(1)/8+9*u(1)/8-u(2)/8;

nu_flux(S.N)=3*du_dh_right*S.dh(S.N-1)/8-9*u(S.N-1)/8+u(S.N-2)/8;

%Fluxes for u equation

u_grad=2*(u(2:S.N-1)-u(1:S.N-2))./(S.dh(2:S.N-1)+S.dh(1:S.N-2));

u_flux=zeros(1,S.N);

u_flux(2:S.N-1)=a_1*P_L(S.K,nu_half(2:S.N-1))+(a_2*zeta(nu_half(2:S.N-1))./nu_half(2:S.N-1)).*u_grad;

% The system of equations

dnu_dt = nu_flux(2:S.N)-nu_flux(1:S.N-1);

du_dt =u_flux(2:S.N)'-u_flux(1:S.N-1)';

% Return the derivatives

dydt = [dnu_dt; du_dt];

end

function dudh = du_dh(y,S)

dudh=NaN(S.N,1);

rho=y(1:S.N-1);

u=y(S.N:end);

dudh(2:S.N-1)=(u(2:S.N-1)-u(1:S.N-2))./(0.5*(S.dh(1:S.N-2)+S.dh(2:S.N-1)));

end

function y=P_L(K,nu)

y=K./nu;

end

function y=zeta(nu)

y=(2/3)*(nu).^(-2)./(nu-1);

end

function x=compute_grid(rho,rho_0,X)

%This function computes the new grid at each time step

%The equation for the new grid is given by: x_X=rho_0/rho, so we solve this

end

6 个评论
显示 4更早的评论隐藏 4更早的评论

Torsten 2024-11-13，14:39

It's hard to understand your discretization.

I'm also not sure whether you need a boundary value for u to be used in the first equation dnu/dt = du/dh.

It's a complicated system.

Mat 2024-11-14，10:11

isothermal_sintering.pdf

It is a complicated system, it's a mixed hyperbolic-parabolic system.

I need a boundary condition for u to get the end points for u, it's a Neumann to Diriclet map basically. I have included the discritisation.

请先登录，再进行评论。

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

ODE15s doesn't seem to be working, one of the equations seem to be unchanged

6 个评论
显示 4更早的评论隐藏 4更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

ODE15s doesn't seem to be working, one of the equations seem to be unchanged

6 个评论 显示 4更早的评论隐藏 4更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

6 个评论
显示 4更早的评论隐藏 4更早的评论