TESTEQUARTAINDEXAÇÃO docx

EM503 - M´etodos Num´ericos aplicados `a Engenharia1◦ semestre de 2024Projeto Computacional 2 - Solu¸c˜ao da equa¸c˜ao de Burgers em duas dimenso˜esEntrega: 23 de maio de 2024.Considere um curso d’água bidimensional sujeito a um campo de velocidades ux,y,v(x,y) dado por: ux,y=αvx,y=α12+sen(π5x)Despeja-se uma substancia química em uma posição (a, b) do domínio bidimensional de análise 0,Lx x [0,Ly] durante um intervalo de tempo T , a uma taxa Q . Utiliza-se a equação de Burgers para descrever o comportamento da concentração Cx, y, t do contaminante. Sendo essa dada por:∂∂tCx,y,t+ux,y∂∂xCx,y,t+vx,y∂∂yCx,y,t-kx,y∂2∂x2Cx,y,t+∂2∂y2Cx,y,t=qcSendokx,y o coeficiente de difusão da sustância com o médio definido pela equaçãokx,y=β|ux,y+v(x,y)|qc o termo de gera¸c˜ao (Q para t ∈0,T )Para este problema pede-se:Considerando T=5s, Q =50kgm3s, α=1.5ms, (a, b) = (6,3), β=0,02m, obtenha a solu¸c˜ao Cx,y,tUtilizando um esquema bidimensional de diferenças finitas, com diferença centrada na discretiza¸c˜ao do domínio e malha temporal segundo o método de Euler.Avalie quantos segundos são necessários para que a concentração da substância química despejada seja considerada desprezível, para tal, considere maxCx, y, t<10-4.Analise a influencia do coeficiente de difusão kx, y na dispersão de contaminante. Para isso, utilize na solução β = 0,2m, β = 0,4m e β = 0,6m. Comente sobre seu impacto na dispersão da substância ao longo do domínio, observando também sua influencia no tempo de purificação do curso d’agua estudado.Observa¸c˜oesPara discretiza¸c˜ao espacial, utilize a malha apresentada na Figura 1, com elementos quadriculares de h=0.5m.Discretize o tempo de forma que γ=∆th2≤0,3 (condição de convergência CFL).Como condi¸c˜ao inicial, utilizeCx,y,0=0 , ∀ 0,Lx x [0, Ly] .Para condi¸c˜oes de contorno, adote a derivada espacial normal ao contorno nula, sendo as condi¸c˜oes descritas como:Cx0,yj,tk= Cx1,yj,tkCxn,yj,tk= Cxn-1,yj,tkCxi,y0,tk= Cxi,y1,tkCxi,ym,tk= Cxi,ym-1,tk.Cheque a cada itera¸c˜ao se Cxi,yj,tk<0, caso ocorra, fa¸ca Cxi,yj,tk=0.y20 mx30 mSoluções: a) Da diferença centrada para as derivadas parciais, temos:∂∂xCx,y,t=Cx+∆x,y,t-Cx-∆x,y,t2∆x∂∂yCx,y,t=Cx,y+∆y,t-Cx,y-∆y,t2∆y∂2∂x2Cx,y,t=Cx+∆x,y,t-2C(x,y,t)+Cx-∆x,y,t(∆x)2∂2∂y2Cx,y,t=Cx,y+∆y,t-2C(x,y,t)+Cx,y-∆y,t(∆y)2 Da equação de Burgers, temos:∂∂tCx,y,t=-ux,y∂∂xCx,y,t-vx,y∂∂yCx,y,t+kx,y∂2∂x2Cx,y,t+∂2∂y2Cx,y,t+qcSubstituindo as equações de diferenças centradas:∂∂tCx,y,t=-ux,yCx+∆x,y,t-Cx-∆x,y,t2∆x-vx,yCx,y+∆y,t-Cx,y-∆y,t2∆y+kx,yCx+∆x,y,t-2Cx,y,t+Cx-∆x,y,t∆x2+Cx,y+∆y,t-2Cx,y,t+Cx,y-∆y,t∆y2+qcE do método de Euler:Cx,y,t+∆t=Cx,y,t+∆t-ux,yCx+∆x,y,t-Cx-∆x,y,t2∆x-vx,yCx,y+∆y,t-Cx,y-∆y,t2∆y+kx,yCx+∆x,y,t-2Cx,y,t+Cx-∆x,y,t∆x2+Cx,y+∆y,t-2Cx,y,t+Cx,y-∆y,t∆y2+qcAgora vamos implementar em MATLAB:Para a discretização espacial, utilizamos dx=dy=0.5. Denotamos por Lx,Ly o número de pontos em x e y respetivamente. Para a discretização do tempo escolhemos dy tal que γ=dydx2≤0,3 (condição de convergência CFL). Logo, escolhemos dt=0.025 e daí, Nt=201(número de passos do tempo).Para a condição inicial, Inicializamos a matriz C como uma matriz nula. E para as condições de contorno, creamos a função function C = aplicarCondicoesDeContorno(C)% Função para aplicar as condições de contornofunction C = aplicarCondicoesDeContorno(C) C(1, :) = C(2, :); % Condição em x_0 C(end, :) = C(end-1, :); % Condição em x_n C(:, 1) = C(:, 2); % Condição em y_0 C(:, end) = C(:, end-1); % Condição em y_mend% Função para aplicar as condições de contornofunction C = aplicarCondicoesDeContorno(C) C(1, :) = C(2, :); % Condição em x_0 C(end, :) = C(end-1, :); % Condição em x_n C(:, 1) = C(:, 2); % Condição em y_0 C(:, end) = C(:, end-1); % Condição em y_mend O código em Matlab:% Inserir o valor de TTmax = input('Insira o valor de T (tempo máximo de simulação em segundos): '); % Parâmetros do problemaLx = 30; % comprimento do domínio em xLy = 20; % comprimento do domínio em yNx = 61; % número de pontos na malha em xNy = 41; % número de pontos na malha em ydt = 0.025; % passo de tempoNt = 1+Tmax/dt; % número de passos de tempo dx = Lx / (Nx - 1);dy = Ly / (Ny - 1);x = linspace(0, Lx, Nx);y = linspace(0, Ly, Ny);[X, Y] = meshgrid(x, y); % Parâmetros físicosalpha = 1.5; % m/sQ_dot = 50; % kg/(m^3 s)a = 6; % posição x da injeçãob = 3; % posição y da injeçãobeta = 0.02; % coeficiente de difusão % Índices da posição (a, b)[~, idx_a] = min(abs(x - a));[~, idx_b] = min(abs(y - b)); % Velocidadesu = alpha * ones(Ny, Nx); v = alpha * (0.5 + sin(pi / 5 * X)); % Coeficiente de difusãok = beta * abs(u + v); % Condição inicialC = zeros(Ny, Nx, Nt); % Variáveis de controle do tempot = 0;n = 1; % índice de tempoconcentracao_maxima = max(C(:));% Simulação temporal com loop whilewhile t <= Tmax % Inicializar nova concentração C_new = C(:, :, n); % Loop sobre a malha interna for i = 2:Ny-1 for j = 2:Nx-1% Inserir o valor de TTmax = input('Insira o valor de T (tempo máximo de simulação em segundos): '); % Parâmetros do problemaLx = 30; % comprimento do domínio em xLy = 20; % comprimento do domínio em yNx = 61; % número de pontos na malha em xNy = 41; % número de pontos na malha em ydt = 0.025; % passo de tempoNt = 1+Tmax/dt; % número de passos de tempo dx = Lx / (Nx - 1);dy = Ly / (Ny - 1);x = linspace(0, Lx, Nx);y = linspace(0, Ly, Ny);[X, Y] = meshgrid(x, y); % Parâmetros físicosalpha = 1.5; % m/sQ_dot = 50; % kg/(m^3 s)a = 6; % posição x da injeçãob = 3; % posição y da injeçãobeta = 0.02; % coeficiente de difusão % Índices da posição (a, b)[~, idx_a] = min(abs(x - a));[~, idx_b] = min(abs(y - b)); % Velocidadesu = alpha * ones(Ny, Nx); v = alpha * (0.5 + sin(pi / 5 * X)); % Coeficiente de difusãok = beta * abs(u + v); % Condição inicialC = zeros(Ny, Nx, Nt); % Variáveis de controle do tempot = 0;n = 1; % índice de tempoconcentracao_maxima = max(C(:));% Simulação temporal com loop whilewhile t <= Tmax % Inicializar nova concentração C_new = C(:, :, n); % Loop sobre a malha interna for i = 2:Ny-1 for j = 2:Nx-1 % Derivadas espaciais dCdx = (C(i+1, j, n) - C(i-1, j, n)) / (2 * dx); dCdy = (C(i, j+1, n) - C(i, j-1, n)) / (2 * dy); d2Cdx2 = (C(i+1, j, n) - 2 * C(i, j, n) + C(i-1, j, n)) / (dx^2); d2Cdy2 = (C(i, j+1, n) - 2 * C(i, j, n) + C(i, j-1, n)) / (dy^2); % Atualização da concentração C_new(i, j) = C(i, j, n) + dt * ( ... - u(i, j) * dCdx ... - v(i, j) * dCdy ... + k(i, j) * (d2Cdx2 + d2Cdy2) ... + Q_dot * (i == idx_a && j == idx_b) ... ); if C_new(i, j)<0 C_new(i, j)=0; end end end % Aplicar condições de contorno C_new = aplicarCondicoesDeContorno(C_new); % Atualizar a concentração para o próximo passo de tempo n = n + 1; t = t + dt; if n <= Nt C(:, :, n) = C_new; else break; end % Atualizar a concentração máxima concentracao_maxima = max(C_new(:)); if concentracao_maxima < concentracao_desprezivel break; endendfprintf(' Cmax %.6f .\n', max(C_new(:))); % Plotar a concentração final % Visualização no espaco[X, Y] = meshgrid(x, y);surf(X, Y, C_new);shading interp;colorbar;xlabel('x (m)');ylabel('y (m)');zlabel('Concentração C(x,y,T)');title('Distribuição de Concentração após T segundos'); % Visualização no espaco no planofigure;contourf(X, Y, C_new);colorbar;title('Concentração de Contaminante');xlabel('x');ylabel('y'); % Derivadas espaciais dCdx = (C(i+1, j, n) - C(i-1, j, n)) / (2 * dx); dCdy = (C(i, j+1, n) - C(i, j-1, n)) / (2 * dy); d2Cdx2 = (C(i+1, j, n) - 2 * C(i, j, n) + C(i-1, j, n)) / (dx^2); d2Cdy2 = (C(i, j+1, n) - 2 * C(i, j, n) + C(i, j-1, n)) / (dy^2); % Atualização da concentração C_new(i, j) = C(i, j, n) + dt * ( ... - u(i, j) * dCdx ... - v(i, j) * dCdy ... + k(i, j) * (d2Cdx2 + d2Cdy2) ... + Q_dot * (i == idx_a && j == idx_b) ... ); if C_new(i, j)<0 C_new(i, j)=0; end end end % Aplicar condições de contorno C_new = aplicarCondicoesDeContorno(C_new); % Atualizar a concentração para o próximo passo de tempo n = n + 1; t = t + dt; if n <= Nt C(:, :, n) = C_new; else break; end % Atualizar a concentração máxima concentracao_maxima = max(C_new(:)); if concentracao_maxima < concentracao_desprezivel break; endendfprintf(' Cmax %.6f .\n', max(C_new(:))); % Plotar a concentração final % Visualização no espaco[X, Y] = meshgrid(x, y);surf(X, Y, C_new);shading interp;colorbar;xlabel('x (m)');ylabel('y (m)');zlabel('Concentração C(x,y,T)');title('Distribuição de Concentração após T segundos'); % Visualização no espaco no planofigure;contourf(X, Y, C_new);colorbar;title('Concentração de Contaminante');xlabel('x');ylabel('y');Execução do programa:% Função para aplicar as condições de contornofunction C = aplicarCondicoesDeContorno(C) C(1, :) = C(2, :); % Condição em x_0 C(end, :) = C(end-1, :); % Condição em x_n C(:, 1) = C(:, 2); % Condição em y_0 C(:, end) = C(:, end-1); % Condição em y_mend% Função para aplicar as condições de contornofunction C = aplicarCondicoesDeContorno(C) C(1, :) = C(2, :); % Condição em x_0 C(end, :) = C(end-1, :); % Condição em x_n C(:, 1) = C(:, 2); % Condição em y_0 C(:, end) = C(:, end-1); % Condição em y_mendAo executar o codigo temos:Da figura 1:Depois de T=5 s, a concentração C(6,3,5) é aproximadamente 28.6979.(que também é a concentração máxima)Da figura 1:Depois de T=5 s, a concentração C(6,3,5) é aproximadamente 28.6979.(que também é a concentração máxima) Da figura 2:Depois de T=5s, o contaminante se dispersa, perto do ponto (6,3) e no resto dos pontos a concentração é 0.Da figura 2:Depois de T=5s, o contaminante se dispersa, perto do ponto (6,3) e no resto dos pontos a concentração é 0. b) Para determinar quantos segundos são necessários para que a concentração da substância química seja considerada desprezível (maxCx, y, t<10-4) precisamos modificar o código do item a) para simular a equação de Burgers até que a concentração máxima no domínio atinja esse valor. Vamos adicionar também uma verificação durante a simulação para parar quando essa condição for satisfeita e registrar o tempo necessário.Vamos a pedir que o usuário insira o Tmax, tempo máximo, e com ajuda de um loop while vamos a procurar o tempo necessário para que a concentração considerada desprezívelPara a verificação da Concentração Máxima, dentro do loop temporal (while), após a atualização da concentração, a condição concentracao_maxima < concentracao_desprezivel verifica se a concentração máxima no domínio é menor que 10-4. Se a condição for satisfeita, tempo é definido como o tempo atual t e o loop é interrompido (break).Após o loop, o código vai verificar se tempo é maior que 0 (indicando que a condição foi satisfeita) e exibirá o tempo necessário. Se a concentração não atingiu o valor desprezível durante o tempo de simulação, uma mensagem será exibida.O código em MATLAB:% Pedir ao usuário para inserir o valor de TmaxTmax = input('Insira o valor de Tmax (tempo máximo de simulação em segundos): '); % Parâmetros do problemaLx = 30; % comprimento do domínio em xLy = 20; % comprimento do domínio em yNx = 61; % número de pontos na malha em xNy = 41; % número de pontos na malha em ydt = 0.025; % passo de tempoNt = 1+Tmax/dt; % número de passos de tempo dx = Lx / (Nx - 1);dy = Ly / (Ny - 1);x = linspace(0, Lx, Nx);y = linspace(0, Ly, Ny);[X, Y] = meshgrid(x, y); % Parâmetros físicosalpha = 1.5; % m/sQ_dot = 50; % kg/(m^3 s)a = 6; % posição x da injeçãob = 3; % posição y da injeçãobeta = 0.02; % coeficiente de difusão% Índices da posição (a, b)[~, idx_a] = min(abs(x - a));[~, idx_b] = min(abs(y - b)); % Velocidadesu = alpha * ones(Ny, Nx); v = alpha * (0.5 + sin(pi / 5 * X)); % Pedir ao usuário para inserir o valor de TmaxTmax = input('Insira o valor de Tmax (tempo máximo de simulação em segundos): '); % Parâmetros do problemaLx = 30; % comprimento do domínio em xLy = 20; % comprimento do domínio em yNx = 61; % número de pontos na malha em xNy = 41; % número de pontos na malha em ydt = 0.025; % passo de tempoNt = 1+Tmax/dt; % número de passos de tempo dx = Lx / (Nx - 1);dy = Ly / (Ny - 1);x = linspace(0, Lx, Nx);y = linspace(0, Ly, Ny);[X, Y] = meshgrid(x, y); % Parâmetros físicosalpha = 1.5; % m/sQ_dot = 50; % kg/(m^3 s)a = 6; % posição x da injeçãob = 3; % posição y da injeçãobeta = 0.02; % coeficiente de difusão% Índices da posição (a, b)[~, idx_a] = min(abs(x - a));[~, idx_b] = min(abs(y - b)); % Velocidadesu = alpha * ones(Ny, Nx); v = alpha * (0.5 + sin(pi / 5 * X)); % Coeficiente de difusãok = beta * abs(u + v); % Condição inicialC = zeros(Ny, Nx, Nt); % Variáveis de controle do tempot = 0;n = 1; % índice de tempoconcentracao_maxima = max(C(:));concentracao_desprezivel = 1e-4; % Simulação temporal com loop whilewhile t <= Tmax % Inicializar nova concentração C_new = C(:, :, n); % Loop sobre a malha interna for i = 2:Ny-1 for j = 2:Nx-1 % Derivadas espaciais dCdx = (C(i+1, j, n) - C(i-1, j, n)) / (2 * dx); dCdy = (C(i, j+1, n) - C(i, j-1, n)) / (2 * dy); d2Cdx2 = (C(i+1, j, n) - 2 * C(i, j, n) + C(i-1, j, n)) / (dx^2); d2Cdy2 = (C(i, j+1, n) - 2 * C(i, j, n) + C(i, j-1, n)) / (dy^2); % Atualização da concentração if t <= 5 C_new(i, j) = C(i, j, n) + dt * ( ... - u(i, j) * dCdx ... - v(i, j) * dCdy ... + k(i, j) * (d2Cdx2 + d2Cdy2) ... + Q_dot * (i == idx_a && j == idx_b) ... ); if C_new(i, j)<0 C_new(i, j)=0; end else C_new(i, j) = C(i, j, n) + dt * ( ... - u(i, j) * dCdx ... - v(i, j) * dCdy ... + k(i, j) * (d2Cdx2 + d2Cdy2)); if C_new(i, j)<0 C_new(i, j)=0; end end end end % Aplicar condições de contorno C_new = aplicarCondicoesDeContorno(C_new); % Coeficiente de difusãok = beta * abs(u + v); % Condição inicialC = zeros(Ny, Nx, Nt); % Variáveis de controle do tempot = 0;n = 1; % índice de tempoconcentracao_maxima = max(C(:));concentracao_desprezivel = 1e-4; % Simulação temporal com loop whilewhile t <= Tmax % Inicializar nova concentração C_new = C(:, :, n); % Loop sobre a malha interna for i = 2:Ny-1 for j = 2:Nx-1 % Derivadas espaciais dCdx = (C(i+1, j, n) - C(i-1, j, n)) / (2 * dx); dCdy = (C(i, j+1, n) - C(i, j-1, n)) / (2 * dy); d2Cdx2 = (C(i+1, j, n) - 2 * C(i, j, n) + C(i-1, j, n)) / (dx^2); d2Cdy2 = (C(i, j+1, n) - 2 * C(i, j, n) + C(i, j-1, n)) / (dy^2); % Atualização da concentração if t <= 5 C_new(i, j) = C(i, j, n) + dt * ( ... - u(i, j) * dCdx ... - v(i, j) * dCdy ... + k(i, j) * (d2Cdx2 + d2Cdy2) ... + Q_dot * (i == idx_a && j == idx_b) ... ); if C_new(i, j)<0 C_new(i, j)=0; end else C_new(i, j) = C(i, j, n) + dt * ( ... - u(i, j) * dCdx ... - v(i, j) * dCdy ... + k(i, j) * (d2Cdx2 + d2Cdy2)); if C_new(i, j)<0 C_new(i, j)=0; end end end end % Aplicar condições de contorno C_new = aplicarCondicoesDeContorno(C_new); % Atualizar a concentração para o próximo passo de tempo n = n + 1; t = t + dt;if n <= Nt C(:, :, n) = C_new; else break; end % Atualizar a concentração máxima concentracao_maxima = max(C_new(:)); if concentracao_maxima < concentracao_desprezivel break; endendfprintf(' Cmax %.6f .\n', max(C_new(:))); % Exibir o tempo necessário para a concentração ser desprezívelif concentracao_maxima < concentracao_desprezivel && t>0 fprintf('O tempo necessário para que a concentração seja desprezível é de %.3f segundos.\n', t);else fprintf('A concentração não atingiu o valor desprezível durante o tempo de simulação.\n');end % Plotar a concentração final se não for constante e t > 0if t > 0 && min(C_new(:)) ~= max(C_new(:)) figure; contourf(X, Y, C_new); colorbar; title('Concentração de Contaminante'); xlabel('x (m)'); ylabel('y (m)');else fprintf('A concentração é constante, não há contornos para plotar.\n');end % Função para aplicar as condições de contornofunction C = aplicarCondicoesDeContorno(C) C(1, :) = C(2, :); % Condição em x_0 C(end, :) = C(end-1, :); % Condição em x_n C(:, 1) = C(:, 2); % Condição em y_0 C(:, end) = C(:, end-1); % Condição em y_mEnd % Atualizar a concentração para o próximo passo de tempo n = n + 1; t = t + dt;if n <= Nt C(:, :, n) = C_new; else break; end % Atualizar a concentração máxima concentracao_maxima = max(C_new(:)); if concentracao_maxima < concentracao_desprezivel break; endendfprintf(' Cmax %.6f .\n', max(C_new(:))); % Exibir o tempo necessário para a concentração ser desprezívelif concentracao_maxima < concentracao_desprezivel && t>0 fprintf('O tempo necessário para que a concentração seja desprezível é de %.3f segundos.\n', t);else fprintf('A concentração não atingiu o valor desprezível durante o tempo de simulação.\n');end % Plotar a concentração final se não for constante e t > 0if t > 0 && min(C_new(:)) ~= max(C_new(:)) figure; contourf(X, Y, C_new); colorbar; title('Concentração de Contaminante'); xlabel('x (m)'); ylabel('y (m)');else fprintf('A concentração é constante, não há contornos para plotar.\n');end % Função para aplicar as condições de contornofunction C = aplicarCondicoesDeContorno(C) C(1, :) = C(2, :); % Condição em x_0 C(end, :) = C(end-1, :); % Condição em x_n C(:, 1) = C(:, 2); % Condição em y_0 C(:, end) = C(:, end-1); % Condição em y_mEndExecução do programa:Para executar o programa, ha que considerar um tempo suficientemente grande para que a concentração da substância química despejada seja considerada desprezível, (isto é, seja muito pequena), ou seja, para satisfazer maxCx, y, t<10-4. Notemos também do item a) o tempo tem que ser maior que Tmax=5. Então, vamos a tentar executar com Tmax = 30. Ao executar o codigo temos: Depois de 21.575 segundos, a concentração máxima é 0.000098 inferior a 10-4, logo a substância química despejada pode ser considerada desprezível. Do gráfico, podemos notar a concentração final da sustância química. c) Analise a influencia do coeficiente de difusão kx, y na dispersão de contaminante. Para isso, utilize na solução β = 0,2m, β = 0,4m e β = 0,6m. Comente sobre seu impacto na dispersão da substância ao longo do domínio, observando também sua influencia no tempo de purificação do curso d’agua estudado.Para analisar a influência do coeficiente de difusão, vamos usar o programa do item b) como uma pequena modificação. Vamos pedir a usuário ingressar o valor da constante β. Para β = 0,2m:Para β = 0,4m:Como a concentração não atingiu o valor desprezível, consideramos: T=22, T=25, T=30, T=35, T=80, T=250 Para β=0.6: Como a concentração não atingiu o valor desprezível, consideramos: T=22, T=25, T=30, T=35, T=80, T=250 Para β = 0,2m:Obtivemos que depois de 18.175 a concentração pode se considerar desprezível, ou seja, tempo de purificação é menor ao caso analisado no item a. Por outro lado, ao passar o tempo, a sustância química vai se dispersando ao longo do domínio e com uma velocidade menor aos casos β = 0,4m e β = 0,6m.Para β = 0,4m:Não obtivemos um tempo para considerar que a concentração seja considerada desprezível. Por outro lado notamos que, ao passar o tempo, a sustância química vai se dispersando ao longo do domínio. Para β = 0,8m:Não obtivemos um tempo para considerar que a concentração seja considerada desprezível. Por outro lado notamos que, ao passar o tempo, a sustância química vai se dispersando ao longo do domínio e com uma velocidade maior ao caso β = 0,2m e β = 0,4m.14