====== SPACAR Light ======
===== Example: Simulating contact =====
This example shows a rudimentary script for the simulation of contact. It can definitely be enhanced and we are open to specific improvements.
Here, we're using a beam that is clamped on one side. The beam is largely rigid, but near the clamping, it is flexible. The flexible part acts like a hinge. We're using an actuation force of 1 N in ''x''-direction on node 3. Without any contact, the tip of the beam would displace 0.0311 m.
==== The concept ====
We're imagining a line at ''x=0.025'', representing a contact surface that cannot be crossed. To simulate this, we're creating a custom force on the node of contact (node number 4), //whenever it crosses// the position of ''x=0.025''. We choose the value of this force to be ''-stiffn*(pos-x0)''. This creates a fictitious force that acts like a spring with stiffness ''stiffn'' and stress-free reference position of ''x0=0.025''. Variable ''pos'' is the current value of the suspected node of contact.
==== In SPACAR ====
This simulation is tricky because it involves iterations to get the system in the right state. We currently have no clear rationale for picking the values of ''stiffn'' of the fictitious force (or the precise mathematical form of this force, for that matter). Also, it's best to do these iterations within one single Spacar call, seeing as each new SPACAR Light call would make the system start from the undeformed configuration again. The current interface for doing a custom force within Spacar is by means of the usersig interface. This interface hasn't been implemented in SPACAR Light properly, but here are some pointers to get a simulation up and running nevertheless. It's currently not elegant and for more details, you're encouraged to read up on the dat-file and the commands ''iterstep'', ''xf'', ''delxf'', as these all come into play under the hood. Also see the [[manuals|manuals]].
* Get the development version of SPACAR Light that supports the ''opt.customdyn'' field. See the [[https://github.com/precision-utwente/spalight/tree/contact|development branch on GitHub]].
* Use ''opt.customdyn = 'usersig PUSHSIG';'' to have Spacar call the m-file PUSHSIG.m (all of this is case-sensitive). This file should be in the same folder as the main SPACAR Light script.
* This m-file will be called repeatedly, depending on the number of ''opt.loadsteps''. This file contains a function with very specific input and output requirements. If you deviate from the template, you'll very quickly run into errors. Within this function, you'll have available the ''x''-vector, which contains the current positions of all nodes. This information can be used to create custom force profiles. Unfortunately, the node numbering in Spacar (full) is different from SPACAR Light. This is actually by design, and something you usually shouldn't have to deal with. Here, though, you do. The conversion is simple and detailed in the file. Note that this simple conversion may not work whenever you use multiple elements per beam with ''eprops(i).nbeams'' > 1. (The .dat-file shows what the actual conversion should be in all cases.)
function [time,sig,f] = PUSHSIG(t,ne,le,e,ep,nx,lnp,x,xp)
x0 = 0.025; %defines contact line
stiffn = 300;
spalightnode = 4;
%conversion between SPACAR Light node and Spacar (full) node.
%(for each SPACAR Light node, there are 2 nodes in Spacar full.
%They are numbered consecutively. First one is the translation node,
%second one is the rotation node. Here, we're taking just the translation node.)
spanode = (spalightnode-1)*2+1;
%get current position
pos = x(lnp(spanode,1:3)); %contains x,y,z coordinate of spalightnode
%define a custom force
node = spanode; %node number (according to Spacar full numbering)
dir = 1; %direction of force: 1 for x, 2 for y, 3 for z
debug_info = 1; %print out some debug info during simulation
if debug_info
fprintf('Current position: %f\n',pos(1))
end
if pos(1) > x0 %if within (forbidden) contact region
val = -stiffn*(pos(1)-x0); %define force to push it back
deriv = stiffn;
f = [node dir val deriv 0]; %apply force
if debug_info
fprintf('Within contact region. Applying force: %f\n',val)
end
else
f = []; %if not within contact region, don't apply contact force
if debug_info
fprintf('Not within contact region.\n')
end
end
sig = []; %can be left empty (is used for specifying generalized stresses)
time = t; %this line needs to be there; don't modify.
end
The script with the SPACAR Light model looks like this.
% EXAMPLE SCRIPT FOR RUNNING SPACAR LIGHT
% This example simulates a simple cross flexure rotating due to an applied moment
clear
clc
close all
% addpath('spacar') %Specify location of spacar folder, if not current
%% NODE POSITIONS
% x y z
nodes = [ 0 0 0; %node 1
0 -0.02 0; %node 2
0 -0.04 0; %node 3
0 -0.06 0]; %node 4
%% ELEMENT CONNECTIVITY
% p q
elements = [ 1 2; %element 1
2 3; %element 2
3 4]; %element 3
%% NODE PROPERTIES
%node 1
nprops(1).fix = true; %Fix node 1
nprops(3).force = [1 0 0];
%% ELEMENT PROPERTIES
eprops(1).elems = [1]; %Add this set of properties to elements 1 and 3
eprops(1).emod = 2e9; %E-modulus [Pa]
eprops(1).smod = 68.9e6; %G-modulus [Pa]
eprops(1).dens = 1150; %Density [kg/m^3]
eprops(1).cshape = 'rect'; %Rectangular cross-section
eprops(1).dim = [5e-3 1e-3]; %Width: 50 mm, thickness: 0.2 mm
eprops(1).orien = [0 0 1]; %Orientation of the cross-section as a vector pointing along "width-direction"
eprops(1).nbeams = 1; %Number of beams used to model this element
eprops(1).flex = 1:6; %Model out-of-plane bending (modes 3 and 4) as flexible
eprops(1).color = 'grey'; %Color
eprops(2).elems = [2 3]; %Add this set of properties to elements 1 and 3
eprops(2).emod = 2e9; %E-modulus [Pa]
eprops(2).smod = 68.9e6; %G-modulus [Pa]
eprops(2).dens = 1150; %Density [kg/m^3]
eprops(2).cshape = 'rect'; %Rectangular cross-section
eprops(2).dim = [5e-3 5e-3]; %Width: 50 mm, thickness: 0.2 mm
eprops(2).orien = [0 0 1]; %Orientation of the cross-section as a vector pointing along "width-direction"
eprops(2).nbeams = 1; %Number of beams used to model this element
eprops(2).flex = 1:6; %Model out-of-plane bending (modes 3 and 4) as flexible
eprops(2).color = 'grey'; %Color
%% OPTIONAL ARGUMENTS
opt.filename = 'contact'; %Filename
opt.loadsteps = 20;
opt.customdyn = {'usersig PUSHSIG'}; %Unfortunately, this is case-sensitive ...
%% CALL SPACAR_LIGHT
out = spacarlight(nodes, elements, nprops, eprops, opt);
pos = out.node(4).p(1,:);
figure
plot(pos)
Some more notes:
* The number of load steps matters, and probably in relation to the value of the fictitious stiffness. The contact force is applied using the ''opt.customdyn'' line, whereas the actuation force is applied using the ''nprops(3).force'' line. Note the differences between ''nprops.force'' and ''nprops.force_initial''. The latter probably won't work because SPACAR Light will ignore the ''opt.loadsteps'' (which you do need when using ''usersig'') when there is no ''nprops.force'' present, it seems.
* When there is some error with the usersig interface or the particular m-file, debugging can be hard, because you only get the general SPACAR Light error. The log-file will provide some details. For more details, it can be helpful to do a manual Spacar call and see what the error then is: in the command window, do ''spacar(-10,'contact')''.
* The input and output variables of the usersig m-file are strict and cannot be deviated from. If you need additional variables in the usersig m-file, a workaround is to use MATLAB's global variables.
* Multiple forces can be specified with a single usersig file. Each row in the ''f''-vector specifies a nodal force. Nodal forces specified with usersig may overwrite any forces specified with SPACAR Light force command, such as ''nprops(i).force''.