NDSolve fails to solve partial differential equation

Question

When solving a partial differential equation, I get

NDSolve::eerri: Warning: estimated initial error on the specified spatial grid in the direction of independent variable x exceeds prescribed error tolerance. And the equation cannot be solved.

My code:

ClearAll
l = 2.5
l0 = 1.5
Eb = 210*1000000000
rb = 10/1000
rhob = 7850
Eg = 35*1000000000
rg = 17.5/1000
vg = 0.25
varphi = 45/180*Pi
Cb = 2*Pi*rb
Ab = Pi*rb^2
R = 10*Eb*rb/((Eg + xE)/2)
Gg = Eg/(2*(1 + vg))
k = Gg*G/((G*Log[rg/rb] + Gg*Log[R/rb])*rb)
YiTa = Tan[varphi]
z1 = b
p[x_, t_] = FullSimplify[YiTa*Sigmaz[r0 + l - x, z1, t]]
PDE1 = D[u1[x, t], {x, 2}] == 
  rhob*D[u1[x, t], {t, 2}]/Eb + k*Cb*u1[x, t]/(Eb*Ab) + 
   p[x, t]*Cb/(Eb*Ab)
initialCondition = {Derivative[0, 1][u1][x, 0] == 
   Vr[r0 + l - x, z1, 0], u1[x, 0] == Ur[r0 + l - x, z1, 0]}
boundaryConditions = {u1[0, t] == Ur[l + r0, z1, t], 
  u1[l0, t] == Ur[r0 + l - l0, z1, t]}
sol = NDSolve[{PDE1, initialCondition, boundaryConditions}, 
  u1, {x, 0, l0}, {t, 0, 0.001}, MaxStepSize -> 0.0000001]
Plot3D[Evaluate[u1[x, t] /. sol], {x, 0, l0}, {t, 0, 0.0001}, 
 PlotRange -> All, AxesLabel -> {"x", "t", "u1(x,t)"}]

result

I tried running the code but got no results. I hope to get results that can be run and plotted.

Answer 1

Too long to put in a comment. As I earlier suggested, I replaced your $MachinePrecision constants with exact rationals.

mu=25/100;rho=2350;xE=18*10^9;G=xE/(2*(1+mu));cp=Sqrt[xE*(1-mu)/(rho*(1+mu)*(1-2*mu))];
r0=25/10;xdb=40/1000;xdp=25/10*40/1000;rho0=1100;xD=3200;sb=4500/1000;r0=25/10;
b=15/10;alpha=4000;
JM=(3*Pi)/(4*r0);m=alpha/cp;j=Sqrt[JM^2+m^2];A0=xdb/(8*sb)*rho0*xD^2/(xdp/xdb)^(22/10);
A=(1-E^(-b*j))*xE((j^2-m^2+(2*j^2-3*m^2)*mu)/Sqrt[r0]*Cos[Sqrt[j^2-m^2]*r0-
  Pi/4]/(1-2*mu)-Sqrt[j^2-m^2]*r0^(-3/2)*Sin[Sqrt[j^2-m^2]*r0-Pi/4]-(3-5*mu)*
  r0^(-5/2)*Cos[Sqrt[j^2-m^2]*r0-Pi/4]/(4*(1-2*mu)))/(b*j*(1+mu));A1=A0/A;
Ur[r_,z_,t_]=A1*(1/2*Cos[Sqrt[j^2-m^2]*r-1/4*Pi]/r^(3/2)+Sqrt[j^2-m^2]*
  Sin[Sqrt[j^2-m^2]*r-1/4*Pi]/Sqrt[r])*E^(-j*z)*E^(-m*cp*t);
Vr[r_,z_,t_]:=-m*cp*A1*(1/2*Cos[Sqrt[j^2-m^2]*r-1/4*Pi]/r^(3/2)+Sqrt[j^2-
  m^2]*Sin[Sqrt[j^2-m^2]*r-1/4*Pi]/Sqrt[r])*E^(-j*z)*E^(-m*cp*t);
Sigmaz[r_,z_,t_]:=xE*A1*((-j^2+(2*j^2-m^2)*mu)*Cos[Sqrt[j^2-m^2]*r-1/4*Pi]/
  ((1-2*mu)*Sqrt[r])-mu*Cos[Sqrt[j^2-m^2]*r-1/4*Pi]/((4-8*mu)*r^(5/2)))*
  E^(-j*z)*E^(-m*cp*t)/(1+mu);
l=25/10;l0=15/10;Eb=210*10^9;rb=10/1000;rhob=7850;Eg=35*10^9;rg=175/10/1000;
vg=025/100;varphi=45/180*Pi;Cb=2*Pi*rb;Ab=Pi*rb^2;R=10*Eb*rb/((Eg+xE)/2);
Gg=Eg/(2*(1+vg));k=Gg*G/((G*Log[rg/rb]+Gg*Log[R/rb])*rb);YiTa=Tan[varphi];z1=b;
p[x_,t_]:=YiTa*Sigmaz[r0+l-x,z1,t];
PDE1={D[u1[x,t],{x,2}]==rhob*D[u1[x,t],t,2}]/Eb+k*Cb*u1[x,t]/(Eb*Ab)+
  p[x,t]*Cb/(Eb*Ab)};
initialCondition={Derivative[0,1][u1][x,0]==Vr[r0+l-x,z1,0],u1[x,0]==Ur[r0+l-x,z1,0]};
boundaryConditions={u1[0,t]==Ur[l+r0,z1,t],u1[l0,t]==Ur[r0+l-l0,z1,t]};
sys=Simplify[N[Join[PDE1,initialCondition,boundaryConditions],64]]

Just before the step of doing the NDSolve that displays (roughly)

{Derivative[2, 0][u1][x, t] == (     
0.00545860389386734035988805739383977751117746615039472979478311832031082160733*
(25.035720659226906140411453104565220866207337272409961236398875988369353985859-
10.*x + 1.*x^2)*Cos[
3.9269908169872415480783042290993786052464617492188822762186807403847705078577- 
.94247779607693797153879301498385086525915081981253174629248337769234492*x])/    
(E^(3999.99999999999999999999999999999999999999999999999999999999999999999*t)* 
(5. - 1.*x)^(5/2)) +    
147.137927874474292837143887032039529266866273227161600896429319953505602294*
u1[x, t] + 3.738095238095238095238095238095238095238095238095238095238095238*^
-8*Derivative[0, 2][u1][x, t], 
Derivative[0, 1][u1][x, 0] == -(((       
(1.6768709894289540203312443887507074674202073320420657991243447623718564+ 0.*I)- 
0.30320332583960679232662864081734027262650464178210400129298177083547+0.*I)*x)*
Cos[0.9424777960769379715387930149838508652591508198125317462924833776923*x] + 
((-1.3551622689671139029350420194226952588448390857789742138054729459828 +0.*I) + 
(0.303203325839606792326628640817340272626504641782104001292981770835+0.*I)*x)*
Sin[0.94247779607693797153879301498385086525915081981253174629248337769*x])/
(5. - 1.*x)^(3/2)), 
u1[x, 0] == -((( 
(-0.00041921774735723850508281109718767686685505183301051644978108619059+0.*I) +
(0.00007580083145990169808165716020433506815662616044552600032324544 +0.*I)*x)*
Cos[0.9424777960769379715387930149838508652591508198125317462924833776923*x] + 
((0.000338790567241778475733760504855673814711209771444743553451368236495+0.*I)- 
(0.0000758008314599016980816571602043350681566261604455260003232454427+0.*I)*x)*
Sin[0.942477796076937971538793014983850865259150819812531746292483377688587*x])/
(5. - 1.*x)^(3/2)), 
u1[0, t] == 
0.000037495975218604724441921501083994167027199906049745772524411502557/
E^(3999.99999999999999999999999999999999999999999999999999999999999999999993*t),
u1[1.5, t] == 
-0.00002665354928541077808491725593466020104224160610364654413085048230532/
E^(3999.9999999999999999999999999999999999999999999999999999999999999999993*t)}

Most of that looks reasonable, except for the E^(3999.99999...*t) terms and all "approximately 0*I" terms.

Can you work through the changes I made, see if I made any mistake doing this, and get rid of those terms that scare me?

I have not removed those, but I have tried the handful of tricks that I've used in the past and it either takes so much time or memory or bails out with the warning about the error being too large.

Is it possible to come up with a substantially simpler system that still correctly represents what you are trying to give NDSolve? If you can do that then that might enable it to solve for you.