Immersed boundary simulations of fluid shear-induced deformation of a cantilever beam

08/05/2019 ∙ by Sudeshna Ghosh, et al. ∙ Simon Fraser University 0

We derive a mathematical model and the corresponding computational scheme to study deflection of a two-dimensional elastic cantilever beam immersed in a channel, where one end of the beam is fixed to the channel wall. The immersed boundary method has been employed to simulate numerically the fluid-structure interaction problem. We investigate how variations in physical and numerical parameters change the effective material properties of the elastic beam and compare the results qualitatively with linear beam theory. We also pay careful attention to "corner effects" -- irregularities in beam shape near the free and fixed ends -- and show how this can be remedied by smoothing out the corners with a "fillet" or rounded shape. Finally, we extend the immersed boundary formulation to include porosity in the beam and investigate the effect that the resultant porous flow has on beam deflection.



There are no comments yet.


This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In this paper, we derive a mathematical model and the corresponding computational scheme to study the deflection of a two-dimensional deformable elastic cantilever beam in response to a surrounding fluid flow. We will study cantilever beam immersed in a channel flow, where one end of the beam is fixed to the channel wall. We will not impose a given shear flow but rather drive the motion via two walls: one fixed and one moving with constant velocity. In the absence of a beam or other channel obstruction, the flow would develop into a simple linear shear flow (also called Couette flow). We also study the effect of introducing a small porosity into the beam structure, and compare the flow over porous and solid (non-porous) beams. The choice of such a porous cantilever is motivated by our interest in modelling biofilm structures in the near future alpkvist-klapper-2007 . A biofilm is a collection of microorganisms (immersed in fluid) that adhere to each other and often also on a nearby surface. Fluid-structure interaction (FSI) plays a key role in several phases of the biofilm life cycle, one being in the deformation of biofilm columns due to fluid forces. We extract from this scenario an idealized 2D model problem in which a rectangular cantilever beam deforms in response to a shear flow. The method we develop may also be useful for simulating other systems arising in applications from biology (e.g., glycocalyx, cilia) and engineering (e.g., MEMS, micropillars).

The study of a cantilever beam is a classical problem in solid mechanics timoshenko-young-1968

that has been applied in studies of a wide range of applications. These studies include the response of a cantilever to a uniformly distributed or varying load

timoshenko-young-1968 or to a surrounding fluid flow kamrin-rycroft-nave-2012 ; pozrikidis-2010 ; pozrikidis-2011 ; farjoun-schaeffer-2011 ; li-schulgasser-cederbaum-1995 . For example, Pozrikidis pozrikidis-2010

studied shear flow over a periodic array of cylindrical rods attached to a substrate in order to determine the macroscopic slip velocity. He also computed the hydrodynamics load exerted along the rod as well as estimating the flow-induced deflection. Shortly afterwards, Pozrikidis 

pozrikidis-2011 studied shear flow past an elastic rod attached to a plane surface (as well as a doubly-periodic array of such rods), motivated by the study of biological flows involving ciliated surfaces. His mathematical framework combined slender-body theory for computing the hydrodynamic load with classical beam theory for computing the rod dynamics. Small vibrations of a thin flexible beam immersed in a laminar flow was studied by Farjoun and Schaeffer farjoun-schaeffer-2011 , where the authors assumed that the dominant restoring force was due to tension caused by shear stress. Goza and Colonius goza-colonius-2017 presented a strongly coupled immersed boundary method for FSI problems involving thin deforming bodies. The method was found to be stable for arbitrary choices of solid to fluid mass ratios.

Because we are motivated by the study of problems involving deforming biofilm layers, where the effects of porosity can be important, we will also be investigating the behaviour of porous flexible cantilevers in response to fluid flow. The closest work we have been able to identify regarding porous beams is a mathematical model for an incompressible poro-elastic beam yang-wang-2007 in which the internal porous fluid flow is confined to the axial direction of the deformed beam. Other related studies of poroelastic beams can be found in li-schulgasser-cederbaum-1995 ; Li-Schulgasser-Cederbaum-1998 ; yang-wang-2007 ; yang-wen-2010 , where again the beam was restricted to be permeable in the axial direction only.

Our modelling approach is based on the immersed boundary or IB method, which has been used in the study of problems in fluid-structure interaction problems from biology in particular, but also in a variety of engineering and other applications (see fauci-dillon-2006 ; peskin-2002 and references therein). The problem under consideration here involves fluid interacting with rigid 1D structures (channel walls) as well as a 2D deformable region (a beam, which can be either porous or solid). Hence, our problem has a mixture of different solid geometries for which the IB method is an ideal choice peskin-1977 ; peskin-2002 . Some researchers have applied the IB method to simulate the dynamics of cilia which individually bear a striking resemblance to a flexible cantilever beam (e.g., dillon-fauci-2000 ; dillon-etal-2007 ; schwartz-etal-1997 ). The issue of porous elastic boundaries has also been addressed in the IB context, for example by Kim and Peskin kim-peskin-2006 who were the first to incorporate porosity within the IB framework in a study of parachute dynamics. Their approach to handling porous air vents at the apex of the chute was to allow allow the normal velocity of the canopy to differ from that of the fluid by an amount proportional to the normal component of the boundary force. Stockie stockie-2009 followed Kim and Peskin’s approach by incorporating porosity directly through Darcy’s law. We extended the ideas in kim-peskin-2006 ; stockie-2009 for a 1D porous membrane with pores directed normal to the surface, to handle the case of solid porous 2D region. We then apply this approach to study deformation of a porous rectangular cantilever beam with no restriction on the direction of the pores. For the general situation we consider here, we have not been able to find any other numerical or experimental results with which we can draw direct comparisons.

The organisation of this paper is as follows. We start in Section 2 by defining the problem geometry and stating the governing equations and numerical scheme. We also define the IB force density that encompasses the influence of the cantilever beam and bounding walls on the flow. In Section 3 we perform numerical simulations of an elastic, cantilever beam and compare the results to experimental and analytical results available in the literature. Finally, in Section 4, we modify the governing equations to introduce porosity into the beam and compare the results from numerical simulations to the non-porous case.

2 Immersed boundary method

The term “immersed boundary method” is commonly used in the literature to refer not only to a numerical method but also to the underlying mathematical formulation peskin-2002 . In this section, we state the governing equations that capture the interaction between a deformable elastic structure and a surrounding fluid. Then we apply a straightforward finite difference discretization and present the resulting numerical algorithm that we employ here to compute approximate solutions to the IB problem.

2.1 Governing equations

The IB formulation consists of a coupled system of nonlinear integro-differential equations that describe both the force generated by a solid, deformable, elastic body and the dynamics of a surrounding incompressible, Newtonian fluid. The method is capable of handling immersed boundaries with a very general shape and configuration, although we restrict ourselves here to the 2D geometry pictured in Figure 1 wherein an initially rectangular cantilever beam is placed inside a fluid-filled channel. The channel is embedded inside a larger, rectangular fluid domain and is represented by two horizontal walls. For simplicity, we impose periodic boundary conditions in both the - and -directions of the fluid domain, and the top and bottom channel walls are horizontal lines separated a small distance from the top and bottom domain boundaries respectively. The lower channel wall is fixed in space along , while the upper wall lies along and moves to the right with a given constant horizontal speed . In the absence of obstructions in the channel, this setup would generate a horizontal linear shear flow with shear rate .

The deformable cantilever beam is initially a rectangle with thickness and length , and has its lower end fixed to the center of the bottom wall. Note that the channel walls and are taken to be idealized one-dimensional structures of zero thickness, while the beam is treated as a solid elastic structure, with a well-defined thickness. The details regarding the specification of the geometry and forces on the walls and beam are delayed until Section 2.3.

Figure 1: A general elastic boundary that consists of several possibly disconnected components is immersed within a doubly-periodic rectangular fluid domain .

In general, we consider an immersed boundary that consists of several disconnected components, , where each is either a one-dimensional elastic membrane (parameterized by a single real dimensionless parameter ) or an elastic solid region (whose specification requires two parameters, and ). We denote the location of the immersed boundary by (with units of []) where is the dimensionless parameterization of points on and represents either a scalar

or a vector

depending on the context. The parameters and are chosen to lie in the interval for simplicity.

The influence of the elastic body on the fluid is to impose a force [] onto the adjacent fluid particle at location , which is incorporated into the incompressible Navier-Stokes equations as a body force term:


Here, is the fluid velocity [], is the pressure [], are the Eulerian coordinates [] for the fluid domain , is density [] and is dynamic viscosity []. The IB force term in the momentum equations (1) is actually specified using an elastic force density [] that is spread onto the surrounding fluid by means of a delta-function convolution


where represents the Cartesian product of two one-dimensional Dirac delta functions. Finally, is required to move with the fluid so that


which is simply the “no-slip” condition for fluid particles located adjacent to the immersed boundary.

In summary, the governing equations consist of (1)–(4), with the IB force density being the only component that remains to be specified. Since it is simplest to write in discrete form, we will first derive the discretized version of the governing equations, after which we will provide a detailed specification for the IB force.

2.2 Numerical algorithm

The algorithm we present next is a semi-implicit scheme that is essentially identical to the method used in stockie-phdthesis-1997 and closely related to the “original” IB method of Peskin peskin-1977 . The fluid domain is divided into an equally-spaced Cartesian grid of points denoted , where , , , and . We also divide the time interval into equally-spaced points with time step and . We may then define discrete approximations of the velocity and pressure variables, and at points . The immersed boundary parameter for each component is likewise taken at equal intervals for (and similarly for , if necessary) and the discrete IB configuration and force density are denoted and respectively.

Using the above notation, we can then introduce finite difference operators that approximate the spatial derivatives appearing in the governing equations. In particular, we define one-sided difference approximations for the –derivative of a grid quantity ,


as well as the centered approximation


Analogous definitions apply for the –derivatives , and , after which the gradient may be replaced by the centered approximation . Finally, the delta functions appearing in the integral terms are replaced by the regularized function


We are now prepared to describe our implementation of the immersed boundary algorithm. In any given time step, we assume that values of the velocity and IB configuration are known from the previous step at time . These quantities are evolved in time to using the following procedure:

  1. Compute the force density based on the IB configuration as described in Section 2.3.

  2. Spread the IB force onto fluid grid points using a discretized version of the integral from (3)


    The geometric scaling factor is inversely proportional to the number of IB points () and has a different definition depending on whether the immersed boundary is a 1D curve (channel wall) or a 2D solid region (cantilever beam). In the case of a wall, is a length, while for a 2D region is an area; in both cases, the factor ensures that the formula (9) scales properly with the number of IB points and that it is a consistent approximation of the corresponding single or double integral. More details on the precise form of (9) and the specification of are provided in Section 2.3.

  3. Integrate the incompressible Navier-Stokes equations using a split-step projection scheme:

    1. Compute an intermediate velocity by applying the elastic force from the IB:

    2. Apply an ADI discretization of the advection and diffusion terms, which splits the spatial derivatives in the and directions:


      where we have linearized the discrete equations by freezing coefficients at the previous time step in the advection terms. As a result, (11)–(12) reduce to a sequence of linear tridiagonal solves for and .

    3. Project the intermediate velocity onto the space of divergence-free vector fields by:

      1. Solving the following Poisson equation for pressure


        Note that represents a wide finite difference stencil for the Laplacian involving the pressure values , , , and . Owing to the periodic boundary conditions on

        , this equation is solved most easily and efficiently by means of the discrete Fourier transform, and making use of the Fast Fourier Transform algorithm.

      2. Updating the velocity via

  4. Evolve the IB in time using the explicit update


This simple semi-implicit time discretization yields a solution that is first-order accurate in time. And although all spatial derivatives are approximated using second-order finite differences, the method is still only first-order accurate in space owing to errors in velocity interpolation near the IB that arise from the use of the regularized delta function 

liu-mori-2012 . It is straightforward to increase the temporal accuracy to second order using an algorithm such as that proposed by Lai and Peskin lai-peskin-2000 , but it is much more difficult to increase the spatial accuracy to second order griffith-etal-2007 . Since the primary aim of the current study is to validate the IB approach for the study of sheared cantilever beams, we have chosen to employ the simple scheme above and leave for future work the implementation of higher order extensions to the algorithm.

2.3 Discrete specification of the IB force density

We next specify the elastic IB forces that are generated by the channel walls and cantilever beam pictured in Figure 1. Keep in mind that because of the two-dimensional geometry, the equivalent three-dimensional flow can be envisioned as extending to infinity in both directions perpendicular to the –plane, and hence the walls are actually horizontal planar surfaces whereas the beam behaves as an infinite cantilever plate. We separate the IB force density into the sum of three terms, , where represents the force generated by the channel walls, is that generated by the beam, and derives from the attachment force between the wall and the bottom edge of the beam. Each of these forces is developed separately in the following three sections.

2.3.1 Force density for channel walls,

Each horizontal wall is discretized using a sequence of IB points that are equally spaced in the fiber parameter . For example, along the stationary bottom wall we define the initial wall point positions by , where and . In the IB framework, the wall is actually permitted to deviate slightly from its target configuration by connecting each wall point to a fixed “tether point” (initially at the same location) using a very stiff spring that exerts a force of the form


where [] is the wall spring stiffness and is the moving IB point location. Any deviation of the wall point from corresponding the tether point location generates a force that drives it back towards the tether point, so that as long as is chosen large enough the wall points can be made to mimic a rigid structure. We emphasize that the tether points do not generate any force themselves but only serve to determine target locations for the moving IB points.

The given horizontal velocity of the top wall points is easily incorporated in the above framework by simply specifying a given motion for the tether points, so that the IB force density becomes


where and the “modulo ” operation enforces the periodic boundary conditions in . The total force density generated by the two channel walls may then be written as


The natural choice of scaling factor for the force spreading step (9) is the wall point spacing, .

2.3.2 Force density for cantilever beam,

The cantilever beam is discretized using a collection of Lagrangian points that lie on its circumference as well as distributed throughout its interior. We employ the unstructured triangular mesh generator DistMesh persson-strang-2004 (implemented in Matlab) to generate an approximately uniform triangulation for the beam such as that depicted in Figure 2a. The nodes of the triangulation are the IB points , for , while the edges of the triangles define a network of springs that acts to maintain the shape of the beam.

(a) Rectangular beam (b) Smoothed beam
(731 points, 1304 triangles) (743 points, 1320 triangles)
Figure 2: Unstructured triangular mesh for the two beam shapes, generated by distmesh2d. Both have the same height and width, and the fillets for the smoothed beam are constructed in such a way that the areas of the two shapes are equal.

We generate the triangulation using the function distmesh2d with the “scaled edge length function” huniform, which is a built-in function that aims to find a mesh that is as uniform as possible. We also set the “initial edge length” parameter equal to , which ensures that in practice the mesh satisfies the following constraint


and hence prevents leakage of fluid through the IB that can occur if IB points become too widely separated peskin-2002 . We will also perform numerical simulations on the “smoothed beam” shape shown in Figure 2b, where the rectangular corners are replaced by fillets or chamfers – the reasons for the choice of this smoothed shape will be explained later in Section 3.4.

We next define the spring forces acting on the network, which we base on the model of Alpkvist and Klapper for viscoelastic biofilm structures alpkvist-klapper-2007 . Let be the vector joining two IB points labeled and , and let be the corresponding distance. We assume that the spring network is initially in equilibrium (i.e., all spring forces are in balance) so that each spring has a resting length equal to its initial length . Let represent an incidence matrix whose entries are either 1 or 0 depending on whether or not points and are connected, respectively. Then the force density acting on the IB point in the network is


where the sum is taken only over those for which is connected to in the network. We have also assumed that the spring stiffness [] is constant for all network connections. The total elastic force density generated by all IB points making up the beam is then given by


The appropriate scaling factor for the force integral (9) is the average area of a triangular grid cell, , where we have employed the asymptotic estimate that for large the number of triangles is approximately to leading order.

2.3.3 Force density for attachment between beam and bottom wall,

In order to hold the lower end of the cantilever beam stationary and coincident with the channel wall, we impose an additional attachment force that connects the points at the base of the beam to the fixed tether points on the bottom wall. The wall and beam points are initially collocated and each pair of points is connected by a very stiff spring. Denote by the attachment point index set which consists of all pairs of IB point indices for which corresponds to the bottom wall tether point that is linked to the beam IB point . Then the force density arising from the spring joining these two points is


where is the vector connecting the wall tether point with fixed location and the beam point at . Note that the springs in this case have a zero resting length, so that the IB forces act to keep the point pairs in the same location. The spring constant is chosen equal to that of the wall–tether connections so that . The total wall–beam attachment force can then be written as


Since this spring force has a zero resting length, there is no obvious length scale for teh attachment forces, and so in this case we simply choose the wall point distance equal to the scaling factor for the force spreading calculation from (9), .

3 Numerical simulations of a solid beam

Our numerical simulations presented in this section consist of a parametric study in which we vary values of the beam spring stiffness , beam length , and top wall velocity , the latter of which is related to shear rate via . We start with a “base case” corresponding to a beam having dimensions and . The fluid domain is a square of size , and the channel walls are located a distance from the boundaries, which means that the channel width is . This choice of geometry is motivated by one of the test cases considered by Alpkvist and Klapper in alpkvist-klapper-2007 , and yields a beam aspect ratio that is large enough to satisfy the thin-beam assumption in the linear theory. This beam length is also small enough (roughly one-third of the channel width) that it doesn’t significantly hinder the bulk flow through the channel. The base value for the top wall velocity is , while the spring stiffness values are and . For the numerical discretization, we take a fluid grid with grid points in each direction, and choose and in order to satisfy the constraint (19) on the spacing between IB points for the beam and walls. The base case values of all parameters are summarized in Table 1.

As parameters are varied from the base case, we focus on the deflection of the beam from its vertical stress-free state and also the changes in the flow structure. In all cases, the fluid velocity is initialized to zero and the speed of the top wall is ramped up linearly from 0 to over the time interval . As the fluid within the channel accelerates and begins to move toward the right, the beam responds by bending downward and to the right as pictured in Figure 3. Because the beam is solid, the flow deflects around the beam, generating a fairly complex flow structure that features recirculation zones in the area near the bottom wall (refer to the streamline patterns in Figure 3a). The beam continues to deform in response to the flow until an equilibrium is attained in which the fluid force acting on the beam and the elastic bending force from the deformed beam are equal. A typical computation requires on the order of 5 to 10  for the beam to reach its steady state.

(a) (b)
Figure 3: Flow streamlines (left) and pressure contours (right) for the base case simulation at steady state.

We remark that during the course of these rectangular beam simulations, especially for cases with large deflection, we sometimes observe large non-physical irregularities in the shape of the beam in the vicinity of the free end and the attachment at the lower wall especially at high shear rates. For this reason, we also perform a series of simulations with a smoothed beam shape for which the corners at the free end and the beam–wall anchor points are rounded using a fillet or chamfer (refer to Fig. 2b). These results are presented in Section 3.5, and we will see that they exhibit more reasonable deformations near the corners.

Parameter Description Value Units
Domain width 0.03
Domain height 0.03
Wall separation distance 0.0033
Beam length 0.0077
Beam width 0.0014
Beam spring stiffness 560
Wall spring stiffness 1000
Attachment spring stiffness 1000
Top wall velocity 0.02
Water density 1.0
Water viscosity 0.01
Re Reynolds number 0.0154
, Number of fluid grid points 64
Number of beam IB points 731
Number of wall IB points 210
Time step
Table 1: “Base case” parameter values for the solid elastic beam problem.

3.1 Euler-Bernoulli beam theory

Before reporting the results of our numerical simulations, we briefly summarize some key results from the Euler-Bernoulli beam theory that pertain to small deflections of a thin cantilever beam. In this case the steady-state configuration of the beam can be desribed by a fourth-order ordinary differential equation 



where is the horizontal deflection from the vertical equilibrium state at location along the beam, and is a given load distribution (in units of force per unit length, or ). The parameter is the flexural rigidity of the beam, which is the product of the Young’s modulus [

] and the second moment of area

[]. The linear theory is derived based on the assumptions that the beam is very thin () and that the deflection is small relative to the length of the beam ().

In particular, we list well-known analytical solutions for the maximum (tip) deflection of the beam in response to two types of applied load: (a) a constant load [] distributed across the entire length of the beam : and (b) a linearly-varying load that decreases from a maximum at the wall to zero at the tip : . The tip deflection in all two cases varies according to , with exponent . These expressions will be compared qualitatively to the simulated deflections in Section 3.2.

3.2 Dependence of deflection on IB spring stiffness,

We begin by investigating the effect of changing the spring stiffness for links in the beam triangulation, which in turn determines the effective flexural rigidity of the overall beam structure. We note that the beam dimensions and corresponding to an aspect ratio of , which is consistent with the thin beam assumption in the linear theory.

Figure 4a depicts the time evolution of the beam tip deflection for values of , where deflection is measured in the -direction from the vertical equilibrium configuration. The steady-state deflection values vary from roughly 10% of the beam length for the stiffest beam up to 65% for the smallest value of . It seems reasonable therefore to expect that the larger simulations will fall within the small deflection regime, while the smaller results might not. Nonetheless, all curves exhibit a similar qualitative behaviour in that the shear flow bends the beam away toward the right and then gradually equilibrates at some maximum deflection. As expected, the stiffer beams have a smaller deflection and also reach their equilibrium state over a shorter time period.

(a) Deflection versus time.

(b) Final deflection versus .

Figure 4: Deflection of the solid cantilever beam as is varied, with and . (a) Tip deflection versus time corresponding to . (b) The final (steady-state) tip deflection is plotted against with diamond points () on top of which is superimposed a solid line corresponding to the linear least squares fit (normalized RMS error 0.055).

Figure 4b plots the steady-state deflection for all simulations versus , along with a number of additional results for in the interval that fill out the range. The data points are approximated very well by a straight line on a log–log scale; indeed, a least squares fit yields a line with with slope which suggests that the deflection is roughly inversely proportional to . We can compare these results to the linear beam theory in Section 3.1 which predicts that the deflection is inversely proportional to the flexural rigidity, , where is Young’s modulus, is the area moment of inertia. Assuming that the linear theory holds here (thin beam, small deflection) our computational results suggest that the effective flexural rigidity of our triangulated spring link structure is directly proportional to the spring stiffness, or in other words that .

3.3 Dependence of deflection on beam length,

Next we perform simulations with different values of the beam length , 0.007, 0.0077, 0.0084, 0.0098, 0.014 (corresponding to aspect ratios of , 5, 5.5, 6, 7, 10 respectively) while holding and constant. The computed results in Figure 5 show that the tip deflection increases with , which is to be expected since a longer beam is clearly more flexible. Furthermore, the dependence of on is again roughly linear on a log-log scale which denotes a power law relationship, and a least squares fit suggests that . This power law exponent is very close to the value 4 predicted by the linear beam theory for the cases when the load force is constant or linearly-varying along the length.

Figure 5: Plot of and final deflection for the non-porous cantilever. Parameters: , . Computed results are shown using diamond points () on top of which is superimposed a solid line corresponding to the linear least squares fit (normalized RMS error 0.0432).

3.4 Dependence of deflection on shear velocity,

In our final sensitivity study, we vary the hydrodynamic force acting on the beam by changing the top wall velocity . In contrast with the previous cases, we consider much stronger shear flows corresponding to values taken from the interval , which at the upper end generates beam deformations that are well outside the linear regime. Figure 6 depicts the deformed beam and corresponding streamlines for values of , computed over long enough times (s) that the beam has essentially reached a steady state configuration.

(c) (base case)
Figure 6: Deflected beam for several values of , showing streamlines (left) and contours of velocity magnitude (right). Other parameter values are , .

If we consider the value of the last streamline lying above the beam tip (from Figure 6), which yields the following values of :

0.0005 0.003 0.17
0.001 0.01 0.1
0.002 0.02 0.1
0.005 0.04 0.125

Consider the ratio of which is approximately 0.1 for all cases, suggesting a clear separation into a boundary layer region where the velocity drops roughly below 10% of the top wall speed.

In each case, the flow separates roughly into two regions: a slow inner flow near the wall that “stalls” ahead of and behind the beam; and a much faster outer shear flow that bypasses the beam. As the outer shear flow increases in strength for higher , the flow is characterized by recirculating eddies ahead of and behind the beam. We remark also that streamlines pass through the beam and the wall, which at first seems counter-intuitive if the immersed boundaries (and the adjacent fluid) are supposed to be stationary at equilibrium. However, it is important to keep in mind that the IB points making up both the wall and beam experience very small-amplitude oscillations even when the flow is near steady state, and so there are still small non-zero velocities in the neighbourhood of the immersed boundaries.

Figure 7: Maximum tip deflection as a function of for the solid beam, with parameters and .

The results of these simulations, as well as for a number of other flow speeds, are summarized in Figure 7 on a plot of tip deflection versus top wall velocity, which clearly shows that the beam deformation can be divided into three regimes according to :

  • For or , the beam deflection is small enough that the linear beam theory applies and the load force varies roughly linearly with .

  • For or , the beam has reached the maximum allowable deflection beyond which any additional increase in has little effect on the shape of the beam at steady-state. This behaviour can be attributed to the fact that at larger values of the shear rate, the beam undergoes a significant vertical deflection as it bends downward and to the right, to the point where a large fraction of the beam is oriented parallel to the bulk flow. This in turn reduces the net fluid load force acting on the beam to the point where increasing any further does not change the deflection (refer to Figure 6d).

  • For intermediate values of or , the beam is in a transitional state between the low and high shear regimes described above.

Another feature of the beam deformations that becomes evident at higher is a large deviation from the original rectangular shape for the IB points near either end of the beam, most notably at the corners. This behaviour is emphasized in Figure 8 which shows two zoomed-in plots of the corner regions near the upper and lower ends of the beam. Near the wall (in the area denoted as region A) large bending forces are generated that cause the beam to “bow out” on both sides. Furthermore, there are noticeable “kinks” at both bottom corners involving the left- and right-most IB points immediately adjacent to the lower wall. At the free end of the beam (zoomed in as region B) the left corner undergoes a significant distortion and the corner point protrudes a significant distance out into the flow owing to the very large shear forces experienced in that region. Both of these behaviours seem non-physical and would not be expected in an actual beam with uniform elastic properties. We attribute the anomalous deformations near both ends of the beam to the existence of sharp corners in the rectangular shape, and the purpose of the next section is to investigate a modified beam shape that aims to eliminate these anomalies.

  B    A
Figure 8: Zoomed view of the large deformations near the corners of the beam when .

3.5 Smoothed cantilever beam

We next consider a modified beam shape that has the corners at both ends smoothed out as shown in Figure 2b using using fillets or chamfers. In particular, all four corners have been replaced using circular arcs that have diameter equal to the beam width . The resulting shape is triangulated as before using distmesh2d, and we choose parameter values that are identical to the simulation from the previous section. The resulting equilibrium beam configuration is depicted in Figure 9, which exhibits a shape that is much more regular than in the rectangular case, even for this high value of the fluid velocity.

  B    A
Figure 9: Zoomed view of the two ends of the smoothed beam when .

Although the change in beam shape due to smoothing out of the corners does influence the solution, these changes are very small as seen in the two plots of tip deflection versus and in Figure 10. Indeed, the relative differences between the two solutions is less than 1%. We will therefore employ the smoothed beam shape in all simulations from this point onward. And with a view to applications, the rounded beam shape also has the advantage that it is much more realistic in the context of biofilm problems.

Figure 10: Maximum tip deflection for the rectangular and smooth shapes, plotted as a function of (a) top wall velocity and (b) beam spring stiffness .

4 Numerical simulations of a porous beam

In this section, we consider a porous deformable structure, so that fluid can flow through the structure in response to transmural pressure gradients. Porous structures abound in biological systems, including such examples as artery walls huang-tarbell-1997 and porous biofilm layers vanloosdrecht-etal-2002 ; nguyen-morgenroth-eberl-2005 ; thullner-baveye-2007 ; davit-etal-2013 . Porous immersed boundaries also appear in engineering applications such as filtration and separation processes.

Several efforts have been made to generalize the IB method to handle porous structures. The first attempt to incorporate porosity within the IB framework was a study of parachute dynamics by Kim and Peskin kim-peskin-2006 , wherein the air vents at the apex on the chute were dealt with by allowing the normal velocity of the canopy to differ from that of the fluid by an amount proportional to the normal component of the boundary force (which according to the jump conditions is also proportional to the pressure jump). Stockie stockie-2009 built on Kim and Peskin’s approach by incorporating porosity directly using Darcy’s law. The IB method has also been used to study flow through granular media at the pore scale by treating the grains making up the medium as immersed boundaries dillon-fauci-2000 , although the grains themselves are rigid and impermeable in these studies. Layton layton-2006 generalized the closely related immersed interface method by introducing a porous slip velocity in the direction normal to the interface that is driven by differences in both transmural water pressure and solute concentration. In this thesis, we will extend the ideas in kim-peskin-2006 ; stockie-2009 for a 1D porous membrane to the case of a solid porous region in 2D. We will use this approach to study the deformation of a porous rectangular cantilever beam, as well as the same smoothed cantilever shape from the previous section.

4.1 Extending the IB formulation to include porosity

The system of equations (1)–(4) is now generalized to include the effect of porosity on the deformable structure. We follow the approaches of Kim and Peskin kim-peskin-2006 and Stockie stockie-2009 by incorporating a porous slip velocity in the fiber evolution equation. The primary assumption made in stockie-2009 and kim-peskin-2006 is that all the pores are directed normal to the fiber. The slip velocity is given by where and are the unit normal and tangential vectors to the fiber. Due to the assumption just mentioned, the tangential component . For a porous 2D region, it is not realistic to continue with their assumption.

Indeed, we assume a more general porous structure with no restriction in the direction of the pores and fluid can flow in any direction in the porous region. The only way to determine the direction of the fluid velocity is to look at the direction of the pressure gradient, , that drives the porous flow. In our approach, we also incorporate porous effects via a porous slip velocity ; however here is related to the pressure gradient (and not the pressure jump) via Darcy’s law:


where represents the structure permeability [].

The porous beam obeys the same governing equations except that the fiber evolution equation (4) is replaced with


The negative sign indicates that when the fluid moves through the porous region outward with slip velocity , the beam in turn moves in the opposite direction, .

In our simulations, we have fixed values of , and , which corresponds to a shear rate of . Permeability lies in the range and Figure 11 shows the streamline plots for different values of . It is evident from the figure that for high value of , the fluid passes through the beam without any obstruction but that is not the case for smaller values of . When is taken as small as , the flow is almost identical to that for the beam in the non-porous case (). Figure 12 shows how the maximum deflection experienced by the cantilever varies with permeability, from which we observe the following:

  • For : the beam is essentially solid and the deflection is the same as in the non-porous case.

  • For : there is a slight increase in maximum deflection. We ascribe this to an increase in porous flow through the beam that acts to enhance the horizontal component of velocity in the region near the beam which in turn increases the shearing force acting on the beam.

  • For : the porous slip velocity is so large that fluid can pass through the solid structure more freely, hence decreasing the forces acting to deform the beam.

Figure 11: Streamline plots at the point of maximum deflection for a smoothed porous beam with different permeability values.
Figure 12: Final deflection () versus permeability () for the smoothed porous beam.

5 Conclusions

The main aim of this study was to derive a mathematical model and the corresponding computational scheme to study the deflection of a two-dimensional deformable elastic cantilever beam in response to a surrounding viscous, incompressible fluid flow. The choice of cantilever was motivated by our efforts to model biofilm structures in the near future. Our study included careful validation of the IB approach so that in the future the model can be easily generalized to irregular and highly deformable 3D structures with porosity, detachment etc. We investigated how variations in physical and numerical parameters change the effective material properties of the elastic beam and also made a qualitative comparison of the results with linear beam theory. We also paid attention to “corner effects” (irregularities in beam shape near free end and near wall connection points) and showed how this can be remedied by smoothing out the corners with a filled or round (a standard technique in solid modelling to reduce stress, and in aerodynamics to reduce interference drag). Finally we extended the previous work done on porous IB membranes, by introducing porosity into our deformable elastic beam. In our study we have not imposed any restrictions on the orientation of the pores. The results obtained are consistent with physical intuition. For example, physical intuition suggests that for higher values of permeability, the fluid should pass through the beam without obstruction as compared to smaller values of permeability. The streamline plots obtained in this study also suggests the same. Along the same lines of argument we can say that as permeability of the beam increases, the deflection of the beam should decrease as the force acting to deform the beam decreases and this argument is in good agreement with our numerical results.

In future, we aim to investigate the relationship by studying two analytical approximations. Our first approach will be by finding a polynomial fit to “typical” computed load, and then deriving the corresponding Euler-Bernoulli beam solution. A second method will apply a simpler spring network and attempt to derive the equilibrium configuration analytically. We will also aim for irregular, deformable shapes corresponding to biofilm columns and streamers. Last but not least, there is a limitation in assuming a 2D geometry and so one of our future aims will be to extend our IB model to 3D cantilevers in the shape of deformed cylinders or parallelipipeds.


  • (1) E. Alpkvist, I. Klapper, Description of mechanical response including detachment using a novel particle model of biofilm/flow interaction, Water Science and Technology 55 (2007) 265–273.
  • (2) S. Timoshenko, D. Young, Elements of Strength of Materials, D. Van Nostrand, New York, 1968.
  • (3) K. Kamrin, C. H. Rycroft, J.-C. Nave, Reference map technique for finite-strain elasticity and fluid-solid interaction, Journal of the Mechanics and Physics of Solids 60 (2012) 1952–1969.
  • (4) C. Pozrikidis, Shear flow over cylindrical rods attached to a substrate, Journal of Fluids and Structures 26 (2010) 393–405.
  • (5) C. Pozrikidis, Shear flow past slender elastic rods attached to a plane, International Journal of Solids and Structures 48 (2011) 137–143.
  • (6) Y. Farjoun, D. G. Schaeffer, A thin cantilever beam in a flow, in: Numerical Analysis and Applied Mathematics ICNAAM 2011, Vol. 1389 of AIP Conference Proceedings, 2011, pp. 1446–1449. doi:10.1063/1.3637895.
  • (7) L. P. Li, K. Schulgasser, G. Cederbaum, Theory of poroelastic beams with axial diffusion, Journal of the Mechanics and Physics of Solids 43 (12) (1995) 2023–2042.
  • (8) A. Goza, T. Colonius, A strongly-coupled immersed-boundary formulation for thin elastic structures, Journal of Computational Physics 336 (2017) 401–411. doi:10.1016/
  • (9) X. Yang, C. Wang, A nonlinear mathematical model for large deflection of incompressible saturated poroelastic beams, Applied Mathematics and Mechanics (English Edition) 28 (12) (2007) 1587–1595.
  • (10) L. P. Li, K. Schulgasser, G. Cederbaum, Large deflection analysis of poroelastic beams, International Journal of Nonlinear Mechanics 33 (1) (1998) 1–14.
  • (11) X. Yang, Q. Wen, Dynamic and quasi-static bending of saturated poroelastic Timoshenko cantilever beam, Applied Mathematics and Mechanics (English Edition) 31 (8) (2010) 995–1008. doi:10.1007/s10483-010-1335-6.
  • (12) L. J. Fauci, R. Dillon, Biofluidmechanics of reproduction, Annual Review of Fluid Mechanics 38 (2006) 371–394. doi:10.1146/annurev.fluid.37.061903.175725.
  • (13) C. S. Peskin, The immersed boundary method, Acta Numerica 11 (2002) 1–39.
  • (14) C. S. Peskin, Numerical analysis of blood flow in the heart, Journal of Computational Physics 25 (1977) 220–252.
  • (15) R. Dillon, L. Fauci, Microscale model of bacterial and biofilm dynamics in porous media, Biotechnology and Bioengineering 68 (5) (2000) 536–547.
  • (16) R. H. Dillon, L. J. Fauci, C. Omoto, X. Yang, Fluid dynamic models of flagellar and ciliary beating, Annals of the New York Academy of Science 1101 (2007) 494–505.
  • (17) E. A. Schwartz, M. L. Leonard, R. Bizios, S. S. Bowser, Analysis and modeling of the primary cilium bending response to fluid shear, American Journal of Physiology: Renal Physiology 272 (1997) F132–F138.
  • (18) Y. Kim, C. S. Peskin, 2-D parachute simulation by the immersed boundary method, SIAM Journal on Scientific Computing 28 (6) (2006) 2294–2312.
  • (19) J. M. Stockie, Modelling and simulation of porous immersed boundaries, Computers and Structures 87 (11-12) (2009) 701–709. doi:10.1016/j.compstruc.2008.11.001.
  • (20) J. M. Stockie, Analysis and computation of immersed boundaries, with application to pulp fibres, Ph.D. thesis, Institute of Applied Mathematics, University of British Columbia, Vancouver, Canada, available at (1997).
  • (21) Y. Liu, Y. Mori, Properties of discrete delta functions and local convergence of the immersed boundary method, SIAM Journal on Numerical Analysis 50 (6) (2012) 2986–3015. doi:10.1137/110836699.
  • (22) M.-C. Lai, C. S. Peskin, An immersed boundary method with formal second-order accuracy and reduced numerical viscosity, Journal of Computational Physics 160 (2) (2000) 705–719.
  • (23) B. E. Griffith, R. D. Hornung, D. M. McQueen, C. S. Peskin, An adaptive, formally second order accurate version of the immersed boundary method, Journal of Computational Physics 223 (1) (2007) 10–49.
  • (24) P. O. Persson, G. Strang, A simple mesh generator in MATLAB, SIAM Review 46 (2) (2004) 329–345.
  • (25) Z. J. Huang, J. M. Tarbell, Numerical simulation of mass transfer in porous media of blood vessel walls, American Journal of Physiology: Heart and Circulatory Physiology 273 (1) (1997) 464–477.
  • (26) M. C. M. van Loosdrecht, J. J. Heijnen, H. Eberl, J. Kreft, C. Picioreanu, Mathematical modelling of biofilm structures, Antonie van Leeuwenhoek 81 (2002) 245–256.
  • (27) V. T. Nguyen, E. Morgenroth, H. J. Eberl, A mesoscale model for hydrodynamics in biofilms that takes microscopic flow effects into account, Water Science and Technology 52 (7) (2005) 167–172.
  • (28) M. Thullner, P. Baveye, Computational pore network modeling of the influence of biofilm permeability on bioclogging in porous media, Biotechnology and Bioengineering 99 (6) (2007) 1337–1351. doi:10.1002/bit.21708.
  • (29) Y. Davit, H. Byrne, J. Osborne, J. Pitt-Francis, D. Gavaghan, M. Quintard, Hydrodynamic dispersion within porous biofilms, Physical Review E 87 (2013) 012718. doi:10.1103/PhysRevE.87.012718.
  • (30) A. T. Layton, Modeling water transport across elastic boundaries using an explicit jump method, SIAM Journal on Scientific Computing 28 (6) (2006) 2189–2207.