Utilização de modelos geométricos simplificados de uma córnea para fins de otimização

Resumo

Otimização por Enxame de Partículas (PSO) é uma técnica de inteligência artificial (AI), que pode ser usada para encontrar soluções aproximadas para problemas numéricos de maximização e minimização extremamente difíceis. Neste trabalho, utilizou-se um algoritmo PSO para comparar os deslocamentos sofridos por uma amostra de córnea humana submetida à uma pressão interna de 45 mmHg com resultados de simulações numéricas e identificar valores otimizados para propriedades hiperelásticas da córnea (µ e α). Por meio dos resultados das simulações via análise inversa pelo Método dos Elementos Finitos (MEF), em conjunto com o algoritmo PSO, foram encontrados valores otimizados de µ = 0,047 e α = 106,7. Quando comparado com resultados otimizados por meio de um software comercial, foram encontrados erros de aproximadamente 0,15%. Por meio dos resultados obtidos, verificou-se ainda que, variando os valores dos coeficientes de inércia da partícula no algoritmo PSO, os resultados podem sofrer ligeira melhoria, o que demonstra potencial uso do PSO em conjunto com análise inversa do MEF para caracterização de materiais hiperelásticos, utilizando modelos geométricos simplificados

Descritores:
Otimização por enxame de partículas; Método dos elementos finitos; Córnea; Humano

Abstract

Particle Swarm Optimization (PSO) is an artificial intelligence technique (AI) that can be used to find approximate solutions to numerical problems of maximization and minimization. In this study, it was used a PSO algorithm to compare displacements from human cornea sample subjected to internal pressure of 45 mmHg with Results of numerical simulations were provided which identified optimized values for hyperelastic properties of the cornea (µ and α). By means of the results from numerical simulations via inverse analysis by the Finite Element Method (FEM), in conjunction with the PSO algorithm, optimized values of µ = 0.047 and α = 106.7 were found. When compared with optimized results from commercial software, errors around 0.15% were found. Results showed that, varying the values of particle inertia coefficients in the PSO algorithm, simulated displacements have improved when compared to experimental data. This demonstrates the potential use of PSO algorithm in conjunction with the FEM inverse analysis for hyperelastic materials characterization, using simplified geometrical models

Keywords:
Particle Swarm optimization; Finite element method; Cornea; Human

Introdução

O desenvolvimento e aprimoramento de técnicas de otimização permite estudar o comportamento em vários materiais e é de suma importância para determinar suas limitações, avaliar e classificar suas propriedades. Assim, há um interesse muito grande na avaliação das propriedades biomecânicas da córnea, sendo considerado um parâmetro importante a ser determinado, uma vez que está relacionado a diversos procedimentos de diagnósticos.

A córnea é a primeira e a mais poderosa superfície de refração do sistema óptico do globo ocular. Pouco se conhece sobre o comportamento da estrutura interna da córnea. Na maior parte das espécies animais, a córnea é constituída por quatro camadas: epitélio, estroma, membrana de Descemet e endotélio.⁽¹⁾ Uma pequena alteração na geometria da córnea notavelmente afeta a sua capacidade óptica. Em decorrência dessa sensibilidade, estudos biomecânicos podem revelar o desempenho e funcionalidade da córnea, a qual possui determinadas características mecânicas que podem ser analisadas com o intuito de aprimorar os métodos para prevenção de doenças e o desenvolvimento de equipamentos mais eficientes na área médica.

O crescente interesse pelo estudo de tecidos biológicos humanos levou a uma necessidade da caracterização de alguns tipos de tecidos com um comportamento mecânico hiperelástico. Experimentalmente, o estudo mecânico desse tipo de tecido tem sido executado por meio de técnicas tradicionais, desenvolvidas para o estudo de materiais com um comportamento linear e isotrópico, o que permite a determinação de valores médios de algumas propriedades mecânicas.

Esses materiais têm comportamentos mecânicos distintos dos materiais convencionais, havendo a necessidade de recorrer-se a algoritmos que viabilizem a solução desses problemas, sendo um deles, a Otimização por enxame de partículas (Particle Swarm Optimization). Soares⁽²⁾ define otimização, como sendo um mecanismo de análise de decisões complexas, envolvendo seleção de valores para variáveis, com o simples objetivo de quantificar o desempenho e medir a qualidade das decisões.

O algoritmo implementado faz uso da deformação sofrida por um corpo-de-prova (córnea humana), tendo como objetivo, identificar as propriedades do material hiperelástico. A caracterização do material testado foi realizada com o auxílio do Método dos Elementos Finitos.

Referencial Teórico

O método de Otimização por enxame de partículas (Particle Swarm Optimization) foi desenvolvido por Kennedy et al.⁽³⁾ a partir da modelagem matemática que implementa uma metáfora do comportamento social de um grupo de pássaros à procura de alimento ou de um lugar para construir o ninho, constituindo uma técnica de computação estocástica.

De acordo com Parsopouloset al.⁽⁴⁾, cada partícula é tratada como um ponto dentro do espaço de busca, que ajusta seu próprio "voo" de acordo com sua própria experiência, bem como a experiência do "voo" de outras partículas.

Já o Método de Elementos Finitos foi utilizado para simular diferentes propriedades mecânicas dos materiais e fornecer dados de deslocamentos em pontos específicos da geometria. Em decorrência do fato de as sub-regiões apresentarem dimensões finitas, estas sub-regiões são chamadas "elementos finitos", em contraste com os elementos infinitesimais utilizados no cálculo diferencial e integral. Advém daí o nome "Método dos Elementos Finitos", estabelecido por Ray Clough.⁽⁵⁾

Ogden⁽⁶⁾ propôs um modelo para descrever o comportamento não linear de materiais complexos tais como: borracha, polímeros e tecidos biológicos. Esse modelo é descrito por meio da Equação 1.

(1) $\Psi =\sum _{\mathit{i}=1}^{\mathit{N}}\frac{\mathit{mi}}{\mathit{ai}}\left({\mathit{I}}_1^{\mathit{ai}}+{\mathit{I}}_2^{\mathit{ai}}+{\mathit{I}}_3^{\mathit{ai}}-3\right)$

onde Ψ é a energia de deformação do material, N é o número total de termos da série e µ_ie α_i são constantes do material.

A função energia de deformação pode ser expandida, por meio de uma série de potências reais e descrita como função dos estiramentos principais demonstrados na Equação 2, com a seguinte condição de estabilidade, µnαn> 0:

(2) $\mathit{m}=\frac{1}{2}\sum _{\mathit{i}=1}^{\mathit{N}}{\mathit{m}}_{\mathit{n}}{\mathit{a}}_{\mathit{n}}$

onde µ é o módulo de cisalhamento e α é o parâmetro de capacidade de endurecimento da deformação (encruamento). µ é uma medida da capacidade de um material para resistir a deformações transversais (medido em KPa), é também conhecido como rigidez e é um parâmetro que quantifica a rigidez do material. α é representativa do comportamento não linear, em razão do endurecimento por deformação do material.

De acordo com Lapper et al.⁽⁷⁾, atualmente têm sido desenvolvidas novas metodologias para a simulação do compor-tamento mecânico hiperelástico de tecidos biológicos. Essas metodologias, geralmente, contemplam a utilização do Método de Elementos Finitos como uma ferramenta de simulações numéricas para a obtenção do comportamento numérico de materiais biológicos.

Métodos

Para os experimentos, buscou-se aliar a capacidade de análise do Método de Elementos Finitos à simplicidade e robustez do PSO. Para a obtenção dos dados de entrada no algoritmo PSO, foram utilizados dados experimentais oriundos da Universidade de Liverpool-UK, sendo possível definir a função objetivo. Nesses experimentos, os corpos-de-prova (córneas) foram fixados na borda, carregamentos foram aplicados com pressão interna de 45 mmHg, a fim de se medir, experimentalmente, os deslocamentos na região central da córnea (pupila) por meio de processamento digital de imagens. Todo o aparato experimental encontra-se ilustrado na figura 1.

Figura 1
Aparato experimental, adaptado de Elsheik.⁽⁸⁾ 1- Região de fixação da córnea; 2: Reservatório p/ aplicação da pressão; 3 – Controlador; 4 – Emissor de Laser; 5 - Câmeras CCD.

Com base nos dados experimentais para córneas humanas com idade de aproximadamente 50 anos, foi obtido o gráfico tempo vs. deslocamento.

A partir da curva gerada experimentalmente, definiu-se a função objetivo por um polinômio de 5º grau, conforme figura 2.

Figura 2
Gráfico da função objetivo (função polinomial de 5º grau).

O algoritmo utilizado, neste trabalho, descreve um problema de "análise inversa", que será solucionado via PSO. Objetivou-se, com este algoritmo gerar curvas de deformação próximas aos dados obtidos experimentalmente, a partir de valores máximos e mínimos de x₀ (primeira incógnita) e x₁ (segunda incógnita), como dados de entrada do algoritmo. Em geral, o PSO é adequado para os problemas numéricos com duas ou mais incógnitas, que é o caso deste trabalho.

O algoritmo PSO trabalha em um laço de repetição. Em cada iteração do laço, a partícula é atualizada com base no espaço de busca e na velocidade. O número de partículas é definido pelo usuário na busca de uma solução ótima em cada iteração. No algoritmo PSO proposto neste trabalho, o usuário entra com quatro valores, aqui denominados por "limites", que são valores máximos e mínimos das incógnitas a serem obtidas, neste caso, mínimo de µ, máximo de µ, mínimo de α e máximo de α, conforme fluxograma apresentado na figura 3.

Figura 3
Fluxograma PSO/FEM⁽⁹⁾.

A soma dos erros ao quadrado (SSE), Equação 3, foi obtida a partir da diferença entre os valores de deslocamento oriundos das simulações e os pontos gerados pela curva da função objetivo (Figura 2):

(3) $\mathrm{SSE}=\sqrt{\sum _1^n\frac{{\left(\mathrm{Pp}-D\right)}^2}{\mathit{n}}}$

onde P_p são os deslocamentos provenientes da equação polinomial da função objetivo, D são os deslocamentos simulados e n é o número de pontos analisados.

A geometria do olho humano foi gerada por um programa na Universidade de Liverpool, chamado de Ocular Model Generation, sendo que algumas geometrias foram modificadas, de acordo com os cenários simulados. Essas modificações se deram por intermédio do software Abaqus^®.

Neste trabalho, o MEF foi utilizado para simular diferentes propriedades mecânicas dos materiais e fornecer dados de deslocamentos em pontos específicos da geometria. Para isso, mediram-se os deslocamentos na região central da córnea (Figura 4b) em função de uma pressão aplicada na região interna da córnea, levando-se em consideração a região de fixação da mesma no globo ocular (Figura 4a). Os resultados de deslocamento foram utilizados pelo algoritmo PSO para comparar com valores da função objetivo em cada iteração.

Figura 4
Exemplo de dados de contorno em córneas para simulação via MEF.¹⁰⁾

Para gerar a geometria, cada um dos elementos foi mode-lado na configuração de 6 nós/elemento e 15 nós/elemento, usando uma forma tetraédrica.

O software Abaqus^® realiza a transformação da geometria C3D6H (6 nós) para a C3 D15H (15 nós), e nesse tipo de configu-ração, deve ser considerada a geometria na forma quadrática.

O arquivo InputFile.inp contém todas as informações que determinam a estrutura da córnea: o número de elementos, o número de nós, pressão, entre outras. Esse arquivo serve como dados de entrada para o software Abaqus^®.

O script Analysis_batch_file.bat é responsável por executar o software Abaqus ^® sem a necessidade de abrir sua interface gráfica, gerar arquivo contendo os dados de saída com extensão *.odb e executar o algoritmo Out_Node.py (desenvolvido na linguagem livre Python^®), que serve para gerar a curva simulada de deslocamento. O algoritmo PSO implementa todo esse processo, trabalhando em sinergia com o Abaqus^®.

O algoritmo PSO (desenvolvido na plataformaVisual Basic^®) executa o arquivo Analysis_batch_file.bat, o qual roda a simulação via Abaqus^® a partir de um arquivo contendo os dados de entrada com extensão *.inp (geralmente denominado por InputFile.inp) e gera um arquivo contendo dados de saída com extensão *.odb (arquivo contendo resultados do Abaqus^®) e executa o arquivo desenvolvido em Python^®, Out_Node.pyque usa o arquivo *.odb para gerar a curva tempo (ou pressão) vs. deslocamento.

Neste trabalho, a curva simulada (tempo vs. deslocamento) foi gerada por meio de um arquivo denominado "displacement.txt", o qual é um arquivo de texto que serve para ser comparado com a função objetivo no programa PSO.

Todo o processo de criação dos modelos geométricos em três dimensões foi gerado via software comercial Abaqus^® (versão 6.14), sendo utilizada sua interface gráfica para configuração e definição dos parâmetros utilizados no decorrer da simulação. Todo esse processo foi baseado no setup utilizado pela Universidade de Liverpool para a coleta dos dados experimentais e na aplicação Ocular Model Generation. O primeiro passo foi realizar a importação do arquivo InputFile.inp, gerado pelo programa Ocular Model Generation. Apesar de a córnea ser um material anisotrópico, a respectiva foi considerada nas simulações como sendo isotrópico, por motivos de simplificações.

Na figura 5, mostra-se a Malha de Elementos Finitos do Globo Ocular obtida pelo software Abaqus^®, composto de 772 nós, 768 elementos e possui regiões que representam toda a estrutura do olho humano.

Figura 5
Malha de elementos finitos do globo ocular.⁽¹⁰⁾

Considerando que, para obter-se a região correspondente da córnea, na figura 6, foi criada a malha de elementos finitos gerada por meio do software Abaqus^®, através do modelo discretizado da córnea de um olho humano a partir de dados geométricos oriundos do programa Ocular Model Generation.

Figura 6
Região discretizada da córnea.⁽¹⁰⁾

O modelo da córnea utilizado foi considerado hiperelástico, sendo o material definido segundo modelo de OGDEN (Equação 1). Foram considerados como referência valores de µ = 0,054 e α = 110,45, obtidos por meio de um software comercial de otimização (HEEDS).

Foi definida uma região interna da córnea, (Figura 7) para aplicação da pressão interna. Além disso, um ponto de monitoramento foi definido para servir como referencial no processo de coleta dos dados do deslocamento da córnea.

Figura 7
Região interna da córnea.⁽¹⁰⁾

Após todo o processo de configurações do cenário simulado, uma pressão com magnitude de 45 mm Hg foi aplicada na região interna da córnea, objetivando, assim, calcular o deslocamento sofrido no ponto de monitoramento quando esta sofre a respectiva pressão interna.

Para cada modelo geométrico gerado foi atribuído um cenário diferente, baseado em valores mínimos e máximos de nós, suportados pela versão do Abaqus^® utilizada, ou seja, 1º cenário ao 5º cenário (Tabela 1).

Tabela 1
Simulações realizadas

Cenários	Nº de Nós	Nº de Elementos	Nó/ elemento
1º cenário	74	54	6 nós/elem
2º cenário	291	54	15 nós/elem
3º cenário	455	108	15 nós/elem
4º cenário	753	150	15 nós/elem
5º cenário	938	864	6 nós/elem

Cada cenário teve seus valores de nós e elementos definidos pela geometria. Valores aleatórios do número de nós, entre o máximo e o mínimo foram atribuídos para os demais cenários.

Resultados e Discussão

Valores limites pré-estabelecidos

Para cada cenário analisado, valores limites pré-definidos para µ₁ e α₁ (Tabela 2) foram estabelecidos para se testar a estabilidade dos resultados em razão da variação dos valores de entrada no algoritmo PSO.

Tabela 2
Simulações realizadas com diferentes limites

PARÂMETROS	µ 1	α1
	LIMITES	LIMITES
	(MIN-MAX)	(MIN-MAX)
1	0.00001-1	20-400
2	0.00001-1	10-200
3	0.02-0.2	20-150
4	0.01-0.1	60-240
5	0.025-0.075	60-150
6	0.001-0.8	40-300

Esses parâmetros foram definidos no intuito de testar a estabilidade do algoritmo, sendo que o mesmo teste foi executado várias vezes, obtendo-se o mesmo resultado, até mesmo em diferentes computadores.

Observando as curvas apresentadas na figura 8, é possível notar que os valores começam a se estabilizar a partir da 40º iteração com valores otimizados próximos ao obtido pelo software comercial HEEDS (µ₁ = 0,054).

Figura 8
Comparação dos valores de μ₁ no Segundo cenário em todos os parâmetros.

Isso ocorre, porque no algoritmo PSO, os valores são gerados aleatoriamente em um espaço de busca, onde cada partícula, ou seja, cada valor corresponde a uma possível solução se desloque até que atinja o melhor resultado. Os valores foram variando e convergiram para valores aproximados dos obtidos pelo software comercial HEEDS.

Observando as curvas apresentadas na figura 9, é possível notar que os valores começam a estabilizar após a 35º iteração, mostrando resultados próximos ao obtido pelo software HEEDS (α₁ = 110,45). Assim como ocorreu com os valores de µ₁, confirmou-se que o segundo cenário foi o mais adequado para essa situação quando comparado com os resultados do software HEEDS.

Figura 9
Comparação dos valores de α₁ no Segundo cenário em todos os parâmetros

Erros dos cenários simulados

O melhor resultado (menor Erro) foi encontrado no 2º cenário / 5º parâmetro, com 291 nós. Na figura 10, apresentam-se os resultados do SSE para o 2º cenário - 291 nós, referente aos parâmetros simulados.

Figura 10
Comparação dos valores de Erro (SSE) no Segundo cenário em todos os parâmetros.

Para cada cenário analisado, os resultados foram obtidos após o programa rodar 50 iterações, notando assim sua estabilidade após 40 iterações. Por meio dessa análise, nota-se que o melhor valor de SSE foi obtido no 2º cenário / 5º parâmetro (aproximadamente 15%), onde o cenário utilizado possui 291 nós.

Conclusão

Comprovou-se que a estabilidade e convergência variando-se limites máximos e mínimos na técnica de análise inversa pelo Método dos Elementos Finitos (MEF), em conjunto com o algoritmo PSO, foram satisfatórios, pois por meio dos resultados obtidos, demonstrou o potencial uso do PSO em conjunto com análise inversa do MEF para caracterização de materiais hiperelásticos.

Valores de µ₁ =0,054 e α₁ = 110,45 foram obtidos por meio de um software comercial (HEEDS), licenciado e utilizado pela Universidade de Liverpool.

Para cenários simulados executados a partir de um número restrito de nós (máximo de 1000 nós por modelo geométrico) em função da versão do software Abaqus utilizada nas simulações.

Em relação ao PSO padrão x Inércia, o valor, mesmo de pouca relevância, leva a concluir que, perante experimentos, foram notados ainda, melhorias no PSO padrão, mas, baseado nesses parâmetros, esses valores foram insuficientes para uma correta avaliação.

Este trabalho conseguiu atingir os objetivos, e a continuação de seus experimentos, será de grande importância, implicando em melhorias tanto no campo de ciências médicas, como em engenharias. Na medicina, poderá contribuir no tratamento de diversas doenças oculares, em equipamentos mais eficientes de diagnósticos, ou mesmo procedimentos cirúrgicos. Poderá ainda, orientar na melhoria e criação de modelos matemáticos, atuando diretamente na melhoria da vida humana.

Para trabalhos futuros, propõe-se utilizar diferentes valores de inércia, ou seja, c1 ≠ c2 e testar valores de µ e α para um número maior de variáveis, por exemplo, modelo de Ogden para N = 3.

Agradecimentos

À Universidade Federal de Lavras (UFLA) e a CAPES pelo apoio financeiro durante a realização desse projeto.

Referências

• ¹
  Samuelson DA. Ophthalmic anatomy. In: Gelatt KN, editor. Veterinary ophthalmology. Iowa: Blackwell. 2007. vol.1. p. 37-148.
• ²
  Soares RC. Otimização de seções transversais de concreto armado sujeitas à flexão: aplicação a pavimentos [dissertação]. São Carlos: Escola de Engenharia de São Carlos, Universidade de São Paulo; 1997.140p.
• ³
  Kennedy J, Eberhart RC. Particle swarm optimization. In: IEEE International Conference on Neural Networks, 1995, Perth. Proceedings. Perth: IEEE; 1995. vol. 4. p. 1942-8.
• ⁴
  Parsopoulos KE, Vrahatis MN. Particle swarm optimization method in multiobjective problems. In: ACM Symposium on Applied Computing, 2002, Madrid. Proceedings. Madrid: ACM; 2002. p. 603-7.
• ⁵
  Clough RW. The finite finite element method in plane stress analysis. In: Conference on Electronic Computation, 2nd, 1960, Pittsburg. Proceedings. Pittsburg: ASCE; 1960.
• ⁶
  Ogden R. Non-linear elastic deformations. New York: Dover; 1984. 560 p.
• ⁷
  Lapper R, Gasson P, Karri V. A hyperelastic finite-element model of human skin for interactive realtime surgical simulation. IEEE Trans Biomed Eng. 2011;58(4):1013-22.
• ⁸
  Elsheikh A. Understanding corneal biomechanics through experimental assessment and numerical simulation. New York: Nova Science; 2010. 63 p.
• ⁹
  Magalhaes RR, Elsheikh A, Buechler P, Whitford C, Wang J. Application of particle swarm optimization in inverse finite element modeling to determine the mechanical behavior of the cornea. Acta Scientiarum Technology. 2017:39(3):325-331.
• ¹⁰
  Almeida PR. Estudo do comportamento mecânico de córneas via simulações numéricas [dissertação]. Lavras:Universidade Federal de Lavras, Lavras; 2016. 80p

Datas de Publicação

• Publicação nesta coleção
  Nov-Dec 2017

Histórico

• Recebido
  13 Dez 2016
• Aceito
  15 Set 2017