# Métodos Numéricos Aplicados à Transferência de Calor

## Introdução

O normal do nosso fluxo de trabalho é iniciar importando as bibliotecas que vamos utilizar no decorrer do material.
Um maior detalhamento sobre elas já foi feito na aula anterior, de modo que agora podemos utilizar diretamente:


```python
import handcalcs.render
import matplotlib.pyplot as plt
import numpy as np
import plotly.express as px
import plotly.figure_factory as ff
import plotly.graph_objects as go
import plotly.io as pio
from tqdm.notebook import tqdm
```

O bloco a seguir é opcional, ele vai alterar o estilo padrão de nossas figuras, e aumentar um pouco o seu tamanho, melhorando a apresentação em nossa aula:


```python
# Definindo um novo estilo para as figuras [opcional]
# Esse bloco modifica alguns dos valores padrões para

plt.rcdefaults()

# https://matplotlib.org/3.1.0/gallery/style_sheets/style_sheets_reference.html
plt.style.use("ggplot")

# https://matplotlib.org/3.1.1/tutorials/introductory/customizing.html
plt.rcParams.update({"figure.dpi": 100, "figure.figsize": (6, 6)})

px.defaults.template = "ggplot2"
px.defaults.height = 600
pio.templates.default = "ggplot2"
```

## Exercícios Resolvidos

### Convecção e difusão transiente unidimensional

Resolver a EDP:

\begin{equation}
    \dfrac{\partial T}{\partial t} = \alpha \dfrac{\partial^2 T}{\partial x^2} - u\dfrac{\partial T}{\partial x}, \quad 0\leq x \leq 1 ; 0\leq t \leq 8
\end{equation}

Condições de contorno:

\begin{equation}
    T(0,t)=T(1,t)=0
\end{equation}

Condição inicial:

\begin{equation}
    T(x,0) =  1 - ( 10 x - 1 )^2 \quad \text{ se $0 \leq x \leq 0,2$}, \quad \text{ senão } T(x,0) = 0
\end{equation}


Discretizando com as derivadas espaciais numa representação por diferença central e a derivada temporal com diferença ascendente:

\begin{equation}
\dfrac{T_{i,n+1}-T_{i,n}}{\Delta t}=\alpha \dfrac{T_{i-1,n}-2T_{i,n}+T_{i+1,n}}{(\Delta x)^2} -u\dfrac{T_{i+1,n}-T_{i-1,n}}{2\Delta x}, \quad 1 \le i \le I - 2, \quad n > 0,
\end{equation}

\begin{equation}
T_{i=0,n} = 0,
\end{equation}
\begin{equation}
T_{i=I-1,n} = 0.
\end{equation}

Agora devemos isolar a incógnita do nosso problema: o termo $T_{i,n+1}$. Perceba que todos os termos à direita são conhecidos, e usamos essa informação para avançar progressivamente no tempo:

\begin{equation}
T_{i,n+1} = T_{i,n} + \Delta t \left( \alpha \dfrac{T_{i-1,n}-2T_{i,n}+T_{i+1,n}}{(\Delta x)^2} -u\dfrac{T_{i+1,n}-T_{i-1,n}}{2\Delta x} \right), \quad 1 \le i \le I - 2, \quad n > 0,
\end{equation}

Veja como podemos programar a solução para o problema, e note, principalmente, como ela se assemelha muito com a notação discretizada da equação acima. Aqui está o código:


```python
def equação_exemplo_1(coord_x, coord_t, alpha, u):
    
    # Condição inicial
    T = np.zeros(shape = (coord_x.size, coord_t.size))
    
    for i, x in enumerate(coord_x):
        if x <= 0.2:
            T[i] = 1. - (10. * x - 1) ** 2.
    
    # Condições de contorno
    T[0, :] = 0.0
    T[-1, :] = 0.0
    
    # Passo de tempo e resolução da malha
    dt = coord_t[1] - coord_t[0]
    dx = coord_x[1] - coord_x[0]

    # Aqui resolve-se a equação
    for n in tqdm(range(0, coord_t.size - 1)):
        for i in range(1, coord_x.size - 1):
                
                T[i, n + 1] = T[i, n] + dt * (
                    alpha * (T[i - 1, n] - 2.0 * T[i, n] + T[i + 1, n]) / dx**2.0
                    - u * (T[i + 1, n] - T[i - 1, n]) / (2.0 * dx)
                )
    
    return T
```

Hora de executar o cálculo utilizando a função que definimos no bloco anterior:


```python
coord_x = np.linspace(0., 1., num=101)
coord_t = np.linspace(0., 8., num=2001)

T = equação_exemplo_1(
    coord_x,
    coord_t,
    alpha = 0.001,
    u = 0.08
)
```

Por fim, exemplificamos a exibição dos resultados com `plotly`:


```python
fig = go.Figure()

for n in range(0, coord_t.size, 500):

    fig.add_trace(
        go.Scatter(
            x=coord_x, y=T[:,n],
            mode='lines',
            name=f't={coord_t[n]}'
        )
    )  

fig.update_layout(
    title='Convecção e difusão transiente',
    xaxis_title='x',
    yaxis_title='Temperatura'
)
    
fig.show()
```


```python
fig = px.imshow(
    T.T,
    x = coord_x,
    y = coord_t,
    color_continuous_scale = 'RdBu_r',
    title = "Temperatura",
    labels = dict(x = "x", y = "tempo"),
    aspect = "auto",
    origin = "lower"
)
fig.show()
```

### Escoamento em Cavidade com Transferência de Calor

Aqui está o sistema de equações diferenciais: a equação da continuidade, duas equações para os componentes de velocidade $u,v$ e uma equação para a temperatura $\Theta$:

$$ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 $$

$$\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y} = - \frac{\partial p}{\partial x}+ \frac{1}{Re} \left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2} \right) $$

$$\frac{\partial v}{\partial t}+u\frac{\partial v}{\partial x}+v\frac{\partial v}{\partial y} = - \frac{\partial p}{\partial y}+ \frac{1}{Re} \left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}\right) + Ri ~ \Theta $$

$$ \frac{\partial \Theta}{\partial t} + u\frac{\partial \Theta}{\partial x} + v\frac{\partial \Theta}{\partial y} = \frac{1}{Re ~ Pr} \left(\frac{\partial^2 \Theta}{\partial x^2}+\frac{\partial^2 \Theta}{\partial y^2} \right) $$

* Equações discretas:

Primeiro, vamos discretizar a equação do momento para $u$, da seguinte maneira:

$$
\begin{split}
& \frac{u_{i,j}^{n+1}-u_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{u_{i+1,j}^{n}-u_{i-1,j}^{n}}{2 \Delta x}+v_{i,j}^{n}\frac{u_{i,j+1}^{n}-u_{i,j-1}^{n}}{2\Delta y} = \\ 
& \qquad -\frac{p_{i+1,j}^{n}-p_{i-1,j}^{n}}{2\Delta x}+\frac{1}{Re}\left(\frac{u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}}{\Delta x^2}+\frac{u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}}{\Delta y^2}\right)
\end{split}
$$

Da mesma forma para a equação do momento para $v$:

$$
\begin{split}
&\frac{v_{i,j}^{n+1}-v_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{v_{i+1,j}^{n}-v_{i-1,j}^{n}}{2\Delta x}+v_{i,j}^{n}\frac{v_{i,j+1}^{n}-v_{i,j-1}^{n}}{2\Delta y} = \\
& \qquad - \frac{p_{i,j+1}^{n}-p_{i,j-1}^{n}}{2\Delta y}
+\frac{1}{Re}\left(\frac{v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}}{\Delta x^2}+\frac{v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}}{\Delta y^2}\right) + Ri ~ \Theta_{i,j}
\end{split}
$$

A obtenção da equação de pressão-Poisson discretizada está fora do escopo dessa aula, por ser parte da disciplina de mecânica dos fluidos computacional, mais detalhes podem ser obtidos nos passos 10 e 11 do curso [CFD com Python](https://github.com/fschuch/CFDPython-BR). De qualquer maneira, a equação pode ser escrita assim:

$$
\begin{split}
& \frac{p_{i+1,j}^{n}-2p_{i,j}^{n}+p_{i-1,j}^{n}}{\Delta x^2}+\frac{p_{i,j+1}^{n}-2p_{i,j}^{n}+p_{i,j-1}^{n}}{\Delta y^2} = \\
& \qquad \left[ \frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right) -\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} - 2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x} - \frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]
\end{split}
$$

Por fim, a equação para temperatura:

$$
\begin{split}
&\frac{\Theta_{i,j}^{n+1}-\Theta_{i,j}^{n}}{\Delta t}+u_{i,j}^{n}\frac{\Theta_{i+1,j}^{n}-\Theta_{i-1,j}^{n}}{2\Delta x}+v_{i,j}^{n}\frac{\Theta_{i,j+1}^{n}-\Theta_{i,j-1}^{n}}{2\Delta y} = \\
& \qquad + \frac{1}{Re ~ Pr}\left(\frac{\Theta_{i+1,j}^{n}-2\Theta_{i,j}^{n}+\Theta_{i-1,j}^{n}}{\Delta x^2}+\frac{\Theta_{i,j+1}^{n}-2\Theta_{i,j}^{n}+\Theta_{i,j-1}^{n}}{\Delta y^2}\right)
\end{split}
$$

**Exercitando:** Você pode escrever essas equações em suas próprias anotações, manualmente, seguindo cada termo mentalmente à medida que as escreve.

Como antes, vamos reorganizar as equações da maneira que as iterações devem proceder no código. Primeiro, as equações de momento para a velocidade no passo de tempo subsequente.

A equação do momento na direção de $u$:

$$
\begin{split}
u_{i,j}^{n+1} = u_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{2\Delta x} \left(u_{i+1,j}^{n}-u_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{2\Delta y} \left(u_{i,j+1}^{n}-u_{i,j-1}^{n}\right) \\
& - \frac{\Delta t}{2\Delta x} \left(p_{i+1,j}^{n}-p_{i-1,j}^{n}\right) \\
& + \frac{1}{Re} \left(\frac{\Delta t}{\Delta x^2} \left(u_{i+1,j}^{n}-2u_{i,j}^{n}+u_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(u_{i,j+1}^{n}-2u_{i,j}^{n}+u_{i,j-1}^{n}\right)\right)
\end{split}
$$


```python
def cavidade_u(u, v, dt, dx, dy, p, re):
    return (
        u[1:-1, 1:-1]
        - u[1:-1, 1:-1] * dt / (2*dx) * (u[2:, 1:-1] - u[0:-2, 1:-1])
        - v[1:-1, 1:-1] * dt / (2*dy) * (u[1:-1, 2:] - u[1:-1, 0:-2])
        - dt / (2 * dx) * (p[2:, 1:-1] - p[0:-2, 1:-1])
        + (1.0 / re)
        * (
            dt / dx ** 2 * (u[2:, 1:-1] - 2 * u[1:-1, 1:-1] + u[0:-2, 1:-1])
            + dt / dy ** 2 * (u[1:-1, 2:] - 2 * u[1:-1, 1:-1] + u[1:-1, 0:-2])
        )
    )
```

A equação do momento na direção de $v$:

$$
\begin{split}
v_{i,j}^{n+1} = v_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{2\Delta x} \left(v_{i+1,j}^{n}-v_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{2\Delta y} \left(v_{i,j+1}^{n}-v_{i,j-1}^{n})\right) \\
& - \frac{\Delta t}{2\Delta y} \left(p_{i,j+1}^{n}-p_{i,j-1}^{n}\right) \\
& + \frac{1}{Re} \left(\frac{\Delta t}{\Delta x^2} \left(v_{i+1,j}^{n}-2v_{i,j}^{n}+v_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(v_{i,j+1}^{n}-2v_{i,j}^{n}+v_{i,j-1}^{n}\right)\right) + \Delta t ~ Ri ~ \Theta_{i,j}
\end{split}
$$


```python
def cavidade_v(u, v, theta, dt, dx, dy, p, re, ri):
    return (
        v[1:-1, 1:-1]
        - u[1:-1, 1:-1] * dt / (2*dx) * (v[2:, 1:-1] - v[0:-2, 1:-1])
        - v[1:-1, 1:-1] * dt / (2*dy) * (v[1:-1, 2:] - v[1:-1, 0:-2])
        - dt / (2 * dy) * (p[1:-1, 2:] - p[1:-1, 0:-2])
        + (1.0 / re)
        * (
            dt / dx ** 2 * (v[2:, 1:-1] - 2 * v[1:-1, 1:-1] + v[0:-2, 1:-1])
            + dt / dy ** 2 * (v[1:-1, 2:] - 2 * v[1:-1, 1:-1] + v[1:-1, 0:-2])
        )
    ) + ri * dt * theta[1:-1, 1:-1]
```

Agora, reorganizamos a equação de pressão-Poisson:

$$
\begin{split}
p_{i,j}^{n} = & \frac{\left(p_{i+1,j}^{n}+p_{i-1,j}^{n}\right) \Delta y^2 + \left(p_{i,j+1}^{n}+p_{i,j-1}^{n}\right) \Delta x^2}{2\left(\Delta x^2+\Delta y^2\right)} \\
& -\frac{\Delta x^2\Delta y^2}{2\left(\Delta x^2+\Delta y^2\right)} \\
& \times \left[\frac{1}{\Delta t}\left(\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}+\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right)-\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x}\frac{u_{i+1,j}-u_{i-1,j}}{2\Delta x} -2\frac{u_{i,j+1}-u_{i,j-1}}{2\Delta y}\frac{v_{i+1,j}-v_{i-1,j}}{2\Delta x}-\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\frac{v_{i,j+1}-v_{i,j-1}}{2\Delta y}\right]
\end{split}
$$


```python
def build_up_b(u, v, dt, dx, dy):
    """Constrói o termo entre chaves na equação acima"""
    
    b = np.zeros_like(u)

    b[1:-1, 1:-1] = (
        1
        / dt
        * (
            (u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dx)
            + (v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dy)
        )
        - ((u[2:, 1:-1] - u[0:-2, 1:-1]) / (2 * dx)) ** 2
        - 2
        * (
            (u[1:-1, 2:] - u[1:-1, 0:-2])
            / (2 * dy)
            * (v[2:, 1:-1] - v[0:-2, 1:-1])
            / (2 * dx)
        )
        - ((v[1:-1, 2:] - v[1:-1, 0:-2]) / (2 * dy)) ** 2
    )

    return b


def pressure_poisson(u, v, dt, dx, dy, p, nit=100):

    b = build_up_b(u, v, dt, dx, dy)

    for q in range(nit):
        pn = p.copy()
        p[1:-1, 1:-1] = (
            (
                (pn[2:, 1:-1] + pn[0:-2, 1:-1]) * dy ** 2
                + (pn[1:-1, 2:] + pn[1:-1, 0:-2]) * dx ** 2
            )
            / (2 * (dx ** 2 + dy ** 2))
            - dx ** 2 * dy ** 2 / (2 * (dx ** 2 + dy ** 2)) * b[1:-1, 1:-1]
        )

        p[-1, :] = p[-2, :]  # dp/dx = 0 at x = 2
        p[:, 0] = p[:, 1]  # dp/dy = 0 at y = 0
        p[0, :] = p[1, :]  # dp/dx = 0 at x = 0
        p[:, -1] = 0  # p = 0 at y = 2

    return p
```

Quase lá! Só falta a equação para temperatura:

$$
\begin{split}
\Theta_{i,j}^{n+1} = \Theta_{i,j}^{n} & - u_{i,j}^{n} \frac{\Delta t}{2\Delta x} \left(\Theta_{i+1,j}^{n}-\Theta_{i-1,j}^{n}\right) - v_{i,j}^{n} \frac{\Delta t}{2\Delta y} \left(\Theta_{i,j+1}^{n}-\Theta_{i,j-1}^{n})\right) \\
& + \frac{1}{Re ~ Pr} \left(\frac{\Delta t}{\Delta x^2} \left(\Theta_{i+1,j}^{n}-2\Theta_{i,j}^{n}+\Theta_{i-1,j}^{n}\right) + \frac{\Delta t}{\Delta y^2} \left(\Theta_{i,j+1}^{n}-2\Theta_{i,j}^{n}+\Theta_{i,j-1}^{n}\right)\right)
\end{split}
$$


```python
def cavidade_theta(u, v, theta, dt, dx, dy, p, re, pr):
    return (
        theta[1:-1, 1:-1]
        - u[1:-1, 1:-1] * dt / (2*dx) * (theta[2:, 1:-1] - theta[0:-2, 1:-1])
        - v[1:-1, 1:-1] * dt / (2*dy) * (theta[1:-1, 2:] - theta[1:-1, 0:-2])
        + (1.0 / re / pr)
        * (
            dt / dx ** 2 * (theta[2:, 1:-1] - 2 * theta[1:-1, 1:-1] + theta[0:-2, 1:-1])
            + dt
            / dy ** 2
            * (theta[1:-1, 2:] - 2 * theta[1:-1, 1:-1] + theta[1:-1, 0:-2])
        )
    )
```

A condição inicial é $u, v, p, \Theta = 0 $ em todos os lugares, e as condições de contorno são:

$u=1$ e $\Theta=0$ em $y=2$ (a "tampa");

$u, v = 0$ e $\Theta=1$ nas fronteiras restantes;

$\frac{\partial p}{\partial y}=0$ em $y=0$;

$p=0$ em $y=2$

$\frac{\partial p}{\partial x}=0$ em $x=0$ e $x=2$

Agora todos juntos em uma nova função para efetivamente resolver o escoamento em uma cavidade:


```python
def cavidade(x, y, t, re, ri, pr):

    # Condição inicial
    u = np.zeros((x.size, y.size))
    v = np.zeros((x.size, y.size))
    p = np.zeros((x.size, y.size))
    theta = np.zeros((x.size, y.size))

    # Passo de tempo e resolução da malha
    dt = t[1] - t[0]
    dx = x[1] - x[0]
    dy = y[1] - y[0]

    # Laço temporal
    for n in tqdm(range(t.size)):
        un = u.copy()
        vn = v.copy()
        thetan = theta.copy()
        
        p = pressure_poisson(un, vn, dt, dx, dy, p)
        u[1:-1, 1:-1] = cavidade_u(un, vn, dt, dx, dy, p, re)
        v[1:-1, 1:-1] = cavidade_v(un, vn, thetan, dt, dx, dy, p, re, ri)
        theta[1:-1, 1:-1] = cavidade_theta(un, vn, thetan, dt, dx, dy, p, re, pr)
        
        # Condições de contorno
        u[:, 0] = 0
        u[0, :] = 0
        u[-1, :] = 0
        u[:, -1] = 1  # Definir velocidade na tampa da cavidade igual a 1
        v[:, 0] = 0
        v[:, -1] = 0
        v[0, :] = 0
        v[-1, :] = 0
        theta[-1, :] = 1
        theta[:, 0] = 1
        theta[0, :] = 1
        theta[:, -1] = 0  # Definir theta na tampa da cavidade igual a 1
    
    return u, v, theta, p
```

Hora da ação:


```python
# Coordenadas
x = np.linspace(start = 0.0, stop = 2.0, num=21)
y = np.linspace(start = 0.0, stop = 2.0, num=21)
t = np.arange(start = 0.0, stop = 5.0, step = 0.001)

u, v, theta, p = cavidade(x, y, t, re = 40.0, ri = 0.0, pr = 1.0)
```

Verificando os resultados:


```python
fig = px.imshow(
    theta.T,
    x=x,
    y=y,
    color_continuous_scale="RdBu_r",
    title="Temperatura",
    labels=dict(x="x", y="y"),
    origin="lower",
)
fig.show()
```


```python
fig = ff.create_streamline(
    x[1:-1], y[1:-1], u.T[1:-1, 1:-1], v.T[1:-1, 1:-1], marker_color="black",
)

fig.add_trace(
    go.Contour(
        z=theta.T,
        x=x,
        y=y,
        colorscale="RdBu_r",
        colorbar=dict(title="Temperatura", titleside="right"),
    )
)

fig.update_layout(
    title="Escoamento em Cavidade com Transferência de Calor",
    xaxis_title="x",
    yaxis_title="y",
    autosize=False,
    width=800,
    height=800,
)

fig.show()
```

> Leitura recomendada:
> * Gostaria de se apronfundar mais no assunto? Veja a aula [Métodos Numéricos com Python](https://github.com/fschuch/metodos-numericos-com-python), que inclui a solução do problema da cavidade com transferência de calor utilizando conceitos Python mais avançados, como a criação e manipulação de classes em programação orientada ao objeto (OOP).

-----

> **Felipe N. Schuch**,<br>
> Pesquisador em Fluidodinâmica Computacional na PUCRS, com interesse em: Escoamentos turbulentos, transferência de calor e massa, e interação fluido-estrutura; Processamento e visualização de dados em Python; Jupyter Notebook como uma ferramenta de colaboração, pesquisa e ensino.<br>
> [felipeschuch@outlook.com](mailto:felipeschuch@outlook.com "Email") [@fschuch](https://twitter.com/fschuch "Twitter") [Aprenda.py](https://fschuch.github.io/aprenda.py "Blog") [@aprenda.py](https://www.instagram.com/aprenda.py/ "Instagram")<br>

-----

4d12b0624fe33a21f9313e12521a05e18056cd63
ipynb
Jupyter Notebook
Aulas/03-Exemplos-Conceccao.ipynb
fschuch/Python-Transferencia-de-Calor
4504639af92a19940f99b5b171739bc960cc9400
MIT
2022-01-14T22:22:34.000Z
2022-01-14T22:22:34.000Z
Aulas/03-Exemplos-Conceccao.ipynb
fschuch/Python-Transferencia-de-Calor
4504639af92a19940f99b5b171739bc960cc9400
MIT
Aulas/03-Exemplos-Conceccao.ipynb
fschuch/Python-Transferencia-de-Calor
4504639af92a19940f99b5b171739bc960cc9400
MIT
Qwen/Qwen-72B
1. YES
2. YES