1
$\begingroup$

I am solving the following coupled system of 3 PDEs modelling a 1D membrane coupled to a 1D fluid flow field underneath.

enter image description here

However, on putting them into NDSolveValue (and trying with FEM), it says the "Matrix is singular". So, I tried the hack suggested by @user21 here which enabled NDSolve to start solving...but ua and p aren't solved at all (remain 0 at all times everywhere) and the solution of v is just that without any coupling. Trying other methods (eg specifying TensorProductGrid) returns back the issue of the matrix being singular. But why is it singular and how do I fix this?

Thanks!

Code:

ClearAll["Global`*"]
Needs["NDSolve`FEM`"]
\[Sigma] = 917.4311927*0.000047;
sf = 1; l = sf 36/100;
lb = 0; rb = 0.36;
vinit[x_] := -0.3 (x + 0.36) (x - 0.36);
eqimp = -0.04 Sin[2 Pi/0.1 t] + vinit[lb];

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

nn = 200; dx = rb/(nn - 1);
grid = Table[dx (i - 1), {i, 1, nn}];
mesh = ToElementMesh[Map[{#} &, grid]]

With[{ua = ua[t, x], p = p[t, x], v = v[t, x]},
\[Rho] = (p + patm) M/(R T);
 aeqs = Simplify[{D[\[Rho] v ua, t] == 0, D[\[Rho] ua, t] + D[\[Rho] ua^2, x] == -D[p, x] + \[Mu] D[ua, x, x](*+NeumannValue[0,x==lb||x==rb]*)}];
 abcs = {(*DirichletCondition[ua==0,x==rb]*){D[ua, x] == 0} /. x -> lb, {D[ua, x] == 0(*ua==0*)} /. x -> rb};
 aics = {ua == 0, p == 0} /. t -> 0;
 
 k = 100;

 \[Phi] = D[v, x] - 1/3 (D[v, x])^3;
 eqs = {\[Sigma] D[v, t, t] == D[v, x, x]+ p Cos[\[Phi]]};
 ics = {v == vinit[x], D[v, t] == 0} /. t -> 0;
 bcs = {v == eqimp /. x -> lb, v == 0 /. x -> rb (*DirichletCondition[v==eqimp,x==lb],DirichletCondition[v==0,x==rb]*)};]

Monitor[{uasol, psol, vsol} = NDSolveValue[Flatten[{aeqs, abcs, aics, eqs, ics, bcs}], {ua, p, v}, {t, 0, tend},x\[Element]mesh, EvaluationMonitor :> (time = t)], time]

$\endgroup$
4
  • 1
    $\begingroup$ Initial condition ua == 0 makes useless equation D[\[Rho] v ua, t] == 0. $\endgroup$ Commented May 14 at 15:34
  • $\begingroup$ Oh my you're right I see the problem with that equation, I've corrected it to the right one, but now I'm facing issues with stabilising the system to solve it at least to 0.04s. By 0.002s the boundaries grow large oscillations...I'll open a new question then? In general when systems of equations are unstable what are the strategies you use to attack the problem? Thanks! $\endgroup$ Commented May 15 at 13:26
  • $\begingroup$ It could be better to use some hand made algorithm like predictor-corrector, see for example this topic mathematica.stackexchange.com/questions/287067/… $\endgroup$ Commented May 15 at 16:12
  • $\begingroup$ Thanks for the suggestion, I'm giving it a try now. Meanwhile I encountered another issue with enforcing the bc such that two variables in my system of equations are equal at all times (NDSolve fails to initialise), I've posted here: mathematica.stackexchange.com/q/303197/88731 Would appreciate your help! Can't quite find a similar question on this site... $\endgroup$ Commented May 16 at 9:42

0