Literature survey of acupuncture study
2.3. SOLVING PROBLEMS IN FREEFEM++
2.3.3 Moving domain problem in computational fluid dynamics
Different methods exist to predict flows in moving domains with the finite element method. Fixed mesh methods use a single mesh covering the domain occupied by the fluid and the domain where there is no fluid to compute the flow. The immersed boundary approach [88] relies on the description of a solid phase by adding a force vector to the governing equations. A similar approach, known as the fictitious domain method, is based on the use of Lagrange multipliers to enforce kinematic condition on the solid phase [89–
91] or alternatively based on a penalty method [92]. Both methods track the solid phase with a characteristic function or a level set function. These methods can be adapted to track moving bodies in fluid. The latest methods have been implemented with FreeFem++
[93].
Moving mesh methods use meshes following the domain occupied by the fluid. The mesh is then given exactly at the boundary. The ALE framework is mathematically rig-orous to describe transport phenomena in time. However, it raises some implementation
2.3. SOLVING PROBLEMS IN FREEFEM++
Figure 2.12: Numerical simulation of flow in idealized 2D aneurysm model, stenosis model, bifurcation model, and bended tube model. Velocity vectors are computed with the FreeFem++ script 2.17.
2.3. SOLVING PROBLEMS IN FREEFEM++
questions on the interface tracking with time discretization. Implementation of the ALE method can be done in FreeFem++. The ALE method of characteristics [92] applied to the mesh moving scheme can be implemented to improve stability properties of the mesh moving scheme. The space-time finite element method is also well-suited to track the moving domain [76, 77]. The space-time finite element method can be implemented in FreeFem++ in 1D and 2D using an extra variable such as the time variable.
Moving mesh methods use meshes following the domain occupied by the fluid. The mesh is then exactly given at the boundary. The ALE framework is mathematically rigor-ous to describe transport phenomena in time and allows some freedom in the description of the mesh motion. However, it raises some implementation questions on the interface tracking with time discretization. Implementation of the ALE method can be done in FreeFem++. The ALE method of characteristics [94] applied to the mesh moving scheme can be implemented to improve stability properties of the mesh moving scheme. The space-time finite element method is also well suited to track the moving domain [78, 79].
The space-time finite element method can be implemented in FreeFem++ in 1D and 2D using and extra variable as the time variable.
A straightforward implementation of the ALE approach is briefly described below.
The implementation in FreeFem++ follows the method given by Decoene A. and Maury B. [95]. For a detailed description of the ALE approach, the reader can refer to [96]. Let Ω(t) be a domain at each time t with regular boundary ∂Ω(t). In the Eulerian representa-tion, the fluid is described by
u(x, t) and p(x, t), ∀x ∈ Ω(t). (2.24) To follow a moving domain, one can define the ALE map as
A : ˜˜ ω × R+ → R2 (˜x, t) → ˜A(˜x, t) := ˜At, (2.25) such that ω(t) = ˜A(˜ω, t), where ˜ω is the reference computational domain. Given an ALE field ˜q : ˜ω × R+ → R, its Eulerian description is given by
∀x ∈ Ω(t), q(x, t) = ˜q( ˜A−1t (x), t) (2.26) In this framework, the computational domain velocity (ALE velocity or grid velocity) is defined as
a(˜˜ x, t) = ∂ ˜A
∂t(˜x, t), ∀˜x ∈ ˜ω, (2.27) so that
a(x, t) = ˜a( ˜A−1t , t). (2.28) The ALE time-derivative is defined as
∂q
∂t ˜
A
= d
dtq( ˜A(˜x, t), t), (2.29)
2.3. SOLVING PROBLEMS IN FREEFEM++
and the following identity holds
∂q
A general method used to construct the mapping, or equivalently the domain velocity a , consists of solving the the following Laplace equation its boundary condition
−∇2a = 0,
a|∂Ω = f , (2.31)
where f is the velocity of the boundary.
A better model for dealing with a computationally more challenging large deformation case can be given by
−∇2a + ∇p = 0,
∇ · a = 0, a|∂Ω = f .
(2.32)
Let the interaction of an oscillating circular cylinder with a fluid at rest be considered.
The problem is to find the velocity vector field u and the pressure p of a flow satisfying the incompressible Navier-Stokes (2.13) in the domain Ω = [−l, l] × [−h, h] with no-slip boundary conditions. no-no-slip boundary conditions on the cylinder and traction-free boundary condition on the outside border.
In the ALE framework, the weak formulation of the incompressible Navier-Stokes is as follows: find (un+1, pn+1) ∈ V such that The horizontal velocity of the cylinder of diameter D is given by uc(t) = −U cos(2π f t), where U = 2π A f . At each time step, the mesh is moved with the command movemesh according to the displacement a ∆t, where a is the solution of (2.31) (see figure 2.13). As soon as the mesh is moved, the computed anand un, that are defined in the previous mesh, are then pushed to the new mesh without interpolation following the scheme proposed in [95]. The solution can be computed with the FreeFem++ script 2.18.
2.3. SOLVING PROBLEMS IN FREEFEM++
1 i n t n = 2 0 , nc = 8 0 , M= 1 0 0 ;
2 r e a l Re = 1 0 0 , KC= 5 . , U= 1 , D= 1 , R=D / 2 . , A=KC∗D / ( 2 ∗p i) , F=U / KC / D ;
3 r e a l L=50∗D , H=30∗D ;
4 r e a l x0 = ( − 1 . / 2 . /p i) ∗s i n( 2 ∗p i∗ 0 ) , y0 = 0 , nu = 1 . / Re , T= 1 , d t =T /M;
5 b o r d e r b1 ( t =−L , L ) {x= t ;y=−H ; l a b e l= 1 ; }
6 b o r d e r b2 ( t =−H , H) {x=L ;y= t ;l a b e l= 1 ; }
7 b o r d e r b3 ( t =L,−L ) {x= t ;y=H ;l a b e l= 1 ; }
8 b o r d e r b4 ( t =H,−H) {x=−L ;y= t ;l a b e l= 1 ; }
9 b o r d e r c ( t = 0 , 2 ∗p i) {x=x0+R∗c o s( t ) ;y=y0+R∗s i n( t ) ; l a b e l= 3 ; }
10 mesh Th=b u i l d m e s h( b1 ( n ) +b2 ( n ∗H / L ) +b3 ( n ) +b4 ( n ∗H / L ) + c (− nc ) ) ;
11 f e s p a c e Vh ( Th ,P2) ;
12 Vh u = 0 , uu , v = 0 , vv , u o l d , v o l d ;
13 Vh a1 , a2 , w , s , f 1 = 0 , f 2 = 0 ;
14 r e a l[ i n t ] tmp ( u [ ] . n ) ;
15 f e s p a c e Ph ( Th ,P1) ;
16 Ph p = 0 , pp ;
17 i n t m= 0 ;
18 p r o b l e m M e s h V e l o c i t y ( [ a1 , a2 ] , [ w , s ] , i n i t=m) / / Lapalce eq : mesh velocity
19 = i n t 2 d( Th ) (dx( a1 ) ∗dx(w) +dy( a1 ) ∗dy(w) +dx( a2 ) ∗dx( s ) +dy( a2 ) ∗dy( s ) )
20 +on( 1 , 2 , a1 = 0 , a2 = 0 )
21 +on( 3 , a1 = f 1 , a2 = f 2 ) ;
22 p r o b l e m a l e N S ( [ u , v , p ] , [ uu , vv , pp ] , i n i t=m) / / ALE Navier-Stokes eq.
23 = i n t 2 d( Th ) ( ( u ∗ uu+v ∗ vv ) / d t
24 + nu ∗ (dx( u ) ∗dx( uu ) +dy( u ) ∗dy( uu ) +dx( v ) ∗dx( vv ) +dy( v ) ∗dy( vv ) )
25 − (dx( u ) +dy( v ) ) ∗ pp
26 + (dx( p ) ∗ uu+dy( p ) ∗ vv )
27 − 1 e −5∗p ∗ pp )
28 −i n t 2 d( Th ) ( (c o n v e c t ( [ u o l d −a1 , v o l d −a2 ] , − d t , u o l d ) ∗ uu
29 + c o n v e c t( [ u o l d −a1 , v o l d −a2 ] , − d t , v o l d ) ∗ vv ) / d t )
30 +on( 3 , u=a1 , v= a2 ) ;
31 f o r (m= 0 ;m<M;m++) {
32 f 1=−U∗c o s( 2 ∗p i∗F∗m∗ d t ) ;
33 f 2 = 0 . ;
34 a l e N S ; / / compute (u,p)
35 M e s h V e l o c i t y ; / / compute the mesh velocity
36 Th = movemesh( Th , [x+ d t ∗ a1 ,y+ d t ∗ a2 ] ) ; / / move the mesh
37 tmp=u [ ] ; u o l d = 0 ; u o l d [ ] = tmp ; / / push to the new mesh
38 tmp=v [ ] ; v o l d = 0 ; v o l d [ ] = tmp ; / / push to the new mesh
39 tmp= a1 [ ] ; a1 = 0 ; a1 [ ] = tmp ; / / push to the new mesh
40 tmp= a2 [ ] ; a2 = 0 ; a2 [ ] = tmp ; / / push to the new mesh
41 p l o t(c o e f= 0 . 1 , [ u , v ] , v a l u e = 1 ,w a i t= 0 ) ; }
Script 2.18: Solving the incompressible Navier-Stokes equations in a moving domain with the ALE method
2.3. SOLVING PROBLEMS IN FREEFEM++
Figure 2.13: Moving meshes in FreeFem++
In figure 2.14, numerical simulation of flow is shown for Re = 100 and D fU = 5, where U = 1, D = 1, T := 1/f = 5, and A = 5/2π. In figures 2.15, 2.3.3, and 2.17, good comparison between the computed solution, solution from [97], and experimen-tal data from [98] is shown at at three different times corresponding to the three phases 2πf t = π, 7π/6, and 11π/6.
Figure 2.14: Numerical simulation of flow accelerated from rest by an oscillating circular cylinder. Velocity magnitude is at t = 92T (left), t = 5512T (middle), and t = 5912T (right).
2.3. SOLVING PROBLEMS IN FREEFEM++ Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
Figure 2.15: Comparison of the computed and referenced solutions along the line x =
−0.6 ∗ D (top) and x = 1.2 ∗ D (bottom) for the values of u = (u, v) at t = 92T . Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
Figure 2.16: Comparison of the computed and referenced solutions along the line x =
−0.6 ∗ D (top) and x = 1.2 ∗ D (bottom) for the values of u = (u, v) at t = 5512T .
2.4. CONCLUDING REMARKS Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
-2 Experimental (Dutsch al., 1998)
Figure 2.17: Comparison of the computed and referenced solutions along the line x =
−0.6 ∗ D (top) and x = 1.2 ∗ D (bottom) for the values of u = (u, v) at t = 5912T .
2.4 Concluding remarks
FreeFem++ makes it easy to generate and manipulate meshes for the resolution of partial differential equations in two and three dimensions. Simple functions presented in this chapter allow the user to perform challenging tasks on a triangulation. External tools and software can be linked with FreeFem++ to deal with more complex tasks especially in three dimensions.
FreeFem++’s major task consists in assembling automatically the stiffness matrix and the right hand side associated with the discretized problem within the framework of the finite element method. Furthermore, FreeFem++ possesses a large panel of linear solvers.
There is no systematic tool to treat evolution and nonlinear problems. Special treat-ments need to be provided by the user to work within the FreeFem++ framework. This gives a certain advantage in making it possible for to user to work efficiently on a given problem. In the specific case of convection problems, FreeFem++ makes it possible to use the method of the characteristics directly.
Finally, while all the tasks presented in this chapter are sequential, they can be paral-lelized by calling different processes simultaneously.
Chapter 3
Modeling and simulation of local physical stress field during needling
Contents
3.1 Introduction . . . 79 3.2 Biological medium . . . 79 3.3 Mathematical modeling . . . 80 3.4 Computational model . . . 81 3.4.1 Scaling and setting for numerical simulations . . . 81 3.4.2 Numerical methods . . . 82 3.5 Results and discussion . . . 84 3.5.1 Effect of needle motion on the interstitial flow . . . 84 3.5.2 Effects of fractional fluid volume and Darcy number on the
in-terstitial flow . . . 85 3.5.3 Shear stress and pressure distributions along the cell membrane 89 3.6 Concluding remarks . . . 92
The work presented in this chapter has association with the following publication.
• Y. Deleuze, M. Thiriet, T.W.H. Sheu. "Modeling and simulation of local physical stress on the mastocytes created by the needle manipulation during acupuncture"
(forthcoming).
Parameter Definition Dimension
Da Darcy number
d cell center/needle tip distance along the x-axis m
fd,s fiber diameter and spacing m
L characteristic length m
n unit outward normal vector
P Darcy permeability m2
pf pressure kg. m−1. s−2
Re Rayleigh number
t time s
¯
u averaged velocity vector m−3. s−1
uf fluid velocity m. s−1
vneedle needle velocity m. s−1
V characteristic velocity m. s−1
x = (x, y) coordinate axes Greek symbols
αf fluid volume fraction
µ dynamic viscosity kg. m−1. s−1
ρ fluid density kg. m−3
θ angle rad
Table 3.1: Nomenclature and parameter dimensions Abstract
In this chapter, the effects of an needle, inserted in the subcutaneous tissue, on the interstitial flow is studied. One of the goals is to describe the physical stress affecting cells during acupuncture treatments. A Brinkman model is considered to describe the flow through a fibrous medium. Numerical studies in FreeFem++ are performed to illustrate the acute interstitial pressure de-veloped by the implanted needle that triggers the physiological reactions of acupuncture.
3.1. INTRODUCTION
3.1 Introduction
The initiation of the effects of acupuncture is due to the local stress generation. The insertion and manipulation of an acupuncture needle cause the deformation of the inter-stitial tissue to occur and in turn makes the corresponding change in interinter-stitial flow. The local stress field is sensed by mastocytes and the sensory transduction is immediately followed by Ca++ entry and subsequently granule exocytosis. After a short period of needle manipulation, the needle remains in the skin until the desired effects have been obtained. This stimulation can thus lead to a cascade of biochemical reactions that drive acupuncture effects.
In the present chapter, focus is given to the effects of interstitial fluid flow during implantation of an acupuncture needle until the tip has reached the desired location within the hypodermis. The objective of this work is to give a description of the physical stress produced by the interstitial fluid and affecting the network of fibers and cells in the tissue.
3.2 Biological medium
The interstitial fluid contains water, ions and other small molecules. Such a fluid is like plasma without macromolecules. It interacts with the ground substance to form a gel-like medium.
A model taking into account individual fibers and cell adhesion complexes is already a falsification of reality. Moreover, it is very costly from the computational sense. When considering an organized homogeneous matrix of fibers, computation of such a model shows the microscopic fluctuations of the fluid shear stress at the protein level [99].
On a microscopic scale, the interstitial tissue is composed of fluid and solid fibers, thereby clearly forming two phases. Homogenized two-phase media coupling with a Newtonian fluid and an elastic matrix have been considered to model soft tissue [100].
This well-posed model exhibits both the fluid and the viscoelastic property of the fluid.
The interstitial tissue can be modeled as a porous medium. Darcy’s law approxi-mates fibers of the media as a continuum and allows us to compute the actual microscopic flow phenomena that occurs in the fibrous media. The phenomenological model cannot give information on unneeded microscopic events but the Darcy equation can describe macroscale flow patterns in porous media.
The Brinkman equation is an extension of the Darcy equation. Introduction of a sec-ond order derivative in the Darcy equation allows the application of no-slip boundary conditions. The Brinkman equation can thus describe the flow field around solid bodies such as the embedded cells in extracellular matrix but the fibrous medium itself is still treated as a continuum.
The interstitial tissue can also be modeled by a poroviscoelastic material [18, 101–