US difraccion limitada

En el VPS se activa FEniCS usando:

conda activate fenics

Se va a la carpeta en la que se quiere trabajar el proyecto y se crea un archivo llamado "Cylinder.geo"

SetFactory("OpenCASCADE");

// Parámetros del cilindro
r = 0.5; // Radio en cm
h = 0.1; // Altura en cm

// Crear un disco y extruirlo para formar el cilindro
Disk(1) = {0, 0, 0, r, r};
Extrude {0, 0, h} { Surface{1}; }

// Asignar tamaño de malla
MeshSize {1} = 0.05;

// Definir volumen físico
Physical Volume("cylinder") = {1};

Se compila por medio de:

gmsh -3 cylinder.geo -format msh2 -o cylinder.msh

Se obtiene el archivo cylinder.xdmf

meshio convert cylinder.msh cylinder.xdmf

Se crea el archivo:

nano piezo_prueba1.py

from dolfin import *
import numpy as np
import meshio

# Definir la geometría: un cilindro (disco piezoeléctrico)
R = 5e-3   # Radio (5 mm)
H = 1e-3   # Espesor (1 mm)
n_mesh = 40  # Resolución de la malla

#mesh = Mesh("cylinder.xml")  # Se reemplaza con malla cilíndrica generada externamente
mesh = Mesh()
with XDMFFile("cylinder.xdmf") as xdmf:
    xdmf.read(mesh)

# Definir el espacio de funciones
V = VectorFunctionSpace(mesh, "P", 1)
vu = TestFunction(V)

# Definir los parámetros del material (PZT-5H)
E = 63e9    # Módulo de Young (Pa)
nu = 0.31   # Coeficiente de Poisson
rho = 7500  # Densidad (kg/m³)

# Tensor de elasticidad de isotrópico
mu = E / (2 * (1 + nu))
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))

def sigma(u):
    return 2 * mu * sym(grad(u)) + lmbda * tr(sym(grad(u))) * Identity(3)

# Definir funciones de prueba y de prueba
u = TestFunction(V)
du = TrialFunction(V)

# Definir la carga mecánica en la parte superior
f = Constant((0, 0, -1e6))  # Fuerza en la dirección -Z (Pa)

# Condiciones de frontera: fijar la base del disco (z = 0)
def boundary_bottom(x, on_boundary):
    return near(x[2], 0)

bc = DirichletBC(V, Constant((0, 0, 0)), boundary_bottom)

# Definir la ecuación de equilibrio
F = inner(sigma(du), grad(vu)) * dx - inner(f, vu) * dx

# Resolver el problema
u_sol = Function(V)
solve(lhs(F) == rhs(F), u_sol, bc)

# Guardar los resultados en formato VTK para ParaView
vtkfile = File("deformation.pvd")
vtkfile << u_sol

print("Simulación completada. Ver resultado en ParaView.")

Se ejecuta como:

python piezo_prueba1.py

Finalmente se abre con paraview

Los archivos se pueden transferir desde el VPS al cliente utilizando el siguiente comando:

scp root@195.35.37.12:~/FENICS/solution* .