def METHOD_NAME(self):
    super().METHOD_NAME()
    # Set up some sets for the space and time indexing
    tset = self.flowsheet().config.time
    # Set up node and face sets and get integer indices for them
    izfaces, iznodes = common._face_initializer(  # pylint: disable=unused-variable
        self, self.config.control_volume_zfaces, "z"
    )
    common._submodel_boilerplate_create_if_none(self)
    common._create_thermal_boundary_conditions_if_none(self, thin=True)
    # Preexponential factor needs to be given in units of ohm*m**2
    self.log_preexponential_factor = pyo.Var(
        initialize=-50, units=pyo.units.dimensionless
    )
    self.thermal_exponent_dividend = pyo.Var(initialize=0, units=pyo.units.K)
    self.contact_fraction = pyo.Var(
        initialize=1, units=pyo.units.dimensionless, bounds=(0, 1)
    )
    @self.Expression(tset, iznodes)
    def contact_resistance(b, t, iz):
        return (
            1
            * pyo.units.ohm
            * pyo.units.m**2
            / b.contact_fraction
            * pyo.exp(
                b.log_preexponential_factor
                + b.thermal_exponent_dividend / (b.temperature[t, iz])
            )
        )
    @self.Expression(tset, iznodes)
    def voltage_drop_total(b, t, iz):
        return b.contact_resistance[t, iz] * b.current_density[t, iz]
    @self.Expression(tset, iznodes)
    def joule_heating_flux(b, t, iz):
        return b.voltage_drop_total[t, iz] * b.current_density[t, iz]
    @self.Constraint(tset, iznodes)
    def heat_flux_x_eqn(b, t, iz):
        return (
            b.heat_flux_x1[t, iz]
            == b.heat_flux_x0[t, iz] + b.joule_heating_flux[t, iz]
        )