COMSOL 6.3 - Fluid–Structure Interaction

The following example demonstrates techniques for modeling fluid–structure interactions in COMSOL Multiphysics. It illustrates how fluid flow can deform structures and how to solve for the flow in a continuously deforming geometry using the arbitrary Lagrangian–Eulerian (ALE) technique.

The model geometry consists of a horizontal flow channel in the middle of which is an obstacle, a narrow vertical structure (Figure 1). The fluid flows from left to right, except where the obstacle forces it into a narrow path in the upper part of the channel, and it imposes a force on the structure’s walls resulting from the viscous drag and fluid pressure. The structure, being made of a deformable material, bends under the applied load. Consequently, the fluid flow also follows a new path, so solving the flow in the original geometry would generate incorrect results.

The ALE method handles the dynamics of the deforming geometry and the moving boundaries with a moving grid. COMSOL Multiphysics computes new mesh coordinates on the channel area based on the movement of the structure’s boundaries and mesh smoothing. The Navier–Stokes equations that solve the flow are formulated for these moving coordinates.

The structural mechanics portion of the model does not require the ALE method, and COMSOL Multiphysics solves it in a fixed coordinate system as usual. However, the strains the model computes in this way are the only source for computing the deformed coordinates with ALE.

Figure 1: Fluid flows into this horizontal flow channel from the left, and it enters with a parabolic velocity profile. A narrow vertical structure in the channel (the straight vertical structure) forces the flow into a narrow path. Due to fluid pressure and viscous drag, the originally vertical structure bends. This simulation models the fluid flow in a deformed, moving mesh that follows the movement of the bending structure.

Model Definition

In this example the flow channel is 100 μm high and 300 μm long. The vertical structure — 5 μm wide, 47.5 μm high, and with a semicircular top — sits 100 μm away from the channel’s left boundary. Assume that the structure is long in the direction perpendicular to the image.

The fluid is a water-like substance with a density ρ = 1000 kg/m3 and dynamic viscosity η = 0.001 Pa·s. To demonstrate the desired techniques, assume the structure consists of a flexible material with a density ρ = 7850 kg/m3 and Young’s modulus E = 200 kPa.

Fluid Flow

The fluid flow in the channel is described by the incompressible Navier–Stokes equations for the velocity field, u = (u, v), and the pressure, p, in the spatial (deformed) moving coordinate system:

In these equations, I denotes the unit diagonal matrix and F is the volume force affecting the fluid. Assume that no gravitation or other volume forces affect the fluid, so that F = 0. The coordinate system velocity is um = (um, vm).

At the channel entrance on the left, the flow has fully developed laminar characteristics with a parabolic velocity profile but its amplitude changes with time. At first flow increases rapidly, reaching its peak value at 0.215 s; thereafter it gradually decreases to a steady-state value of 5 cm/s. The centerline velocity in the x direction, uin (see Figure 4), with the steady-state amplitude U comes from the equation

where t must be expressed in seconds.

At the outflow (right-hand boundary), the condition is p = 0. On the solid (nondeforming) walls, no slip conditions are imposed, u = 0, v = 0, while on the deforming interface the velocities equal the deformation rate, u0 = ut and v0 = vt (the default condition; note that u and v on the right-hand sides refer to the displacement components).

Structural Mechanics

The structural deformations are solved for using an elastic formulation and a nonlinear geometry formulation to allow large deformations.

The obstacle is fixed to the bottom of the fluid channel. All other object boundaries experience a load from the fluid, given by

where n is the normal vector to the boundary. This load represents a sum of pressure and viscous forces.

Moving Mesh

The Navier–Stokes equations are solved on a freely moving deformed mesh, which constitutes the fluid domain. The deformation of this mesh relative to the initial shape of the domain is computed using Yeoh smoothing. Inside the obstacle, the moving mesh follows the deformations of the obstacle. At the exterior boundaries of the flow domain, the deformation is zero in all directions.

Results and Discussion

Figure 2 shows the geometry deformation and flow at t = 4 s when the system is close to its steady state. Due to the channel’s small dimensions, the Reynolds number of the flow is small (Re << 100), and the flow stays laminar in most of the area. The swirls are restricted to a small area behind the structure. The amount of deformation as well as the size and location of the swirls depend on the magnitude of the inflow velocity.

Figure 2: Flow velocity and geometry deformation at t = 4 s. The streamlines indicate the flow direction and the color indicates flow-velocity magnitude.

Figure 3 shows the mesh velocity at t = 0.15 s. The boundaries of the narrow structure are the only moving boundaries of the flow channel. Therefore the mesh velocity also has its largest values near the structure. Depending on the current state of the deformation — whether it is increasing, decreasing or stationary — the mesh velocity can have a very different distribution. Figure 4 further illustrates this point; it compares the average inflow velocity to the horizontal mesh velocity and the horizontal mesh displacement just beside the top of the structure. Most of the time the deformation follows the inflow velocity quite closely. Whenever the inflow velocity starts to decrease, the deformation also decreases, which you can observe as the negative values on the horizontal mesh velocity. Toward the end of the simulation, when inflow and structure deformation approach their steady-state values, the mesh velocity also decreases to zero.

Figure 3: Mesh velocity (arrows) and mesh and geometry deformation at t = 0.15 s.

Figure 4: Inflow velocity, horizontal mesh velocity, and mesh deformation. The blue curve shows the average x direction velocity at the inflow boundary (m/s); the green shows 104× mesh displacement in the x direction (dx_ale; m) at the geometry point (1.05·10-4, 0.5·10-4); and the red curve shows 103× mesh velocity in the x direction (xt; m/s), also at the point (1.05·10-4, 0.5·10-4).

Figure 5 compares the meshes at different times. The first image shows the initial mesh, which you generate prior to solving the model. This mesh is equally distributed around the top of the structure. The second image shows the mesh in its deformed form. Because the structure deforms more in the horizontal direction, the mesh also changes more in this direction: On the left, the mesh elements are stretched; on the right, they are compressed in the x direction.

Figure 5: Geometry and mesh near the top of the structure at t = 0 s and 2 s.

Notes About the COMSOL Implementation

This example implements the model using Fluid–Solid Interaction multiphysics interface. The interface includes Laminar Flow and Solid Mechanics interfaces to model fluid and solid domains. The Deforming Domain feature is added to represent the domain and boundaries where the mesh can deform. By default, the Deforming Domain node has an empty selection. In this model, the fluid domain deforms significantly due to the deformation of the solid. Therefore, add a fluid domain to the selection list in the Deforming Domain node.

To model the interaction at the fluid–solid boundaries, The Fluid–Structure Interaction coupling feature is included in the interface. The feature automatically identifies the fluid–solid interaction boundaries and assigns the boundary condition to those boundaries.

To get a more accurate computation of the large strains, large deformation analysis is the default setting.

MEMS_Module/Fluid-Structure_Interaction/fluid_structure_interaction

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Parameters 1

In the window, under click .

In the window for , locate the section.

In the table, enter the following settings:


Name	Expression	Value	Description
U	3.33[cm/s]	0.0333 m/s	Inlet mean velocity at steady state
H	100[um]	1E-4 m	Channel height

Definitions

Analytic 1 (an1)

In the toolbar, click

In the window for , locate the section.

In the text field, type t.

In the text field, type t^2/sqrt((0.04[s^2]-t^2)^2+(0.1[s]*t)^2).

Locate the section. In the table, enter the following settings:


Argument	Unit
t	s

Geometry 1

In the window, under click .

In the window for , locate the section.

From the list, choose .

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 300.

In the text field, type H.

Click

Rectangle 2 (r2)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5.

In the text field, type 47.5.

Locate the section. In the text field, type 100.

Click

Fillet 1 (fil1)

In the toolbar, click

On the object , select Points 3 and 4 only.

It might be easier to select the points by using the window. To open this window, in the toolbar click and choose . (If you are running the cross-platform desktop, you find in the main menu.)

In the window for , locate the section.

In the text field, type 2.5.

Click

The geometry should look like that in the figure below.

Modify the domain selections for the and interface.

Laminar Flow (spf)

In the window, under click .

Select Domain 1 only.

Solid Mechanics (solid)

In the window, under click .

Select Domain 2 only.

Moving Mesh

Deforming Domain 1

In the window, under > click .

Select Domain 1 only.

Laminar Flow (spf)

In the window, under click .

Inlet 1

In the toolbar, click

and choose .

Select Boundary 1 only.

In the window for , locate the section.

From the list, choose .

Locate the section. In the Uav text field, type U*an1(t).

This gives a parabolic velocity profile with the specified mean velocity.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 8 only.

Solid Mechanics (solid)

In the window, under click .

Fixed Constraint 1

In the toolbar, click

and choose .

Select Boundary 5 only.

Materials

Material 1 (mat1)

In the window, under right-click and choose .

Select Domain 1 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	1e3	kg/m³	Basic
Dynamic viscosity	mu	1e-3	Pa·s	Basic

Material 2 (mat2)

Right-click and choose .

Select Domain 2 only.

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Young’s modulus	E	2e5	Pa	Young’s modulus and Poisson’s ratio
Poisson’s ratio	nu	0.33	1	Young’s modulus and Poisson’s ratio
Density	rho	7850	kg/m³	Basic