Newton BC (06-bc-newton)

Git reference: Tutorial example 06-bc-newton.

This example uses the same domain \Omega as the previous one, and the same equation

(1)-\mbox{div}(\lambda \nabla u) - C_{src} = 0,

but the Neumann condition on the “Outer” boundary is replaced with a Newton condition of the form

(2)\lambda \frac{\partial u}{\partial n} = \alpha (T_{ext} - u)

where \alpha > 0 is a heat transfer coefficient and T_{ext} exterior temperature. On the boundary parts “Bottom”, “Left” and “Inner” we keep the nonconstant Dirichlet condition

u(x, y) = Ax + By + C.

The weak formulation is

(3)\int_\Omega \lambda \nabla u \cdot \nabla v \;\mbox{d\bfx} - \int_{\Gamma_{Outer}} \lambda \frac{\partial u}{\partial n}v \;\mbox{dS} - \int_\Omega C_{src} v\;\mbox{d\bfx} = 0.

Using the Newton condition (2), we obtain

(4)\int_\Omega \lambda \nabla u \cdot \nabla v \;\mbox{d\bfx} + \int_{\Gamma_{Outer}} \alpha u v \;\mbox{dS} - \int_{\Gamma_{Outer}} \alpha T_{ext} v \;\mbox{dS} - \int_\Omega C_{src} v\;\mbox{d\bfx} = 0.

Thus in addition to all types of weak forms from the previous example 05-bc-neumann we now also have a surface matrix form

- \int_{\Gamma_{Outer}} \alpha u v \;\mbox{dS}

Default surface matrix form

This integral can be added using the default H^1 form DefaultMatrixFormSurf that has the constructors

DefaultMatrixFormSurf(int i, int j, std::string area = HERMES_ANY,
                      HermesFunction* coeff = HERMES_ONE,
                      GeomType gt = HERMES_PLANAR);

DefaultMatrixFormSurf(int i, int j, Hermes::vector<std::string> areas,
                      HermesFunction* coeff = HERMES_ONE,
                      GeomType gt = HERMES_PLANAR);

The header of the custom form CustomWeakFormPoissonNewton is defined in definitions.h:

class CustomWeakFormPoissonNewton : public WeakForm
{
public:
  CustomWeakFormPoissonNewton(std::string mat_al, HermesFunction* lambda_al,
                              std::string mat_cu, HermesFunction* lambda_cu,
                              HermesFunction* vol_src_term, std::string bdy_heat_flux,
                              double alpha, double t_exterior);
};

and its constructor in definitions.cpp:

CustomWeakFormPoissonNewton::CustomWeakFormPoissonNewton(std::string mat_al, HermesFunction* lambda_al,
                                                         std::string mat_cu, HermesFunction* lambda_cu,
                                                         HermesFunction* vol_src_term, std::string bdy_heat_flux,
                                                         double alpha, double t_exterior) : WeakForm(1)
{
  // Jacobian forms - volumetric.
  add_matrix_form(new WeakFormsH1::DefaultJacobianDiffusion(0, 0, mat_al, lambda_al));
  add_matrix_form(new WeakFormsH1::DefaultJacobianDiffusion(0, 0, mat_cu, lambda_cu));

  // Jacobian forms - surface.
  add_matrix_form_surf(new WeakFormsH1::DefaultMatrixFormSurf(0, 0, bdy_heat_flux, new HermesFunction(alpha)));

  // Residual forms - volumetric.
  add_vector_form(new WeakFormsH1::DefaultResidualDiffusion(0, mat_al, lambda_al));
  add_vector_form(new WeakFormsH1::DefaultResidualDiffusion(0, mat_cu, lambda_cu));
  add_vector_form(new WeakFormsH1::DefaultVectorFormVol(0, HERMES_ANY, vol_src_term));

  // Residual forms - surface.
  add_vector_form_surf(new WeakFormsH1::DefaultResidualSurf(0, bdy_heat_flux, new HermesFunction(alpha)));
  add_vector_form_surf(new WeakFormsH1::DefaultVectorFormSurf(0, bdy_heat_flux, new HermesFunction(-alpha * t_exterior)));
};

Sample results

The output for the parameters C_{src} = 0, \lambda_{Al} = 236, \lambda_{Cu} = 386, \alpha = 5, T_{ext} = 50, A = 0, B = 0 and C = 20 is shown below:

Solution of the Newton problem.

Table Of Contents

Previous topic

Neumann BC (05-bc-neumann)

Next topic

General 2nd-Order Linear Equation (07-general)