next up previous contents
Next: System Matrix Assembling and the Up: Biomechanical Fish Model and Locomotion Previous: Muscles and Hydrodynamics

Numerical Solution

To simulate the dynamics of the fish model, the differential equations of motion (gif) 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. (gif) may be written in the standard matrix form for nonlinear, damped elasto-dynamic systems [Bathe and Wilson1976] as follows:

  equation554

where tex2html_wrap_inline2514 , tex2html_wrap_inline2516 and tex2html_wrap_inline2518 are the mass, damping and stiffness matrices, respectively. Matrix tex2html_wrap_inline2514 is a constant diagonal matrix with tex2html_wrap_inline2330 on the main diagonal, while tex2html_wrap_inline2516 and tex2html_wrap_inline2518 are time-dependent nondiagonal matrices. These three matrices are of dimension tex2html_wrap_inline2528 where n is the number of nodes in the biomechanical model. Matrices tex2html_wrap_inline2532 , tex2html_wrap_inline2534 , tex2html_wrap_inline2536 consist of nodal position, velocity and acceleration vectors, respectively, and tex2html_wrap_inline2538 is the matrix of external force vectors. These four matrices are of dimension tex2html_wrap_inline2540 , 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 tex2html_wrap_inline2542 and the stiffness term tex2html_wrap_inline2544 in the above equation are most conveniently collected into the form of an effective stiffness matrix tex2html_wrap_inline2546 :

displaymath2512

and hence Eq. (gif) can be re-written as follows:

  equation580

where tex2html_wrap_inline2548 is time-dependent since the viscoelastic forces depend nonlinearly on the nodal position and velocity variables tex2html_wrap_inline2550 and tex2html_wrap_inline2552 .

Our method computes the left-hand side of Eq. (gif) implicitly and the right-hand side explicitly. In particular, we discretize continuous time into time steps, tex2html_wrap_inline2554 and solve for tex2html_wrap_inline2556 in

  equation591

where t represents simulation time and tex2html_wrap_inline2560 is the simulation time step. We do not solve for the position vectors tex2html_wrap_inline2562 directly because, first, the velocity vectors tex2html_wrap_inline2564 is readily available since it has to be stored at each time step for computing tex2html_wrap_inline2566 (more details can be found in the next section) and, second, using tex2html_wrap_inline2564 saves us from having to store tex2html_wrap_inline2570 at each time step. We use the following finite difference approximations for tex2html_wrap_inline2562 and tex2html_wrap_inline2574 :

  equation608

and

  equation616

Substituting the above equations into Eq. (gif) and collecting all known vectors on the right-hand side, we obtain a linear system for tex2html_wrap_inline2556 :

  equation629

where

  equation634

is the effective system matrix and

  equation640

is the effective load matrix. As we will explain in the next section, the structure the system matrix tex2html_wrap_inline2578 is such that an efficient solution to Eq. (gif) can be obtained by assembling tex2html_wrap_inline2578 and factorizing it as follows:

equation647

where tex2html_wrap_inline2582 is diagonal and tex2html_wrap_inline2584 is lower triangular. The new nodal velocities tex2html_wrap_inline2556 are then computed by first solving the lower triangular system

equation652

for tex2html_wrap_inline2588 by forward-substitution and finally solving the upper triangular system

equation655

for tex2html_wrap_inline2556 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. gif for tex2html_wrap_inline2592 .


next up previous contents
Next: System Matrix Assembling and the Up: Biomechanical Fish Model and Locomotion Previous: Muscles and Hydrodynamics
Xiaoyuan TuJanuary 1996