Biorobotic Waterfowl Flipper with Skeletal Skins in a Computational Framework: Kinematic Conformation and Hydrodynamic Analysis

Cormorants (Phalacrocoraxe), types of aquatic birds, utilize the compliance/flexibility of the flippers and exploit hydrodynamic/biomechanic processes to accomplish diverse operations. Particularly, the flipper‐propelled locomotion exhibits traits such as super‐redundancy and large deformations, necessitating depiction of both movements of the rigid skeletons as well as local deformations of the soft tissues. However, there are few well‐established kinematic/hydrodynamic framework models and constitutive equations for such rigid–flexible intrinsically coupled biosystems. Herein, combined with a skeletal skinning algorithm to handle the deformation of a flexible body attached to a rigid body, a numerical computation framework for an in‐depth fluid–structure interaction is presented, which enables the capture of viscoelastic and anisotropic characteristics of a highly compliant 3D rigid–flexible coupled model in a low‐Reynolds‐number flow. Considering the biorobotic cormorant flipper with a nonuniformly distributed stiffness as a representative, the challenging issue of controlling a biomechanically compliant flipper to synthesize realistic locomotion sequences, including rigid skeleton movements and soft tissue deformations, is addressed. Furthermore, a numerical computational hydrodynamic analysis is performed to demonstrate that the cormorant flipper can generate 5 N fluid force and 0.45 N m fluid moment during the turning operation in 0.8 s, which is consistent with the former experimental results.

a reaction that drives the fluid vortex further. In this context, the hindlimbs reversibly generate rapid and stable structural deformations due to the fluid force. [14] In addition, the fluid force magnitude depends on the size and shape of the hindlimbs (particularly the flippers), orientation of the flippers relative to the locomotion, and speed of the flipper propulsion. Recently, we reported that the flipper has a slight wave component in the low-speed oscillatory motion during the aquatic posture adjustment of the cormorant, as shown in Movie S1, Supporting Information. In other words, the flexible flipper realizes the wave movement through a small swing of the toes and passive coordination of the pad surface, as shown in Figure 1.
Conceptually, the flexibility and potency of the cormorant flippers are considered the foremost features for an optimal swimming performance. We previously studied the high-speed flapping motion of the cormorant flippers during the water surface takeoff, [15][16][17] in which the flippers are modeled as solid bodies or piece models. To handle the oscillatory motions of the biomimetic flippers with a strong resemblance, we need to investigate the dynamic musculoskeletal architecture modeling with compliant mechanics. The unsteady locomotion in liquid or gaseous fluids is a vital aspect of the shape deformation/morphing of animals and robots. The underwater locomotion of the cormorant involves an intricate fluid-structure interaction (FSI), with the flexible flipper characterized by compliance and prone to substantial displacement deformations. Additionally, the flipper ray stiffness varies from proximal to distal end. [18] Studies indicate that the use of appropriate structural stiffness values in biomimetic robot locomotion may increase the propulsion performance. [19] Furthermore, the cormorant may actively manipulate the flipper rays to form specific propulsive surface conformations to efficiently cope with the complex and unstructured flow field. Consequently, due to their intrinsic distributedness, suppleness, and compliance, biological musculoskeletal architectures exhibit an ability to outsource control tasks to their embodiments, an emerging paradigm denoted as morphological computation or mechanical intelligence. Hence, the underwater locomotion of flippers is primarily affected by the dynamic issue of the bidirectional mutual coupling between the flexible structure and fluid. Thus, it is required to address the flow field evolution, structural motion, and interaction between the two as accurately as possible.
As the flow field comprises vortices rotating at constantly varying directions and speeds, the fluid forces and moments may be obtained based on disordered velocity/vorticity fields, which are subtly determined either experimentally or computationally. For experimental implementations, the two principal means are hydrodynamic force platform (HFP) and particle image velocimetry (PIV). HFP requires a direct mechanical integration of the pressure and shear fields via instrumentation, whereas PIV utilizes an indirect numerical integration of the measured flow field. The experimental approaches are appropriate for a constrained locomotion, but are ineffective for an analysis of deforming bodies in freely locomoting animals and robots. [20] In bionics, aside from animal welfare concerns, imposing constraints may also impede or disrupt the natural movement, particularly when considerable body deformations are required to generate propulsion. [21] For computational implementations, several new nonintrusive approaches have been applied to study biomimetic aquatic locomotion, which integrate measured kinematics with computational fluid dynamics (CFD). [22][23][24] Among them, the cardinal is to synthesize the www.advancedsciencenews.com www.advintellsyst.com anatomically biomechanical characteristics and comprehensive biomimetic locomotion sequences with an impeccable reality. [25] In our recent research, we quantitatively analyzed the 3D maneuvering trajectories of the flipper with underwater stereo-videography. [26] However, current computational frameworks generally require certain simplification assumptions of model rigidization, neglecting the real complex curved surfaces, relative internal locomotion, and active volumetric variable deformations. Conversely, natural flippers have uneven shapes and intricate biological components, leading to a problematic acquisition of adequate deformation parameters. Moreover, the hydrodynamic flipper locomotion is accompanied by elastic and volumetric deformation caused by muscle expansion/ contraction. Such inadequacies are addressed concurrently in this study. Musculoskeletal systems consist of bones, muscles, tendons, ligaments, and other connective tissues that together provide function and structure to flippers. To synthesize realistic anatomically detailed flipper locomotion in a physics-based manner, we construct a comprehensive flipper model with synthetic hard (skeleton) and soft (skin) tissues that are properly coupled and simulated to produce natural life-like flipper motions in its environment while accounting for the mentioned accuracy issues and compensating for them with relevant algorithms. The skeleton transformations, i.e., the aforementioned rigid kinematics maneuvering trajectories, characterize the position and orientation of the articulations. The skin deformation is computed by skinning techniques, in which the core idea is to establish a hierarchical skeleton for the articulated character to attach a mesh model (skin) to the frame, and then drive the mesh model movement using the animation data of the skeleton. [27,28] Regarding skinning, more attention was devoted to skeleton-based linear blend skinning (LBS) methods, [29] in which the surface deformation is restricted to the skeletal posture that entirely specifies the skin deformation. However, LBS generally relies on optimization in a vast parameter space with numerous degrees of freedom, which is computationally costly and frequently leads to the issue of local minima. Therefore, LBS is typically confined to basic characters/scenarios [30] and cannot be anticipated to capture complex deformations. Moreover, the deformation is restricted to the subspace of the affine transformation of the skeletons, so that such a simplex method has challenges in deforming the skin near articulations due to volume loss or well-known "candy-wrapper" artifacts, which can be eliminated by nonlinear blending methods such as dual quaternion skinning (DQS). [31] In general, it is not feasible to rule out any particular mode of deformation a priori, particularly when compatible architectures interact with unstructured and dynamic surroundings, hence necessitating efficient, robust, and more broadly applicable solvers. Skeleton-skin-driven motion generation has attracted considerable attention as the demand for similarity and realism in motion increases. Researchers have been developing increasingly sophisticated biomechanical models. [32,33] Nevertheless, they seldom consider heterogeneous structures capable of undergoing all modes of deformation.
In this study, we consider the biorobotic flippers of the cormorant as a representative to explore the hydrodynamics process in the flipper-propelled locomotion, attempting a general research solver approach and performing a related verification. We constructed a morphological 3D model reconstruction of the hindlimb skeleton, whose rigid-body kinematics in 3D space can be adequately expressed mathematically. We then demonstrate a methodology in which the skeleton skinning approach is used to represent bionic heterogeneous bone-muscle layouts and biomimetic rigid-flexible kinematics realistically. Furthermore, biorobotic systems with various material characteristics are combined with the flow field through suitable boundary conditions. Finally, we propose a versatile numerical computational framework of the musculoskeletal architectures' locomotion/ deformation and extend the validation of CFD in which fluid flows around an active deformation by merging the newly developed skeleton-skin-driven modeling approach.

Musculoskeletal Skinning
In general, an intricate locomotion can be decomposed into the overall rigid-body motion of individual nodes and flexible deformation inside the surface formed by the nodes. For a cormorant, as shown in Note S1, Supporting Information, the webbing surface deformation is caused by the skeletal displacement/rotation and muscle expansion/contraction; hence, it is vital to achieve a realistic and accurate musculoskeletal skin-deep deformation to improve the quality of the numerical computation. We investigated morphological feature data of the cormorant flipper and performed a relevant modeling, as shown in Note S2 and Figure S1, Supporting Information. To facilitate the processing of the surface skinning mesh and subsequent numerical calculations, the model was smoothed without affecting the motion characteristics, as shown in Figure S2, Supporting Information. Owing to the difficulties in directly capturing the muscular skin surface mesh deformations from the continuous locomotion of the hindlimb, we proceeded with the standing posture, which is readiest to validate surface mesh details, to develop the skeleton skin mapping algorithm. The modified standing state composing the positions of 13 articulations is described by the minimax optimization algorithm Moreover, the positions of V 11 -V 13 can be obtained by the morphological formulations in the preceding subsection. In particular, the lengths of bones and initial angles between toes can be obtained by the standing state articulation positions, as shown in Table S1 and S2, Supporting Information. Hence, the resulting modified parameters according to Equation ( The skeletal skinning approaches the correlation between the skeletons and 3D muscular meshes via the former standing posture, and then executes the mapping transformation based on the relative positioning relationship of the skeletons between time series. For the muscles covering the skeleton, real-time and intuitive interactive skinning deformation algorithms are indispensable to alter the articulated musculoskeletal models. We utilized a skeletal skinning framework to drive the geometric and behavioral features of the cormorant's flippers, consisting of three parts: the skeleton, skin (mesh), and bonding between them. The skeleton is a collection of bones representing geometric features, the skin refers to the body surface mesh, and various weight functions and interpolation methods perform the binding. An updated position of the skin vertex of the flexible deformation is acquired by weighing and summing the multiple bones affecting the vertex, where the weighting factor is chosen depending on the relative position of the bones, physical and biological constraints, etc. Algorithmically, we select a series of control nodes on the geometric model, and calculate the influence weighting factor of the model surface by the above control nodes. We then drag the control node, and the geometric model deforms correspondingly. The weighting factor in the first stage impacts the effectiveness of the skeletal skinning algorithm. Therefore, if the flesh and muscles of the skeletal flipper are to undergo natural and high-quality deformation, there must be an accurate weight factor calculation method to build the connectivity of the skeletal nodes and surface meshes. The utilized strategy is the approach of bounded biharmonic weights [34] to address the skeleton weight distribution and interpolation method of DQS [35] to implement the skeletal skinning process, which can achieve comparable speeds to those of LBS while increasing the visual quality of the animation. The bounded biharmonic weight formulation is described by arg min W j, j¼1, : : : , m where the boundary conditions are The resulting framework is intuitive for existing physics-based animation methods as we alter both surface and entire volumetric model. The obtained distribution of skeletal weight is presented in Movie S2, Supporting Information, and more details can be found in Supporting Information Visualization.

Mesh Reconstruction
Based on the quaternion formulation of each bone and locomotion fitting of each articulation node, we can derive the rotation quaternion and dual quaternion for each bone of the cormorant's flipper at each instant. Combined with the above skeletal weighting distribution, the skin meshes of the flipper can be constructed to match the skeleton locomotion at each instant by DQS interpolation.
where b q B i, stand ðtÞ is the corresponding dual quaternion indicating the conversion from the standing posture to the current condition of the bone B i , x stand is the mesh location at the standing posture, x t ð Þ is the mesh location at the current condition, and w i x stand Þ ð is the weighting factor of the bone B i at x stand . Surface skin meshes of the flipper at each instant based on the standing posture can be constructed from the above algorithm through mapping transformation of Equation (4), as shown in Figure 2a. The reconstructed meshes retain the characteristics of the flipper at the standing posture: the toes are "separating" (η toe is relatively large) so that an obvious distinction is observed between the toes and webs, among which more meshes describe the webs. The metatarsophalangeal (MTP) articulation is "extending" (θ mtp < 0) while standing and "flexing" (θ mtp > 0) during the majority of the power stroke, which changes largely. Therefore, the reconstructed meshes seem unnatural at MTP articulation during flexion. These properties pose the following challenges for the hydrodynamic computations: 1) When the toes are "gathering" (η toe is relatively small) (e.g., t ¼ 0 ms), the meshes of the web tend to overaggregate and lead to overlapping. Nevertheless, a clear distinction between the webs and toes is not conducive to the mesh formation required for hydrodynamics; and 2) When the MTP articulation is "flexing" (e.g., t ¼ 700 ms), the proximal meshes surrounding the MTP articulation tend to overaggregate and overlap. In contrast, the distal meshes tend to be overdispersed, with sharp edges.
Therefore, although the above algorithm can reconstruct meshes that are faithful to the standing posture with distinct characteristics at the toes, the morphology and details of the flipper have significantly changed from the standing posture to the motion process. Thus, the mesh distribution set for the standing posture is inappropriate for the entire motion process. Consequently, the DQS mapping from the standing posture (not t ¼ 0 ms) is not favorable to the hydrodynamic computation of the whole motion chain, necessitating an additional smoothing optimization.

Mesh Optimization
To address the stated challenges on the originally reconstructed meshes, we apply the following approaches for optimization: 1) Shrink-wrapping approach [36] is applied to reestablish uniform meshes at t ¼ 0 ms (not the standing posture) above. The overlapping meshes are removed, and restricted meshes are set to describe the toes and webs, which is the foundation for mapping reconstruction; and 2) The skeleton weight functions and DQS quaternion are recalculated for the altered meshes at t ¼ 0 ms. The mapping reconstruction is carried out by Figure 2b illustrates the web meshes that are no longer sufficiently aggregated to cause overlapping. The mesh size distribution is relatively uniform at each place for various times. The drawback comprises the unobvious characteristics of toes and webs, so that the dimensional changes between the toe and webs are blurred at t ¼ 400 ms. However, the impact of such characteristics on the hydrodynamic force condition is limited, and the skins still follow the locomotion of the skeletons accurately in the displacement and deformation variation. As the most relevant hydrodynamic characteristics are reserved, the process of smoothing optimization is reasonable.

Hydrodynamic Model
We provide a multiphysics computation framework for biological locomotion based on the finite-volume approach, within which we develop a thorough biomechanical model of the flipper. As illustrated above, the smoothing optimized reconstructed meshes are utilized as surface meshes of the flipper. In ANSYS FLUENT 19, [37] the fluid region is established between the red cuboid boundary and flipper surface, as shown in Figure 2c. We assume that the fluid is nonviscous and that the boundary layer has no transition, which inevitably diminishes the accuracy of the calculation results but enhances the efficiency and feasibility of the computational processing. The proposed CFD solver has been partially validated in our previous biomimetic numerical computations. Relevant details are provided in Note S3, Supporting Information. In the coordinate system, z ¼ 0 represents the water surface. At t ¼ 0 ms, the x and y coordinates of the knee articulation are 0, while the þx direction is directly in front of the body. The hydrostatic pressure at z ¼ 0 equals 0 Pa. The distances from the origin to the cuboid boundaries are shown in Figure 2d. Assuming that the flipper is sufficiently far away from the cuboid boundary, we set the cuboid boundary as the given static pressure boundary, whose surface pressure minus the static water pressure equals 0 Pa. The surface boundary of the flipper is a nonslip wall.
Three million nonuniform hexahedral meshes are constructed in the fluid area for a fast processing and compatible accuracy. The meshes near the flipper surface are progressively intensified to polyhedrons for an improved adaptability in the complex www.advancedsciencenews.com www.advintellsyst.com structure-fluid area. The flippers' solid portions are constructed with tetrahedral meshes to obtain the relatively best computation accuracy. Considering the mesh-driven incompressible singlephase turbulence in hydrodynamics, we develop Reynolds turbulence governing equations, as shown in Note S4, Supporting Information. The provided detailed information on hydrodynamic variables such as the velocity, flow field, and pressure of biomimetic locomotion at macro-and microlevels can be obtained by solving the mathematical and physical models iteratively, in which the computation time step is 1 ms, as shown in Figure 3.

Skeleton Morphology
Our snapshots of the cormorant's lower limbs indicate that the "knees" seem to be bent backward in a manner that differs from that of the human bones and articulations, as shown in Figure 4a.  Figure 4b. For cormorants, the majority of their hindlimbs' musculature is located very high on the thigh bone and hip, whereas the lower swinging elements of the legs are comparatively light, moved by long mass-reducing tendons. Analogously, the arrangement optimizes the cormorant's legs and flippers for a high-velocity locomotion, providing them with both large step length and high step frequency. Below we list the nomenclature of the bones and articulations, as shown in Figure 4c. Finally, we unify and provide a concise nodal description of the bones and articulations for a better mathematical and physical analysis, as shown in Figure 4d.
To facilitate the commissioning of hydrodynamic numerical computations and comply with the physical laws of anatomical bones and articulation nodes, we propose the following morphological definitions and formulations: 1) The knee, ankle, MTP, and second toe (P kne , P ank , P mtp , P t21 , P t22 , and P t23 ) are coplanar on the leg plane with its normal vector z leg pointing inward (left side). The relative coordinates of the leg use P kne as an origin, the y leg axis is collinear with the TIB bone, and the leg plane can be constituted by x leg and y leg axes, as shown in Figure 5a; 2) The MTP and first articulation of each toe (P mtp , P t21 , P t31 , and P t41 ) are coplanar on the flipper plane with its normal vector E toe1 www.advancedsciencenews.com www.advintellsyst.com pointing to the sole of the flipper, which is perpendicular to the leg plane (E fli ⊥ E leg ), as shown in Figure 5b; 3) The MTP and each toe are coplanar on the appropriate toe plane with its normal vector pointing inward (left side); e.g., P mtp , P t11 , P t12 , P t13 , and P t14 are coplanar on the first toe plane with its normal vector E toe1 pointing inward (left side). Likewise, P mtp , P t31 , and P t32 are coplanar on the third toe plane with its normal vector E toe3 pointing inward (left side). Each toe plane is perpendicular to the flipper plane (E fli ⊥ E toe1 , E fli ⊥ E toe3 ), as shown in Figure 5c; and 4) The angle between the toe planes changes proportionally. The initial angles in the relaxed standing state are θ t12,0 (negative), θ t32,0 (positive), and θ t42,0 (positive), as shown in Figure 5d. The toe narrowing factor is defined as 5) The angle between the adjacent articulation bones is same, indicated by the toe flexion angle θ toe , as shown in Figure 5e. Therefore, the reference frame of leg P kne xyz leg contains the origin at the knee articulation P kne (x kne , y kne , z kne ) and rotation vector α leg (α x , α y , α z ) of xyz leg relative to the world coordinate. α leg (α x , α y , α z ) is described in terms of interaxial angles, which represent the rotation about the axis of the unit vector (α x , α y , α z )/ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi Compared to the quaternions, the interaxial angle can be determined by the least three scalars, while the quaternion formulation requires four. The interaxial angle can uniquely predict the shortest rotation path, while the Euler angle cannot. Along with the ankle angle θ ank , MTP angle θ mtp , toe narrowing factor η toe , and bone flexion angle θ toe , we propose a mathematical model that simplifies the analysis via ten state parameters of the cormorant's flipper position ξ ¼ x kne , y kne , z kne , α x , α y , α z , θ ank , θ mtp , θ toe , η toe Â Ã Based on the preceding assumptions, when the bone lengths and initial angles between the toes (i.e., θ t12,0 , θ t32,0 , and θ t42,0 ) are provided, the position state of the cormorant flipper at any instant may be established, granted that the above ten sequential parameters are derived. Accordingly, the following detailed kinematic transformations can also be determined. Leg relative rotation  [38] from John Wiley and Sons. In cormorants, the bone between the ankle and toes, the tarsometatarsus, is considerably longer than in humans and serves as a functional equivalent to the human's shin bone. c) Morphological classification based on ref. [38]. www.advancedsciencenews.com www.advintellsyst.com ¼ cos α þ w 2 x ð1 À cos αÞ w x w y ð1 À cos αÞ À w z sin α w x w z ð1 À cos αÞ À w y sin α w y w x ð1 À cos αÞ À w z sin α cos α þ w 2 y ð1 À cos αÞ w y w z ð1 À cos αÞ À w x sin α w z w x ð1 À cos αÞ À w y sin α w z w y ð1 À cos αÞ À w x sin α cos α þ w 2 z ð1 À cos αÞ where y leg is the unit direction vector of bone B 1 . Bone vectors. Unit direction vector of B 2 E leg2 ¼ X leg sin θ ank þ Y leg cos θ ank (9) Vertical component of the B 2 vector (E leg2 ) on the leg plane E leg2v ¼ X leg cos θ ank À Y leg sin θ ank (10) Unit vector in the direction of B 7 www.advancedsciencenews.com www.advintellsyst.com E t21 ¼ E leg2v sin θ mtp þ E leg2 cos θ mtp (11) Unit normal vector of the flipper plane E fli ¼ ÀE leg2v cos θ mtp þ E leg2 sin θ mtp (12) Unit vectors in the directions of toes Articulation positions Bone rotations.
To facilitate the computations, we integrated the representation of a bone rotation relative to the parent bone coordinate in both terms of interaxial angle and quaternion. The former has been mentioned previously and is described as (a, b, c), which is the unit vector of the rotation axis. The latter is described as Q (w, (a, b, c)) corresponding to the interaxial angle where w is the rotation angle. The relationship between the bones is known from Figure 5. The relative rotation quaternion solutions for all bones are presented in Table 1.
Based on the parent-child relationship between the bones, as illustrated in Figure 5, the representations of the absolute rotation quaternions can be solved by the chain rule (parent bone quaternion Â child bone relative quaternion), as shown in Table 1. For example, the absolute rotation quaternion of B 5 of a case is Qðθ ank , ð0, 0, À 1ÞÞ ⋅ Qðη toe θ t12,0 , ð0, À 1, 0ÞÞ⋅ Qðθ toe , ð0, 0, 1ÞÞ ⋅ Qðθ toe , ð0, 0, 1ÞÞ In summary, the position vector of the root node, i.e., the knee articulation P kne (x kne , y kne , z kne ) and absolute rotation of each bone, can be described by the above ten sequential parameters, as shown in Equation (7), to fully convey the skeleton locomotion of the flipper.

Skeleton Kinematics
The propulsion performance of an oscillating flipper is governed by numerous characteristics such as the motion amplitude, oscillation frequency, fin geometry, and fin stiffness. This study intends to capture the locomotion of the main articulations of the leg and flipper through two high-speed cameras, acquire the image information of the spatial feature points, and reconstruct the image in a real 3D scene using the image projection on the 2D plane, as shown in Note S5, Supporting Information. The visual detection algorithm calculates the critical spatial positions of the articulations when the cormorant is at a particular posture and obtains the real-time 3D coordinate data of each articulation point through data processing, as shown in Figure S3, Supporting Information. The measured position data of V 1 -V 10, shown in Figure 4b, can be obtained from an experiment. [26] Preliminary processing of experimental data is required, and a transformation is carried out: the height direction is z direction, the plane z ¼ 0 is coincident with the water surface, while x is pointing in the direction of the knee-head vector. To describe the locomotion of the cormorant by fitting the experimental data of its flipper, the fitting expression ξ(t) needs to be optimized to minimize the objective function. Considering the multiple fitting objectives in this study, we use the minimax method [39] min The parameter i is the serial number of target data at an arbitrary time point and j is the time serial number. Thus, the total number of objective functions equals i Â j. In the application of the above method, two main challenges are encountered: choice of the form of ξ (t) and objectives that should be listed as the objective function F i (ξ (t)). This study aims to achieve hydrodynamic calculations, for which the cormorant's flipper trajectory needs to be remarkably smoothed. The physical locomotion of the cormorant's flipper should also be continuous and smooth so that the tenth-order function is chosen to describe the ξ (t) variation with time Objective function Table 1. Bone variables and their relative rotation quaternion parameters.

Bone
Relative rotation quaternion Bone Relative rotation quaternion The constraints are θ ank > 0 and η toe > 0.2, where the fitting parameter ξ m contains 11 Â 10 scalars. The variations in the components of ξ m with time are shown in Table S3, Supporting Information.
The comparisons of the minimax-based skeleton generation results to the original articulation positions are shown in Figure 6. We also list the calculated deviations of the minimum locomotion deviation results of the ten articulation (V 1 -V 10 ) Figure 6. Multimax skeletal articulation optimization results based on the minimax method, illustrated by Equation (21). Each plot represents the optimized flipper locomotion (solid black line) and observed articulation position (red square frame) at given time. The light green shading groups quantitatively show deviation effects of the optimization results of each articulation relative to the observed values, whose detailed values are shown in Figure S5, Supporting Information. The fitting deviations are infinitesimal, which indicates that the optimized locomotion can reflect the natural movement of the flipper and can be applied for a hydrodynamic calculation analysis.
www.advancedsciencenews.com www.advintellsyst.com points, marked at the top of each subplot, as shown in Figure S4, Supporting Information.

Discussion
As shown in Figure 7a, the displacement of the knee articulation in the z axis is restricted and disregarded. It mainly varies along the x and y axes, particularly on the y axis. The key characteristics of the knee articulation locomotion include almost absence of change in the first 450 ms, slight increase in the x coordinate, and continuous decrease in the y coordinate after 450 ms. When the knee articulation has a large displacement in the -y direction, the cormorant's body rotates counterclockwise around the z axis (i.e., the angular velocity component of the cormorant's body in the z direction is larger than 0). The counterclockwise rotation is related to the force on the cormorant's flipper, so that the following analysis focuses on the forces and moments before and after 350 ms. The variation curve of the TIB rotation vector α leg (α x , α y , α z ) is presented in Figure 7b. Although the rotational locomotion of TIB can only be collectively characterized by the three components of the rotation vector rather than by any one of them, a qualitative analysis can still be carried out with each component.
α y can qualitatively describe the pitch swing (up and down) of the TIB. The more considerable value represents that the TIB is closer to the horizontal posture. When the value is approximately equal to 0, the TIB is closer to the bottom. During the period of 790 ms, α y is initially relatively large, dropping to near-zero rapidly in the end, which indicates that the TIB sways in a higher position at the early stage before the last swing driving it near the yz plane. α y exhibits a significant downward swing of the TIB from 150 to 350 ms and substantial uplift from 350 to 500 ms.
During the remainder, approximately 100 ms, it swings downward.
α z can qualitatively represent the yaw swing of the TIB. α z oscillates around 0, which indicates a limited yaw vacillation, so that the TIB locomotion is basically parallel to the xy plane. α x qualitatively describes the rolling swing of the TIB and TMT articulation (hereinafter referred to as "TIB-TMT"). α x was negative during the period of 790 ms, which indicates that the direction of the TIB-TMT has a nonzero rightward component (-y direction) during the whole period. α x decreased significantly from 200 to 400 ms, which implies that, simultaneously as the TIB-TMT downswings, the articulation and instep are turning rightward (-y direction). At 350 ms, the TIB movement transits from the end of the downswing to the rising phase, and the TIB-TMT turns rightward. After 400 ms, the TIB-TMT maintains a significant tilt rightward.
The locomotion of the TMT and flipper relative to the TIB can be described by the variation in fitting angles at the ankle, MTP, and toes presented in Figure 7c. θ ank as 0 implies that the articulation is in a straightened state. It is reorganized, larger than 0, only in the first 100 ms and 300-600 ms, with the remainder of the period swinging around 0. At 300 ms, the TIB is swaying downward and rightward when a sharp increase in θ ank implies that the TMT accelerates rearward, carrying the MTP articulation. At 420 ms, the TIB-TMT is obviously tilted rightward with the TIB rising. θ ank begins to fall, which indicates that the TIB is driving the ankle and MTP leftward.
At 300 ms, θ mtp is still at a flexion angle close to 90°and tends to further flexion. Thus, it can be speculated that the hydrodynamic force has a negligible effect on the flipper. At 400 ms, θ mtp also increases, which indicates that the MTP stretches the flipper. However, θ toe is at its peak, so that the toes are in www.advancedsciencenews.com www.advintellsyst.com an excessive flexion. Hence, similarly, the impact of the hydrodynamic force on the flipper can be regarded as negligible. However, after 420 ms, due to the ankle driving the flipper leftward, the toes begin to straighten with the decrease in θ toe . In contrast, the flipper begins to stretch with the increase in η toe , and the MTP articulation continues to straighten with the increase in θ mtp . These actions represent that, after 420 ms, the flipper exerts a significant force on the water in the þy direction (i.e., the fluid exerts a considerable force on the flipper in the -y direction) by increasing the stressed area and exerting a leftward force. Until 500 ms, both θ toe and θ ank are close to 0, and η toe is relatively large, which implies that the flipper and toes are straightened, and the webs are stretched to a rather large extent. Therefore, it can be speculated that the flipper is under a heavy load at the moment.
Later, from 500 to 570 ms, θ mtp decreases again, and the MTP starts to bend although θ ank still decreases. The acting force pulling the MTP leftward by the ankle still exists. θ toe is relatively small, with a large η toe representing a rather sizeable stressed area. Nevertheless, the leftward force decreases with the decrease in θ mtp .
The results of the kinematic analysis can be summarized as follows. From approximately 200 ms, the TIB begins to swing downward. The TIB-TMT is in a state of chained wavy transmission. Owing to the flexion and retraction of the MTP and toes, it is estimated that the flipper is subjected to a smaller hydrodynamic force during the period. Until 420 ms, with the locomotions such as the TIB-TMT articulation straightening and leftward lifting, MTP and toe straightening, and web stretching, there is a leftward force against the water at the bottom of the flipper until 670 ms. At 500 ms, θ mtp reaches a maximum; meanwhile, the flipper's force on the water should also reach its peak. After 500 ms, the force of the flipper on the water decreases with the decreases in θ mtp and η toe . The force has a crucial role in the cormorant's anticlockwise turning mentioned above, consistent with the time point. The description in this section corroborates the force analysis results below.
Fluid forces can be recovered by measuring the control surface's velocity, pressure, and shear stress fields. The moments acting on a deforming body determine its stability and maneuverability. If the cormorant turns anticlockwise without initial rotation, it needs to be subjected to a counterclockwise external M z, tot in the z direction on the cormorant's center of mass (COM); it is not the COM of the leg. The moment can be decomposed into two parts: M z, tot ¼ M z, leg þ F y, leg L x : 1) The cormorant's leg is subjected to the external M z, leg , and the center of M z, leg is chosen as the origin in this study; 2) and The cormorant's leg is subjected to F y, leg L x þ F y, leg L y generated from the external force F leg , using the distance vector from the cormorant's COM to the equivalent center of F y, leg as the moment arm. Therefore, either the -y direction or þx direction of F leg positively contributes to the counterclockwise moment. However, as the arm L x of F y, leg (the projection distance from the center of M z, leg to the cormorant's COM in the x direction) is considerably larger than the arm L y of F x, leg (the projection distance from the center of M z, leg to the cormorant's COM in the y direction), F y, leg is the most critical determinant of F y, leg L x þ F y, leg L y .
Before the analysis of the force details, it is recommended to clarify that we primarily focus on the data between 100 and 700 ms owing to two reasons: 1) The initial effect affects the calculation of the initial 100 ms, which leads to an exaggeration of the hydrodynamic force. The initial effect refers to the phenomenon when the hydrodynamic analysis is performed, assuming that the fluid velocity is 0 while the moving mesh wall velocity does not equal 0. It will lead to an excessive impact effect and result in a more considerable calculated fluid stress than the actual; and 2) The low fitting accuracy of the tracking data of the flipper during the first and last 100 ms makes the results of these two stages unreliable because the time-continuous method is used for data fitting in this study. In the fitting of the data at a particular time, a high fitting accuracy can be guaranteed only when there are more data before and after the present time. Figure 8a shows the curve of M leg with time, including M z, leg . Between 100 and 700 ms, there are two significant bulges in the M z, leg curve: one at approximately 170 ms and other around 500 ms; that at 500 ms is more significant. Both of them have positive values, which have a positive role in the anticlockwise   Figure 8b shows the curve of F leg with time, including F y, leg . F y, leg also has two prominent bulges between 100 and 700 ms: one at approximately 170 ms and other around 500 ms; that at 500 ms is more significant. Both of them are negative, and have a positive effect on the anticlockwise rotation of the cormorant's body. F y, leg remains positive for a relatively long period between 200 and 400 ms, but the value is small. The negative interval around 170 ms almost offsets the momentum. The behavior of F y, leg between 200 and 400 ms corresponds to the description from the above kinematic analysis. The cormorant's leg is a wavy transmission during this period, driving the cormorant rightward. However, due to the MTP flexion, toe retraction, and web closing, the rightward force against the water exerted by the flipper exists but is limited. Further information such as the contours of fluid pressure is presented in Movie S3, Supporting Information. Figure S6, Supporting Information shows the distribution of pressure and flow field of the cormorant's leg surface viewed from the þy direction. The hydrostatic pressure is the water pressure formed by gravity when the water is still, which has a limited influence on the single-phase incompressible fluid flow. Thus, it is not displayed when the flow field is shown. At 500 ms, the flipper is under the maximum pressure; meanwhile, the instep is loaded by a relatively sizeable negative pressure. In addition, during the locomotion of the cormorant's flipper, the water flows continuously along the þy direction in the area on the left side of the flipper bottom, which is another proof of the leftward force exerted on the flow. The locomotion and stress condition indicate a one-to-one correspondence between kinematic and dynamic analysis results, as shown in Movie S4, Supporting Information.

Conclusion
In this article, we presented a newly developed multiphysics framework to address bioinspired FSI swimming issues, synthesizing the locomotion of a biomechanical flipper along with the aquatic environment in which it is situated, as well as the flexible physical interaction. Besides, the skeletal subspace deformation algorithm was applied to establish the relationship between the skeletons and skins (meshes) for an active deformation and passive compliant locomotion of the surface muscles. First, the morphological and kinematic analysis of the skeleton rigid body is performed. Then, the skeleton control points are set, and the weight of each skeleton segment is calculated. Next, the parameters of the skinning algorithm for each bone are adjusted to obtain the locomotion and deformation characteristics of the skin soft body. Finally, by entering them into our FSI computational framework, we can extract the related hydrodynamics results in fluid locomotion. The cormorant can complete a turn in 0.8 s by a lateral drive of the feet. A single flipper can produce a peak fluid force of 5 N and a peak fluid moment of 0.45 N m, which are basically consistent with the experimental results. [40,41] The results help better understand the characteristics of the locomotion mechanism and hydrodynamics in the unsteady locomotion of the compliant flippers, which are significant in exploring the energy accumulation and kinematic adjustments of the cormorant's pretakeoff phase. Our approach to the biological prototype is rather general and may be exploited, thus providing a perspective direction to extend various fluid applications, such as aquatic soft biorobots inspired by the morphology and swimming kinematics. Biorobotic prototypes with rigid skeleton-soft skin materials can be given kinematic locomotion parameters in our multiphysics framework to obtain specific deformation solutions for a computational hydrodynamic analysis, in which the surface soft skins can deform regularly with the movement of the central rigid skeletons. The numerical results can inform and improve rigid-flexible coupling structures, biomimetic soft composites, and integrated technologies (flexible actuators, sensors, controllers). In the future, with technological developments in bionic drives and materials, higher fidelity bionic robot validation can be attempted, which, in turn, will facilitate the optimization of the numerical computational framework.