Created Samstag 07 Mai 2016
A dformable object is simulated as mass particles that are interconnected with spring-damper systems.
Three particles are connected in series:
- The first paricle is fixed and can not be moved.
- The spring that connects the first and the second particle is at its rest/natural length.
- The spring that connects the second and the third particle is initially streched.
We introduce a connection matrix A.
1 → connection exists
0 → connection does not exist
- particle-1 is connected to particle-2
- particle-2 is connected to particle-1 and particle-3
- particle-3 is connected to particle-2
Nex, we introduce stiffnes matrix (KS) and damping matrix (BS)
KS = A*K
BS = A*B
where K and B are stiffness and damping coefficient, respectively.
function simple_msd close all clear all clc M = 10; r = 0.1; % Natural length of the spring x = [0 0; 0 0.1; 0 0.3]; % position, particle-3 is streched by 0.1 unit length! v = [0 0; 0 0 ; 0 0 ]; % velocity a = [0 0; 0 0 ; 0 0 ]; % acceleration N = length(x); % Stiffness Matrix ks = 10; KS = [0 ks 0; ks 0 ks; 0 ks 0]; % Damping Matrix bs = 10; BS = [0 bs 0; bs 0 ks; 0 bs 0]; % Time vector ts = 1e-3; tsim = 50; timespan = 0:ts:tsim; xk = x; % Perform the simulation for k = 1:length(timespan) FS = [0 0; 0 0; 0 0]; FD = [0 0; 0 0; 0 0]; for i = 1:N for j = 1:N FS(i,:) = FS(i,:) + fs(x(i,:), x(j,:), r, KS(i,j)); FD(i,:) = FD(i,:) + fd(v(i,:), v(j,:), BS(i,j)); end end % Particle-1 is fixed, don't modify its states! a(2:end, :) = (FS(2:end, :)+FD(2:end,:))./M; v(2:end, :) = v(2:end,:) + a(2:end,:) .* 0.001; x(2:end, :) = x(2:end,:) + v(2:end,:) .* 0.001; xk = [xk;x]; end figure hold on % plot positions of all particles for i = 1:N for j = 1:length(timespan) xk_(j, :) = xk((j-1)*N+i, :); end plot(timespan,xk_(:,2)); end legend('Particle-1', 'Particle-2', 'Particle-3') end function output = fs(xi,xj,r,k) if k > 0 % avoid division by zero xij = xi-xj; xij_hat = xij/norm(xij); output = -k*(norm(xij)-r)*xij_hat; else output = [0 0]; end end function output = fd(vi,vj,b) if b > 0 % avoid division by zero vij = vi-vj; output = -b*vij; else output = 0; end end
Download the code: ./simple_msd.m
- http://hugi.scene.org/online/hugi28/hugi 28 - coding corner uttumuttu implementing the implicit euler method for mass-spring systems.htm