TESTEDOCX docx

Projeto Computacional 1 - Dimensionamento de um eixo submetido á torção puraConsidere o eixo circular maciço apresentado na Figura 1.Figura 1: Eixo circular de diâmetro constante.Item AddxG Jxdϕxdx =-t0 , x∈(0,L)s.aϕ0=0TL=TResolva analiticamente a equação diferencial de equilíbrio para o eixo considerando que o momento polar de inércia J e o módulo de elasticidade transversal G são constantes. Solução da equação diferencial de equilíbrio: Como G e J são constantes temos a equação: ddxdϕxdx =-t0G JIntegrando na equação, temos: ddϕxdx =-t0G Jdx ⇒ dϕxdx=-t0G Jx+CIntegrando novamente, temos: dϕx= -t0G Jx+C dx Portanto, a solução geral da equação diferencial de equilíbrio é: ϕx= -t02G Jx2+Cx+K Como Tx=GJdϕxdx=GJ-t0G Jx+C=-t0x+GJCUsando a condição TL=T, então: T =-t0L+GJC ⇒ C=1GJT+t0LUsando a condição ϕ0=0, chegamos a: K=0. Portanto,ϕx= -t02G Jx2+1GJT+t0LxItem B Utilizando um esquema de diferenças centradas, implemente computacionalmente a solução da equação diferencial de equilíbrio, determinando o diâmetro mínimo do eixo para suportar uma tensão máxima τadm=65 MPa e ao mesmo tempo um ângulo de torção máximo ϕadm=1.5°. Utilize um método iterativo, considerando o diâmetro inicial do eixo d0=25 mm com incrementos de 1 mm, e um número máximo de iterações de 50. Para solução utilize: • G = 77, 2 GPa • t0 = 100 Nm/m; • T= 200 Nm; • L = 1, 5 mImplementação computacionalProgramemos o código em MATLAB que resolva a equação diferencial de equilíbrio para determinar o diâmetro mínimo de um eixo que deve suportar uma tensão máxima admissível e um ângulo de torção máximo admissível. Comecemos definindo os dados do problema:% Dados do problematau_adm = 65e6; % Tensão máxima admissível em Paphi_adm = deg2rad(1.5); % Ângulo de torção máximo admissível em radianosL = 1.5; % comprimento do eixo em metrosT = 200; % Condição de contorno no extremo direito em Nmd_0 = 25e-3; % Diâmetro inicial do eixo em metrosincremento = 1e-3; % Incremento no diâmetro do eixo em metrosnum_iter = 50; % Número máximo de iteraçõesG = 77.2e9; % Módulo de cizalhamento do material (Pa)Agora vamos inicializar vetores e matrizes:% Inicialização de matrizes e vetoresx = linspace(0, L, N)'; % Vetor de posições normalizado (de 0 a L)phi = zeros(N, 1); % Vetor de ângulos de torção inicializados a zerophi(1) = 0; % Condição de contorno no extremo esquerdo % Inicialização de matriz bidimensional para armazenar todos os valores de phi em cada iteraçãophi_values = zeros(N, num_iter);Agora vamos com a implementação que nos proporcionara o diâmetro mínimo:O código começa com um loop for para iterar sobre diferentes diâmetros do eixo.for j = 1:num_iterPara cada iteração calculemos o momento de inércia polar J para cada diâmetro. Usamos j-1, já que temos que começar com o diâmetro d0=65 mm. d = d_0 + (j-1)*incremento;Também para cada iteração calculemos o valor do momento de inércia polar J. Estamos considerando um eixo circular maciço então o momento de inércia polar pode-se calcular como J=πd432.% Momento de inércia polar do eixo na iteração atualJ = pi*(d^4)/32;Por outro lado, para programar usando o método das diferenças centradasϕx+h=2ϕx-ϕx-h-h2D2ϕxPrecisamos de duas condições iniciais, ϕ0=0 e ϕ(1). Para determinar ϕ(1) vamos usar o fato dϕxdx=-t0G Jx+1GJT+t0L ⇒ dϕ0dx=1GJT+t0LE o método das diferenças para a primeira derivada:ϕx+h=ϕx+hDϕxAssim, no MATLAB temos o seguinte código:% Calcular dphi_dx dphi_dx = (T+150) / (G * J); % Primeira derivada em x = 0 % Calcular phi(2) usando a aproximação de diferenças finitas para frente phi(2) = phi(1) + dphi_dx * h; % Aproximação de phi(2)Usando for i = 2:N-1,determinamos os ângulo de torção ϕ para cada ponto usando o método das diferenças centradas.phi(i+1) = 2*phi(i) - phi(i-1) - h^2 * (10^(-7) / (77.2 * J));Também salvamos os valores de ϕ de cada iteração na matriz phi_values.% Armazenar os valores de phi da iteração atual na matriz bidimensional phi_values(:, j) = phi;Antes de verificar se as condições de tensão e ângulo de torção máxima se satisfazem, se calcula o máximo valor de ϕ em cada iteração. Se as condições forem atendidas, imprime o diâmetro mínimo do eixo.% Calcular o valor máximo de phi na iteração atual max_phi = max(abs(phi));% Verificar se as condições de término são atendidas if max_phi < phi_adm && T*d/(2*J) < tau_adm fprintf('O diâmetro mínimo do eixo é: %.2f mm\n', d*1e3); break; % Sair do loop se as condições de término forem atendidas endFinalmente, se o loop terminar sem encontrar uma solução satisfatória, imprime uma mensagem informando que o número máximo de iterações foi alcançado sem encontrar uma solução.% Se o loop terminar sem encontrar uma solução satisfatóriaif j == num_iter fprintf(' Não se achou solução.\n');endExecução do programa: Para N=100, o diâmetro mínimo do eixo é: 38 mm. Utilizando um esquema de diferença centrada com um ponto auxiliar, para aproximar o ponto T(0), apresente os diagramas de torque, distribuição de tensão cisalhante e ângulo de torção ao longo do eixo, comparando a solução exata e aproximada. Avalie e apresente o erro relativo para malhas de 25, 50 e 100 pontos igualmente esparsados. Para solução utilize: • G = 77, 2 GPa; • t0 = 100 Nm/m; • T = 200 Nm; • L = 1, 5 m.Código em MATLAB Considera-se d=0.038m para a solução aproximada. O código solicita ingressar um ponto auxiliar, para fazer o calculo das derivadas. Sugestão para o ponto: 201.5. Os gráficos correspondem a malha de 25. Para visualizar os gráficos para as malhas de 50, 100 ha que mudar o vetor x = linspace(0, L, num_points(1)), por 2 e 3 respetivamente e também no plot. % Dados do problemaG = 77.2e9; % Módulo de cisalhamento em Pat0 = 100; % Torque por unidade de comprimento em Nm/mT = 200; % Torque total em NmL = 1.5; % Comprimento do eixo em metrosd = 0.038; J=pi*d^4/32;pontoaux = input('Por favor, digite um ponto auxiliar: '); % Número de pontos para cada malhanum_points = [25, 50, 100]; % Inicializar vetores para armazenar os resultadosder = cell(1, numel(num_points));s = cell(1, numel(num_points));torque_exact = cell(1, numel(num_points));torque_approx = cell(1, numel(num_points));shear_stress_exact = cell(1, numel(num_points));shear_stress_approx = cell(1, numel(num_points));angle_exact = cell(1, numel(num_points));angle_approx = cell(1, numel(num_points));error_torque = cell(1, numel(num_points));error_shear_stress = cell(1, numel(num_points));error_angle = cell(1, numel(num_points));x = cell(1, numel(num_points)); for i = 1:numel(num_points) % Passo de discretização h = L / (num_points(i) - 1); % Vetor de posições normalizadas x{i} = linspace(0, L, num_points(i)); s{i}(1) = T/(J*G); s{i}(2) = pontoaux/(J*G); for j = 3: num_points(i) s{i} (j) = s{i}(j-2) + (2*h*t0)/(J*G); end for k = 1:num_points(i) der{i}(k) = s{i}( num_points(i)-k+1); end torque_approx{i}(num_points(i)) =T; shear_stress_approx{i}(num_points(i)) =T*d/(2*J); for k = 1:num_points(i)-1 torque_approx{i}(k) = J*G*der{i}(k); shear_stress_approx{i}(k) = 0.5*d*G*der{i}(k); end angle_approx{i}(1)=0; % Calcular phi(2) usando a aproximação de diferenças finitas angle_approx{i}(2) = angle_approx{i}(1) + der{i}(1) * h; % Aproximação de phi(2) for k = 2:num_points(i)-1 % Atualizar o ângulo de torção usando diferenças centradas angle_approx{i}(k+1) = 2*angle_approx{i}(k) - angle_approx{i}(k-1) - h^2 * (10^(-7) / (77.2 * J)); end % Soluçoes exatas angle_exact{i} = (-t0 * x{i}.^2 / (2 * G*J)) + ((T + t0 * L) * x{i} / (J*G)); torque_exact{i} = T + t0*(L- x{i}); shear_stress_exact{i} = torque_exact{i} * 0.5 * d/J; %Calcular erros relativos error_torque{i} = abs(torque_exact{i} - torque_approx{i}) ./ torque_exact{i}; error_shear_stress{i} = abs(shear_stress_exact{i} - shear_stress_approx{i}) ./ shear_stress_exact{i}; error_angle{i} = abs(angle_exact{i} - angle_approx{i}) ./ angle_exact{i}; end disp(sprintf(' 25 Pontos: T(0) = %f\n 50 Pontos: T(0) = %f\n 100 Pontos: T(0) = %f', torque_approx{1}(1), torque_approx{2}(1), torque_approx{3}(1))); % Gráficos e análise de errofigure;x = linspace(0, L, num_points(1));subplot(2, 2, 1);plot(x, torque_exact{1}, 'r', 'LineWidth', 2);hold on;plot(x, torque_approx{1}, 'b--', 'LineWidth', 2);xlabel('Posição (m)');ylabel('Torque (Nm)');legend('Solução Exata', 'Solução Aproximada');title('Comparação do Torque'); subplot(2, 2, 2);plot(x, shear_stress_exact{1}, 'r', 'LineWidth', 2);hold on;plot(x, shear_stress_approx{1}, 'b--', 'LineWidth', 2);xlabel('Posição (m)');ylabel('Tensão Cisalhante (Pa)');legend('Solução Exata', 'Solução Aproximada');title('Comparação da Tensão Cisalhante'); subplot(2, 2, 3);plot(x, angle_exact{1}, 'r', 'LineWidth', 2);hold on;plot(x, angle_exact{1}, 'b--', 'LineWidth', 2);xlabel('Posição (m)');ylabel('Ângulo de Torção (rad)');legend('Solução Exata', 'Solução Aproximada');title('Comparação do Ângulo de Torção'); subplot(2, 2, 4);plot(x, error_angle{1}, 'k', 'LineWidth', 2);xlabel('Posição (m)');ylabel('Erro Relativo (%)');title('Erro Relativo da Solução Aproximada');Execução do programa: Consideramos o ponto auxiliar 201.5 e obtemos o seguinte dados:T(0) para as diferentes malhas:25 Pontos: T(0) = 350.000000 50 Pontos: T(0) = 348.438776 100 Pontos: T(0) = 349.984848Gráficos para malha de 25 pontos: Gráficos para malha de 50 pontos: Gráficos para malha de 50 pontos: