MODELAGEM SIGMOIDAL COMPLETA VERSUS ABORDAGEM LINEAR NA MECÂNICA DA FRATURA: ANÁLISE DE SENSIBILIDADE COMPUTACIONAL DOS MODELOS DE PARIS E NASGRO

REGISTRO DOI: 10.70773/revistatopicos/789070601

RESUMO
A previsão da vida útil remanescente de estruturas e dispositivos submetidos a carregamentos cíclicos é um pilar crítico para o gerenciamento de integridade baseado na filosofia da tolerância ao dano. Este artigo apresenta um estudo comparativo e uma análise de sensibilidade computacional entre duas abordagens matemáticas de propagação de trincas por fadiga: a clássica formulação linear da Lei de Paris e o modelo sigmoidal completo da equação NASGRO. Para avaliar a divergência numérica entre os critérios sob diferentes perfis de tensões, foram desenvolvidos e implementados dois estudos de caso em ambiente MATLAB através de um algoritmo de integração ciclo a ciclo com passo adaptativo. O primeiro cenário modela uma tubulação de gás natural pressurizada (aço API 5L X70) submetida a tensões de tração uniforme (hoop stress), enquanto o segundo avalia uma viga estrutural sob gradiente linear de flexão. Os resultados da varredura paramétrica demonstraram que, para o duto de gás, o modelo NASGRO prolonga significativamente a vida útil calculada devido à incorporação do efeito protetivo do fechamento de trinca induzido pela plasticidade (crack closure), gerando predições até 55% menos conservadoras que a Lei de Paris e evitando eventuais paradas de manutenção prematuras. Inversamente, no caso da viga sob flexão, a equação NASGRO revelou-se indispensável para capturar a desaceleração geométrica do fator de intensidade de tensão, identificando regimes de parada de trinca (crack arrest) induzidos pelo limiar de fadiga (ΔKth), de forma a mitigar o risco de condenações prematuras de estruturas com defeitos estacionários. Conclui-se que o algoritmo desenvolvido em MATLAB oferece uma ferramenta de baixo custo computacional e alta fidelidade, fornecendo suporte científico robusto para a transição de critérios tradicionais de projeto para metodologias avançadas de tolerância ao dano.
Palavras-chave: Mecânica da Fratura Elástica Linear; Lei de Paris; Equação NASGRO; Sensibilidade Computacional; Tolerância ao Dano.

ABSTRACT
The prediction of the remaining useful life of structures and devices subjected to cyclic loadings is a critical pillar for integrity management based on the damage tolerance philosophy. This paper presents a comparative study and a computational sensitivity analysis between two mathematical approaches for fatigue crack propagation: the classic linear formulation of the Paris Law and the full sigmoidal model of the NASGRO equation. To evaluate the numerical divergence between the criteria under different stress profiles, two case studies were developed and implemented in a MATLAB environment using a cycle-by-cycle integration algorithm with an adaptive step. The first scenario models a pressurized natural gas pipeline (API 5L X70 steel) subjected to uniform tensile stresses (hoop stress), while the second evaluates a structural beam under a linear bending gradient. The parametric sweep results demonstrated that, for the gas pipeline, the NASGRO model significantly extends the calculated useful life due to the incorporation of the protective effect of plasticity-induced crack closure, generating predictions up to 55% less conservative than the Paris Law and avoiding premature maintenance shutdowns. Conversely, in the case of the beam under bending, the NASGRO equation proved to be indispensable for capturing the geometric deceleration of the stress intensity factor, identifying crack arrest regimes induced by the fatigue threshold (ΔKth), thereby mitigating the risk of premature condemnation of structures with stationary defects. It is concluded that the algorithm developed in MATLAB offers a low computational cost and high-fidelity tool, providing robust scientific support for the transition from traditional design criteria to advanced damage tolerance methodologies.
Keywords: Linear Elastic Fracture Mechanics; Paris Law; NASGRO Equation; Computational Sensitivity; Damage Tolerance.

1. INTRODUÇÃO

O fenômeno da fadiga mecânica em componentes estruturais e sistemas de tubulações industriais representa uma das principais causas de falhas catastróficas em serviço, exigindo abordagens rigorosas de engenharia orientadas pela filosofia de tolerância ao dano. Tradicionalmente, o projeto de componentes submetidos a carregamentos cíclicos, a exemplos de redes de transporte e distribuição de gás natural ou vigas submetidas a carregamentos dinâmicos diversos, fundamenta-se em abordagens conservadoras baseadas na vida total ou na mecânica da fratura elástica linear (MFLE). No entanto, a integridade estrutural dessas instalações é constantemente desafiada pela nucleação e propagação de trincas superficiais, cuja evolução sob flutuações de pressão ou tensões mecânicas alternadas necessita de ferramentas de modelagem computacional precisas para mitigar falhas e estender a vida útil remanescente

Historicamente, a formulação proposta por Paris e Erdogan (1963) — amplamente conhecida como a Lei de Paris — consolidou-se como o modelo pioneiro e mais difundido para quantificar a taxa de crescimento de trincas por fadiga (da/dN) em função da variação do fator de intensidade de tensão (ΔK). Apesar de sua reconhecida simplicidade, a abordagem de Paris adota um perfil estritamente linear em escala bi logarítmica, limitando sua validade matemática à Região II do regime de propagação [5]. Por ser uma aproximação linear simplificada, a Lei de Paris falha em descrever a evolução da trinca nas extremidades da resposta mecânica do material, conforme apontado em revisões recentes da literatura sobre a evolução do dano estrutural [Mourad, 2025]: a Região I, onde o crescimento da trinca cessa ao se aproximar do limiar de fadiga (ΔK_th), e a Região III, associada à aceleração instável que precede a ruptura catastrófica do material.

Para superar tais limitações e possibilitar uma modelagem sigmoidal completa, a agência espacial norte-americana (NASA) desenvolveu a equação NASGRO [Forman, 1998]. Este modelo avançado introduz fatores de correção não-lineares que incorporam o fenômeno do fechamento de trinca induzido pela plasticidade (crack closure) e o efeito da razão de tensões (R), parametrizando com fidelidade todo o histórico de propagação do defeito. Pesquisas recentes demonstram que em estruturas críticas de dutos de aço carbono (como os aços de classe API X70), por exemplo, as flutuações de pressão promovem gradientes de tensão complexos que alteram severamente a resposta ao fenômeno da fadiga, exigindo modelos matemáticos capazes de mapear o dano sob condições complexas de contorno [Maachou, 2024].

Contudo, apesar dos benefícios teóricos do modelo NASGRO, a engenharia prática frequentemente opta pela utilização da Lei de Paris devido à escassez de dados experimentais necessários para calibrar os múltiplos expoentes empíricos exigidos pelo algoritmo da NASA. No caso de vigas estruturais expostas a flexões cíclicas ou tubulações de gás natural operando sob pressões pulsantes, a divergência matemática entre usar uma extrapolação linear simples e uma modelagem sigmoidal completa pode resultar em resultados subdimensionados, comprometendo a segurança operacional. Dessa forma, quantificar a sensibilidade computacional e mapear os desvios acumulados entre esses dois critérios constitui uma lacuna de conhecimento vital para o avanço dos algoritmos de gerenciamento da integridade estrutural.

Diante desse cenário, o presente artigo apresenta um estudo comparativo focado na análise de sensibilidade computacional dos modelos de Paris e NASGRO. Utilizando duas frentes de aplicação como estudo de caso — uma tubulação de gás natural sob regime cíclico e uma viga estrutural submetida a carregamentos transversais alternados —, este trabalho desenvolve uma rotina numérica integrada em ambiente MATLAB para avaliar o impacto das variáveis do fechamento de trinca e proximidade dos limites físicos (KI_c) e (ΔK_th) na curva de evolução do dano. Ao final, os resultados numéricos são confrontados para estabelecer diretrizes operacionais sobre os limites de segurança de cada abordagem, servindo como ferramenta de suporte para tomadas de decisão em projetos de engenharia baseados na tolerância ao dano.

2. FUNDAMENTOS DA MECÂNICA DA FRATURA TOLERANTE AO DANO

O gerenciamento da integridade estrutural e a previsão da vida útil de componentes baseiam-se na filosofia de tolerância ao dano. Esta abordagem assume que as estruturas contêm defeitos ou trincas pré-existentes oriundos de processos de fabricação, soldagem ou fadiga mecânica, estabelecendo que tais descontinuidades podem propagar-se de forma segura em serviço, desde que o tamanho da trinca seja rigorosamente monitorado e mantido abaixo de um limite crítico admissível.

2.1. 2.1 Mecânica da Fratura Elástica Linear (MFEL) Vs. Representação Sigmoidal Completa

A formulação matemática clássica para quantificar a taxa de crescimento de trincas por fadiga (da/dN), expressa em metros por ciclo (m/ciclo), baseia-se na Mecânica da Fratura Elástica Linear (MFLE). A premissa central estabelece que o campo de tensões e deformações na vizinhança da ponta da trinca é governado pelo fator de intensidade de tensão (K). Sob carregamentos cíclicos, o parâmetro que impulsiona o dano corresponde à variação do fator de intensidade de tensão (ΔK), calculado conforme a Equação (1):$$∆K=K_{max}-K_{min}$$

A modelagem proposta por Paris e Erdogan (1963) assume uma dependência linear e direta em escala bi logarítmica entre a taxa de crescimento e a variação do fator de intensidade de tensão, sendo expressa por:$$\frac{da}{dN}=C.{\left(∆K\right)}^m$$

Onde C e m são constantes empíricas obtidas por ensaios de laboratório sob condições controladas (ASTM E647). Embora a Lei de Paris forneça uma representação robusta para a Região II (regime linear estável), ela impõe um erro físico ao omitir os limites termodinâmicos e mecânicos do material, extrapolando o crescimento linear indefinidamente para regiões onde as restrições microestruturais alteram o regime de fluxo.

A transição para uma representação sigmoidal completa torna-se, dessa forma, necessária para capturar o comportamento físico em toda a curva sigmoidal, abrangendo as Regiões I, II e III. Na Região I, o crescimento da trinca desacelera de forma abrupta à medida que (ΔK) se aproxima do limiar de intensidade de tensão de fadiga (ΔKth), abaixo do qual o material é considerado imune à propagação macroscópica por fadiga. Na Região III, ocorre uma aceleração descontrolada do dano decorrente da proximidade da fratura instável iminente, que se inicia quando o fator de intensidade máximo (K_max) converge para a tenacidade à fratura por deformação plana do material (K_Ic). O modelo NASGRO equaciona matematicamente essas não-linearidades integrando os limites assintóticos por meio da seguinte formulação:$$\frac{da}{dN}=C.{\left(∆K\right)}^m{\left(\frac{1-F}{1-R}\right)}^n.\frac{{\left(1-\frac{∆K_{th}}{∆K}\right)}^p}{{\left(1-\frac{K_{max}}{K_{Ic}}\right)}^q}$$

A robustez da equação NASGRO em relação à Lei de Paris fundamenta-se na introdução do termo ((1-F)/(1-R))ⁿ, o qual atua como um fator de escala responsável pela calibração da sensibilidade da velocidade de propagação de acordo com a razão de tensões (R = K_min/K_max), bem como pelo fator de fechamento de trinca de Newman (F). O parâmetro F quantifica o nível de tensão necessário para abrir as superfícies da trinca que permanecem em contato em uma porção do ciclo de descarregamento devido à plasticidade residual acumulada atrás da ponta do defeito (crack closure), permitindo que o modelo preveja com precisão, do atraso na propagação sob flutuações complexas de carga.

2.2. 2.2 Caracterização dos Parâmetros de Entrada e de Material nos Estudos de Caso

Para quantificar a sensibilidade e os desvios dessas formulações, as constantes constitutivas das equações de Paris e NASGRO devem ser mapeadas de acordo com as propriedades microestruturais e as condições de contorno geométricas de dois estudos de caso distintos: nomeadamente uma tubulação de gás natural pressurizada e uma viga sob flexão alternada.

Estudo de Caso I – Tubulação de Gás Natural Pressurizada

Para o duto de gás natural de aço carbono API 5L X70 (material com tensão de escoamento = 485 MPa), com diâmetro externo de 323,8 mm e espessura de parede de 12,7 mm, assume-se a presença de uma trinca longitudinal interna semielíptica na parede do tubo, conforme representado pela Figura 1. A modelagem desta estrutura justifica-se pelo fato de que tubulações industriais enterradas operam sob regimes de pressões altamente oscilantes e pulsantes, reguladas por normas estritas de segurança como a ASME B31.8 [ASME, 2022].

Sob este arranjo geométrico, o fator de intensidade de tensão (K) varia ao longo da profundidade da trinca (a) e depende diretamente do diâmetro externo da tubulação (D_e), da espessura de parede da tubulação (t) e da pressão interna dinâmica cíclica (ΔP), sendo formulado conforme as diretrizes do código API 579-1/ASME FFS-1, por:$$∆K_I=\left[M_m{M}_1+M_2{\left(\frac{a}{t}\right)}^2+M_3{\left(\frac{a}{t}\right)}^4\right]\frac{∆P\left(De-t\right)}{2t}$$

Onde M_m, M₁, M₂ e M₃ são fatores de forma geométricos tabelados dependentes da razão entre a profundidade e o comprimento total da trinca (a/c). Como o duto opera sob o regime de pressão pulsante estabelecido pelas flutuações de demanda da rede de gás, a variação de pressão (ΔP = P_max – P_min) gera uma variação correspondente (ΔK_I), ativando as constantes do material em alta pressão. O comportamento de aços de classe API X70 sob carregamentos de espectro variável tem sido amplamente investigado por Maachou et al. (2024), cujos dados experimentais validam que as transições de amplitude alteram a plasticidade na ponta do defeito. Adicionalmente, o parâmetro F na equação NASGRO simula, nesse cenário, a compressão do gás dentro das paredes da trinca durante o ciclo de alívio, modificando a tensão de abertura real.

Estudo de Caso II – Viga Estrutural sob Flexão Alternada

Para a viga de secção transversal retangular (b x h = 50 x 100 mm) submetida a um carregamento transversal cíclico gerando um momento fletor alternado (ΔM), modela-se uma trinca de borda transversal reta de profundidade (a). A solução de (ΔK_I) para essa geometria segue a formulação clássica elástica dada por:$${∆K}_I=Y\frac{6∆M}{b{h}^2}\sqrt{\pi a}$$

Onde Y é o fator geométrico adimensional de intensidade para corpos sob flexão pura, calculado polinomialmente como uma função da profundidade relativa da trinca (α = a/h), por:$$Y=1,122-1,40\left(\frac{a}{h}\right)+7,33{\left(\frac{a}{h}\right)}^2-13,08{\left(\frac{a}{h}\right)}^3+14,0{\left(\frac{a}{h}\right)}^4$$

Diferente do caso da tubulação, onde o carregamento de pressão gera uma tensão de tração hidrostática uniforme (hoop stress), a viga sob flexão experimenta um gradiente linear de tensões que passa de tração na fibra externa para compressão na zona oposta. No escoamento plástico da ponta da trinca sob flexão, as forças compressivas cíclicas aceleram a taxa de fechamento de trinca, tornando os coeficientes de assíntota da região III (q) e do limiar (p) do modelo NASGRO altamente sensíveis ao decaimento do momento à medida que o defeito migra para regiões internas da viga.

3. IMPLEMENTAÇÃO NUMÉRICA E ALGORITMO

A avaliação da sensibilidade computacional entre uma abordagem estritamente linear (Paris) e uma formulação sigmoidal completa (NASGRO) requer a discretização das equações diferenciais constitutivas em um modelo discreto integrado. Como o fator de intensidade de tensão (K_I) é uma função transcendental dependente da geometria instantânea da trinca (a), a evolução temporal do dano ciclo a ciclo (ΔN) não possui solução analítica exata para geometrias complexas. Por esse motivo, desenvolveu-se uma rotina de integração numérica em ambiente MATLAB baseada no Método de Euler Unidirecional Modificado, otimizada para monitorar o acúmulo de dano sem incorrer em instabilidades numéricas nas proximidades das assíntotas físicas.

3.1. 3.1 Arquitetura Lógica do Algoritmo e Critério de Convergência

A arquitetura do algoritmo executa uma varredura incremental baseada no histórico de carregamento. O fluxo lógico do programa foi estruturado em quatro blocos principais de processamento sequencial, conforme apresentado na Tabela 1:

Tabela 1. Fluxo Lógico da Arquitetura de Processamento.

1. Bloco de Inicialização e Entrada de Dados - o programa realiza a leitura das constantes elásticas e geométricas dos materiais (E, σ_esc, K_Ic, ΔK_th), além de configurar os parâmetros geométricos iniciais dos estudos de caso, tais como o tamanho da trinca inicial (a₀), a razão de tensões (R), a variação de pressão (ΔP) na tubulação de gás natural e a variação de momento (ΔM) na viga estrutural.

2. Mecânica do Loop de Integração: o avanço do tamanho da trinca é calculado discretamente atualizando a profundidade instantânea (a_j+1) a cada incremento de blocos de ciclos (ΔN) através da seguinte relação recursiva:$$a_{j+1}=a_j+{\left(\frac{da}{dN}\right)}_j∆N$$

Para garantir a convergência e evitar erros de truncamento numérico, o algoritmo adota um critério de passo adaptativo onde (ΔN = 1) ciclo para regiões de alto gradiente (Região III) e (ΔN = 100) ciclos na zona de estabilidade linear (Região II).

3. Mapeamento de Restrições e Tratamento de Singularidades: a principal contribuição do código está no monitoramento em tempo real do termo assintótico do modelo NASGRO. Nas proximidades da falha imediata, quando (K_max → K_Ic), o denominador da equação tende a zero, fazendo com que a derivada (da/dN) tenda para o infinito. O algoritmo previne a quebra de execução do código inserindo uma condicional lógica de interrupção baseada no limiar: se a razão (K_max/K_Ic ≥ 0,98), a integração numérica é cessada, etapa na qual se consolida o registro do tamanho final da trinca e se computa o número de ciclos totais até a falha (N_f).

4. Bloco de Saída e Pós-Processamento: após a finalização do laço (seja pelo atingimento de a_critico ou pela condição de trava descrita na etapa anterior), o algoritmo compila os dados históricos de (ΔK_I) versus (da/dN) e exporta as curvas de evolução espacial (a vs. N) em gráficos vetoriais com escala adequada.

3.2. Fluxograma de Execução do Código MATLAB

O comportamento dinâmico e o controle de fluxo das variáveis dentro do ambiente de simulação são regidos pelo algoritmo detalhado na Tabela 2.

Tabela 2. O Algoritmo em MATLAB.

%=========================================================================%
        %          ANÁLISE COMPARATIVA - LEI DE PARIS vs. CRITÉRIO NASGRO         %
        %=========================================================================%
        clear; clc; close all;
        %% 1. PARÂMETROS DO MATERIAL E CONFIGURAÇÃO
        max_iter = 50000;
        a_inicial = 0.002; a_critica = 0.012;
        KIc = 34.0; DK_th = 3.2;
        % Parâmetros dos Modelos
        C_paris = 1.5e-11; m_paris = 3.2;
        C_nasgro = 1.3e-11; n = 3.2; p = 0.5;
        q = 0.5; F = 0.22;
        % Geometrias básicas
        De = 0.3238; t = 0.0127; % Caso I (Duto de Gás)
        b = 0.05; h = 0.1;       % Caso II (Viga Estrutural)
        % Definição da matriz de carregamentos para os 3 cenários de R
        % Linha 1: R = 0.0 (Crítico) | Linha 2: R = 0.1 (Padrão) | Linha 3: R = 0.4 (Suave)
        cenarios_duto = [5.0, 0.0; 5.0, 0.5; 5.0, 2.0];
        % Momentos fletores de entrada em N.m
        cenarios_viga = [15000, 0; 15000, 1500; 15000, 6000];
        num_cenarios = size(cenarios_duto, 1);
        %% 2. PROCESSAMENTO AUTOMÁTICO DO CASO I: TUBULAÇÃO DE GÁS NATURAL
        fprintf('=======================================================\n');
        fprintf(' PROCESSANDO CASO I: TUBULAÇÃO DE GÁS NATURAL \n');
        fprintf('=======================================================\n');
        for k = 1:num_cenarios
            Pmax = cenarios_duto(k, 1);
            Pmin = cenarios_duto(k, 2);
            R_local = Pmin / Pmax;
            
            % --- Simulação Paris
            a_paris = zeros(1, max_iter); N_paris = zeros(1, max_iter);
            a_paris(1) = a_inicial; j_p = 1;
            
            while (a_paris(j_p) < a_critica) && (j_p < max_iter)
                Delta_K = calcular_K_duto(a_paris(j_p), Pmax, Pmin, De, t);
                Kmax = Delta_K / (1 - R_local);
                if Delta_K <= DK_th || Kmax >= KIc, break; end
                
                da_dN_paris = C_paris * (Delta_K)^m_paris;
                dN = (da_dN_paris > 1e-6) * 1 + (da_dN_paris <= 1e-6) * 100;
                
                a_paris(j_p+1) = a_paris(j_p) + da_dN_paris * dN;
                N_paris(j_p+1) = N_paris(j_p) + dN;
                j_p = j_p + 1;
            end
            a_paris = a_paris(1:j_p); N_paris = N_paris(1:j_p);
            
            % --- Simulação NASGRO
            a_nasgro = zeros(1, max_iter); N_nasgro = zeros(1, max_iter);
            a_nasgro(1) = a_inicial; j_n = 1;
            
            while (a_nasgro(j_n) < a_critica) && (j_n < max_iter)
                Delta_K = calcular_K_duto(a_nasgro(j_n), Pmax, Pmin, De, t);
                Kmax = Delta_K / (1 - R_local);
                if Delta_K <= DK_th || Kmax >= KIc, break; end
                
                termo_escala = C_nasgro * ((1 - F)/(1 - R_local))^n;
                termo_limiar = (1 - (DK_th / Delta_K))^p;
                termo_ruptura = (1 - (Kmax / KIc))^q;
                
                da_dN_nasgro = termo_escala * (Delta_K^n) * (termo_limiar / termo_ruptura);
                dN = (da_dN_nasgro > 1e-6) * 1 + (da_dN_nasgro <= 1e-6) * 100;
                
                a_nasgro(j_n+1) = a_nasgro(j_n) + da_dN_nasgro * dN;
                N_nasgro(j_n+1) = N_nasgro(j_n) + dN;
                j_n = j_n + 1;
            end
            a_nasgro = a_nasgro(1:j_n); N_nasgro = N_nasgro(1:j_n);
            
            % --- Geração da Janela Gráfica para este R ---
            figure('Color', 'w', 'Units', 'inches', 'Position', [1, 1, 6.0, 4.0]);
            plot(N_paris, a_paris*1000, '--', 'Color', '#d9534f', 'LineWidth', 2, 'DisplayName', 'Modelo de Paris');
            hold on;
            plot(N_nasgro, a_nasgro*1000, 'Color', '#1f4e79', 'LineWidth', 2, 'DisplayName', 'Modelo NASGRO');
            grid on; gca.GridAlpha = 0.25;
            xlabel('Número de Ciclos Acumulados (N)', 'FontSize', 10);
            ylabel('Profundidade da Trinca a (mm)', 'FontSize', 10, 'FontWeight', 'bold');
            title({['Caso I: Duto de Gás - Curva de Resposta para R = ' num2str(R_local, '%.2f')], ...
                ['[P_{max} = ' num2str(Pmax) ' MPa | P_{min} = ' num2str(Pmin) ' MPa]']}, 'FontSize', 10, 'FontWeight', 'bold');
            legend('Location', 'southeast', 'FontSize', 9);
            hold off;
            
            fprintf('Cenário R = %.2f calculado. Vida Paris: %7d | Vida NASGRO: %7d\n', R_local, N_paris(end), N_nasgro(end));
        end
        %% 3. PROCESSAMENTO AUTOMÁTICO DO CASO II: VIGA ESTRUTURAL SOB FLEXÃO
        fprintf('\n=======================================================\n');
        fprintf(' PROCESSANDO CASO II: VIGA ESTRUTURAL SOB FLEXÃO \n');
        fprintf('=======================================================\n');
        for k = 1:num_cenarios
            Mmax = cenarios_viga(k, 1);
            Mmin = cenarios_viga(k, 2);
            R_local = Mmin / Mmax;
            
            % --- Simulação Paris
            a_paris = zeros(1, max_iter); N_paris = zeros(1, max_iter);
            a_paris(1) = a_inicial; j_p = 1;
            
            while (a_paris(j_p) < a_critica) && (j_p < max_iter)
                Delta_K = calcular_K_viga(a_paris(j_p), Mmax, Mmin, b, h);
                Kmax = Delta_K / (1 - R_local);
                if Delta_K <= DK_th || Kmax >= KIc, break; end
                
                da_dN_paris = C_paris * (Delta_K)^m_paris;
                dN = (da_dN_paris > 1e-6) * 1 + (da_dN_paris <= 1e-6) * 100;
                
                a_paris(j_p+1) = a_paris(j_p) + da_dN_paris * dN;
                N_paris(j_p+1) = N_paris(j_p) + dN;
                j_p = j_p + 1;
            end
            a_paris = a_paris(1:j_p); N_paris = N_paris(1:j_p);
            
            % --- Simulação NASGRO
            a_nasgro = zeros(1, max_iter); N_nasgro = zeros(1, max_iter);
            a_nasgro(1) = a_inicial; j_n = 1;
            
            while (a_nasgro(j_n) < a_critica) && (j_n < max_iter)
                Delta_K = calcular_K_viga(a_nasgro(j_n), Mmax, Mmin, b, h);
                Kmax = Delta_K / (1 - R_local);
                if Delta_K <= DK_th || Kmax >= KIc, break; end
                
                termo_escala = C_nasgro * ((1 - F)/(1 - R_local))^n;
                termo_limiar = (1 - (DK_th / Delta_K))^p;
                termo_ruptura = (1 - (Kmax / KIc))^q;
                
                da_dN_nasgro = termo_escala * (Delta_K^n) * (termo_limiar / termo_ruptura);
                dN = (da_dN_nasgro > 1e-6) * 1 + (da_dN_nasgro <= 1e-6) * 100;
                
                a_nasgro(j_n+1) = a_nasgro(j_n) + da_dN_nasgro * dN;
                N_nasgro(j_n+1) = N_nasgro(j_n) + dN;
                j_n = j_n + 1;
            end
            a_nasgro = a_nasgro(1:j_n); N_nasgro = N_nasgro(1:j_n);
            
            % --- Geração da Janela Gráfica para este R ---
            figure('Color', 'w', 'Units', 'inches', 'Position', [1, 1, 6.0, 4.0]);
            plot(N_paris, a_paris*1000, '--', 'Color', '#d9534f', 'LineWidth', 2, 'DisplayName', 'Modelo de Paris');
            hold on;
            plot(N_nasgro, a_nasgro*1000, 'Color', '#1f4e79', 'LineWidth', 2, 'DisplayName', 'Modelo NASGRO');
            grid on; gca.GridAlpha = 0.25;
            xlabel('Número de Ciclos Acumulados (N)', 'FontSize', 10);
            ylabel('Profundidade da Trinca a (mm)', 'FontSize', 10, 'FontWeight', 'bold');
            % AJUSTE: Exibição das unidades corrigidas nos títulos para N.m
            title({['Caso II: Viga - Curva de Resposta para R = ' num2str(R_local, '%.2f')], ...
                ['[M_{max} = ' num2str(Mmax) ' N\cdotm | M_{min} = ' num2str(Mmin) ' N\cdotm]']}, 'FontSize', 10, 'FontWeight', 'bold');
            legend('Location', 'southeast', 'FontSize', 9);
            hold off;
            
            fprintf('Cenário R = %.2f calculado. Vida Paris: %7d | Vida NASGRO: %7d\n', R_local, N_paris(end), N_nasgro(end));
        end
        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % FUNÇÕES LOCAIS GEOMÉTRICAS (Fatores de Intensidade de Tensão K)
        %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        function DK = calcular_K_duto(a_inst, Pmax, Pmin, De, t)
            DP = Pmax - Pmin;
            hoop_stress = (DP * (De - t)) / (2 * t);
            F_geo = 1.12;
            DK = F_geo * hoop_stress * sqrt(pi * a_inst);
        end
        function DK = calcular_K_viga(a_inst, Mmax, Mmin, b, h)
            DM = Mmax - Mmin;
            sigma_flexao = ((6 * DM) / (b * h^2)) / 1e6;
            alpha = a_inst / h;
            Y = 1.122 - 1.40*alpha + 7.33*alpha^2 - 13.08*alpha^3 + 14.0*alpha^4;
            DK = Y * sigma_flexao * sqrt(pi * a_inst);
        end
        

Os resultados obtidos são representados graficamente por meio das Figuras 2 a 4, para o caso da tubulação de gás natural pressurizada e das Figuras 5 a 7, para o caso da viga submetida a carregamentos cíclicos. Recomenda-se observar, nessas circunstâncias, as variações impostas para a razão de tensões (R), no transcorrer das simulações.

Tabela 3. Resumo Comparativo da Vida Útil Residual Predita (Nf) – Tubulação de Gás Natural.

4. RESULTADOS E DISCUSSÕES: EFEITO DA RAZÃO DE TENSÕES NA SENSIBILIDADE DOS MODELOS

A expansão do modelo numérico para múltiplos cenários operacionais da razão de tensões (R = 0,00; R = 0,10 e R = 0,40) consolidou a necessidade de uma abordagem não-linear, revelando fenômenos físicos distintos e opostos para cada frente de análise.

No Caso I (Tubulação de Gás Natural), observou-se que sob severa ciclicidade (R = 0,00), o modelo NASGRO projeta uma vida útil significativamente estendida em comparação à Lei de Paris (1.926.500 contra 649.400 ciclos). Esse comportamento quantifica a magnitude do efeito protetivo de fechamento de trinca induzido pela plasticidade (crack closure), que atua reduzindo o fator de intensidade de tensão efetivo (ΔK_eff).

Contudo, ao elevar a razão de tensões para R = 0,40, ocorre uma convergência quase exata entre as predições (uma diferença residual inferior a 2%, com 3.329.100 ciclos para Paris e 3.278.100 para NASGRO). Fisicamente, esse fenômeno explica-se pelo fato de que o aumento da contrapressão (ou pressão mínima de operação) mantém a trinca permanentemente aberta durante todo o ciclo de carregamento. Sob alta razão R, o efeito protetivo de fechamento é mitigado, e o dano passa a responder predominantemente ao regime linear da Região II, onde ambas as leis constitutivas compartilham parâmetros equivalentes.

Em contrapartida, os resultados do Caso II (Viga Estrutural sob Flexão) revelam uma inversão crítica de sensibilidade matemática. Sob baixas razões de tensões (R = 0,00) e (R = 0,10), o comportamento assemelha-se ao do duto, com o modelo NASGRO fornecendo predições de vida ligeiramente superiores devido à sutil atuação do fechamento mecânico. No entanto, ao analisar o cenário de flutuação suave sob torque e flexão moderada (R = 0,40), a divergência entre os modelos atinge sua divergência máxima: a Lei de Paris projeta uma vida útil de 115.200 ciclos, enquanto o modelo NASGRO aponta o colapso ou interrupção do regime estável em apenas 37.722 ciclos — uma discrepância drástica de mais de 200%.

Investigações numéricas semelhantes conduzidas por Carpinteri et al. (2013) sobre tubos e barras submetidos a carregamentos combinados de pressão e flexão confirmam que o gradiente de tensões altera significativamente os fatores de forma geométricos ao longo do avanço da trinca. Essa mudança de comportamento sob flexão pura foi validada recentemente por Maachou e Mechab (2026), que demonstraram que a aplicação isolada de modelos estritamente lineares falha em prever a transição elasto-plástica quando a frente da trinca migra através de regiões com gradientes de tensão não-uniformes.

Esse desvio massivo no regime de (R = 0,40) na viga é justificado pela interação entre o gradiente linear de tensões e as assíntotas de projeto da Região I e III. Conforme a trinca se aprofunda na secção retangular, caminhando em direção à linha neutra, a força motriz mecânica original (ΔK) decai de forma acentuada. No modelo NASGRO, essa redução conduz o sistema para a zona de influência direta do limiar de fadiga (ΔK_th). Paralelamente, como o momento máximo aplicado gera níveis de tensão localizados elevados na fibra extrema, o fator máximo (K_max) aproxima-se rapidamente do patamar de tenacidade à fratura do material (K_Ic).

O termo de ruptura da equação NASGRO, expresso por (1 - K_max/K_Ic)^q), assume o controle do balanço diferencial e acelera a taxa de crescimento (da/dN) de forma assintótica na iminência da falha. A Lei de Paris, por ser rigidamente linear e indiferente à proximidade do colapso estático (K_Ic) ou à proteção do limiar (ΔK_th), falha em perceber a degradação acelerada do ligamento remanescente da viga. Essa divergência paramétrica prova que a omissão das assíntotas elastoplásticas pela Lei de Paris pode induzir a erros de subdimensionamento não-conservativos em componentes operando sob gradientes de flexão combinados com altas razões de carregamento estático.

A justificativa física para o comportamento discrepante observado no Caso II (Viga) reside na transição geométrica do campo de forças. Enquanto a tubulação de gás sofre uma tensão de tração pura e uniforme em toda a espessura da parede (hoop stress), a viga sob flexão transversal experimenta um gradiente linear de tensões: tração máxima na fibra externa, decaindo até zero na linha neutra e passando para compressão na metade oposta. Essa transição geométrica altera drasticamente a sensibilidade dos modelos computacionais por dois motivos físicos e numéricos principais:

1. A desaceleração geométrica do (ΔK): à medida que a trinca se aprofunda na viga, ela caminha em direção a zonas de menor tensão axial de tração. O polinômio geométrico (Y) tenta compensar a profundidade, mas o decaimento da tensão real desacelera a taxa de crescimento primária.

2. A Ativação precoce da Região III (ruptura instável): devido à alta componente de tensão estática gerada pela combinação do momento máximo com a alta razão de tensões (R = 0,40), o fator máximo (K_max) aproxima-se rapidamente do patamar de tenacidade à fratura do material (K_Ic). O modelo NASGRO captura essa aproximação assintótica através do termo (1 - K_max/K_Ic)^q), acelerando drasticamente a taxa de crescimento nas camadas intermediárias da viga, o que resultou na falha prematura com apenas 37.722 ciclos. A Lei de Paris, por ser rigidamente linear, ignora esse limiar superior de colapso estrutural e prevê erroneamente que a trinca continuará propagando de forma lenta e estável até 115.200 ciclos, gerando um erro de previsão perigoso contra a segurança.

4.1. Análise de Sensibilidade dos Casos I e II

Diferente do comportamento observado na tubulação de gás natural, a avaliação da sensibilidade computacional no Caso II (viga estrutural sob flexão) revela uma inversão na tendência dos desvios, evidenciando o impacto do gradiente linear de tensões sobre as leis de propagação. Na viga, à medida que a profundidade da trinca (a) avança em direção ao núcleo geométrico da seção, ela migra da região de máxima tração fibrosa na fibra extrema para zonas de menor intensidade de tensão próximas à linha neutra.

Esse decaimento espacial do campo de tensões força o fator de intensidade de tensão induzido (ΔK) a decrescer de forma acentuada, posicionando a força motriz do sistema na zona de influência das assíntotas elasto-plásticas do material. O impacto dessa transição geométrica altera drasticamente a sensibilidade dos modelos computacionais por dois motivos físicos e numéricos principais:

1. A ativação do limiar (ΔK_th) em baixas razões de R: sob condições de fadiga de alto ciclo e amplitudes moderadas (R = 0,00 e R = 0,10), o decaimento do campo de tensões conduz o ΔK local para a zona de transição da Região I. O modelo NASGRO assimila imediatamente a assíntota do limiar através do termo (1 - ΔK_th/ΔK)^p). À medida que o ΔK aproxima-se de (ΔK_th), a taxa de crescimento (da/dN) converge para zero, caracterizando o fenômeno físico de aprisionamento ou parada da trinca (crack arrest). Por outro lado, a rigidez matemática da aproximação linear de Paris introduz uma simplificação acentuada neste cenário. Como a Lei de Paris desconsidera a existência de um limiar inferior protetivo, ela assume que mesmo variações infinitesimais de ΔK produzem dano contínuo, superestimando a velocidade de propagação nas fases intermediárias e resultando em predições de vida mais conservadoras (22.550 ciclos contra 35.255 do NASGRO para R = 0,00).

2. A ativação da região III (ruptura instável) em altas razões de R: uma inversão matemática drástica ocorre ao analisar o cenário de flutuação suave sob alta componente estática (R = 0,40). Em virtude do elevado momento fletor máximo combinado (15.000 N·m), o fator de intensidade máximo (K_max) aproxima-se rapidamente do patamar de tenacidade à fratura do material (K_Ic). O modelo NASGRO captura essa aproximação assintótica através do termo (1 -  K_max/K_Ic)^q), acelerando violentamente a taxa de crescimento da/dN na iminência do colapso, o que resulta em uma falha prematura com apenas 37.722 ciclos. A Lei de Paris, por sua linearidade rígida, ignora esse limiar superior de colapso estrutural e assume erroneamente que a trinca continuará propagando de forma lenta e estável até 115.200 ciclos, gerando um erro de predição não-conservativo que compromete a segurança estrutural.

Essa divergência paramétrica bidirecional prova que a omissão das assíntotas elastoplásticas superiores e inferiores pela Lei de Paris pode induzir tanto a penalizações econômicas por conservadorismo excessivo quanto a erros de subdimensionamento catastróficos, a depender estritamente da interação entre a razão de tensões (R) e o gradiente de forças do componente.

5. CONCLUSÕES

A análise comparativa e de sensibilidade computacional realizada neste trabalho permitiu quantificar as divergências estruturais entre a abordagem linear clássica da Lei de Paris e a modelagem sigmoidal avançada da equação NASGRO. A partir dos ensaios numéricos aplicados aos dois estudos de caso práticos, extraem-se as seguintes conclusões fundamentais:

1. Efeito protetivo do fechamento de trinca no duto de gás: para sistemas mecânicos governados por campos de tração uniforme e contrapressão, como a tubulação de gás natural pressurizada sob regime cíclico (hoop stress), o modelo NASGRO demonstrou prolongar significativamente a estimativa de vida útil residual do componente. Esse comportamento decorre da incorporação matemática do fenômeno de fechamento de trinca induzido pela plasticidade (crack closure). Ao assimilar o atraso na taxa de propagação gerado pelo contato das faces do defeito nas fases de alívio de pressão, o modelo da NASA evita condenações prematuras de engenharia e paradas de manutenção antecipadas de linhas operacionais que ainda guardam margem de segurança estática.

2. Efeito assintótico do limiar de fadiga na viga sob flexão: no cenário avaliado para a viga estrutural, o gradiente linear de tensões impõe uma desaceleração geométrica do fator de intensidade de tensão (ΔK) à medida que a trinca migra para a linha neutra. Nesse contexto, a equação NASGRO revelou-se um instrumento indispensável para identificar regimes de parada de trinca (crack arrest), impulsionada pelo termo assintótico do limiar de fadiga (ΔK_th) em baixas razões de R, bem como para prever a aceleração catastrófica por ruptura instável (região III) sob altas componentes estáticas de carregamento. Enquanto a rigidez linear de Paris assume erroneamente uma propagação contínua e destrutiva mesmo sob forças motrizes infinitesimais, o modelo sigmoidal completo mapeia a estabilização real do dano, impedindo o sobredimensionamento de intervenções e custos de manutenção em estruturas com defeitos estacionários, ao mesmo tempo em que mitiga riscos catastróficos ao identificar falhas prematuras induzidas pela aproximação da tenacidade à fratura (K_Ic).

3. Viabilidade da implementação computacional proposta: a rotina integrada desenvolvida em ambiente MATLAB consolidou-se como uma plataforma robusta, de baixo custo computacional e de alta fidelidade analítica. O algoritmo adaptativo provou ser capaz de gerenciar com precisão as singularidades matemáticas das Regiões I e III sem incorrer em instabilidades numéricas, oferecendo uma solução viável e acessível para a transição de critérios tradicionais de projeto para ferramentas avançadas de gerenciamento de integridade por tolerância ao dano.

REFERÊNCIAS BIBLIOGRÁFICAS

AMERICAN SOCIETY OF MECHANICAL ENGINEERS (ASME), ASME B31.8: Gas Transmission and Distribution Piping Systems, ASME, New York, 2022.

CARPINTERI, A., et al., Fatigue crack propagation analysis of round pipes with external surface flaws under cyclic pressure and bending, Journal of Pressure Vessel Technology 135 (3) 031402, 2013.

FORMAN, R.G., et al., Fatigue crack growth computer program NASGRO Version 3.0 - Reference Manual, NASA Johnson Space Center, Houston, TX, Report No. JSC-22267, 1998.

MAACHOU, S., et al., Fatigue crack growth behavior of API X70 steel pipelines under variable amplitude and spectrum loading, Materials Science and Engineering: A 892 145020, 2024.

MAACHOU, S., MECHAB, B. A new investigation of three fatigue models for assessing crack propagation in pipelines under bending loading, Journal of the Serbian Society for Computational Mechanics 19 (2) 48-61, 2026.

MOURAD, A.H.I., et al., Fatigue life and crack growth prediction of metallic structures: A review, Structures 76 109031, 2025.

PARIS, P., ERDOGAN, F. A critical analysis of crack propagation laws, Journal of Basic Engineering 85 (4) 528-533, 1963.

¹ Docente do Programa de Pós-Graduação em Métodos Numéricos em Engenharia – UFPR. E-mail: [clique para visualizar o e-mail]acesse o artigo original para visualizar o e-mail

² Docente do Curso de Engenharia Mecânica - UFPR. E-mail: [clique para visualizar o e-mail]acesse o artigo original para visualizar o e-mail

³ Discente do Programa de Pós-Graduação em Métodos Numéricos em Engenharia – UFPR. E-mail: [clique para visualizar o e-mail]acesse o artigo original para visualizar o e-mail