home

Scientific Computing (Psychology 9040a)

Fall, 2021

Hanging Chain: Dynamical Simulation

In a previous topic we saw how to optimize the potential energy of a hanging chain, and end up with a catenary curve. Here we will extend this idea to produce a dynamical simulation of the motion of a hanging chain. The idea will be to define the initial locations (and velocities) of the chain nodes, the masses of the nodes, the stiffness (and damping) of the wires connecting the nodes, and then let the chain fall with gravity, and describe the dynamical behaviour as it falls.

Here is the main script that runs a simulation, starting with a horizontal chain, and letting it fall:

t = 0:0.01:10;
n = 7;
x0 = 0; y0 = 0;
xn = 5; yn = 0;
m = 0.5/n; k = 50; b = 0.06;
x_init = linspace(x0,xn,n+2)'; x_init = x_init(2:end-1);
y_init = zeros(n,1);
Lr = (xn-x0)/(n+1);
state0 = [x_init; y_init; zeros(n,1); zeros(n,1)]; % x;y;xd;yd
myfun = @(t,state) spring2Dn(t,state,n,m,k,b,x0,y0,xn,yn,Lr);
[t,state] = ode45(myfun, t, state0);
figure
AnimateChain(t,state,x0,y0,xn,yn);

The variable n defines the number of inner nodes of the chain (not including the anchor nodes on the left and right). We define a time vector from 0 to 10 seconds, in steps of 10 ms. The m vector defines the masses of each node. The k value defines the stiffness of the wires between each node (assuming they are all equal) and the b value defines the damping. Without damping the chain will oscillate forever (try it by setting b=0).

The spring2Dn() function defines the dynamical behaviour of the hanging chain:

function stated = spring2Dn(t,state,n,m,k,b,x0,y0,xn,yn,Lr)

% (0,0)---(x1,y1)-...-(xm,ym)---(1,0)
% +x to the right, -y down
%
% state is [x; y; xd; yd]

% unpack state vector into x, y, xd, yd
x = state(1:n);
y = state(n+1:(2*n));
xd = state((2*n)+1:(3*n));
yd = state((3*n)+1:(4*n));

% add the anchors (x0,y0) and (xn,yn)
xx = [x0; x(:); xn];
yy = [y0; y(:); yn];

% compute lengths of each chain segment
L = zeros(n+1,1);
for i=1:n+1
   L(i) = sqrt((xx(i+1)-xx(i)).^2 + (yy(i+1)-yy(i)).^2); 
end
% compute stretch of each segment relative to resting length
S = (L-Lr).*((L-Lr)>0); % elastic bands; only pull, no push

g = 9.81; % the gravitational constant

% compute accelerations xdd and ydd of each state
xdd = zeros(n,1);
ydd = zeros(n,1);
for i=1:n
    cos1 = (xx(i+1)-xx(i))/L(i);
    cos2 = (xx(i+2)-xx(i+1))/L(i+1);
    sin1 = (yy(i+1)-yy(i))/L(i);
    sin2 = (yy(i+2)-yy(i+1))/L(i+1);
    xdd(i) = (-k*S(i)*cos1 + k*S(i+1)*cos2 - b*xd(i)      ) / m;
    ydd(i) = (-k*S(i)*sin1 + k*S(i+1)*sin2 - b*yd(i) - m*g) / m;
end

stated = [xd; yd; xdd; ydd];

end

I won’t explain the details of how these equations of motion were derived, it’s not so important for our purposes, which is to play with the various parameters of the simulation, and the system, and see how the behaviour changes. The basic idea is to characterize the sum of forces acting on each node, and then take advantage of Newton’s second law of motion, \(F=ma\), to arrive at an expression that gives the acceleration in \(x\) and \(y\) for each node.

Here is the MATLAB function that will animate the simulation of the hanging chain:

function out = AnimateChain(t,state,x0,y0,xn,yn)

n = size(state,2)/4;
p0=plot(x0,y0,'rs','markersize',10,'MarkerFaceColor',[1 0 0]);
hold on
pn=plot(xn,yn,'rs','markersize',10,'MarkerFaceColor',[1 0 0]);
for i=1:n
    xi = i;
    yi = n+i;
    p(i)=plot(state(1,xi),state(1,yi),'bs','markersize',10,'MarkerFaceColor',[0 0 1]);
end
l(1) = line([x0,state(1,1)],[y0,state(1,n+1)],'color','b');
for i=2:n
    xi = i;
    yi = n+i;
    l(i) = line([state(1,xi),state(1,xi+2)],[state(1,yi),state(1,yi+2)],'color','b'); 
end
l(n+1) = line([xn,state(1,xi)],[yn,state(1,yi)],'color','b'); 
tl = title(sprintf('%.3', t(1)));
min_x = min([min(state(:,1:n)),x0,xn])-0.05;
max_x = max([max(state(:,1:n)),x0,xn])+0.05;
min_y = min([min(state(:,n+1:n*2)),y0,yn])-0.05;
max_y = max([max(state(:,n+1:n*2)),y0,yn])+0.05;
axis([min_x max_x min_y max_y])
axis square
drawnow
pause
for i=1:size(state,1)
    for j=1:n
        xj = j;
        yj = n+j;
        p(j).XData = state(i,xj); p(j).YData = state(i,yj);
    end
    l(1).XData = [x0,state(i,1)]; l(1).YData = [y0,state(i,n+1)];
    for j=2:n
        xj = j;
        yj = n+j;
        l(j).XData = [state(i,xj),state(i,xj-1)];
        l(j).YData = [state(i,yj),state(i,yj-1)];
    end
    l(n+1).XData = [xn,state(i,xj)];
    l(n+1).YData = [yn,state(i,yj)];
    tl.String = sprintf('%.3f', t(i));
    drawnow
end
end

Things to try