I would like to solve the 3 coupled PDEs describing a damped, nonlinear (i.e displacements in the $x$ direction along the beam need to be considered along with the $y$ displacements normally considered), inextensible, euler-bernoulli cantilever beam's motion. The equations are as follows:
1: where $u(t,s)$ is the displacement in the $x$ direction $$\mu\frac{\partial^2 u}{\partial t^2}+k\frac{\partial u}{\partial t}-\frac{\partial}{\partial s}\left(\lambda\left(1+\frac{\partial u}{\partial s}\right)\right)=0$$ 2: where $v(t,s)$ is the displacement in the $y$ direction
EDIT: To make things simpler, we expand the elasticity term to first order only. $$\mu\frac{\partial^2 v}{\partial t^2}+k\frac{\partial v}{\partial t}-\frac{\partial}{\partial s}\left(\lambda\frac{\partial v}{\partial s}\right)+EI\frac{\partial^4 v}{\partial s^4}=0$$ 3: inextensibility of beam condition $$\left(1+\frac{\partial u}{\partial s}\right)^2+\left(\frac{\partial v}{\partial s}\right)^2-1=0$$ $\mu$, $EI$ and $k$ are constants representing mass density, stiffness and drag coefficient, while EDIT: $\lambda(t,s)$ is the lagrange multiplier needed to maintain the inextensibility condition.
My boundary conditions are as follows (for a beam of unit length):
EDIT: New b.c.s. that include $\lambda(t,s)$
2 for $u$
$u(t,0)=0$
$\lambda(t,1)\left(1+u^{(0,1)}(t,1)\right)=0$
4 for $v$
$v(t,0)=0$
$v^{(0,1)}(t,0)=0$
$v^{(0,2)}(t,1)=0$
$EI\left(v^{(0,3)}(t,1)\right)-\lambda(t,1)v^{(0,1)}(t,1)=1$
EDIT: I use a point load of 1N at the end of the beam to perturb it instead so the beam is completely stationary in the i.c.s
Next, my initial conditions are of a fully stationary, straight beam:
$u(0,x)=0$
$u^{(1,0)}(0,x)=0$
$v(0,x)=0$
$v^{(1,0)}(0,x)=0$
My approach is to use the pdetoode function by @xzczd. I had no problems following and implementing the solution for the linearized Euler-Bernoulli equation as in this earlier question (that is, only small deflections in the $y$ direction are considered and inextensibility is ignored), but I have issues with generating the correct number of ODEs for NDSolve to solve for the nonlinear system.
Here is my implementation within Mathematica:
EDIT: Here is the code with the updated functions and b.c.s.
μ = 1;
EI = 10;
k = 1;
eqn1 = μ*D[u[t, s], {t, 2}] + k*D[u[t, s], t] -
D[λ[t, s]*(1 + D[u[t, s], s]), s] == 0;
eqn2 = μ*D[v[t, s], {t, 2}] + k*D[v[t, s], t] -
D[λ[t, s]*D[v[t, s], s], s] + EI*D[v[t, s], {s, 4}] == 0;
eqn3 = (1 + D[u[t, s], s])^2 + D[v[t, s], s]^2 - 1 == 0;
bc1 = {u[t, 0] == 0, λ[t, 1]*(1 + Derivative[0, 1][u][t, 1]) ==
0};
bc2 = {v[t, 0] == 0, Derivative[0, 1][v][t, 0] == 0,
Derivative[0, 2][v][t, 1] == 0,
EI*Derivative[0, 3][v][t, 1] - λ[t, 1]*
Derivative[0, 1][v][t, 1] == 1};
ic1 = {u[0, s] == 0, Derivative[1, 0][u][0, s] == 0};
ic2 = {v[0, s] == 0, Derivative[1, 0][v][0, s] == 0};
Generate the finite-difference grid (i use 9 points for now to speed things up but more would be needed for accuracy):
lb = 0;
rb = 1;
torder = 2;
sdifforder = 2;
points = 9;
grid = Array[# &, points, {lb, rb}];
Generate ODEs using pdetoode:
removeredundant1 = #[[2 ;; -2]] &;
removeredundant2 = #[[3 ;; -3]] &;
ptoofunc = pdetoode[{u, v, λ}[t, s], t, grid, sdifforder];
odeqn1 = eqn1 // ptoofunc // removeredundant1;
odeqn2 = eqn2 // ptoofunc // removeredundant2;
odeqn3 = eqn3 // ptoofunc;
odeic1 = Flatten[removeredundant1 /@ ptoofunc@ic1];
odeic2 = Flatten[removeredundant2 /@ ptoofunc@ic2];
odebc1 = bc1 // ptoofunc;
odebc2 = bc2 // ptoofunc;
EDIT: For the ODEs generated from eqn 1 and i.c. 1, I remove 2 equations to make space for 2 b.c.s. For the ODEs generated from eqn 2 and i.c. 2, I remove 4 equations to make space for 4 b.c.s. I remove no equations from the ODEs generated from eqn 3. The system is now correctly determined and I no longer get the error NDSolve::overdet:
tEnd = 1;
sollst = NDSolveValue[
Join[odebc1, odebc2, odeic1, odeic2, odeqn1, odeqn2, odeqn3],
Join[u /@ grid, v /@ grid, λ /@ grid], {t, 0, tEnd},
MaxSteps -> Infinity];
Now, the issue is as follows:
NDSolveValue::icfail: Unable to find initial conditions that satisfy the residual function within specified tolerances. Try giving initial conditions for both values and derivatives of the functions.
WhenEvent(these methods are already working for me for the lineaerized equtions). $\endgroup$