I would like to simulate the motion of an inextensible membrane in 1D as an inextensible string subject to spatially-dependent body forces. For now, I would like it to be constrained at both ends. I use a Lagrange multiplier to enforce inextensibility, it is effectively tension in the string. We consider large displacements (ie both the horizontal and vertical displacements).
I first used NDSolveValue, but I was not quite sure what bcs to set for the multiplier so I set them to zero at the ends. However, inextensibility clearly isn't being enforced in my simulation (yellow being a later time and blue being the initial condition, here I set an initial velocity to a few points along the string). It is also not sensitive to the direction of loads applied.
Then, I borrowed code from here, removing the terms for bending moment, and including code for the spatially-dependent Lagrange multiplier. However, for a vertical stress larger than 2 N / m^2, the code slips into complex values within several timesteps and decreasing the timestep or increasing the points does not seem to resolve the issue. Furthermore, FindRoot takes far too long.
How should I solve my equations with reasonable computation time? I would like to extend this model much further with free ends and coupling to a pressure field beneath it. Thanks!
Code modified from @Alex Trounev :
ClearAll["Global`*"]
sf = 1(*0.98*);
(*Parameters*)l = sf 36/100; W0 = 23/100; H0 =
0.047*10^-3; n = 50; dx = l/(n - 1);
σ = 917.4311927*0.000047;(*mass per unit area*)
dt = 0.1 10^-2; tmax = 20 dt;
sgrid = Range[0, l, dx];tgrid = Range[0, tmax, dt]; jmax = nt = Length[tgrid];
nn = Length[sgrid];
ds = Table[NDSolve`FiniteDifferenceDerivative[i, sgrid, "DifferenceOrder" -> 4]["DifferentiationMatrix"], {i, 4}];
U = Array[uu, {nn, nt}]; V = Array[vv, {nn, nt}]; lambda = Array[ll, {nn, nt}];
Parray = Array[pe, {nt}];
fv = Table[If[0.05 < i dx && i dx < 0.08, -1, 0], {i, 1, nn}];
us = Sqrt[1 - (ds[[1]] . V)^2] - 1;(*inextensibility constraint*)
Do[uu[i, j] = If[i == 1, 0, dx/2 Sum[(us[[k - 1, j]] + us[[k, j]]), {k, 2, i}]];, {i, nn}, {j, nt}];(*u[t,s]=Integrate[us,{x,0,s}]*)
Do[ll[i, j] = If[i == 1, 0, (0 + dx/2 (Sum[σ (uu[k, j] + uu[k, j - 2] - 2 uu[k, j - 1])/dt^2, {k, 1, i - 1}] + Sum[σ (uu[k, j] + uu[k, j - 2] - 2 uu[k, j - 1])/dt^2, {k, 2, i}]))/(1 + us[[i, j]])];, {i, nn}, {j, 3, nt}]; (*Lagrange multiplier*)
(*ics*)
f[x_] :=(*1.1*)1.1 x (sf l - x);
Do[vv[i, 1] = f[(i - 1) dx]; vv[i, 2] = f[(i - 1) dx];, {i, nn + 1}]
Do[ll[i, j] = 0, {i, nn}, {j, 2}];
test = Table[.1, {nn}];
Monitor[Do[(*beam implicit step*)
r2 = ds[[1]] . lambda[[All, j + 1]] (ds[[1]] . V[[All, j + 1]]) +
lambda[[All, j + 1]] (ds[[2]] . V[[All, j + 1]]) + fv(*lambda[[
All,j+1]] (ds[[2]].V[[All,
j+1]])*);
eq2 = σ (-V[[All, j + 1]] - V[[All, j - 1]] + 2 V[[All, j]]) + dt^2 r2;
eqs = Table[Drop[Drop[eq2, 2], -2][[i]] == 0, {i, nn - 4}];
varAll = V[[All, j + 1]];
bC0 = {(ds[[2]] . V[[All, j + 1]])[[1]] == 0, V[[1, j + 1]] == 0, (ds[[2]] . V[[All, j + 1]])[[-1]] == 0, V[[-1, j + 1]] == 0};(*simply supported end*)
eqAll = Join[eqs, bC0];
ini = Table[{varAll[[i]], test[[i]]}, {i, Length[varAll]}];
sol =(*Quiet@*)FindRoot[eqAll, ini, Method -> {"Newton", "StepControl" -> "TrustRegion"}, MaxIterations -> 30];
test = varAll /. sol;
Do[vv[i, j + 1] = test[[i]], {i, nn}];, {j, 2, jmax - 1}],
j] // AbsoluteTiming
NDSolve code:
ClearAll["Global`*"]
sf = 1(*0.98*);
(*Parameters*)l = sf 36/100; W0 = 23/100; H0 = 0.047*10^-3;
σ = 917.4311927*0.000047;(*mass per unit area*)
lb = 0; rb = l;
Clear[fv, vinit, uinit, λinit]
fv[s_] := Piecewise[{{-1, 0.05 < s && s < 0.1}}, 0]
vinit[x_] :=(*1.1*)1.1 x (sf l - x);
uinit[x_] := 0;
λinit[x_] := 0;
With[{u = u[s, t], v = v[s, t], λ = λ[s, t]},
eqs = {σ D[u, t, t] - D[λ (1 + D[u, s]), s] == 0,
σ D[v, t, t] - D[λ (D[v, s]), s] == 0(*10000fv[s]*)(*fv[s]*),
(1 + D[u, s])^2 + D[v, s]^2 == 1};
bcs = Flatten[{{u == 0, v == 0, λ == 0(*,D[u,s,s]==0,D[v,s,s]==0*)} /. s -> lb, {u == 0, v == 0, λ == 0(*,D[u,s,s]==0,D[v,s,s]==0*)} /. s -> rb}];
ics = {u == uinit[s], D[u, t] == 0, v == vinit[s], D[v, t] == fv[s], λ == λinit[s]} /. t -> 0;
];
tend = 0.1;
Monitor[{usol, vsol, lsol} = NDSolveValue[{eqs, ics, bcs}, {u, v, λ}, {t, 0, tend}, {s,lb, rb}, Method -> {"IndexReduction" -> Automatic}, EvaluationMonitor :> (time = t)], time]
