1
$\begingroup$

I have a 4 pdes for 4 variables: v, v2, ua, and p. I would like to specify that v==v2 on the left and right boundaries at all times. Here are the equations:

enter image description here

I discretise the system using PDEtoODE. However, NDSolve fails to initialise. I thought it could be a problem with the bc being an algebraic equation so I applied diffbc, but to no avail. Aborting the computation crashes the kernel.

Would appreciate any guidance on this! Thanks! I can't quite find anything on this site for such boundary conditions...

For the same system, I also have another question: how could I couple a vertical grid to solve the last 2 equations, while the first 2 remain only functions of x? The first 2 equations describe the displacement of the top and bottom of the vertical boundaries. A scaling of the vertical coordinate by the distance between v and v2 at each x position? But at the boundaries the distance between them would go to zero...

Code: (pdetoode is here)

\[Sigma] = 917.4311927*0.000047;
sf = 1; l = sf 36/100;
lb = 0; rb = 0.36;
Clear[vinit]
hinit = 0.02;
vinit[x_] := -1 x (x - 0.36) + hinit;
v2init[x_] :=-((0.03*1)/2) (Tanh[30 x])*(Tanh[30 (-x + 0.36)]) + hinit;
glevel = v2init[0.18]


T = 298; M = 29*10^-3; R = 8.31; patm = 101000; \[Mu] = 10^-5;
g = 9.81; W0 = 0.23;

mi = 41*10^-3; zvinit = -3.6; ri = 0.021; droppos = 0.075;

impfunc = 1/2 1/2 (Tanh[100 (x - (droppos - ri))] + 1)*(Tanh[100 (-x + (droppos + ri))] + 1)(*1/21/2(Tanh[10000(x-(W0/2-ri))]+1)*(Tanh[10000(-x+(W0/2+ri))]+1)*);(*position of load*)

Nforc = (mi (-D[v[t, x], t, t] + -g))/(Pi ri^2)*impfunc;(*load on top*)

Ntable[a_] := Simplify`PWToUnitStep@Piecewise[{{0.2*10^6 (glevel - a), a < glevel}}, 0]; (*normal force on bottom*)

\[Lambda] = 0.001;

With[{ua = ua[t, x], p = p[t, x], v = v[t, x], v2 = v2[t, x]},
\[Rho] = (p + patm) M/(R T);
 aeqs = Simplify[{D[\[Rho] (v - v2), t] + D[\[Rho] (v - v2) ua, x] == 0, D[\[Rho] ua, t] + ua D[\[Rho] ua, x] == -D[p, x] + \[Mu] D[ua, x, x]}];
 abcs = {{D[ua, x] == 0, D[p, x] == 0} /. x -> lb, {D[ua, x] == 0, D[p, x] == 0} /. x -> rb};
 aics = {ua == 0, p == 0} /. t -> 0;
 
 \[Phi] = D[v, x] - 1/3 (D[v, x])^3; \[Phi]2 = D[v2, x] - 1/3 (D[v2, x])^3;
 eqs = {\[Sigma] D[v, t, t] == \[Lambda]   D[v, x, x]+ p Cos[\[Phi]] + Nforc,
\[Sigma] D[v2, t, t] ==\[Lambda] D[v2, x, x]- p Cos[\[Phi]2] + Ntable[v2]};
 ics = {v == vinit[x], D[v, t] == zvinit*impfunc, v2 == v2init[x], D[v2, t] == 0} /. t -> 0;
 bcs = {{D[v, x] == 0, D[v2, x] == 0} /.  x -> lb, {D[v, x] == 0, D[v2, x] == 0)} /. x -> rb};
equality = {v == v2 /. x -> lb, v == v2 /. x -> rb};]

xdomain = {0, 36/100};
xpoints = 36*4;
xgrid = Array[# &, xpoints, xdomain];
difforder = 4;
varua = Outer[ua, xgrid];
varp = Outer[p, xgrid];
varv = Outer[v, xgrid]; varv2 = Outer[v2, xgrid];

ptoo = pdetoode[{v[t, x], v2[t, x], ua[t, x], p[t, x]}, t, xgrid, difforder];

del = #[[2 ;; -2]] &; del3 = #[[3 ;; -3]] &;

ode@1 = ptoo@aeqs[[1]] // del; ode@2 = (ptoo@aeqs[[2]]) // del; 
ode@3 = ptoo@(eqs[[1]]) // del; ode@4 = ptoo@(eqs[[2]]) // del3;

odeic@1 = ptoo@aics; odeic@2 = ptoo@ics;
With[{sf = 1},
 odebc@1 = DeleteDuplicates[Flatten[abcs // ptoo]]; 
 odebc@2 = DeleteDuplicates[Flatten[ptoo@bcs]];
 odebc@3 = DeleteDuplicates[Flatten[diffbc[t, sf]@equality // ptoo]];]

tend = 0.004;
Monitor[sollst = NDSolveValue[{ode[1], ode[2], ode[3], ode[4], odeic[1], odeic[2], odebc[1], odebc[2], odebc[3]}, {varua, varp, varv, varv2}, {t, 0, tend}, Method -> {"EquationSimplification" -> "Residual", "IndexReduction" -> Automatic}, EvaluationMonitor :> (time = t)], time];

{uasol, psol, vsol, v2sol} = rebuild[#, xgrid, 1] & /@ sollst;
$\endgroup$
2
  • $\begingroup$ If we add boundary conditions in a form of odebc[3] then we need to remove 2 boundary conditions, for example from odebc[2] . What is you suggestion? $\endgroup$ Commented May 16 at 18:55
  • $\begingroup$ My apologies for the delay, there was a typo in my code, I applied del3 to one of the equations so that there are the correct number of equations. The code I was running on my end already included it but the error still stands. I've updated the question, turns out the solution is to apply EquationSimplification and IndexReduction but I don't understand why? $\endgroup$ Commented May 22 at 4:07

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.