-1
$\begingroup$

I have to solve $$\left\{\begin{matrix} \ddot{R}=-\frac{3}{2}\frac{\dot{R}^{2}}{R}+\frac{1}{R\rho}\left [P_{v}(T_{\infty})-P_{0}-a_{1}t^{2}-a_{2}t+\frac{(L\rho_{g})^{2}}{hT_{\infty}}\dot{R}-4\mu \frac{\dot{R}}{R}-\frac{2S}{R} \right ]\\ h=\frac{2\hat{\sigma}}{2-\hat{\sigma}}\frac{L^{2}\rho_{g}}{T_{B}}\left ( \frac{\bar{M}}{2\pi \bar{R}T_{B}} \right )^{1/2}\left [ 1-\frac{1}{2L\rho_{g}}\left ( P_{v}(T_{\infty})+\frac{(L\rho_{g})^{2}}{hT_{\infty}}\dot{R} \right ) \right ]\\ \ln\left ( \frac{\rho_{g}}{\rho_{c}}\right )=c_{1}\theta^{1/3}+c_{2}\theta^{2/3}+c_{3}\theta^{4/3}+c_{4}\theta^{3}+c_{5}\theta^{37/6}+c_{6}\theta^{71/6} \end{matrix}\right.$$ where $$\left\{\begin{matrix} R(0)=R_{0}\\ \dot{R}(0)=0\\ T_{B}=T_{\infty}+\frac{L\rho_{g}}{h}\dot{R} \end{matrix}\right.$$

(I do not know if this may be helpful, but note the order of magnitude of $a_{1}$ ($10^{14}$) and of time span ($10^{-6}$). The problem is that I receive a warning: "Warning: Failure at t=5.142410e-005. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.826951e-019) at time t." How can I fix it? I tried switching to ode23t but that did not solve the problem. enter image description here

function out=densita_20_latente_christ(t,y,R0)
S = 0.073; %N/m
rho = 998;    
mi = 1.005e-3;
P0 = 101325;
L = 2454000;
Tinf = 293;   
lambda = 0.5982;
cp = 4184;   
Pv_Tinf = 2329.6;
sigma_hat=0.03;
M_bar=18;
R_bar=8314;
c1 = -2.03150240;
c2 = -2.6830294;
c3 = -5.38626492;
c4 = -17.2991605;
c5 = -44.7586581;
c6 = -63.9201063;
Tc = 647.096;
rhoc = 322;
a1 = 2.488125e14;
a2 = -9.9525e9;
%
out=[y(2)
    -1.5 * y(2)^2 / y(1) + 1 / rho / y(1) * (Pv_Tinf - a1 * t^2 - a2 * t - P0 +...
    (L * y(4))^2 / Tinf * y(2) / y(3) -...
    4 * mi * y(2) / y(1) - 2 * S / y(1))
     -y(3) + 2 * sigma_hat / (2 - sigma_hat) *...
     L^2 * y(4) / (Tinf + L * y(4) * y(2) / y(3))*...
    (M_bar / (2 * pi * R_bar * (Tinf+ L * y(4) * y(2) / y(3))))^0.5 *...
    (1 - 1 / (2 * L * y(4)) *...
    (Pv_Tinf + (L * y(4))^2 * y(2) / y(3) / Tinf))
    -y(4) + rhoc * exp(c1*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(1/3)...
    + c2*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(2/3)...
    + c3*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(4/3)...
    + c4*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(3)...
    + c5*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(37/6)...
    + c6*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(71/6)) ];

run file

R0 = 510e-6;
y0 = [510e-6; 0; 11326.8525; 0.01716];
M = [1 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0];
options = odeset('Mass',M);
[t,y] = ode15s(@(t,y) densita_20_latente_christ(t,y,R0),[0 52.424e-6],y0,options);
t=t*1000000;
y=y*1000000;
[t,y(:,1)];
plot(t,y(:,1))
$\endgroup$
0
$\begingroup$

Did you try scaling your problem? If $a_1$ is very large and $t$ is very small, why not scale $t$ (and all the other parameters) accordingly? If $t$ is now expressed in seconds, try to write down the problem in microseconds... $a_1$ and $a_2$ are both large so for $t$ it will work out. I don't know how it will go with the other parameters.

$\endgroup$
1
  • $\begingroup$ I scaled his ode such that all coefficients in front of r'^2/r, 1/r, t^2/r, t/r, r'/r, r'/r^2,1/r^2 are all more or less around 1, e.g. 1.5, 1., 1., -2.0053, 0.0049, 0.0020, 0.0074. This does not help. I think there must be something inherently wrong (e.g. a pole as the plot also indicates) when approaching this time. $\endgroup$
    – Bort
    Nov 24 '16 at 13:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.