## 1 Introduction  
The immersed boundary (IB) method was first proposed by [PERSON] in the 1970s to simulate the fully fluid-flexible structure interaction problem. This method uses a fixed Eulerian mesh to simulate the flow field and a set of Lagrangian points to represent the boundary immersed in the fluid. The interaction between the fluid and the immersed boundary is implemented by a discrete Dirac delta function, which spreads the force of the deformable elastic boundary to the nearby fluid grid nodes and interpolates the boundary velocity from the local flow velocity to update the boundary position. In the most early studies [1; 2; 3; 4], the immersed boundary is considered as a tensile structure, and hardly has any bending resistance. However, most real flexible materials are nearly incompressible (cannot deform under stretching or compression) and have bending rigidity (can resist bending and torsion). Therefore, the bending rigidity cannot be ignored, particularly when the deformation curvature is relatively large. The bending rigidity is ignored in the early literature, mostly because of the stability restraint of the IB method, where a small time step is required to maintain stability (much smaller than on the imposed time step by explicitly differencing the advective or diffusive term) [5; 6]. The restraint may be particularly severe for a larger bending rigidity or a low fluid viscosity. Much effort has been expended to alleviate this severe restriction. For example, the fluid viscosity has been artificially increased by two orders of magnitude [7]. However, this solution can only be used in special problems. A three-step approach to update the IB boundary velocity was presented by [PERSON] [8], but it could not obtain an obvious effect [9; 10]. A subspace iteration method was presented by [PERSON] et al. [11], which could remedy the conventional IB scheme to a limited extent but was difficult to expand to three-dimensional cases. A subgrid resolution IB method that improved the stability against the slender fibers issue, was proposed by [PERSON] [12]. Compared with the above methods, an unconditional stability IB scheme presented by [PERSON] et al. [13] has more advantages in terms of effectiveness and applicability. However, due to the need to iteratively solve large matrices, the scheme may consume more computational resources and time, so it has limited practicability. The lattice Boltzmann (LB) method is an alternative simulation technique for complex fluid systems. Its easy implementation, intrinsic parallelism, and good suitability for numerous fluid flow problems have been demonstrated by many works. The immersed boundary-lattice Boltzmann (IB-LB) integrates the advantages of the two methods, and has extensive prospects in research and application [2; 4; 7]. However, the IB-LB coupling scheme also suffers from the stability restriction when it simulates a flexible structure. Some scholars have attempted to alleviate this instability. [PERSON] et al. [14] integrated the three-step approach [8] with the IB-LB scheme, but this approach did not give the expected effect. [PERSON] and [PERSON] [15] introduced an unconditional stable scheme, where the time step could be increased up to 30 times compared with that of the explicit method. [PERSON] et al. [16] presented an implicit IB-LB coupling method based on the unconditional stability IB scheme [13]. This method can alleviate problems such as instability, restrictive choices of time step and boundary rigidity, and easily incorporate the boundary mass. [PERSON] et al. [17] proposed a robust IB-LB method employing the fractional step technique to simulate fluid-structure interaction problems. Unlike the methods mentioned above, which require iteratively solving large matrices, [PERSON] et al. [18] improved the stability of the IB-LB coupling scheme by using an iterative external forcing method, which is the basis of the present study, but no thorough investigation was reported in their paper.