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 :
Substituting the above equations into Eq. () and collecting all known vectors on the right-hand side, we obtain a linear system for :
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|