The effect of variable stiffness of tuna-like fish body and fin on swimming performance

The work in this paper focuses on the examination of the effect of variable stiffness distributions on the kinematics and propulsion performance of a tuna-like swimmer. This is performed with the use of a recently developed fully coupled fluid-structure interaction solver. The two different scenarios considered in the present study are the stiffness varied along the fish body and the caudal fin, respectively. Our results show that it is feasible to replicate the similar kinematics and propulsive capability to that of the real fish via purely passive structural deformations. In addition, propulsion performance improvement is mainly dependent on the better orientation of the force near the posterior part of swimmers towards the thrust direction. Specifically, when a variable body stiffness scenario is considered, the bionic body stiffness profile results in better performance in most cases studied herein compared with a uniform stiffness commonly investigated in previous studies. Given the second scenario, where the stiffness is varied only in the spanwise direction of the tail, similar tail kinematics to that of the live scombrid fish only occurs in association with the heterocercal flexural rigidity profile. The resulting asymmetric tail conformation also yields performance improvement at intermediate stiffness in comparison to the cupping and uniform stiffness.


Introduction
Tuna fish, as one of the most derived members of the family Scombridae, has long been thought as an efficient swimmer during high-speed swimming (Donley and Dickson 2000, Fierstine and Walters 1968, Mariel-Luisa et al 2017. It is streamlined with a tear-drop-shaped body and a narrow caudal peduncle. One crucial feature of tuna fishes is their highlyforked semilunate caudal fin with dorsal-ventrally symmetric external and intrinsic morphology. The tail has a high aspect-ratio shape (defined as the ratio of the square of the span relative to the surface area) and mostly composes of bone and collagen fibres. The stiffened fin rays, the hypural plate, and the collagen fibres together form the main structure of the caudal fin, which withstands the majority of resistance fish experiences when it swims. These intrinsic configurations indicate that tuna tail is a composite structure, similar to the fish body. Anatomy studies revealed that the internal formulation of a tuna body includes anterior pointing arms, backbone, red muscle, main horizontal septum, myosepta, etc (Westneat et al 1993). The non-uniform distribution of vertebrae and caudal fin rays impart anisotropic structural flexibility and contribute significantly to the fish's swimming behaviour (Affleck 1950 andMcHenry et al 1995).
Attributed to the above mentioned morphological feature and associated high swimming efficiency of the scombrid fish, it is prudent to study the underlying swimming mechanism, and thus provide insight into the design of bio-inspired artificial underwater vehicles. Unfortunately, the aforementioned material stiffness and its distribution were either excluded in previous studies by assuming that the physical models are rigid (Triantafyllou et al 1993, Triantafyllou andTriantafyllou 1995); or simplified by using predefined body/fin deformation [see a series of work at University of Virginia, e.g., (Han et al 2020 and Zhong et al 2019)]. Indeed, stiffness was considered in many previous studies. However, those models were composed of uniformly distributed bending stiffness (Dai et al 2012b andHeathcote andGursul 2007).
The study of the impact of stiffness on the swimming behaviour of fish-like models can be found from a series of experimental work at Harvard University (Feilich and Lauder 2015, Lucas et al 2015and Mariel-Luisa et al 2017. In the early study by Lucas et al (2015), four rectangular plates were tested with the variation of stiffness along the length of the foils. Their results indicated that models with high stiffness anteriorly and low stiffness posteriorly outperformed the others with uniform stiffness profiles in terms of thrust generation and self-propelled speed. In another investigation, Feilich and Lauder (2015) modified the physical model to include an anterior part mimicking the fish body and a posterior part representing the caudal fin. Detailed examinations on the shapes, which varied from a forked tuna-like tail to an unforked tail with a deep peduncle, covered three different stiffness variations with a heave motion imposed at the leading edge of the model. It was found that there was no single 'optimal' tail shape exhibiting the best performance in all metrics studied, highlighting a complex interaction among body, tail shape and its material stiffness. Following the work by Feilich and Lauder (2015), Mariel-Luisa et al (2017) extended the tuna-like foil models to a broader parametric span including the variation of structural flexibility, heave amplitude and frequency. They compared the foil models' kinematics with the counterparts from a live tuna fish, and they found that stiffness and kinematics interacted subtly in effect on hydrodynamic performance, with no one stiffness maximizing both the thrust and efficiency.
Although some progress has been made in the study of stiffness effect on tuna-like swimmers, the main limitations in the above work are three folds. Firstly, the availability of realistic materials to build up the physical models is restricted for an experimental study, resulting in the limitation to select the parameters considered. Secondly, most of the foil models were constructed from plastic uniform material stiffness, while in reality, tuna fish body is characterized by their stiffness variation along the body length (Mariel-Luisa et al 2017). Moreover, the experiment is restricted to observe some detailed flow field information, e.g., the surface force and wake vorticity, which can be compensated by numerical modelling.
On the other hand, the conformation of the tuna tail is not extensively explored in the study of (Mariel-Luisa et al 2017), as the caudal fin is one of the significant factors affecting the swimmer propulsion performance. Experimental observation of live Scombrid fish revealed that the locomotion of tail is asymmetric during steady swimming within a wide range of swimming speeds, i.e., 1.2-3.0 L/s (L is the fish body length) (Gibb et al 1999). However, this dorsal-ventrally asymmetry is also frequently observed for bluegill sunfish during its braking manoeuvring, resulting in the loss of thrust and efficiency as compared with symmetric tail locomotion, e.g., the cupping tail movement (Esposito et al 2012, Flammang and Lauder 2009, Luo et al 2020b. Given the differences in morphology between these two fish species, i.e., scombrid and bluegill sunfish, what is the actual role of biologically asymmetric tail locomotion on the swimming behaviour for tuna-like swimmers remains unknown [see examples in (Feilich and Lauder 2015) and (Krishnadas et al 2018)].  Inspired by the above studies (Lucas et al 2015 andMariel-Luisa et al 2017), we systematically investigate the effects of non-uniform distributions of flexural stiffness on the kinematics and propulsion performance of a tuna-like swimmer (figure 1) using our recently developed fully coupled fluid-structure interaction (FSI) solver (Luo et al 2020b). Distinguish from previous studies, both the fish body and fin stiffness variation is considered. Specifically, through a non-uniform distribution of stiffness along the body (figure 2), the main objective is to examine whether a bio-inspired stiffness profile improves performance and yields more fish-like kinematics. Additionally, with the non-uniform stiffness distribution in the spanwise direction of the caudal fin (figure 3), we aim to explore the possibility to passively control the fin deformation and replicate some features of tail kinematics observed from live fish (Gibb et al 1999), and understand the effect of these tail conformations on the hydrodynamic force production. Particularly, we are curious about the role of the aforementioned asymmetrical tail conformation during steady swimming of scombrid fish and the comparison with the experimental and numerical findings of bluegill sunfish (Esposito et al 2012 andLuo et al 2020b).
The remainder of this paper is organized as follows: the geometry, locomotion kinematics, structural properties of the tuna-inspired flexible swimmer are introduced in section 2. The metrics to evaluate swimming performance is also defined in this section. In the next section, the governing equations of fluids and solids, as well as implemented numerical techniques, are presented. Section 4 provides numerical results including structural deformations, propulsion performance and flow field etc. The discussion is then presented in section 5. Finally, the conclusions are drawn in section 6.

Problem formulation
As mentioned earlier, the present tuna-like model in figure 1(b) is inspired by the experimental studies of (Feilich and Lauder 2015) and (Mariel-Luisa et al 2017) and has the same dimension and size as that in the experiments. In this study, we consider the stiffness variations of the body (from the leading edge to the peduncle, as shown in figure 1(b)) and the caudal fin separately, with an attempt to shed light upon their respective effect on propulsive performance and kinematics.
However, it is worthwhile to note that the present study does not attempt to reproduce the real fish in terms of its lifelike geometry or material features in-vivo. Instead, following the studies by Esposito et al (2012) and Zhu and Bi (2017), we focus on some key characteristics, e.g., anisotropic flexural rigidities and associated FSI extracted from a real fish.
The length of the model L is defined as the characteristic length in this problem, and its thickness h is 0.139 cm. The leading edge of this model matches 30% total body length point on a real live fish (Mariel-Luisa et al 2017). All the edges of the model are chamfered to ease our fluid solver mesh generation. In accordance with the experimental studies in (Feilich andLauder 2015 andMariel-Luisa et al 2017), the swimmer performs heave motion in the y-direction, i.e., the leading edge moves laterally in heave without the pitch motion, in a uniform flow along the positive x-direction with a velocity of U ∞ . The timedependent heave motion of the model is described by y = y 0 sin 2πft , where y 0 is the maximum heave amplitude and f denotes the oscillation frequency.
The dimensionless parameters are defined in this study as the Reynolds number Re = U ∞ L/ν with ν is the kinematic viscosity of the fluid; the mass ratio m * = ρ s h/ρ f L, with ρ s and ρ f represents the density of the solid and fluid, respectively; the reduced frequency where E denotes Young's modulus and I = b h 3 /12 is the area moment of inertia of the cross-section per unit height. It is worth noting that the height of the model b is variable along the body length. We take the unit height b as a reference here for simplification (Dai et al 2012b).

Structural models of the tuna-like swimmer
As aforementioned, the stiffness distributions of the body and caudal fin are considered separately to avoid the interactive effect. Specifically, two scenarios are considered here, i.e., one is varying stiffness along the body while the other has different spanwise stiffness distributions in the tail.

Stiffness variation along the body length
The body of this model which is considered as 70% total length of a real fish is chosen to be composed of 21 segments in our study [the vertebral number for various scombrid fishes ranges from 22 to 66 (Fierstine and Walters 1968)], as shown in figure 2. Each body segment is assigned with a unique K, i.e., for the ith segment, the normalized flexural rigidity is . . . , N, where N = 21). To our knowledge, there are no direct stiffness measurements of tuna body in literature. Instead, inspired by the stiffness distributions measured by McHenry et al (1995) of a pumpkinseed sunfish, the variation pattern of K bi is written in the form of part of an elliptic equation as where K c is a constant and denotes the normalized flexural rigidity of the first segment. Meanwhile, a uniform profile K bi = K c is also used for comparison. These two variation patters are denoted by NU (non-uniform) and UB (uniform along the body length) modes for simplification. The variation patterns of the flexural rigidities of the body segments are depicted in figure 2(b). In all these cases, the stiffness of the first segment is 10 times that of the segment near the peduncle, which is within the range of stiffness variation naturally observed in live fishes (McHenry et al 1995). When the stiffness of the body segments is varied, the caudal fin has a uniform stiffness with a value of 0.04 K c which is derived from the estimation from the measurement in (McHenry et al 1995).

Stiffness variation along the tail
Regarding the flexural rigidities variation of the fin surface, inspired by the studies by Zhu and Bi (2017), the following stiffness profiles are used with an attempt to replicate some deformation patterns of the caudal fin observed from scombrid fishes in (Fierstine and Walters 1968) and (Gibb et al 1999), i.e., a cupping and heterocercal fashion. The total number of principal caudal-fin rays of a tuna fish is almost always 17 according to Fierstine and Walters (1968), and therefore, there are 17 segments in the fin part of our model, as shown in figure 3(a). Following the study by Zhu and Bi (2017) and our previous studies (Luo et al 2020b andShi et al 2019), the variation styles of K fi corresponding to different deformation fashions can be described as: (a) Cupping distribution: Here, R = 1 N N fi=1 R fi , and K m is the mean value of the stiffness of all the fin segments. The parameter γ is used to determine the ratio of the stiffness between the least flexible segment and the most flexible one. Like the study in (Luo et al 2020b), γ = 10 is selected in this work. The uniform distribution is also introduced here for a comparison. The three stiffness fashions are notated by CF (cupping style fin segments), HF and UF, and their corresponding stiffness profiles of the fin segments are presented in figure 3(b). Like the above practice, when the stiffness of the fin segments is varied, the stiffness of the body is uniform and has a value of 25K m . Following an experimental measurement of the kinematics of Sombrid fishes by Gibb et al (1999), we also placed seven maker points, as shown in figure 1(b), on the fin surface to monitor the fin deformation during the locomotion.

Performance metrics
The propulsion performance of the tuna-like swimmer is characterised by the mean thrust coefficient C T , the mean energy expenditure coefficient C P , the mean lateral force coefficient C y in the y-direction and the mean vertical force coefficient C z in the z-direction. These mean values can be obtained by averaging the instantaneous values over one locomotion period T. The instantaneous thrust generated by the model is defined as where S is the reference area, i.e., the area of the model in xz plane, and F x is the component of total hydrodynamic force in the x-direction. Similarly, the lateral and vertical force coefficients are written as where F y and F z are the components of hydrodynamic force in the y-and z-direction, respectively. For the present tethered model in the flow, the power expenditure coefficient can be evaluated as (Olivier and Dumas 2016) With the mean value of C T and C P over one period, the propulsion efficiency is given by

Governing equations and numerical approach
The present FSI model involves a finite volume method based fluid dynamics solver, a finite element method based solid dynamics solver and the coupler between the two. The fluid solver solves the unsteady, viscous and compressible flow which is governed by the laws of the conservation of mass, momentum and energy, and it can be written in an integral form as where W = {ρ f , ρ f v, ρ f E} T is the conservative variable vector with v is the velocity vector and E denotes the total energy, the vector F c represents the convective and pressure fluxes and F d is the fluxes arising from the viscous shear stress and thermal diffusion. Ω f denotes the fluid control volume with the boundary Γ f . n is the unit vector in the outward direction.
The fluid governing equation is discretized using a cell-centred finite volume method based on a multiblock structured grid system. The fluid domain Ω f is divided into an array of hexahedral grid cells. For each of these hexahedral cells indexed by i, j, k , the conservation laws are applied and reformulated in the semi-discrete form given by where ΔΩ f is the volume of the cell i, j, k , W i,j,k denotes the average flow variables of the cell, and R i,j,k is the residual which measures the net fluxes entering the hexahedral cell through all the six cell faces. An artificial viscosity term is introduced in R i,j,k to stabilize the scheme and eliminate the spurious numerical oscillations (Jameson et al 1981). The dual-time stepping scheme (Jameson 1991) is employed for time-dependent simulations, where equation (8) is reformulated as a steady-state flow problem with a pseudo-time t * where where the solution vectors of two previous time levels denoted by n and n − 1 are used here to yield a second-order accuracy. Equation (10) is integrated by a hybrid multistage Runge-Kutta scheme.
In this study, message passing interface based parallelization is achieved by domain decomposition to enable large-scale computation. The different grid blocks are automatically distributed over a number of processors by the block size with the application of a load-balancing algorithm. Furthermore, the local time-stepping and multigrid method are implemented to accelerate the convergence.
It should be noted that the present fluid solver resolves compressible Navier-Stokes equations. To ensure that the compressibility is small enough to be negligible, the freestream Mach number, defined as Ma = U ∞ /a ∞ with a ∞ denotes the speed of sound of the freestream, is chosen as 0.06, which is far below the critical value 0.3 where compressibility effect becomes pronounced, but meanwhile, still sufficiently large to ensure numerical stability. Besides, the local Mach numbers in the whole computational domain are under monitor during computation to guarantee that they are below the critical value. This compressible fluid solver has been successfully applied to different incompressible flow simulations in our previous biomimetic studies (Liu et  Regarding the structural dynamics, the basic equation is the weak form of the balance of momentum which is written in the differential form as where the acceleration of the material point is obtained by the second derivatives of displacement vector U of the structure, and surface forces are modelled by the second Piola-Kirchoff stress tensor P and body force of per unit mass such as gravity is represented by f. A constitutive equation describing the relation between the stress and the strain is used to close up equation (11). Specifically, for a Saint Venant-Kirchhoff material, the second Piola-Kirchoff stress tensor P is obtained by where C is the elasticity tensor, E represents the Green-Lagrange strain tensor, the deformation gradient is characterized by F and δ is the unit tensor.
The general governing equation of the solid dynamics, i.e., equation (11), is discretized using the finite element method. With the application of the virtual work method, we obtain a linear algebraic equation system by the discretization in the complete solid domain: The time domain is discretized using the αmethod here (Dhondt 2004). Denoting the velocity vector {V} := U and acceleration vector {A} := Ü , the solution at time level n + 1 can be obtained by where Ṽ n+1 and Ũ n+1 can be considered as the predictor at time level n + 1 which are only dependent on the values at time level n.
In this work, the finite element method based solid solver is CalculiX written by Dhondt (2004), in which a variety of element types are used to discretize the solid domain and define the shape functions.
Another main ingredient of an FSI solver is the coupling between the fluid and structure solver. To reduce the effort to adapt the original computation codes and preserve both the advanced features of the fluid and structure solvers, the two solvers are coupled via a partitioned framework based on preCICE (Bungartz et al 2016 and Luo et al 2020b). It is a challenge to simulate strongly coupled FSI problems (Tian et al 2014), like the current flexible swimmer propulsion, in which numerical instabilities may cause divergence when the densities of fluid and solid are comparable (Causin et al 2005). Therefore, within this framework, an implicit scheme is designed in which sub-iterations are introduced during each time step to ensure numerical stability and convergence.
The details of the numerical methods and the validation tests are provided in (Luo et al 2020a andLuo et al 2020b). This FSI solver has also been applied to the simulation of flexible swimmer propulsion in our previous work (Luo et al 2020b andLuo et al 2019).

Results
For all our simulations, the Reynolds number Re = 8000, mass ratio m * = 0.0089, the heaving amplitude y 0 = 1 cm, the Poison ratio ν s = 0.25. Most of the parameters are chosen to match with that in the experiment by Mariel-Luisa et al (2017). It is worth noting that the flow in our study is assumed to be laminar. At this Reynolds number regime (below or in the order of 10 3 ), the turbulence effect may play an insignificant role on the flow field, which was proved from some previous studies (Bozkurttas et al 2009 andBuchholz andSmits 2006  A self-consistency study is performed to assess the appropriate mesh and time-step resolution for f * = 2.5 when the stiffness of the body segments is varied in the UB pattern with K c = 0.1. Three grids are generated: a coarse grid with 2628 096 cells of which the minimum grid spacing is 1.48 × 10 −3 L, a medium gird with 4056 000 cells of which the minimum grid spacing is 9.73 × 10 −4 L, and a fine grid with the cell number is 5679 360 of which the minimum grid spacing is 5.95 × 10 −4 L. The computational domain and the medium fluid mesh around the tuna-like locomotor are shown in figure 4. On the model surface, the no-slip condition is applied, while for the other boundaries, the non-reflective farfield boundary condition is imposed. The structural mesh contains 4937 quadratic tetrahedral elements. For the three fluid meshes, different non-dimensional time-step sizes defined as Δt = Δt/T are used, i.e., Δt = 0.0087 corresponding to the coarse mesh, Δt = 0.0069 corresponding to the medium mesh and Δt = 0.0056 for the fine mesh. The time variation of C T within one locomotion period, when three meshes with different time-step sizes are used, are compared in figure 5. As seen, the results yielded by the medium and fine mesh are quite close. Therefore, the medium mesh and Δt = 0.0069 are used for our following simulations to reduce computational cost while retaining sufficient accuracy.

Midline kinematics
The midline kinematics envelopes of the tuna-like models with uniform and non-uniform stiffness variations along the body length when f * = 2.5 are depicted in figure 6. The variation pattern of the model when K c = 2 corresponds to the first mode defined by Michelin and Llewellyn Smith (2009). They characterized the vibration modes of the flexible wings according to the number of necks in the enclosing envelope. Therefore, other patterns correspond to a second mode. We also find that material flexural rigidity governs the number of waves  and an increase of wavelengths along the model with increasing inflexibility is observed here which was also demonstrated in (Dai et al 2012a, Feilich and Lauder 2015and Mariel-Luisa et al 2017.
With a closer inspection of the quantitative results of the lateral displacements, as shown in figure 7, the difference in the excursion in the lateral direction between the two stiffness distributions is mainly seen between 40%-80% of the model length. The maximum tail tip displacement is experienced by the intermediate-stiffness model when K c = 0.2 of the model with NU stiffness variation pattern, followed by that of the model with UB profiles when K c = 0.2 and 0.06.
The total lateral displacements presented by live tuna kinematics measured by Donley and Dickson (2000), the S3 foil model in the experimental study (Mariel-Luisa et al 2017) and the current FSI tunalike swimmer model are compared in figure 8. It can be found that the latter two show different kinematics variation patterns compared with that observed from a real kawakawa fish. The most evident difference between the tuna-like models and tuna is that the minimum lateral displacement of the tuna fish is located at around 20%-40% of its 'thrust producing' body length which corresponds to 44%-58% of its total length, while for most of the other models they are observed around 50%-80% of the models' length, as shown in figures 7 and 8.

Propulsion performance
The results of C T , C P and η when K c is varied under different heave frequencies are depicted in figure 9. As can be observed from figure 9(a), the generated thrusts generally increase with a larger locomotion frequency. The effects of body stiffness distribution patterns on the thrust production are not monotonous. Under the parameters studied here, the NU mode generally creates larger thrust than the UB mode by a small majority (15 out of the total 26 cases). Nevertheless, the advantages of the NU mode are only presented in low frequencies, i.e., f * = 2 and f * = 2.5, where larger thrust force is yielded by the NU style in more than 70% of the cases (13 out of 18 cases). In contrast, for a higher frequency when f * = 3.7, the UB mode produces larger thrust for most cases (6 out of the 8 cases). Based on an experimental study of flexible rectangular foils in (Shelton et al 2014), one would expect the maximum thrust occurs at K c = 8 which has the largest stiffness when the frequency is fixed. However, our results indicate that at this point, thrust is not the largest when f * = 3.7 even there is no thrust generated under a small frequency. This finding corroborates a similar conclusion by Mariel-Luisa et al (2017) that the stiffer models do not always produce more massive thrust.
The kinematic patterns are not reliable indicators in predicting the swimming performance, as suggested by Mariel-Luisa et al (2017). For example, the models with the NU stiffness variation have the most 'fish-like' kinematic curvatures at K c = 0.2 (see figure 6), and therefore they are expected to present high performance. However, their thrust production is poor in some circumstances, e.g., at f * = 2.5 and K c = 0.2, compared with that of the cases when K c = 0.06 and K c = 0.5 at the same frequency. On the other hand, based on the experimental results of rectangular foils by Lucas et al (2015), they suggested that a larger lateral displacement, especially the tail tip displacement, leads to a larger thrust. By the comparison between figures 7 and 9(a), we find that for the UB mode when K c = 0.06 and f * = 2.5, the tip displacement is 2.48 cm and C T = 0.12. Nevertheless, for the NU mode, it yields a tip displacement as 2.41 cm and a larger thrust coefficient C T = 0.16 at the same stiffness and frequency. Consistent with the experimental results in (Mariel-Luisa et al 2017), this indicates that the tail tip displacement does not necessarily predicate the propulsion performance in isolation.
In figure 9(b), significant distinctions in C P between the two modes are observed for very flexible and stiff cases when f * = 3.7, and at this frequency, the NU mode is more energy-saving (7 out of the 8 cases) compared with the UB mode. When the frequency is smaller, i.e., f * = 2 and 2.5, the difference of C P between the two stiffness styles are less noticeable, except for a few stiffness values, e.g., K c = 2. A quantitative comparison between these two modes reveals that a smaller C P is seen for the NU mode in a majority of the total cases (18 out of the total 26), although some of the differences are marginal, e.g., when K c is near 0.2.
Regarding the variations of propulsion efficiency, the effect of frequency on η is not so dominant as that on thrust and power expenditure. Namely, a high frequency does not always yield high efficiency, especially when f * = 3.7. Generally, the more flexible models are almost always more efficient, which is in line with the experimental results in (Mariel-Luisa et al 2017). This may not be applicable when the flexibility is sufficiently high, as indicated by FSI studies in (Dai et al 2012b). Given the same locomotion frequency, the NU mode performs more efficiently than the UB mode in most cases (19 out of the 26 cases).
It is interesting to compare the present numerical results with that of flexible flapping wings. For example, in the numerical study of a 2D flexible flapping wing in forward flight, Tian et al (2013) found that the thrust value always peaked at certain wing flexibility as the flexibility was varied (see figures 3(a) and 4(a) in their paper). However, in this work, a global and several local thrust peaks are reached as the bending stiffness is varied at a fixed frequency, as shown in figure 9(a). The difference of variation patterns of power expenditure coefficient and efficiency between the current tuna-like model and the flapping wing can also be found in contrast with the results in (Tian et al 2013). This may be attributed to the different model shapes and kinematics imposed on the models. The present model comprises a 3D body and a forked tail while a 2D flexible plate model was used in the study of Tian et al (2013). Besides, in (Tian et al 2013), the flapping wing performed asymmetrical combined translational and rotational locomotion, while only a heave motion is applied to the tuna-like swimmer here. When pure heave locomotion was applied to the models, the occurrence of several local thrust peaks as the variation of stiffness was also reported in (Dai et al 2016, Ryu et al 2019and Zhu et al 2014. As shown in figures 9(a) and (b), the data plotted using conventional dimensionless parameters defined in section 2.2 is not organized concisely. Thus, it may be interesting to investigate a new scaling parameter to present the data. Inspired by the scaling parameter study by Kang et al (2011), two nondimensional parameters are defined here, i.e., the effective stiffness 1 = Eh 3 / 12 1 − ν s 2 ρ f U 2 ∞ with h = h/L is the thickness ratio, and the relative tip deformation λ = w tip − w root /y 0 where w tip and w root is the displacement of the tail tip and the root of the model. The resulting scaling plotted in the logscale for the UB and NU stiffness styles is presented in figure 10. Two linear fits are used to approximate the correlation between log 10 (C T / 1 ) and λ when the value of log 10 (C T / 1 ) is positive and negative, respectively. When the frequency is small, i.e., f * = 2, the values of log 10 (C T / 1 ) are all smaller than zero and their relation to λ is well represented by the linear fit with a coefficient of determination (R 2 ) as 0.93. Under higher frequency, especially when f * = 3.7, the variation of log 10 (C T / 1 ) with λ is less regular, indicating a more complex interaction between the structure and the fluid. With an inspection of figure 10(b), we find that more than half of the points with high frequency, i.e., f * = 3.7, lie above the linear fitted line. In contrast, for other smaller frequencies, they are more likely to be seen below this line, which indicates that frequency has a significant effect on the power expenditure. The time histories of the thrust and power input within one locomotion period of the two stiffness distribution modes along the body length are depicted in figure 11. As can be observed, the non-uniform stiffness profile only slightly changes the phase positions of the peak and valley values of C T and C P compared with uniform distribution. Meanwhile, it significantly increases the amplitude of the instantaneous thrust, e.g., a 23% increase of the peak thrust from the UB to the NU mode. Moreover, the case with non-uniform stiffness distribution generates no drag throughout the entire motion period, which is reminiscent of a previous study on a flexible pectoral fin which suggested that fish can avoid the creation of drag force by complex 3D conformations (Mittal et al 2006). In contrast, the power expenditure only shows a minor difference, thus leading to a significant increase of propulsion efficiency (45 % from the UB mode) at K c = 0.2 and f * = 3.7, as observed in figure 9(c).

Near-body flow field
The wake structure around the tuna-inspired models is visualized in figure 12. Remarkable cow-horn shaped posterior body vortices (PBVs) are generated near the dorsal edge (PBV(D)) and the ventral edge (PBV(V)) in the wake of swimmers with the UB and NU mode. Similar PBVs were also reported in the numerical simulations of the swimming a Crevalle Jackfish by Liu et al (2017) and a bluegill sunfish by (Han et al 2020). The dorsal and ventral PBVs are compressed towards the root of the caudal fin, as shown in figures 12(e) and (f ), which has also been presented in (Zhu et al 2002) (see figure 8 in their paper) and (Liu et al 2017) (see figure 11 in their paper). This vortex compression is due to the narrowing peduncle at the posterior fish body. Leading-edge vortices (LEVs) and trailing-edge vortices (TEVs) are seen near the caudal fin whose strength is weaker than the PBVs. In comparison, the previous shed TEVs of the swimmer with the NU stiffness profile is stronger than that of the UB mode (see figures 12(e) and (f)). Additionally, the tooth-shaped vortices are seen near the first quarter model length which is covered by high-pressure when the swimmer flaps at the right-most position and is about to stroke reversal.
The Z-vorticity formulation within the xy plane around the locomotors with the NU and UB stiffness fashion is presented in figure 13. As can be seen, the vorticity of the two cases is qualitatively similar. At this instant, a pair of LEVs of the body and TEVs of the tail can be observed clearly. With the illustration of the streamline, we can observe remarkable vortex flow, especially for the NU mode, near the left surface of the leading edge part of the propulsor. By comparison, the clockwise vortices near the trailing edge of the tail with the NU mode are slightly larger than that of the UB mode. The dense distribution of the streamline near the wake of the trailing edge also indicates the suction effect on the flow velocity, which may contribute to the thrust production.
To visualize the pressure distribution along the model surface, we depict the pressure coefficient contours on both sides of the model in the xz plane in figure 14. The area and magnitude of the high (left side surface, figure 14(a)) and low (right side surface, figure 14(b)) pressure regions of the UB mode are larger than that of the NU mode (figures 14(d) and (e)). These high and low-pressure regions may correspond to the counterclockwise and clockwise LEVs of the body, respectively, as shown in figure 13. For instance, the counterclockwise LEVs dominate at the left side surface of the body. Likewise, high-pressure distribution is observed near the anterior body part on the same side surface as shown in figures 14(a) and (d). This is reminiscent of a numerical simulation of the mackerel-like swimmers by Borazjani and Daghooghi (2013) which suggested that the LEVs could alter the pressure distribution on the tail. By comparison, the pressure difference in the anterior part of the UB mode is more pronounced than that of the NU through direct observation. However, it appears hard to apply this method to evaluate the pressure difference in the posterior part. Despite this, the configurations of the propulsors in the xy plane present very different bending patterns, i.e., the tail of the model with the NU mode as shown in figure 14(f) flexes to a more considerable extent and thus showing a much larger pitch angle. This leads to a better orientation of the hydrodynamic forces along the negative x-axis direction, which benefits the thrust generation directly. The force vector, the magnitude of C T and C y for the two stiffness distribution fashions within one motion period are presented in figure 15. An inspection of figure 15(a) reveals that the force generated by the model with the NU mode is almost always better oriented in the thrust direction, although in some cases the magnitude of the force is smaller than that of the UB mode. This leads the NU mode to have a larger thrust than the UB mode in the entire motion cycle, as demonstrated in figure 15(b). As shown in figure 9(c), the values of C P for the two stiffness styles are quite close at K c = 0.2 and f * = 3.7, which indicates that this orientation of forces attributed to the flexing patterns of the tail does not require additional power expenditure. As a result, higher propulsion efficiency is obtained by the NU mode. The larger pressure difference at the anterior part of the model with the  UB mode as aforementioned and shown in figure 14, however, only leads to a greater lateral force in the ydirection as depicted in figures 15(a) and (c), which may be adverse for the straight cruising of swimmers.

Results when stiffness distribution of the tail is varied 4.2.1. Tail kinematics
The instantaneous deformation patterns of the swimmers with two different stiffness profiles assigned on the fin segments are presented in figure 16. The dorsal and ventral lobes of the tail with the CF mode are symmetrical with respect to the middle horizontal plane, and similar conformation is observed from the UF mode and thus not shown. The tail with the HF stiffness profile yields an asymmetry of movement, i.e., the dorsal lobe leads the ventral lobe during the flapping. To quantitatively analyse the tail kinematics, the movement of the dorsal tail tip, i.e., point 7 in figure 1(b), in x, y and z-direction is plotted in figure 17. As can be seen, the present numerical tuna-like propulsors present relatively little tail movement in the vertical (z) and horizontal (x) dimensions but the majority of the locomotion occurs in the lateral (y) direction during one tail-beat cycle, which aligns with the experimental measurement of scombrid fishes by Gibb et al (1999).
By observing the time histories of the dorsal tailtip lateral displacement and tail height in figure 18, we find that the amplitude of the tip displacement which is around 2.6 cm is close to that range around 3 cm observed from live Scomber japonicus fishes with a  similar total body length (around 25 cm) (Gibb et al 1999). Besides, the variation range of the tail height is about 0.5 cm and its variation period is around half of that of the varied tail-tip displacement, which is in line with the measurements by Gibb et al (1999).
The maximum displacements in x, y and zdirection of the seven points 1-7 on the fin are shown in figure 19. The amplitude of the tail excursions tends to be smallest on the peduncle and mid-tail regions and larger at the tail-tips in all the three dimensions, which agrees with the measurements of live S. japonicus fishes by Gibb et al (1999). The excursions in x (horizontal) and z (vertical) direction are quite small compared with that in y (lateral) dimension whose magnitude is almost an order larger than the former.
The wave of the lateral displacement is propagated posteriorly as shown in figure 20. For instance, the ventral and dorsal peduncle points reach the maximum excursions in lateral dimension around 10% of the flapping period time ahead of the tail tip. However, the ventral tail-tip reaches its maximum displacement approximately 7% of the cycle time behind the dorsal tail-tip.

Propulsive capabilities
The results of C T , C y , C z , C P and η when the stiffness is varied for the three different stiffness profiles  are presented in figure 21. In general, the effects of different flexural rigidity distribution patterns on propulsive performance are mainly noticeable in the vertical-forces generation when the stiffness is the same, while others present close results. It also applies to the scenarios when the locomotion frequency is varied, and therefore they are not shown in this study.
With an inspection of the curves of thrust coefficient, we find that the tuna-like swimmers generate quite close thrusts unless at the very flexible and moderate stiff cases. At intermediate stiffness, models with the HF stiffness profile produce larger thrust compared with the others. For example, at K m = 0.005 and 0.02, the thrust generated by the swimmer with the HF stiffness patterns increases by 4.8% and 4.0%, respectively, from that of the CF mode. Regarding the lateral (y) force production, the magnitudes of C y for both the CF and UF fashions firstly decrease as stiffness is increased under highly flexible cases, and then experience a general increase with a larger K m . In comparison, the lateral force by the HF mode almost always increases as the models become stiffer under the parameters studied. Only the swimmer with the HF stiffness profile yields non-negligible vertical forces, as presented in figure 21(c). The values of lift force (C z ) by the HF stiffness model are generally smaller by an order of the magnitude of the thrust forces, which is in line with the experimental measurements of the forces produced by chub mackerel fishes (Nauen and Lauder 2002). In terms of power  . Phase lag measured as the percentage of tail-beat cycle period illustrating the effect of location on the fin on the timing of lateral (y) locomotion of the fin. The dorsal tail-tip is defined as the reference location and therefore, has a zero phase shift. A negative value indicates that the point reaches its maximum lateral displacement before the dorsal tail-tip. expenditure, they all experience a continuous increase with larger inflexibility. On the contrary, the propulsion efficiency monotonously decreases as the swimmers become stiffer after the peak at K m = 0.001. The performance drop presented in figures 21(a) and (e) at very flexible cases, i.e., K m = 0.0008, may be attributed to the declining ability of a highly flexible swimmer to communicate momentum to the flow to induce thrust production, which has been demonstrated in (Michelin andLlewellyn Smith 2009 andOlivier andDumas 2016).
To accurately distinguish the effects of body and fin on the thrust production, we split the total thrust of the swimmer into two parts as shown in figure 22. The cases when the body is rigid are also considered for comparison. As can be seen, the difference in thrust among different stiffness distributions is indeed derived from the tail. Especially when the body Figure 21. The time-averaged coefficient of thrust C T , lateral forces C y , vertical forces C z , power input C P , and efficiency η when the stiffness along the fin is varied when f * = 2.5.  is rigid, the thrust of the HF stiffness profile doubles compared with the others, which indicates a remarkable interaction between the body and tail during flapping.
The instantaneous variations of thrust and power expenditure over one locomotion cycle are depicted in figure 23. The total value of the model including the body and fin part and the partial value of the  fin alone are both presented for comparison. With a closer inspection of figure 23(b), we find that the stiffness distribution along the fin appears to have little effect on the variations of C P both of the entire model and the tail alone. Nevertheless, it has a more significant impact on the thrust production of the fin part, although this difference is almost entirely eliminated when combined into the overall consideration. Figure 24 demonstrates the wake structure near the tail of the tuna-like swimmers with the CF and HF inflexibility patterns. The wake structures near the body are similar to that when the body stiffness is varied, as shown in figure 12. Therefore, only the vortex formation near the caudal fin is presented here. With an inspection of figure 24, we can find that the TEV of the CF mode generally presents a good symmetry relative to the middle line in the z-direction. A closed vortex ring is observed at the wake of the trailingedge of the swimmer with CF mode. In contrast, the counterpart near the tail of the swimmer with the HF stiffness style has an opening at the dorsal lobe, which indicates the symmetry is broken here.

Flow field near the swimmer
The Y-vorticity contour and near-body streamline by the locomotor with the CF and HF profile are depicted in figure 25. As seen, the height of the vortices near the tail tip, i.e., the secondary trailingedge vortices, is approximately equal to the caudal fin height, which is consistent with the wake structure of chub mackerel fish obtained by using digital particle image velocimetry techniques (Nauen and Lauder 2002). Their results also showed that the vortex jet was oriented at a slightly negative angle around −3 degree relative to the horizontal x-axis. This is also presented in figure 25(b) where the perpendicular line to the streamline tilts from the vertical direction, indicating that the flow is slightly pushed downward along the negative z-direction as pointed by the black line.
The pressure distribution along the swimmer surface is presented in figure 26. The main difference between the CF and HF mode is that much lower pressure (marked by a black circle in figure 26(d)) is located at the ventral lobe of the tail at the right side surface for the HF stiffness fashion. Therefore, a more considerable pressure difference between the left (high pressure) and the right (low pressure) is generated by the ventral lobe of the swimmer with the HF stiffness pattern. With an observation of the tail conformation at this instant in figure 16, one may find that the ventral half of the tail rolls upward. This orientation of the tail gives rise to a positive vertical component of forces resulted from the pressure difference, and it cannot be balanced by the dorsal lobe which is almost upright in vertical. This may explain the production of lift by the HF stiffness profile. In contrast, although there is also pressure difference generated both on ventral and dorsal lobes of the tail of the swimmer with the CF mode, these two forces are counteracted by the symmetrical distribution dorsally and ventrally.

Passive control via non-uniform stiffness distribution
A few previous studies have indicated that it is possible to imitate some morphological features of fish using passive control via imposing appropriate stiffness distribution (McHenry et al 1995). A study by Videler (1993) revealed that the rigidities of the pectoral fish fin are enhanced near leading-edge as rays are bonded together. This phenomenon was reinforced with numerical modelling of a flexible fin ray by Shoele and Zhu (2012). The mechanism behind was the lessening of the effective angle of attack in the vicinity of the leading edge, which was reflected by the mitigation of LEVs separation (Shoele and Zhu 2012). The experimental work from (Lucas et al 2015) also provided evidence that the model with a biologically relevant stiffness, i.e., the stiffer anterior, presented more fish-like kinematics than uniform foils.
Similar fish kinematic features are also revealed in the present study when the attention is paid to fishtail. For example, the symmetrical and asymmetrical rigidity styles lead to rather different tail conformations, as shown in figure 16. Generally, the HF profile causes similar features of scombrid fishtail previously observed in the experiment of Gibb et al (1999). It is noted that such asymmetry only presents by HF profile, where the lateral excursion of the dorsal tail-tip is 9.6% larger than that of the ventral tail-tip as depicted in figure 19. The time-dependent dorsal-ventral asymmetry tail movements are also noticeable for the HF stiffness fashion (see figure 20). Although previous studies by Fierstine and Walters (1968) and Gibb et al (1999) reported similar trends with a caudal fin of Skipjack tuna fish and a chub mackerel fish, to our knowledge, for the first time, such asymmetry trends in lateral displacement magnitude and phase shift of the tail-tip is well replicated in this numerical FSI study covering a variety of flow and structure parameters.
However, it is challenging to replicate true-tonature fish kinematics entirely relying on pure passive control via non-uniform stiffness. For instance, the location where the valley value of the lateral displacement occurs by the NU stiffness profile does not match with that of a live tuna fish (figure 8). Additionally, the variation pattern of the lateral displacement of the numerical swimmer differs much with that of the real fish data at the first 60% model length. This may suggest that active muscle contractions/extractions play a dominant role in the formulation of body waveform. Regarding the tail kinematics, the present models show a difference in the phase shift of the curves of tail-tip displacement and tail height with the experimental observations (Gibb et al 1999), as depicted in figure 18. Particularly, their results suggested that the tail is maximally compressed at the maximum lateral displacement of the tail tip. Instead, in our results, the maximum lateral displacement almost appears at the same instant as the maximum abduction of the tail. Such difference is much likely the consequence of the sophisticate caudal fin muscle active control, which is hard to be achieved by purely passive deformations. This finding is also reminiscent of the speculation that the cyclical vertical compression of the Scombrid fishtail is a result of the action by interradialis muscle which is positioned such that contraction of this muscle could draw the dorsal and ventral rays towards one another by Gibb et al (1999). Our numerical results corroborate this opinion.
In addition to the above locomotion kinematics, passive control on the flow field can also be achieved indirectly by the non-uniform stiffness distributions as presented in (Shoele and Zhu 2012), where the strengthened leading-edge caused the mitigation of LEVs separation. In the present study, the bio-inspired non-uniform body stiffness profiles yield slightly stronger trailing edge vortices (figure 13) and alter the pressure distribution, reflected by reduced pressure at the anterior body surface and near the peduncle as shown in figure 14. Collectively, the nonuniform stiffness causes reorientation of the fluid force so that it points more towards the swimming direction and thus increases thrust as presented in figure 15. This applies to the thrust augmentation by the heterocercal stiffness profile along the fin surface. HF profile does not change the pressure magnitude very much but results in the rolling motion of the ventral lobe of the caudal fin, yielding a resultant force along the thrust direction (see figures 26 and 16). As a consequence of that, a lift force is also generated (figure 21(c)).

The function of heterocercal conformation of mackerel fishtail
Scombrid fish has a homocercal tail with dorsal-ventrally symmetrical external and internal morphology. Due to this symmetric feature, it has usually been thought to function as a homocercal model (Gibb et al 1999). Kinematics measurement suggested that its tail flaps asymmetrically so as to provide lift force during steady swimming (Gibb et al 1999), which is presented by the heterocercal stiffness profile in our numerical models. Indeed, previous research has indicated that the fish body is negatively buoyant, tending to push fish moving towards the substratum (Magnuson 1973). To prevent this sinking trend, the body is tipped up to deliver additional upward lift near the anterior part of the fish which is balanced by the vertical lift produced posteriorly by the tail (Aleev 1969). This theory leads to a general prediction that the neutrally buoyant fish do not show asymmetrical tail conformation during swimming.
Interestingly, biological observations have suggested that significant dorsal-ventral asymmetry and tilting may also appear in some teleost fishes with near-neutral buoyancy. (Gibb et al 1999, Lauder 1989, 2000and Webb 1993. Subsequent research revealed that this dorsal-ventrally asymmetrical tail deformation of bluegill sunfish is closely related to their manoeuvring behaviours (Flammang and Lauder 2009). Inspired by this finding, a biomimetic bluegill sunfish tail whose shape is plump without a sharp fork was numerically studied in our previous study (Luo et al 2020b). Our previous results suggested that the asymmetrical heterocercal tail conformation continuously yielded the smallest thrust and lowest efficiency during steady swimming among all the deformation patterns, including cupping and uniform styles. Similar conclusions have been drawn from the experiment on a robotic fish caudal fin with imposed tail motion derived from biological observation of the bluegill sunfish (Esposito et al 2012). However, in the present study on the swimmer with a tuna-like tail, a heterocercal stiffness profile and its resultant asymmetrical deformation do not cause the deterioration of propulsion performance. Instead, in some instances, when the stiffness is at an intermediate level, the locomotors with the HF pattern even outperform the others in terms of thrust generation and propulsion efficiency (see figures 21(a) and (e)).
By comparing the present results with the other studies, we may propose one additional explanation/hypothesis to the asymmetric scombrid fishtail conformation: the role of it may play is not only to contribute to the lift force balance but also the thrust generation and propulsion efficiency. This is suggested at all speeds during steady swimming of mackerel fish. However, this does not apply to bluegill sunfish with a different tail shape. Most of the bluegill fish are neutrally buoyant, and thus there is no need to balance the gravity and buoyancy force when they are swimming. In this situation, they adopt the asymmetric tail movement to offer additional lift force along with the thrust reduction during manoeuvring.
It is also found that the stiffness profiles adopted in this study do not yield as remarkable impact on the propulsion performance as that in the results from (Zhu and Bi 2017) and (Luo et al 2020b) where a bluegill sunfish inspired tail model was used. For example, the largest relative thrust difference was 29.3% seen between the models with cupping and heterocercal stiffness profile in (Luo et al 2020b). In contrast, the largest thrust distinction is seen between the HF and CF stiffness styles with a difference of 4.8% in this study. This may be related to the different tail shapes tested and the different intrinsic musculature conformations of the tails of their biological prototypes. Morphologically, the scombrid fishtail has a larger aspect ratio with a highly forked trailing edge while the bluegill sunfish has an unforked tail with a smaller aspect ratio. Hydrodynamically, the different tail deformations due to variable stiffness are more likely to induce different force production for the bluegill sunfish thanks to their large control surface. On the other hand, biologically, the intrinsic tail myology also determines the role of the tail plays on swimming behaviour. Anatomical studies revealed that there is an extensive complement of intrinsic musculature of bluegill sunfish (Lauder 2015). It is believed that the complex conformations can control adduction and abduction of individual fin rays, the movement of fin rays and the relative motion of upper and lower tail lobes. The utilization of these intrinsic muscles activities enables excellent control of the tail surface during different locomotor behaviours (Flammang and Lauder 2009). In comparison, the intrinsic caudal fin musculature significantly reduces for scombrid fishes (Nursall 1963), and there is even no intrinsic tail musculature as to black skipjack tuna fish (Fierstine and Walters 1968).
In summary of the above results and discussions, the tuna-like tail may not be a favourable prototype when the manoeuvrability purpose is the focus. This is suggested from the results that the change of tail shape has little influence on the force generation, although their semilunate conformation is believed to offer high propulsion efficiency during high-speed (Nauen and Lauder 2002).

Conclusions
By the utilization of a fully coupled three-dimensional FSI solver, we have numerically studied a tunainspired swimmer. Specifically, we investigate the effects of variable stiffness distributions along the body and tail on the kinematics and dynamics of the locomotors separately. Firstly, a bio-inspired nonuniform rigidity profile of the body is compared with a uniform mode through systematic simulations. The numerical results indicate that given the parameters studied in this work, the larger thrust produced by the model with the bionic stiffness fashion is mainly seen under low frequencies. Instead, when the frequency is high, the swimmer with uniform stiffness produced larger thrust in most cases. The enhanced performance by the non-uniform stiffness mode is more noticeable in terms of propulsion efficiency where more than 73% of the total cases saw an increased efficiency, and this improvement is seen for all the three frequencies. Secondly, among the three different distributions of tail stiffness, i.e., heterocercal, cupping and uniform, the swimmer with a heterocercal pattern shows resemblance to that of real scombrid fishes in terms of tail kinematics. Additionally, those with the heterocercal stiffness profile also outperform that with other inflexibility distributions at an intermediate stiffness. The lift force produced by them is absent for the other two stiffness patterns. These findings suggest that the asymmetrical tail conformation does not only provides additional lift to balance swimming body but may also contribute to efficient propulsion during steady swimming of scombrid fish. This heterocercal tail deformation has distinctive functions compared to that of a bluegill sunfish whose caudal fin has superior abilities in the manoeuvre.
Throughout our results, we also find that it is impossible to achieve entirely real fish-like kinematics if only the passive control, via the variable body and fin stiffness proposed here, is adopted. These are reflected by the comparison between our results with the experiment in (Donley andDickson 2000 andGibb et al 1999) as we discussed in section 5.
It is reasonable to conjecture that these discrepancies are induced by the subtle and advanced active control of vertebrae and tail muscular activities by fish, which are unable to be considered in this study. More work needs to be done in the future in order to fully explore the complex interactions between the swimmer and its surrounded environment, which are driven by muscular actuation morphology and structural properties, and the resultant swimming kinematics and performance.