FreeFEM Documentation on GitHub

stars - forks

Free boundary problems

The domain \(\Omega\) is defined with:

 1 // Parameters
 2 real L = 10; //length
 3 real hl = 2.1; //left height
 4 real hr = 0.35; //right height
 5 int n = 4;
 6 
 7 // Mesh
 8 border a(t=0, L){x=t; y=0;}; //bottom: Gamma_a
 9 border b(t=0, hr){x=L; y=t;}; //right: Gamma_b
10 border f(t=L, 0){x=t; y=t*(hr-hl)/L+hl;}; //free surface: Gamma_f
11 border d(t=hl, 0){x=0; y=t;}; // left: Gamma_d
12 mesh Th = buildmesh(a(10*n) + b(6*n) + f(8*n) + d(3*n));
13 plot(Th);
../_images/FreeBoundary_Mesh1.png

Fig. 177 The mesh of the domain \(\Omega\)

The free boundary problem is:

Find \(u\) and \(\Omega\) such that:

\[\begin{split}\left\{ \begin{array}{rcll} -\Delta u &=& 0 & \mbox{ in }\Omega\\ u &=& y & \mbox{ on }\Gamma_b\\ \partial u \over \partial n &=& 0 & \mbox{ on }\Gamma_d\cup\Gamma_a\\ \partial u \over \partial n &=& {q \over K} n_x & \mbox{ on }\Gamma_f\\ u &=& y & \mbox{ on }\Gamma_f \end{array} \right.\end{split}\]

We use a fixed point method;

\(\Omega^0 = \Omega\)

In two step, fist we solve the classical following problem:

\[\begin{split}\left\{ \begin{array}{rcll} -\Delta u &=& 0 & \mbox{ in }\Omega^n\\ u &=& y & \mbox{ on }\Gamma^n_b\\ \partial u \over \partial n &=& 0 & \mbox{ on }\Gamma^n_d\cup\Gamma^n_a\\ u &=& y & \mbox{ on }\Gamma^n_f \end{array} \right.\end{split}\]

The variational formulation is:

Find \(u\) on \(V=H^1(\Omega^n)\), such than \(u=y\) on \(\Gamma^n_b\) and \(\Gamma^n_f\)

\[\int_{\Omega^n}\nabla u \nabla u' = 0,\ \forall u' \in V \mbox{ with } u' =0 \mbox{ on }\Gamma^n_b \cup \Gamma^n_f\]

And secondly to construct a domain deformation \(\mathcal{F}(x,y)=[x,y-v(x,y)]\) where \(v\) is solution of the following problem:

\[\begin{split}\left\{ \begin{array}{rcll} -\Delta v &=& 0 & \mbox{ in }\Omega^n\\ v &=& 0 & \mbox{ on }\Gamma^n_a\\ \partial v \over \partial n &=& 0 & \mbox{ on }\Gamma^n_b\cup\Gamma^n_d\\ \partial v \over \partial n &=& {\partial u \over \partial n} - {q\over K} n_x & \mbox{ on }\Gamma^n_f \end{array} \right.\end{split}\]

The variational formulation is:

Find \(v\) on \(V\), such than \(v=0\) on \(\Gamma^n_a\):

\[\int_{\Omega^n} \nabla v \nabla v' = \int_{\Gamma_f^n}({\partial u \over \partial n} - { q\over K} n_x )v',\ \quad \forall v' \in V \mbox{ with } v' =0 \mbox{ on }\Gamma^n_a\]

Finally the new domain \(\Omega^{n+1} = \mathcal{F}(\Omega^n)\)

Tip

Free boundary

The FreeFEM implementation is:

 1 // Parameters
 2 real L = 10; //length
 3 real hr = 2.1; //left height
 4 real hl = 0.35; //right height
 5 int n = 4;
 6 
 7 real q = 0.02; //incoming flow
 8 real K = 0.5; //permeability
 9 
10 // Mesh
11 border a(t=0, L){x=t; y=0;}; //bottom: Gamma_a
12 border b(t=0, hr){x=L; y=t;}; //right: Gamma_b
13 border f(t=L, 0){x=t; y=t*(hr-hl)/L+hl;}; //free surface: Gamma_f
14 border d(t=hl, 0){x=0; y=t;}; // left: Gamma_d
15 mesh Th = buildmesh(a(10*n) + b(6*n) + f(8*n) + d(3*n));
16 plot(Th);
17 
18 // Fespace
19 fespace Vh(Th, P1);
20 Vh u, v;
21 Vh uu, vv;
22 
23 // Problem
24 problem Pu (u, uu, solver=CG)
25     = int2d(Th)(
26           dx(u)*dx(uu)
27         + dy(u)*dy(uu)
28     )
29     + on(b, f, u=y)
30     ;
31 
32 problem Pv (v, vv, solver=CG)
33     = int2d(Th)(
34           dx(v)*dx(vv)
35         + dy(v)*dy(vv)
36     )
37     + on(a, v=0)
38     + int1d(Th, f)(
39           vv*((q/K)*N.y - (dx(u)*N.x + dy(u)*N.y))
40     )
41     ;
42 
43 // Loop
44 int j = 0;
45 real errv = 1.;
46 real erradap = 0.001;
47 while (errv > 1e-6){
48     // Update
49     j++;
50 
51     // Solve
52     Pu;
53     Pv;
54 
55     // Plot
56     plot(Th, u, v);
57 
58     // Error
59     errv = int1d(Th, f)(v*v);
60 
61     // Movemesh
62     real coef = 1.;
63     real mintcc = checkmovemesh(Th, [x, y])/5.;
64     real mint = checkmovemesh(Th, [x, y-v*coef]);
65 
66     if (mint < mintcc || j%10==0){ //mesh too bad => remeshing
67         Th = adaptmesh(Th, u, err=erradap);
68         mintcc = checkmovemesh(Th, [x, y])/5.;
69     }
70 
71     while (1){
72         real mint = checkmovemesh(Th, [x, y-v*coef]);
73 
74         if (mint > mintcc) break;
75 
76         cout << "min |T| = " << mint << endl;
77         coef /= 1.5;
78     }
79 
80     Th=movemesh(Th, [x, y-coef*v]);
81 
82     // Display
83     cout << endl << j << " - errv = " << errv << endl;
84 }
85 
86 // Plot
87 plot(Th);
88 plot(u, wait=true);
FreeBoundary_Sol

Fig. 178 The final solution on the new domain \(\Omega^{72}\)

FreeBoundary_Mesh2

Fig. 179 The adapted mesh of the domain \(\Omega^{72}\)

Table of content