A Very Simple Procedure for Constructing the Solution of Nonlinear Elliptic PDEs (with Nonlinear Boundary Conditions) Arising from Heat Transfer Problems in Solids

Abstract

This work presents a procedure for constructing the solution of a second order nonlinear partial differential equation subjected to nonlinear boundary conditions in a very simple and systematic way. The class of equations to be considered here includes, but is not limited to, the partial differential equations which govern the steady-state heat transfer process in bodies with temperature-dependent thermal conductivity and temperature-dependent internal heat source. In addition, it will be considered a large class of nonlinear boundary conditions, especially those arising from the description of the heat exchange process from/to bodies at high temperature levels. Proofs of the solution’s existence and uniqueness are presented.

1. Introduction

Conduction heat transfer in rigid bodies is usually described under severe approximations that, many times, are not physically realistic, such as Dirichlet boundary conditions [1,2], Neumann boundary conditions [3,4], linear boundary conditions [4,5], constant source term [1,6], and constant thermal conductivity [2,5]. Such approximations have, as their main goal, the purpose of simplifying the mathematical description, allowing the use of more available and low-cost methods for simulations.

Nevertheless, when the heat transfer process involves great temperature variations and, in particular, involves high temperature levels, these “linear approaches” may give rise to considerably inaccuracies [7,8,9]. The dependence of the thermal conductivity on the temperature has been the main subject of many works in the last twenty years, especially those which employ the Kirchhoff Transform [10,11,12,13,14,15].

The solution of nonlinear second order differential equations is the subject of many papers as, for instance, references [16,17,18,19,20]. In reference [16], the steady state heat transfer is treated by means of the local multi-quadratic radial basis function collocation method. In reference [17], a double iteration process is employed, while in [18], an enthalpy-based finite element method is employed, taking into account a phase change process. In [19], a nonstandard finite difference method is used for a heat transfer problem. Meshless methods, that represent the solution as a superposition of shape functions, have also been used for heat transfer problems [20].

The goal of this work is to provide a reliable and simple procedure for describing and simulating heat transfer problems, without the need for the usual simplifications, carried out for minimizing mathematical complexities.

In order to illustrate, in an explicit way, the main objective of this work, let us consider the following problems:

Problem 1. Find $v$, solution of the linear problem below [21]

$$\begin{array}{l}\Delta v-\overline{c}v+s=0\hspace{1em}in \mathsf{\Omega },\hspace{1em}s=\widehat{s}\left(x\right),\hspace{1em}\overline{c}=constant\ge 0\\ -\nabla v·n=\alpha v-\gamma \hspace{1em}on \partial \mathsf{\Omega },\hspace{1em}\gamma =\widehat{\gamma }\left(x\right),\hspace{1em}\alpha =constant>0\end{array}$$

(1)

Problem 2. Find  $T$, solution of

$$\begin{array}{l}\nabla ·\left(k \nabla T\right)+\dot{q}=0\hspace{1em}in \mathsf{\Omega },\hspace{1em}k=\tilde{k}\left(T\right),\hspace{1em}\dot{q}=\widehat{q}\left(x,T\right)\\ -k\nabla T·n=\widehat{f}\left(x,T\right)\hspace{1em}on \partial \mathsf{\Omega }\end{array}$$

(2)

in which $\mathsf{\Omega }$ is a bounded open set with piecewise smooth boundary $\partial \mathsf{\Omega }$ (see Figure 1).

Clearly, Problem 1 is one of the most studied and well known mathematical problems. On the other hand, Problem 2 is a nonlinear problem that requires much more effort to be simulated than what is necessary for Problem 1.

In this work, we represent the solution of Problem 2 by means of the solution of problems such as Problem 1, provided the following conditions hold:

• There exists a positive constant $\delta$ such that, for any $T$, $\tilde{k}\left(T\right)\ge \delta >0$;
• $\dot{q}=\widehat{q}\left(x,T\right)$ is a non-increasing function of $T$, for any $x$;
• $\widehat{f}\left(x,T\right)$ is a strictly increasing function of $T$, for any $x$;
• $\underset{T\to \infty }{\lim }\widehat{f}\left(x,T\right)=+ \infty$, for any $x$;
• $\underset{T\to -\infty }{\lim }\widehat{f}\left(x,T\right)=- \infty$, for any $x$;
• $\widehat{f}\left(x,0\right)\le 0$, for any $x\in \partial \mathsf{\Omega }$;
• $\widehat{q}\left(x,a\right)\ge 0$, for any $x\in \mathsf{\Omega }$, for any $a\in \left(-\infty ,\infty \right)$.

In other words, the solution of Problem 2 will be represented by the limit of the sequence $\left[p_1,p_{2,}{p}_3\dots \right]$, in which

$$\begin{array}{l}\Delta v_i-\overline{c}{v}_i+s_{i-1}=0\hspace{1em}in \mathsf{\Omega },\hspace{1em}{s}_{i-1}={\widehat{s}}_{i-1}\left(x\right),\hspace{1em}\hspace{1em}\overline{c}=non negative constant\\ -\nabla v_i·n=\alpha v_i-{\gamma }_{i-1}\hspace{1em}on \partial \mathsf{\Omega },\hspace{1em}{\gamma }_{i-1}={\widehat{\gamma }}_{i-1}\left(x\right),\hspace{1em}\alpha =positive constant\\ v_i=\tilde{\mathsf{\Theta }}\left(p_i\right),\hspace{1em}\frac{d}{dz}\tilde{\mathsf{\Theta }}\left(z\right)=\tilde{k}\left(z\right)\ge \delta >0,\hspace{1em}\delta =constant\end{array}$$

(3)

The (known) functions ${\widehat{s}}_{i-1}\left(x\right)$ and ${\widehat{\gamma }}_{i-1}\left(x\right)$ are obtained from the knowledge of $v_{i-1}$.

It is to be noticed that almost all steady-state heat transfer problems can be described by Problem 2 and satisfy the six above conditions. In addition, it is remarkable that the function $\mathsf{\Theta }$ always admits an inverse. In fact, when $k$ represents the thermal conductivity (or any diffusion coefficient), the existence of a positive lower bound is always ensured.

When Equation (2) represents a heat transfer process, $\dot{q}$ is a heat source (per unit time and per unit volume), $k$ is the thermal conductivity, and $T$ is the temperature. The function $f$ represents the heat loss (per unit time and per unit area) from the body boundary.

The main subject of this work is concerned with real physical applications arising from heat transfer problems. In this way, the domain $\mathsf{\Omega }$ is assumed to have the cone property and that its boundary is composed by a finite number of smooth parts. In addition, it is assumed that $\dot{q}$ and $f$ are bounded and are piecewise continuous functions of $x$ for each given $T$. As physically expected, the function $T$ (temperature) will be a continuous function in $\mathsf{\Omega }$ and will satisfy Problem 2 almost everywhere.

It is to be noticed that, in general, real heat transfer problems are nonlinear and usually are simulated under a severe hypothesis in order to avoid complex mathematical tools for their solutions. With the procedure proposed in this work, a large class of nonlinear heat transfer problems can be treated with simple basic tools, such as those used for solving the linear problems. The proposed procedure allows a low-cost simulation for problems that, usually, require sophisticated methods.

In addition, since the solution will be obtained as the limit of a non-decreasing sequence (with a known upper bound), it is easy to verify the consistency of a given simulation. In other words, if the sequence does not behave as expected, the simulation must be stopped, because something must be corrected (probably $\alpha$ or $\overline{c}$ is not large enough).

2. Construction of the Non-Decreasing Sequence $\left[v_1,v_2,v_3\dots \right]$

Let us consider the problem below:

$$\begin{array}{l}\Delta \mathsf{\Phi }+Q=0\hspace{1em}in \mathsf{\Omega },\hspace{1em}Q=\widehat{Q}\left(x,\mathsf{\Phi }\right)\\ -\nabla \mathsf{\Phi }·n=\widehat{F}\left(x,\mathsf{\Phi }\right)\hspace{1em}on \partial \mathsf{\Omega }\end{array}$$

(4)

in which the functions $\widehat{F}\left(x,\mathsf{\Phi }\right)$ and $\widehat{Q}\left(x,\mathsf{\Phi }\right)$ have the following properties

• $\widehat{Q}\left(x,\mathsf{\Phi }\right)$ is a non-increasing function of $\mathsf{\Phi }$, for any $x\in \mathsf{\Omega }$;
• $\widehat{F}\left(x,\mathsf{\Phi }\right)$ is a strictly increasing function of $\mathsf{\Phi }$, for any $x\in \partial \mathsf{\Omega }$;
• $\underset{\mathsf{\Phi }\to +\infty }{\lim }\widehat{F}\left(x,\mathsf{\Phi }\right)=+\infty$, for any $x\in \partial \mathsf{\Omega }$;
• $\underset{\mathsf{\Phi }\to -\infty }{\lim }\widehat{F}\left(x,\mathsf{\Phi }\right)=- \infty$, for any $x\in \partial \mathsf{\Omega }$;
• $\widehat{F}\left(x,0\right)\le 0$, for any $x\in \partial \mathsf{\Omega }$;
• $\widehat{Q}\left(x,a\right)\ge 0$, for any $x\in \mathsf{\Omega }$, for any $a\in \left(-\infty ,\infty \right)$.

Thus, the sequence $\left[v_1,v_2,v_3\dots \right]$, whose elements are obtained from the solution of

$$\begin{array}{l}\Delta v_i-\overline{c}{v}_i+s_{i-1}=0\hspace{1em}in \mathsf{\Omega },\hspace{1em}i=1,2,\dots \\ -\nabla v_i·n=\alpha v_i-{\gamma }_{i-1}\hspace{1em}on \partial \mathsf{\Omega },\hspace{1em}i=1,2,\dots \\ s_{i-1}={\widehat{s}}_{i-1}\left(x\right)=\overline{c}{v}_{i-1}+\widehat{Q}\left(x,v_{i-1}\right)\\ {\gamma }_{i-1}={\widehat{\gamma }}_{i-1}\left(x\right)=\alpha v_{i-1}-\widehat{F}\left(x,v_{i-1}\right)\\ v_0\equiv 0\end{array}$$

(5)

is a convergent non-increasing sequence, provided $\overline{c}$ and $\alpha$ are sufficiently large constants in order to ensure that

$$\begin{array}{l}\overline{c}\left(v_i-v_{i-1}\right)\ge \widehat{Q}\left(x,v_{i-1}\right)-\widehat{Q}\left(x,v_i\right),\hspace{1em}i=1,2,\dots \\ \alpha \left(v_i-v_{i-1}\right)\ge \widehat{F}\left(x,v_i\right)-\widehat{F}\left(x,v_{i-1}\right),\hspace{1em}i=1,2,\dots \end{array}$$

(6)

As a first step, let us prove that $v_1$ is non-negative everywhere. Aiming for this, we write, from Equation (5), for $i=1$

$$\begin{array}{l}\Delta v_1-\overline{c}{v}_1+s_0=0\hspace{1em}in \mathsf{\Omega }\\ -\nabla v_1·n=\alpha v_1-{\gamma }_0\hspace{1em}on \partial \mathsf{\Omega }\\ s_0=\overline{c}{v}_0+\widehat{Q}\left(x,v_0\right)\\ {\gamma }_0=\alpha v_0-\widehat{F}\left(x,v_0\right)\\ v_0\equiv 0\end{array}$$

(7)

Now, let us assume that $v_1$ assumes its minimum at a point $x_{\ast }\in \mathsf{\Omega }$. Thus, there exists a small neighborhood of $x_{\ast }$ in which $\Delta v_1\ge 0$. Therefore, in this neighborhood, $-\overline{c}{v}_1+s_0\le 0$. Since $v_0\equiv 0$ and $\widehat{Q}\left(x,0\right)\ge 0$, we are able to conclude that, if $v_1$ has a local minimum in $\mathsf{\Omega }$, then this minimum is non-negative. On the other hand, if this minimum does not occur in $\mathsf{\Omega }$, we have that there exists a subset $\partial {\mathsf{\Omega }}_1^{-}\subseteq \partial \mathsf{\Omega }$, such that [21,22]

$$\underset{\partial \mathsf{\Omega }}{\inf }{v}_1=\underset{\partial {\mathsf{\Omega }}_1^{-}}{\inf }{v}_1=\underset{\mathsf{\Omega }}{\inf }{v}_1$$

(8)

On the subset $\partial {\mathsf{\Omega }}_1^{-}\subseteq \partial \mathsf{\Omega }$, from the boundary condition, we have

$$-\nabla v_1·n=\alpha v_1-\alpha v_0+\widehat{F}\left(x,v_0\right)\ge 0\hspace{1em}on \partial {\mathsf{\Omega }}_1^{-}$$

(9)

Therefore, since $\widehat{F}\left(x,0\right)\le 0$, we conclude that $v_1$ is non-negative everywhere. Therefore, $v_1\ge v_0\equiv 0$.

From Equation (5), we may write

$$\begin{array}{l}\Delta \left(v_{i+1}-v_i\right)-\overline{c}\left(v_{i+1}-v_i\right)+s_i-s_{i-1}=0\hspace{1em}in \mathsf{\Omega },\hspace{1em}i=1,2,\dots \\ -\nabla \left(v_{i+1}-v_i\right)·n=\alpha \left(v_{i+1}-v_i\right)-\left({\gamma }_i-{\gamma }_{i-1}\right)\hspace{1em}on \partial \mathsf{\Omega },\hspace{1em}i=1,2,\dots \\ s_i-s_{i-1}=\overline{c}\left(v_i-v_{i-1}\right)+\left(\widehat{Q}\left(x,v_i\right)-\widehat{Q}\left(x,v_{i-1}\right)\right)\\ {\gamma }_i-{\gamma }_{i-1}=\alpha \left(v_i-v_{i-1}\right)-\left(\widehat{F}\left(x,v_i\right)-\widehat{F}\left(x,v_{i-1}\right)\right)\\ v_0=0\hspace{1em}on \partial \mathsf{\Omega }\end{array}$$

(10)

Combining Equation (10) with Equation (6), we conclude that, if $v_i\ge v_{i-1}$ everywhere, then

$$\Delta \left(v_{i+1}-v_i\right)\le 0\hspace{1em}in \mathsf{\Omega }$$

(11)

and $v_{i+1}-v_i$ does not reach a minimum in $\mathsf{\Omega }$. Again, using Equation (6) and taking into account the boundary conditions, we conclude that

$$v_i\ge v_{i-1}\hspace{1em}in \mathsf{\Omega }\hspace{1em}\Rightarrow \hspace{1em}\underset{\partial \mathsf{\Omega }}{\inf }\left(v_{i+1}-v_i\right)=\underset{\mathsf{\Omega }}{\inf }\left(v_{i+1}-v_i\right)\ge 0$$

(12)

Since $v_1\ge v_0\hspace{1em}in \mathsf{\Omega }$, the above result enables us to conclude (non-decreasing sequence) that

$$\dots v_i\ge v_{i-1}\ge \dots v_2\ge v_1\ge 0\hspace{1em}in \mathsf{\Omega }$$

(13)

3. An Upper Bound for the Sequence

The solution of Equations (4), denoted by $\mathsf{\Phi }$, is an upper bound for the elements of the sequence $\left[v_1,v_2,v_3\dots \right]$. In order to prove this affirmation, let us combine Equations (4) and (5) to obtain

$$\begin{array}{l}\Delta \left(\mathsf{\Phi }-v_i\right)+Q\left(x,\mathsf{\Phi }\right)+\overline{c}{v}_i-\overline{c}{v}_{i-1}-\widehat{Q}\left(x,v_{i-1}\right)=0\hspace{1em}in \mathsf{\Omega }\\ -\nabla \left(\mathsf{\Phi }-v_i\right)·n=\widehat{F}\left(x,\mathsf{\Phi }\right)-\alpha v_i+\alpha v_{i-1}-\widehat{F}\left(x,v_{i-1}\right)\hspace{1em}on \partial \mathsf{\Omega }\end{array}$$

(14)

or, in a more convenient form,

$$\begin{array}{l}\Delta \left(\mathsf{\Phi }-v_i\right)+Q\left(x,\mathsf{\Phi }\right)-\widehat{Q}\left(x,v_i\right)+\overline{c}\left(v_i-v_{i-1}\right)+\widehat{Q}\left(x,v_i\right)-\widehat{Q}\left(x,v_{i-1}\right)=0\hspace{1em}in \mathsf{\Omega }\\ -\nabla \left(\mathsf{\Phi }-v_i\right)·n=\widehat{F}\left(x,\mathsf{\Phi }\right)-\widehat{F}\left(x,v_i\right)-\alpha \left(v_i-v_{i-1}\right)+\widehat{F}\left(x,v_i\right)-\widehat{F}\left(x,v_{i-1}\right)\hspace{1em}on \partial \mathsf{\Omega }\end{array}$$

(15)

Now, let us suppose that the difference $\mathsf{\Phi }-v_i$ assumes its minimum in $\mathsf{\Omega }$. In this case we must have $\Delta \left(\mathsf{\Phi }-v_i\right)\ge 0$ in a small neighborhood of the local minimum. Therefore, in this neighborhood, taking into account Equation (6), we have

$$Q\left(x,\mathsf{\Phi }\right)-\widehat{Q}\left(x,v_i\right)+\overline{c}\left(v_i-v_{i-1}\right)+\widehat{Q}\left(x,v_i\right)-\widehat{Q}\left(x,v_{i-1}\right)\le 0\hspace{1em}\Rightarrow \hspace{1em}\mathsf{\Phi }\ge v_i$$

(16)

On the other hand, if $\mathsf{\Phi }-v_i$ does not assume a minimum in $\mathsf{\Omega }$, we must have a subset $\partial {\mathsf{\Omega }}_i^{-}$, such that

$$\underset{\partial \mathsf{\Omega }}{\inf }\left(\mathsf{\Phi }-v_i\right)=\underset{\partial {\mathsf{\Omega }}_i^{-}}{\inf }\left(\mathsf{\Phi }-v_i\right)=\underset{\mathsf{\Omega }}{\inf }\left(\mathsf{\Phi }-v_i\right)$$

(17)

and, on this subset, taking into account the boundary conditions and Equation (6),

$$\widehat{F}\left(x,\mathsf{\Phi }\right)-\widehat{F}\left(x,v_i\right)\ge 0\hspace{1em}\Rightarrow \hspace{1em}\mathsf{\Phi }\ge v_i\hspace{1em}on \partial {\mathsf{\Omega }}_i^{-}$$

(18)

Thus, it is proven that $\mathsf{\Phi }\ge v_i$ in $\mathsf{\Omega }$. Clearly, it must be proven that (4) admits a solution and that $\mathsf{\Phi }\equiv v_{\infty }$. In the next section, with the aid of a variational formulation, we prove that the solution $\mathsf{\Phi }$ exists and is unique.

Since (from Equation (5))

$$\begin{array}{l}\Delta v_i-\overline{c}\left(v_i-v_{i-1}\right)+\widehat{Q}\left(x,v_{i-1}\right)=0\hspace{1em}in \mathsf{\Omega },\hspace{1em}i=1,2,\dots \\ -\nabla v_i·n=\alpha \left(v_i-v_{i-1}\right)+\widehat{F}\left(x,v_{i-1}\right)\hspace{1em}on \partial \mathsf{\Omega },\hspace{1em}i=1,2,\dots \\ v_0\equiv 0\end{array}$$

(19)

we have that

$$\begin{array}{l}\Delta v_{\infty }+\widehat{Q}\left(x,v_{\infty }\right)=0\hspace{1em}in \mathsf{\Omega }\\ -\nabla v_{\infty }·n=\widehat{F}\left(x,v_{\infty }\right)\hspace{1em}on \partial \mathsf{\Omega }\end{array}$$

(20)

Therefore, $\mathsf{\Phi }\equiv v_{\infty }$ (the limit of the sequence $\left[v_1,v_2,v_3\dots \right]$ consists of a solution of Problem 4, provided it admits a solution).

4. Existence and Uniqueness of the Solution of Equation (4)

The solution of Equation (4) is the field which minimizes the functional $I\left[\eta \right]$, defined by

$$\begin{array}{l}I\left[\eta \right]=\frac{1}{2}{\displaystyle \underset{\mathsf{\Omega }}{\int }{\left|\nabla \eta \right|}^{2}dV}-{\displaystyle \underset{\mathsf{\Omega }}{\int }\widehat{A}\left(x,\eta \right)dV}+{\displaystyle \underset{\partial \mathsf{\Omega }}{\int }\widehat{B}\left(x,\eta \right)dA},\hspace{1em}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\frac{\partial }{\partial \mathsf{\Phi }}\widehat{A}\left(x,\mathsf{\Phi }\right)=\widehat{Q}\left(x,\mathsf{\Phi }\right),\hspace{1em}\frac{\partial }{\partial \mathsf{\Phi }}\widehat{B}\left(x,\mathsf{\Phi }\right)=\widehat{F}\left(x,\mathsf{\Phi }\right)\end{array}$$

(21)

Since, for any $\theta \in \left(0,1\right)$, the following inequality holds:

$$\begin{array}{l}\left(1-\theta \right)\theta {\left|\nabla \left({\eta }_1-{\eta }_2\right)\right|}^2=\theta \nabla {\eta }_1·\nabla {\eta }_1\\ +\left(1-\theta \right)\nabla {\eta }_2·\nabla {\eta }_2-{\theta }^2\nabla {\eta }_1·\nabla {\eta }_1-{\left(1-\theta \right)}^2\nabla {\eta }_2·\nabla {\eta }_2-2\theta \left(1-\theta \right)\nabla {\eta }_1·\nabla {\eta }_2\ge 0\end{array}$$

(22)

we ensure that

$$\theta {\displaystyle \underset{\mathsf{\Omega }}{\int }{\left|\nabla {\eta }_1\right|}^{2}dV}+\left(1-\theta \right){\displaystyle \underset{\mathsf{\Omega }}{\int }{\left|\nabla {\eta }_2\right|}^{2}dV}-{\displaystyle \underset{\mathsf{\Omega }}{\int }{\left|\nabla \left\{\theta {\eta }_1+\left(1-\theta \right){\eta }_2\right\}\right|}^{2}dV}\ge 0,\hspace{1em}\theta \in \left(0,1\right)$$

(23)

In addition, since the derivative (with respect to $\mathsf{\Phi }$) of $-\widehat{A}\left(x,\mathsf{\Phi }\right)$ is non-negative and the derivative (with respect to $\mathsf{\Phi }$) of $\widehat{B}\left(x,\mathsf{\Phi }\right)$ is positive, we may write (convexity of the functional $I\left[\eta \right]$)

$$\theta I\left[{\eta }_1\right]+\left(1-\theta \right)I\left[{\eta }_2\right]-I\left[\theta {\eta }_1+\left(1-\theta \right){\eta }_2\right]\ge 0,\hspace{1em}\theta \in \left(0,1\right)$$

(24)

A necessary condition for ensuring the equality in Equation (22), for two different functions ${\eta }_1$ and ${\eta }_2$, is ${\eta }_1-{\eta }_2=constant$.

Nevertheless, from the definition of $\widehat{B}\left(x,\mathsf{\Phi }\right)$, this constant must be zero for ensuring the equality in (24). Therefore, we conclude that the functional $I\left[\eta \right]$ is strictly convex. This ensures the uniqueness of the minimum, provided it exists.

In order to prove the existence of the minimum, it is sufficient to show that the convex functional $I\left[\eta \right]$ is such that (coerciveness) [23]

$$\underset{\gamma \to \infty }{\lim }\frac{I\left[\gamma \eta \right]}{\gamma }=\infty ,\hspace{1em}\Vert \eta \Vert =1$$

(25)

in which $\Vert \eta \Vert$ is the usual norm of $H^1\left(\mathsf{\Omega }\right)$ [24].

Taking into account that $\Vert \eta \Vert =1$ and considering the behavior of $\widehat{Q}\left(x,\mathsf{\Phi }\right)$ and of $\widehat{F}\left(x,\mathsf{\Phi }\right)$, the above limit may be represented as follows:

$$\begin{array}{l}\underset{\gamma \to \infty }{\lim }\frac{I\left[\gamma \eta \right]}{\gamma }=\underset{\gamma \to \infty }{\lim }\left\{\frac{1}{2\gamma }{\displaystyle \underset{\mathsf{\Omega }}{\int }{\left|\nabla \gamma \eta \right|}^{2}dV}-{\displaystyle \underset{\mathsf{\Omega }}{\int }\frac{\widehat{A}\left(x,\gamma \eta \right)}{\gamma }dV}+{\displaystyle \underset{\partial \mathsf{\Omega }}{\int }\frac{\widehat{B}\left(x,\gamma \eta \right)}{\gamma }dA}\right\}\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}=\underset{\gamma \to \infty }{\lim }\left\{\frac{\gamma }{2}{\displaystyle \underset{\mathsf{\Omega }}{\int }{\left|\nabla \eta \right|}^{2}dV}-{\displaystyle \underset{\mathsf{\Omega }}{\int }\eta \widehat{Q}\left(x,\gamma \eta \right)dV}+{\displaystyle \underset{\partial \mathsf{\Omega }}{\int }\eta \widehat{F}\left(x,\gamma \eta \right)dA}\right\}=\infty \end{array}$$

(26)

ensuring the coerciveness (solution existence).

At this point, it is convenient to show that the minimization of $I\left[\eta \right]$ corresponds, in a weak sense, to the solution of (4). Aiming to this, based on Equation (21), let us calculate the first variation of $I\left[\eta \right]$ as follows

$$\delta I\left[\eta \right]={\displaystyle \underset{\mathsf{\Omega }}{\int }\nabla \eta ·\nabla \delta \eta dV}-{\displaystyle \underset{\mathsf{\Omega }}{\int }\left\{\frac{\partial }{\partial \eta }\widehat{A}\left(x,\eta \right)\right\}\delta \eta dV}+{\displaystyle \underset{\partial \mathsf{\Omega }}{\int }\left\{\frac{\partial }{\partial \eta }\widehat{B}\left(x,\eta \right)\right\}\delta \eta dA}$$

(27)

or, in a more convenient way, as

$$\delta I\left[\eta \right]={\displaystyle \underset{\mathsf{\Omega }}{\int }\nabla \eta ·\nabla \delta \eta dV}-{\displaystyle \underset{\mathsf{\Omega }}{\int }\widehat{Q}\left(x,\mathsf{\Phi }\right)\delta \eta dV}+{\displaystyle \underset{\partial \mathsf{\Omega }}{\int }\widehat{f}\left(x,\mathsf{\Phi }\right)\delta \eta dA}$$

(28)

Since [25]

$${\displaystyle \underset{\mathsf{\Omega }}{\int }\nabla \eta ·\nabla \delta \eta dV}={\displaystyle \underset{\mathsf{\Omega }}{\int }\left(\nabla ·\left(\delta \eta \nabla \eta \right)-\nabla ·\left(\nabla \eta \right)\delta \eta \right)dV}=-{\displaystyle \underset{\mathsf{\Omega }}{\int }\nabla ·\left(\nabla \eta \right)\delta \eta dV}+{\displaystyle \underset{\partial \mathsf{\Omega }}{\int }\delta \eta \nabla \eta ·ndA}$$

(29)

we may write

$$\delta I\left[\eta \right]=-{\displaystyle \underset{\mathsf{\Omega }}{\int }\nabla ·\left(\nabla \eta \right)\delta \eta dV}+{\displaystyle \underset{\partial \mathsf{\Omega }}{\int }\delta \eta \nabla \eta ·ndA}-{\displaystyle \underset{\mathsf{\Omega }}{\int }\widehat{Q}\left(x,\eta \right)\delta \eta dV}+{\displaystyle \underset{\partial \mathsf{\Omega }}{\int }\widehat{f}\left(x,\eta \right)\delta \eta dA}$$

(30)

Thus, taking into account that $\delta \eta$ is arbitrary, $\delta I\left[\eta \right]=0$ gives rise to (Euler–Lagrange equation and natural boundary condition)

$$\begin{array}{l}\nabla ·\left(\nabla \eta \right)+\widehat{Q}\left(x,\eta \right)=0\hspace{1em}in \mathsf{\Omega }\\ \nabla \eta ·n+\widehat{F}\left(x,\eta \right)=0\hspace{1em}on \partial \mathsf{\Omega }\end{array}$$

(31)

proving the equivalence between (4) and the minimization of the functional $I\left[\eta \right]$ (note that $\nabla ·\nabla \equiv \Delta \equiv div\left(grad\right)$ represents the Laplacian).

In other words, we have that, in a weak sense, (4) admits one, and only one, solution. In fact, this minimum will belong to $H^2\left(\mathsf{\Omega }\right)$. Since the domain $\mathsf{\Omega }$ belongs to ${ℝ}^3$ and has the cone property, this minimum will be continuous in $\mathsf{\Omega }$ (Sobolev imbedding theorem [24]).

It is natural to ask about the use of $I\left[\eta \right]$ for reaching the solution $\mathsf{\Phi }$. Surely, the minimization of $I\left[\eta \right]$ is a valid and reliable procedure, but, the procedure proposed in this work is also reliable and involves much more simple tools.

5. Construction of the Solution of Problem 2

Since $\tilde{k}\left(T\right)\ge \delta >0$, the knowledge of the elements of the sequence $\left[v_1,v_{2,}{v}_3\dots \right]$, as well as its limit, enable us to determine, in a unique way, the elements of the sequence $\left[p_1,p_{2,}{p}_3\dots \right]$, as well as its limit. In other words, the knowledge of $\mathsf{\Phi }$ determines, in a unique way, the solution of Problem 2.

The functional relationship between $\mathsf{\Phi }$ (solution of (4)) and $T$(solution of (2)) depends on the choice of the function $\tilde{\mathsf{\Theta }}\left(z\right)$. There are infinite possible choices. For instance, when the physically considered range for the temperature is the interval $T\in \left[a,b\right]$, the function $\mathsf{\Theta }$ may be defined as [26]

$$\mathsf{\Theta }=\tilde{\mathsf{\Theta }}\left(T\right)=\{\begin{array}{c}\left\{\widehat{k}\left(a\right)\right\}T,\hspace{1em}if\hspace{1em}T<a\\ {\displaystyle \underset{a}{\overset{T}{\int }}\widehat{k}\left(\xi \right)d\xi }+\left\{\widehat{k}\left(a\right)\right\}a,\hspace{1em}if\hspace{1em}a\le T\le b\\ {\displaystyle \underset{a}{\overset{b}{\int }}\widehat{k}\left(\xi \right)d\xi }+\left\{\widehat{k}\left(a\right)\right\}a+\left\{\widehat{k}\left(b\right)\right\}\left(T-b\right),\hspace{1em}if\hspace{1em}b<T\end{array}$$

(32)

It is to be noticed that, for any real heat transfer process, $a$ and $b$ are positive valued constants, provided $T$ represents an absolute temperature.

In order to illustrate the calculation of $T$ from the knowledge of $\mathsf{\Phi }$, let us consider a particular case in which, within the physically admissible range,

$$k=\widehat{k}\left(T\right)=k_a+\frac{T-a}{b-a}\left(k_b-k_a\right),\hspace{1em}{k}_a\ge \delta >0,\hspace{1em}{k}_b\ge \delta >0$$

(33)

In this case

$$\mathsf{\Phi }=\tilde{\mathsf{\Phi }}\left(T\right)=\{\begin{array}{c}{k}_{a}T,\hspace{1em}if\hspace{1em}T<a\\ k_a\left(T-a\right)+\frac{k_b-k_a}{b-a}\frac{{\left(T-a\right)}^2}{2}+k_{a}a,\hspace{1em}if\hspace{1em}a\le T\le b\\ \left(k_a-k_b\right)\frac{a+b}{2}+k_{b}T,\hspace{1em}if\hspace{1em}b<T\end{array}$$

(34)

and the temperature is given by

$$T=\tilde{T}\left(\mathsf{\Phi }\right)=\{\begin{array}{c}\underset{}{\frac{\mathsf{\Phi }}{k_a},\hspace{1em}if\hspace{1em}\mathsf{\Phi }<k_{a}a}\\ \begin{array}{l}\frac{a{k}_b-b{k}_a}{k_b-k_a}+\left(\frac{b-a}{k_b-k_a}\right)\sqrt{{\left(k_a\right)}^2-2\left(\frac{k_b-k_a}{b-a}\right)\left(k_{a}a-\mathsf{\Phi }\right)},\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}if\hspace{1em}{k}_{a}a\le \mathsf{\Phi }\le k_{a}a+\left(k_a+k_b\right)\frac{b-a}{2}\end{array}\\ \stackrel{}{\frac{\mathsf{\Phi }}{k_b}-\frac{\left(k_a-k_b\right)\left(a+b\right)}{2k_b},\hspace{1em}if\hspace{1em}{k}_{a}a+\left(k_a+k_b\right)\frac{b-a}{2}<\mathsf{\Phi }}\end{array}$$

(35)

The linear relationship between the thermal conductivity $k$ and the temperature $T$ is the most usual for describing heat transfer processes. Other relationships may be used as, for instance, the one proposed in [26] and used in [27].

Figure 2 illustrates a particular case of Equation (34).

6. The Constants $\overline{c}$ and $\alpha$

The constants $\overline{c}$ and $\alpha$ must satisfy Equation (6). Since (4) admits a unique solution and $\mathsf{\Phi }\ge v_i$, the existence of the constants $\overline{c}$ and $\alpha$ is ensured.

In order to provide a sufficiently (not necessarily) large value for each of these two constants, it is convenient to establish an upper bound for $\mathsf{\Phi }$ before solving the problem.

Aiming at this, let us define the function $w$, (any) solution of the following partial differential equation:

$$\Delta w+\underset{x\in \mathsf{\Omega }}{\sup }\left[\widehat{Q}\left(x,0\right)\right]=0\hspace{1em}in \mathsf{\Omega }$$

(36)

For instance, an always possible choice (not a unique choice and, probably, often not the best) for $w$ is

$$w=-\frac{1}{6}\left(x·x\right)\left\{\underset{x\in \mathsf{\Omega }}{\sup }\left[\widehat{Q}\left(x,0\right)\right]\right\},\hspace{1em}x\in \mathsf{\Omega }$$

(37)

Combining Equation (4) with Equation (36), we have

$$\Delta \left(\mathsf{\Phi }-w\right)+\widehat{Q}\left(x,\mathsf{\Phi }\right)-\underset{x\in \mathsf{\Omega }}{\sup }\left[\widehat{Q}\left(x,0\right)\right]=0\hspace{1em}\Rightarrow \hspace{1em}\Delta \left(\mathsf{\Phi }-w\right)\ge 0\hspace{1em}in \mathsf{\Omega }$$

(38)

Therefore, there exists a subset $\partial {\mathsf{\Omega }}^{+}\subseteq \partial \mathsf{\Omega }$ in which $\nabla \left(\mathsf{\Phi }-w\right)·n\ge 0$, such that

$$\underset{\partial {\mathsf{\Omega }}^{+}}{\sup }\left(\mathsf{\Phi }-w\right)=\underset{\partial \mathsf{\Omega }}{\sup }\left(\mathsf{\Phi }-w\right)=\underset{\mathsf{\Omega }}{\sup }\left(\mathsf{\Phi }-w\right)$$

(39)

Taking into account the boundary condition, we may write

$$-\nabla \left(\mathsf{\Phi }-w\right)·n=\widehat{F}\left(x,\mathsf{\Phi }\right)+\nabla w·n\hspace{1em}on \partial \mathsf{\Omega }\hspace{1em}\Rightarrow \hspace{1em}\widehat{F}\left(x,\mathsf{\Phi }\right)+\nabla w·n\le 0\hspace{1em}on \partial {\mathsf{\Omega }}^{+}$$

(40)

Since,

$$\underset{\partial {\mathsf{\Omega }}^{+}}{\sup }\left(\mathsf{\Phi }-w\right)=\underset{\partial \mathsf{\Omega }}{\sup }\left(\mathsf{\Phi }-w\right)=\underset{\mathsf{\Omega }}{\sup }\left(\mathsf{\Phi }-w\right)$$

(41)

we are able to write

$$\underset{\mathsf{\Omega }}{\sup }\mathsf{\Phi }\le \underset{\partial {\mathsf{\Omega }}^{+}}{\sup }\left(\mathsf{\Phi }-w\right)+\underset{\mathsf{\Omega }}{\sup }w\le \underset{\partial {\mathsf{\Omega }}^{+}}{\sup }\mathsf{\Phi }-\underset{\partial {\mathsf{\Omega }}^{+}}{\inf }w+\underset{\mathsf{\Omega }}{\sup }w$$

(42)

in which $w$ is a known function.

Therefore, any constant ${\mathsf{\Phi }}^{+}$ satisfying the inequality

$$\widehat{F}\left(x,{\mathsf{\Phi }}^{+}\right)+\nabla w·n\le 0,\hspace{1em}x\in \partial {\mathsf{\Omega }}^{+}$$

(43)

will be greater than the supremum of $\mathsf{\Phi }$ on $\partial {\mathsf{\Omega }}^{+}$. Therefore,

$$\underset{\mathsf{\Omega }}{\sup }\mathsf{\Phi }\le {\mathsf{\Phi }}^{+}-\underset{\partial {\mathsf{\Omega }}^{+}}{\inf }w+\underset{\mathsf{\Omega }}{\sup }w$$

(44)

Once known, an upper bound for $\mathsf{\Phi }$, the constants $\overline{c}$ and $\alpha$ may be chosen from the following inequalities (sufficient conditions, but not necessary):

$$\begin{array}{l}\overline{c}\left(b-a\right)\ge \widehat{Q}\left(x,a\right)-\widehat{Q}\left(x,b\right),\hspace{1em}\underset{\mathsf{\Omega }}{\sup }\mathsf{\Phi }>b>a>0\\ \alpha \left(b-a\right)\ge \widehat{F}\left(x,b\right)-\widehat{F}\left(x,a\right),\hspace{1em}\underset{\mathsf{\Omega }}{\sup }\mathsf{\Phi }>b>a>0\end{array}$$

(45)

It must be highlighted that $\overline{c}$ and $\alpha$ satisfying (45) ensure the convergence of the non-decreasing sequence $\left[v_1,v_2,v_3\dots \right]$. Nevertheless, in general, this non-decreasing convergent sequence may be reached with lower values for $\overline{c}$ and $\alpha$.

When $\widehat{Q}\left(x,\mathsf{\Phi }\right)$ and $\widehat{F}\left(x,\mathsf{\Phi }\right)$ are differentiable with respect to $\mathsf{\Phi }$, within the interval of interest, (45) may be rewritten as follows:

$$\begin{array}{l}\overline{c}\ge -\frac{\partial }{\partial \xi }\left\{\widehat{Q}\left(x,\xi \right)\right\},\hspace{1em}\underset{\mathsf{\Omega }}{\sup }\mathsf{\Phi }>\xi >0,\hspace{1em}for any x\in \mathsf{\Omega }\\ \alpha \ge \frac{\partial }{\partial \xi }\left\{\widehat{F}\left(x,\xi \right)\right\},\hspace{1em}\underset{\mathsf{\Omega }}{\sup }\mathsf{\Phi }>\xi >0,\hspace{1em}for any x\in \partial \mathsf{\Omega }\end{array}$$

(46)

Clearly, from Equation (45), if $\widehat{Q}\left(x,\mathsf{\Phi }\right)$ does not depend on $\mathsf{\Phi }$, then $\overline{c}$ may be zero.

7. Conclusions

This work provides a powerful and simple tool for simulating any steady-state heat transfer phenomenon and any other process, described by a nonlinear second order partial differential equation such as Equation (2), even when subjected to nonlinear boundary conditions (for instance, arising from thermal radiation heat transfer from/to the body) and involving temperature-dependent source terms (even when this dependence is not linear).

The results presented here generalize common situations in which the source term is assumed to be known (i.e., did not depend on the temperature) and where the boundary condition is the classical linear Newton law of cooling.

The proposed procedure for constructing the solution may be used together with any numerical scheme available for a linear problem such as Problem 1, making a large class of nonlinear problems accessible to any student (of engineering, physics, or mathematics). An example is presented in Appendix A.

The conditions imposed for ensuring the convergence are sufficient conditions, but are not necessary. In other words, the procedure may work even outside of the ensured range.

The disadvantage of the proposed procedure is the fact of not knowing the optimal value of the constant $\alpha$ (the lower one which ensures convergence). On the other hand, if the convergence is reached for a given $\alpha =\overline{\alpha }$, it will be reached for any $\alpha$ greater than $\overline{\alpha }$.

It is to be noticed that linear problems, in particular those involving linear boundary conditions, admit an equivalent minimum principle represented by a quadratic functional. Nonlinear problems are not equivalent to the minimization of quadratic functionals.