En mecánica de fluidos, el operador bi-Laplaciano con condiciones de contorno de Neumann homogéneas surge al transformar las ecuaciones de Navier-Stokes a la formulación de vorticidad-velocidad. En el caso de problemas con una dirección periódica, el problema puede transformarse en múltiples problemas elípticos independientes de cuarto orden bidimensionales. Se diseñó y validó un método eficiente para resolver estos operadores bi-Laplacianos bidimensionales con condiciones de contorno de Neumann homogéneas utilizando esquemas de diferencias finitas compactas en 2D. La solución se formula como una combinación lineal de soluciones auxiliares, tantas como el número de puntos en el contorno, un método que era prohibitivo hace algunos años debido a los grandes requisitos de memoria para almacenar todas estas funciones auxiliares. La validación se ha realizado para diferentes configuraciones de campo, tamaños de rejilla y estarcidos del esquema numérico, mostrando su potencial para abordar campos de alto gradiente como los que se pueden encontrar en flujos turbulentos.