I am solving the following coupled system of 3 PDEs modelling a 1D membrane coupled to a 1D fluid flow field underneath.
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]
ua == 0
makes useless equationD[\[Rho] v ua, t] == 0
. $\endgroup$