To simulate the dynamics of the fish model, the differential equations
of motion (
) must be integrated over time. The system
is intrinsically stiff and can become unstable in a variety of
situations that may occur in the simulated dynamic world
(e.g. executing a left turn to avoid a sudden dynamic obstacle before
completing a right turn). Therefore, the stability of the numerical
solution is an important consideration. We employ a simple,
numerically stable, semi-implicit Euler method.
The system of equations of motion Eq. () may be written
in the standard matrix form for nonlinear, damped elasto-dynamic systems
[Bathe and Wilson1976] as follows:
where 
, 
and 
are the mass, damping and
stiffness matrices, respectively. Matrix is a constant
diagonal matrix with 
on the main diagonal, while and
are time-dependent nondiagonal matrices. These three
matrices are of dimension 
where n is the number of
nodes in the biomechanical model. Matrices 
, 
,
consist of nodal position, velocity and acceleration
vectors, respectively, and 
is the matrix of external force
vectors. These four matrices are of dimension 
, such that
each row contains the three components of each vector.
Because of the special properties of the viscoelastic units in our
dynamic fish system, the damping term 
and the stiffness term 
in the above equation
are most conveniently collected into the form of an effective
stiffness matrix 
:
and hence Eq. () can be re-written as follows:
where 
is time-dependent since the viscoelastic forces depend
nonlinearly on the nodal position and velocity variables 
and 
.
Our method computes the left-hand side of Eq. ()
implicitly and the right-hand side explicitly. In particular, we
discretize continuous time into time steps, 
and solve for 
in
where t represents simulation time and 
is the simulation
time step. We do not solve for the position vectors 
directly because, first, the velocity vectors 
is
readily available since it has to be stored at each time step for
computing 
(more details can be found in the next section)
and, second, using saves us from having to store 
at each time step. We use the following finite
difference approximations for and 
:
and
Substituting the above equations into Eq. () and
collecting all known vectors on the right-hand side, we obtain a
linear system for :
where
is the effective system matrix and
is the effective load matrix. As we will explain in the next section,
the structure the system matrix 
is such that an efficient
solution to Eq. () can be obtained by assembling and factorizing it as follows:
where 
is diagonal and 
is lower triangular. The new
nodal velocities are then computed by
first solving the lower triangular system
for 
by forward-substitution and finally solving the upper
triangular system
for by back-substitution. Hence, the
numerical simulation of the fish biomechanical model requires the
assembly and factorization solution of a sequence of linear systems of
the form of Eq. for 
.
| Xiaoyuan Tu | January 1996 |
