# 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$

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}$
Fig. 35: The adapted mesh of the domain $\Omega^{72}$