Free boundary problem

The domain \Omega is defined with:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
// Parameters
real L = 10; //length
real hr = 2.1; //left height
real hl = 0.35; //right height
int n = 4;

// Mesh
border a(t=0, L){x=t; y=0;}; //bottom: Gamma_a
border b(t=0, hr){x=L; y=t;}; //right: Gamma_b
border f(t=L, 0){x=t; y=t*(hr-hl)/L+hl;}; //free surface: Gamma_f
border d(t=hl, 0){x=0; y=t;}; // left: Gamma_d
mesh Th = buildmesh(a(10*n) + b(6*n) + f(8*n) + d(3*n));
plot(Th);
Fig. 33: The mesh of the domain \Omega
Mesh1

The free boundary problem is:

Find u and \Omega such that:

We use a fixed point method;

\Omega^0 = \Omega

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

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

The variational formulation is:

Find u on V=H^1(\Omega^n), such than u=y on \Gamma^n_b and \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{equation*} \left\{ \begin{array}{rcll} -\Delta v &=& 0 & \mbox{ in }\Omega^n\\ v &=& 0 & \mbox{ on }\Gamma^n_a\\ \p v \over \p n &=& 0 & \mbox{ on }\Gamma^n_b\cup\Gamma^n_d\\ \p v \over \p n &=& {\p u \over \p n} - {q\over K} n_x & \mbox{ on }\Gamma^n_f \end{array} \right. \end{equation*}

The variational formulation is:

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

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

Free boundary

The FreeFem++ implementation is:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
// Parameters
real L = 10; //length
real hr = 2.1; //left height
real hl = 0.35; //right height
int n = 4;

real q = 0.02; //incoming flow
real K = 0.5; //permeability

// Mesh
border a(t=0, L){x=t; y=0;}; //bottom: Gamma_a
border b(t=0, hr){x=L; y=t;}; //right: Gamma_b
border f(t=L, 0){x=t; y=t*(hr-hl)/L+hl;}; //free surface: Gamma_f
border d(t=hl, 0){x=0; y=t;}; // left: Gamma_d
mesh Th = buildmesh(a(10*n) + b(6*n) + f(8*n) + d(3*n));
plot(Th);

// Fespace
fespace Vh(Th, P1);
Vh u, v;
Vh uu, vv;

// Problem
problem Pu (u, uu, solver=CG)
    = int2d(Th)(
          dx(u)*dx(uu)
        + dy(u)*dy(uu)
    )
    + on(b, f, u=y)
    ;

problem Pv (v, vv, solver=CG)
    = int2d(Th)(
          dx(v)*dx(vv)
        + dy(v)*dy(vv)
    )
    + on(a, v=0)
    + int1d(Th, f)(
          vv*((q/K)*N.y - (dx(u)*N.x + dy(u)*N.y))
    )
    ;

// Loop
int j = 0;
real errv = 1.;
real erradap = 0.001;
while (errv > 1e-6){
    // Update
    j++;

    // Solve
    Pu;
    Pv;

    // Plot
    plot(Th, u, v);

    // Error
    errv = int1d(Th, f)(v*v);

    // Movemesh
    real coef = 1.;
    real mintcc = checkmovemesh(Th, [x, y])/5.;
    real mint = checkmovemesh(Th, [x, y-v*coef]);

    if (mint < mintcc || j%10==0){ //mesh too bad => remeshing
        Th = adaptmesh(Th, u, err=erradap);
        mintcc = checkmovemesh(Th, [x, y])/5.;
    }

    while (1){
        real mint = checkmovemesh(Th, [x, y-v*coef]);

        if (mint > mintcc) break;

        cout << "min |T| = " << mint << endl;
        coef /= 1.5;
    }

    Th=movemesh(Th, [x, y-coef*v]);

    // Display
    cout << endl << j << " - errv = " << errv << endl;
}

// Plot
plot(Th);
plot(u, wait=true);
Fig. 34: The final solution on the new domain \Omega^{72}
Sol
Fig. 35: The adapted mesh of the domain \Omega^{72}
Mesh2