FreeFEM Documentation on GitHub

stars - forks

Fluid-structure coupled problem

This problem involves the Lamé system of elasticity and the Stokes system for viscous fluids with velocity \(\mathbf{u}\) and pressure \(p\):

\[\begin{split}\begin{array}{rcll} -\Delta\mathbf{u} + \mathbf{\nabla} p &=& 0 & \hbox{ in }\Omega\\ \nabla\cdot\mathbf{u} &=& 0 & \hbox{ in }\Omega\\ \mathbf{u} &=& \mathbf{u}_\Gamma & \hbox{ on }\Gamma=\partial\Omega \end{array}\end{split}\]

where \(u_\Gamma\) is the velocity of the boundaries. The force that the fluid applies to the boundaries is the normal stress

\[\mathbf{h} =(\nabla\mathbf{u} +\nabla\mathbf{u}^T)\mathbf{n} -p\mathbf{n}\]

Elastic solids subject to forces deform: a point in the solid at (x,y) goes to (X,Y) after. When the displacement vector \(\mathbf{v}=(v_1,v_2) = (X-x, Y-y)\) is small, Hooke’s law relates the stress tensor \(\sigma\) inside the solid to the deformation tensor \(\epsilon\):

\[\sigma_{ij} = \lambda \delta_{ij} \nabla.\mathbf{v} + 2\mu\epsilon_{ij}, \, \epsilon_{ij} = {1\over 2}({\partial v_i\over\partial x_j} + {\partial v_j\over\partial x_i} )\]

where \(\delta\) is the Kronecker symbol and where \(\lambda\), \(\mu\) are two constants describing the material mechanical properties in terms of the modulus of elasticity, and Young’s modulus.

The equations of elasticity are naturally written in variational form for the displacement vector \(v(x)\in V\) as:

\[\int_\Omega \left[2\mu\epsilon_{ij}(\mathbf{v})\epsilon_{ij}(\mathbf{w}) +\lambda \epsilon_{ii}(v)\epsilon_{jj}(\mathbf{w})\right] =\int_\Omega \mathbf{g}\cdot \mathbf{w} +\int_\Gamma \mathbf{h}\cdot \mathbf{w},%\`{u} \forall \mathbf{w}\in V\]

The data are the gravity force \(\mathbf{g}\) and the boundary stress \(\mathbf{h}\).

Tip

Fluide-structure In our example, the Lamé system and the Stokes system are coupled by a common boundary on which the fluid stress creates a displacement of the boundary and hence changes the shape of the domain where the Stokes problem is integrated. The geometry is that of a vertical driven cavity with an elastic lid. The lid is a beam with weight so it will be deformed by its own weight and by the normal stress due to the fluid reaction. The cavity is the \(10 \times 10\) square and the lid is a rectangle of height \(l=2\).

A beam sits on a box full of fluid rotating because the left vertical side has velocity one. The beam is bent by its own weight, but the pressure of the fluid modifies the bending.

The bending displacement of the beam is given by \((uu, vv)\) whose solution is given as follows.

 1// Parameters
 2int bottombeam = 2; //label of bottombeam
 3real E = 21.5;
 4real sigma = 0.29;
 5real gravity = -0.05;
 6real coef = 0.2;
 7
 8// Mesh (solid)
 9border a(t=2, 0){x=0; y=t; label=1;}
10border b(t=0, 10){x=t; y=0; label=bottombeam;}
11border c(t=0, 2){x=10; y=t; label=1;}
12border d(t=0, 10){x=10-t; y=2; label=3;}
13mesh th = buildmesh(b(20) + c(5) + d(20) + a(5));
14
15// Fespace (solid)
16fespace Vh(th, P1);
17Vh uu, w, vv, s;
18
19// Macro
20real sqrt2 = sqrt(2.);
21macro epsilon(u1, u2) [dx(u1), dy(u2), (dy(u1)+dx(u2))/sqrt2] //
22macro div(u1, u2) (dx(u1) + dy(u2)) //
23
24// Problem (solid)
25real mu = E/(2*(1+sigma));
26real lambda = E*sigma/((1+sigma)*(1-2*sigma));
27solve Elasticity([uu, vv], [w, s])
28    = int2d(th)(
29          lambda*div(w,s)*div(uu,vv)
30        + 2.*mu*(epsilon(w,s)'*epsilon(uu,vv))
31    )
32    + int2d(th)(
33        - gravity*s
34    )
35    + on(1, uu=0, vv=0)
36    ;
37
38plot([uu, vv], wait=true);
39mesh th1 = movemesh(th, [x+uu, y+vv]);
40plot(th1, wait=true);

Then Stokes equation for fluids at low speed are solved in the box below the beam, but the beam has deformed the box (see border h):

1// Mesh (fluid)
2border e(t=0, 10){x=t; y=-10; label= 1;}
3border f(t=0, 10){x=10; y=-10+t ; label= 1;}
4border g(t=0, 10){x=0; y=-t; label= 2;}
5border h(t=0, 10){x=t; y=vv(t,0)*( t>=0.001 )*(t <= 9.999); label=3;}
6mesh sh = buildmesh(h(-20) + f(10) + e(10) + g(10));
7plot(sh, wait=true);

We use the Uzawa conjugate gradient to solve the Stokes problem like in Navier-Stokes equations.

 1// Fespace (fluid)
 2fespace Xh(sh, P2);
 3Xh u1, u2;
 4Xh bc1;
 5Xh brhs;
 6Xh bcx=0, bcy=1;
 7
 8fespace Mh(sh, P1);
 9Mh p, ppp;
10
11// Problem (fluid)
12varf bx (u1, q) = int2d(sh)(-(dx(u1)*q));
13varf by (u1, q) = int2d(sh)(-(dy(u1)*q));
14varf Lap (u1, u2)
15    = int2d(sh)(
16          dx(u1)*dx(u2)
17        + dy(u1)*dy(u2)
18    )
19    + on(2, u1=1)
20    + on(1, 3, u1=0)
21    ;
22
23bc1[] = Lap(0, Xh);
24
25matrix A = Lap(Xh, Xh, solver=CG);
26matrix Bx = bx(Xh, Mh);
27matrix By = by(Xh, Mh);
28
29
30func real[int] divup (real[int] & pp){
31    int verb = verbosity;
32    verbosity = 0;
33    brhs[] = Bx'*pp;
34    brhs[] += bc1[] .*bcx[];
35    u1[] = A^-1*brhs[];
36    brhs[] = By'*pp;
37    brhs[] += bc1[] .*bcy[];
38    u2[] = A^-1*brhs[];
39    ppp[] = Bx*u1[];
40    ppp[] += By*u2[];
41    verbosity = verb;
42    return ppp[];
43}

do a loop on the two problems

1// Coupling loop
2for(int step = 0; step < 10; ++step){
3    // Solve (fluid)
4    LinearCG(divup, p[], eps=1.e-3, nbiter=50);
5    divup(p[]);

Now the beam will feel the stress constraint from the fluid:

1// Forces
2Vh sigma11, sigma22, sigma12;
3Vh uu1=uu, vv1=vv;
4
5sigma11([x+uu, y+vv]) = (2*dx(u1) - p);
6sigma22([x+uu, y+vv]) = (2*dy(u2) - p);
7sigma12([x+uu, y+vv]) = (dx(u1) + dy(u2));

which comes as a boundary condition to the PDE of the beam:

 1// Solve (solid)
 2solve Elasticity2 ([uu, vv], [w, s], init=step)
 3= int2d(th)(
 4      lambda*div(w,s)*div(uu,vv)
 5    + 2.*mu*(epsilon(w,s)'*epsilon(uu,vv))
 6)
 7+ int2d(th)(
 8    - gravity*s
 9)
10+ int1d(th, bottombeam)(
11    - coef*(sigma11*N.x*w + sigma22*N.y*s + sigma12*(N.y*w+N.x*s))
12)
13+ on(1, uu=0, vv=0)
14;
15
16// Plot
17plot([uu, vv], wait=1);
18
19// Error
20real err = sqrt(int2d(th)((uu-uu1)^2 + (vv-vv1)^2));
21cout << "Erreur L2 = " << err << endl;

Notice that the matrix generated by Elasticity2 is reused (see init=i). Finally we deform the beam:

1// Movemesh
2th1 = movemesh(th, [x+0.2*uu, y+0.2*vv]);
3plot(th1, wait=true);

Fluid velocity and pressure, displacement vector of the structure and displaced geometry in the fluid-structure interaction of a soft side and a driven cavity are shown Fig. 173, Fig. 174 and Fig. 175

../_images/FluidStructure1.png

Fig. 173 Velocity and pressure

FluidStructure2

Fig. 174 Displacement

FluidStructure3

Fig. 175 Moved mesh

Table of content