n = 6;4;
variables = Table[ToExpression["φ" <> ToString[i]], {i, 1, n}];
xComponents = Table[Cos[variables[[i]][t]], {i, 1, n}];
yComponents = Table[-Sin[variables[[i]][t]], {i, 1, n}];
endsOfRods = Table[{xComponents[[i]], yComponents[[i]]}, {i, 1, n - 1}];
xTotal = Table[
l (Sum[endsOfRods[[i, 1]], {i, 1, j - 1, 1}] + 1/2 endsOfRods[[j, 1]]), {j,
2, n - 1, 1}];
yTotal = Table[
l (Sum[If[i <= n/2, endsOfRods[[i, 2]], (-1)*endsOfRods[[i, 2]]], {i, 1,
j - 1, 1}] +
If[j > n/2, (-1)*1/2 endsOfRods[[j, 2]], 1/2 endsOfRods[[j, 2]]]), {j,
2, n - 1, 1}];
centreOfMass = Table[{xTotal[[i]], yTotal[[i]]}, {i, 1, Length[xTotal], 1}];
velocityCentreOfMass =
Table[D[centreOfMass[[i]], t].D[centreOfMass[[i]], t], {i, 1,
Length[xTotal], 1}];
kineticEnergy =
Sum[If[i == 1 || i == n, 1/2*(1/3 m l^2)*D[variables[[i]][t], t]^2,
1/2 m velocityCentreOfMass[[i - 1]] +
1/2*(1/12 m l^2)*D[variables[[i]][t], t]^2], {i, 1, n, 1}];
potentialEnergy =
Sum[If[i == 1 || i == n, -m g l 1/2 Sin[variables[[i]][t]] ,
m g yTotal[[i - 1]]], {i, 1, n, 1}](*gravity*)+
Sum[If[i < n/2, 1/2 k (variables[[i]][t] - variables[[i + 1]][t])^2,
If[i == n/2, 1/2 k (variables[[i]][t] + variables[[i + 1]][t])^2,
1/2 k (variables[[i + 1]][t] - variables[[i]][t])^2]], {i, 1, n - 1,
1}](*rotational springs*);
lagrange = kineticEnergy - potentialEnergy;
Nonconservative part
F = {0, d D[(-1)*(yTotal[[n/2 - 1]] - l/2 Sin[variables[[n/2]][t]]), t]};
r = {(xTotal[[n/2 - 1]] + l/2 Cos[variables[[n/2]][t]]), (yTotal[[n/2 - 1]] -
l/2 Sin[variables[[n/2]][t]])};
generalizedForces = Table[F.D[r, variables[[i]][t]], {i, 1, n/2, 1}];
Numerical solution
equations =
Table[If[i <= n/2,
D[D[lagrange, D[variables[[i]][t], t]], t] -
D[lagrange, variables[[i]][t]] - generalizedForces[[i]] == 0,
D[D[lagrange, D[variables[[i]][t], t]], t] -
D[lagrange, variables[[i]][t]] == 0], {i, 1, n - 1, 1}];
AppendTo[equations, (-1)*(yTotal[[n/2]] + l/2 Sin[variables[[n/2 + 1]][t]] +
l Sin[variables[[n]][t]]) == 0];
m = 250(*kg*);
l = 0.75(*m*);
d = 4000(*Ns/m*);
EE = 2.1*10^(11)(*Pa*);
II = 1.5*10^(-6)(*m^4*);
g = 9.81(*m/s^2*);
k = EE II/l;
AppendTo[equations, Table[variables[[i]][t] == 0 /. t -> 0, {i, 1, n, 1}]];
AppendTo[equations,
Table[D[variables[[i]][t], t] == 0 /. t -> 0, {i, 1, n, 1}]];
solution =
NDSolve[equations, Table[variables[[i]][t], {i, 1, n, 1}], {t, 0, 5}];
