iTranslated by AI

The content below is an AI-generated translation. This is an experimental feature, and may contain errors. View original article
🐥

[Mesh-free] Solving Complex Geometries (Holes & Airfoils) with PINNs | A Practical Guide for CAE Engineers Vol. 5

に公開

Introduction

To all CAE engineers, what percentage of your life have you spent on "mesh generation"?

Tiny fillets, complex channels, stress concentration around holes... The process of re-meshing every time the geometry changes slightly is the primary bottleneck in analysis workflows.

In Chapter 5, we will implement the true essence of PINNs: "meshless geometry definition." Specifically, we will solve the physical field around a "plate with a circular hole." This shape serves as a foundation for fluid analysis (flow around a cylinder) and structural analysis (stress concentration).


1. Objective (Physics)

We will simulate the flow (potential flow) in a "channel with a circular obstacle in the center."

In traditional CAE, it is necessary to finely mesh the area around the circle, but with PINNs, this is completed simply by performing Python operations: "discard points (data) inside the circle and place points (boundary conditions) on the surface of the circle."

  • Phenomenon: 2D incompressible, irrotational flow (potential flow)
  • Domain: A rectangular region with a circular hole in the center

2. Theory

Geometry by Boolean Operation

"Geometry definition" in PINNs is essentially a filtering operation that determines whether a coordinate point (x, y) is included in the domain.

  • Overall Domain: x \in [-2, 2], \quad y \in [-1, 1]
  • Hole (Obstacle): A circle centered at the origin with radius R=0.5
  • Computational Domain: x^2 + y^2 > R^2 (outside the circle)

Governing Equation (Laplace Equation)

For an ideal fluid where viscosity is ignored, the stream function \psi follows the Laplace equation.

\frac{\partial^2 \psi}{\partial x^2} + \frac{\partial^2 \psi}{\partial y^2} = 0

Flow velocity can be obtained by u = \frac{\partial \psi}{\partial y}, \quad v = -\frac{\partial \psi}{\partial x}.

Boundary Conditions

  • Inlet/Outlet (Left/Right): Uniform flow u=1, so \psi = y
  • Top/Bottom Walls: Slip walls \psi = y_{wall}
  • Cylinder Surface: Since the flow does not penetrate the object, the surface matches the streamline: \psi = 0

3. Implementation (Code)

The highlight of this code is the CSG_Geometry class. We use a technique called "Rejection Sampling" to remove points that fall inside the hole from the randomly scattered points using a mask.

import torch
import torch.nn as nn
import numpy as np
import matplotlib.pyplot as plt

# ==========================================
# 0. Config
# ==========================================
torch.manual_seed(1234)
np.random.seed(1234)
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")

# Geometry parameters
R = 0.5  # Radius of the cylinder

# ==========================================
# 1. Geometry Class (The "Meshless" Magic)
# ==========================================
class CSG_Geometry:
    def __init__(self, x_range, y_range, radius):
        self.x_min, self.x_max = x_range
        self.y_min, self.y_max = y_range
        self.radius = radius

    def sample_interior(self, n_samples):
        """Sample points inside the domain (Rejection Sampling)"""
        # 1. First, scatter points generously across the entire rectangular area
        x = (self.x_max - self.x_min) * torch.rand(int(n_samples * 1.5), 1) + self.x_min
        y = (self.y_max - self.y_min) * torch.rand(int(n_samples * 1.5), 1) + self.y_min
        
        # 2. Filtering (keep only points outside the hole)
        mask = (x**2 + y**2) > self.radius**2
        x_valid = x[mask]
        y_valid = y[mask]
        
        return x_valid[:n_samples].to(device), y_valid[:n_samples].to(device)

    def sample_hole_boundary(self, n_samples):
        """Generate points on the surface of the hole (on the circumference)"""
        theta = 2 * np.pi * torch.rand(n_samples, 1)
        x = self.radius * torch.cos(theta)
        y = self.radius * torch.sin(theta)
        return x.to(device), y.to(device)

# ==========================================
# 2. Model & Loss Function
# ==========================================
class PINN(nn.Module):
    def __init__(self):
        super().__init__()
        self.net = nn.Sequential(
            nn.Linear(2, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 1)
        )
    def forward(self, x, y):
        return self.net(torch.cat([x, y], dim=1))

def gradients(outputs, inputs):
    return torch.autograd.grad(outputs, inputs, torch.ones_like(outputs), create_graph=True)[0]

# ==========================================
# 3. Training Loop
# ==========================================
geom = CSG_Geometry(x_range=(-2, 2), y_range=(-1, 1), radius=R)
model = PINN().to(device)
optimizer = torch.optim.Adam(model.parameters(), lr=0.001)

print("Start Training for Complex Geometry...")
for epoch in range(5001):
    optimizer.zero_grad()
    
    # --- Loss 1: Governing Equation (Laplace Eq) ---
    x_int, y_int = geom.sample_interior(2000)
    x_int.requires_grad = True
    y_int.requires_grad = True
    psi = model(x_int, y_int)
    
    psi_x = gradients(psi, x_int)
    psi_y = gradients(psi, y_int)
    psi_xx = gradients(psi_x, x_int)
    psi_yy = gradients(psi_y, y_int)
    
    loss_pde = torch.mean((psi_xx + psi_yy)**2)
    
    # --- Loss 2: Cylinder Surface (psi = 0) ---
    x_cyl, y_cyl = geom.sample_hole_boundary(200)
    psi_cyl = model(x_cyl, y_cyl)
    loss_cyl = torch.mean(psi_cyl**2)
    
    # --- Loss 3: Outer Boundaries (Inlet/Outlet/Walls) ---
    x_box = (torch.rand(500, 1, device=device) * 4 - 2)
    y_box = (torch.rand(500, 1, device=device) * 2 - 1)
    mask_box = (torch.abs(x_box) > 1.9) | (torch.abs(y_box) > 0.9)
    x_bc, y_bc = x_box[mask_box], y_box[mask_box]
    
    psi_bc_pred = model(x_bc, y_bc)
    psi_bc_true = y_bc # psi = U * y (U=1)
    loss_bc = torch.mean((psi_bc_pred - psi_bc_true)**2)
    
    loss = loss_pde + loss_cyl * 10.0 + loss_bc
    loss.backward()
    optimizer.step()
    
    if epoch % 1000 == 0:
        print(f"Epoch {epoch}: Loss {loss.item():.5f}")

# ==========================================
# 4. Visualization
# ==========================================
x_test = np.linspace(-2, 2, 200)
y_test = np.linspace(-1, 1, 100)
X, Y = np.meshgrid(x_test, y_test)
mask_hole = (X**2 + Y**2) < R**2

x_flat = torch.tensor(X.flatten(), dtype=torch.float32).view(-1, 1).to(device)
y_flat = torch.tensor(Y.flatten(), dtype=torch.float32).view(-1, 1).to(device)

model.eval()
with torch.no_grad():
    psi_pred = model(x_flat, y_flat).cpu().numpy().reshape(100, 200)

psi_pred[mask_hole] = np.nan

plt.figure(figsize=(10, 5))
plt.contourf(X, Y, psi_pred, levels=50, cmap='viridis')
plt.colorbar(label='Stream Function $\psi$')
plt.contour(X, Y, psi_pred, levels=20, colors='black', linewidths=0.5)
circle = plt.Circle((0, 0), R, color='white')
plt.gca().add_patch(circle)
plt.title("Flow around a Cylinder (Mesh-free PINNs)")
plt.xlabel("x"); plt.ylabel("y")
plt.axis('equal')
plt.show()

4. Validation

Check the output figure (contours and streamlines).

Streamlines: You should see the flow entering from the left, smoothly avoiding the central circle (white blank area), splitting above and below, and merging again on the right side.

Boundary Alignment: You can confirm that the streamlines follow the circle near the surface (the \psi=0 contour line matches the circumference).

Meshless: There is absolutely no "mesh generation" or "node connection" processing in the code. We simply wrote the conditional expression:
x2 + y2 > R**2


5. Discussion

Application to Airfoils: This time it was a simple circle, but the process is the same even for airfoils (NACA profiles, etc.).
Prepare a function y = f(x) that represents the surface of the wing.
Use the sample_interior function to discard random points that fall inside the wing.
Use the sample_boundary function to generate point coordinates on the wing surface and include the boundary conditions at those points in the Loss.

If you were to perform airfoil analysis using traditional CFD,
you would need specialized mesh generation technologies such as "C-grids" or "O-grids."

With PINNs, shape definition is completed simply by "coordinate filtering."
This is the greatest福音 (blessing) for CAE engineers.

(This article is Chapter 5 of the book currently being written, "PINNs Implementation Bible for CAE Engineers")

GitHubで編集を提案

Discussion