neumann_bc.mac (1223B)
1 load("generic_equation.mac")$ 2 3 /* ethalpy equation for the positive and negative w[k] cases */ 4 eq[pm] := E[k-1]*l[pm] + E[k]*d[pm] + E[k+1]*u[pm] = b[pm]; 5 6 /* Additional equation implementing the Neumann B.C. */ 7 neumann[k] := (E[k+1] - E[k-1]) / (2 * dz) = G[k]; 8 9 /* Basal equation */ 10 E_base : solve(neumann[0], E[-1])[1]; 11 eq_base[pm] := facsum(subst([k=0, E_base], eq[pm]), E[0], E[1]); 12 13 /* Matrix entries for the basal equation */ 14 d_base[pm] := facsum(coeff(lhs(eq_base[pm]), E[0]), w[0]); 15 u_base[pm] := facsum(coeff(lhs(eq_base[pm]), E[1]), w[0]); 16 17 b_base[pm] := 18 expandwrt( 19 facsum( 20 expand( 21 rhs(eq_base[pm]) - (lhs(eq_base[pm]) - E[0]*d_base[pm] - E[1]*u_base[pm])), 22 dz), 23 rho); 24 25 /* Surface equation */ 26 eq_surface[pm] := facsum(subst([k=k_s, E_surface], eq[pm]), E[k_s-1], E[k_s]); 27 E_surface : solve(neumann[k_s], E[k_s+1])[1]; 28 29 /* Matrix entries for the surface equation */ 30 l_surface[pm] := facsum(coeff(lhs(eq_surface[pm]), E[k_s-1]), w[k_s]); 31 d_surface[pm] := facsum(coeff(lhs(eq_surface[pm]), E[k_s]), w[k_s]); 32 33 b_surface[pm] := 34 expandwrt( 35 facsum( 36 expand( 37 rhs(eq_surface[pm]) - (lhs(eq_surface[pm]) - E[k_s-1]*l_surface[pm] - E[k_s]*d_surface[pm])), 38 dz), 39 rho); 40 41 G_eq : G = phi / K;