Exploring Relaxation Techniques in Physics Applications
Written on
Chapter 1: Introduction to Relaxation Methods
This article presents an insightful overview of numerical analysis techniques known as relaxation methods, particularly in the realms of electrostatics and fluid dynamics.
Let’s denote ϕ as a function of x and y, while A, B, ..., G are also functions of these variables. Any second-order partial differential equation (PDE) with ϕ as its solution can be articulated in a specific format. The subscripts indicate partial derivatives, such as ϕₓ=∂ϕ/∂x. The presence of 2B instead of just B arises from the equality ∂²ϕ/∂x∂y=∂²ϕ/∂y∂x. Our focus will be on PDEs where B²-AC is strictly negative, classifying them as elliptic equations.
Among the simplest elliptic equations is Poisson's equation expressed as ∇²ϕ=f(x,y). The specific case where ∇²ϕ=0 is recognized as Laplace's equation. The operator ∇² is referred to as the Laplacian:
Given that B=0 and A=C=1, it is evident that both Laplace's and Poisson's equations are elliptic. Poisson's equation plays a pivotal role in electrostatics, emerging from Gauss's law: since ∇⋅E=ρ/ε₀ and E=-∇ϕ, we derive ∇⋅E=∇⋅(-∇ϕ), leading to ∇²ϕ=-ρ/ε₀, where ρ denotes the charge density function. Consequently, Laplace's equation is applicable in charge-free regions, such as the space between capacitor plates.
Laplace's equation is also relevant in fluid dynamics. If a fluid’s velocity field V is irrotational (i.e., ∇×V=0), it implies that V is conservative. Therefore, a function ϕ exists such that V=∇ϕ. Taking the divergence of both sides yields ∇⋅V=∇²ϕ. According to the continuity equation -∂ρ/∂t=∇²ϕ, where ρ represents the fluid density, and if the fluid is incompressible, then ρ remains constant, leading to ∇²ϕ=0. This potential is measured in m²/s.
In this discussion, we will solve Laplace's equation in regions confined by surfaces where ϕ is predetermined. Such problems are termed Dirichlet problems for Laplace's equation.
Analyzing Dirichlet problems analytically can be tedious, especially for straightforward boundary conditions, and for many practical scenarios, pursuing an analytical solution may not be worthwhile. Thus, numerous numerical approximation methods exist for addressing Dirichlet problems, with relaxation techniques forming a crucial subset of these methods. Essentially, a relaxation method approximates the solution at each point by treating it as a simple function of neighboring points.
To devise a relaxation algorithm for Laplace's equation, we can impose a square grid over the region of interest, with points spaced by a distance h. We then approximate ϕ at (x,y) by averaging the values at (x+h,y), (x,y+h), (x-h,y), and (x,y-h). The accuracy of this approach can be confirmed using Taylor's theorem. Start by expanding ϕ for the nearest neighbors on the grid:
The notation 𝓞(h⁴) denotes terms where h appears to the fourth power or higher. By summing these four equations and rearranging, we arrive at:
The term associated with h² corresponds to the Laplacian, transforming the equation into h²f as ∇²ϕ=f. Furthermore, as h approaches 0, the higher-order terms diminish significantly compared to the h² term. Thus, when h is considerably small in relation to the problem's length scale, the following approximation holds:
Having covered the mathematical aspect, let's discuss the relaxation algorithm itself. In this article, we will only consider boundary lines parallel to the x and y axes, as illustrated in the example below, where black dots represent grid points spaced 1 centimeter apart.
Example 1: Basic Grid Setup
More intricate geometries that involve smooth curves or non-perpendicular angles are beyond the scope of this introductory piece. For simplicity, we will also assume f=0.
For demonstration purposes, I will utilize MATLAB, which is particularly suited for these applications. If you lack MATLAB but wish to experiment with the code, Scilab serves as a free and open-source alternative compatible with basic MATLAB syntax. It's important to remember that in MATLAB, (1,1) refers to the top left element of an array, while (i,j) indicates the element located at row i and column j.
The steps are as follows:
- Initialize an (M+1)×(N+1) rectangular array named phi(i,j), where M and N correspond to the width and height of the modeled area divided by the grid spacing. For instance, in the illustration above, the width and height are 4 centimeters, and the grid spacing is 1 centimeter, which results in M=N=4.
- Assign boundary values to the appropriate points.
- For all points (i,j) not located on the boundary, apply the formula: PHI(i,j)=0.25×[PHI(i+1,j)+PHI(i-1,j)+PHI(i,j+1)+PHI(i,j-1)].
- Conclude after reaching the desired number of iterations or a specific level of accuracy.
If you wish, both MATLAB and Scilab feature built-in routines for plotting the results.
The next two images illustrate the outcomes of applying relaxation to Example 1, with the grid spacing adjusted to a more practical 0.125 centimeters. To prevent clutter in the article, I have placed the code on Pastebin, which you may find helpful if you choose to experiment with it. The code is designed to work in Scilab.
These images also highlight the screening effect: A space enclosed by a grounded conductor is entirely shielded from external potentials. Nonetheless, complete enclosure is not a prerequisite for effective interference protection, which embodies the principle of a Faraday cage.
Example 2: Faraday Cage Simulation
Consider a one-meter-long plate maintained at 10kV, with a Faraday cage consisting of centimeter-thick bars positioned 50 centimeters away from the plate, spaced 10 centimeters apart. The code I created for this example can be found here.
Let's first visualize the scenario without the Faraday cage.
As anticipated, the potential diminishes approximately linearly with increasing distance from the plate, converging to an exact linearity if the plate were infinitely long in the y-direction. Now, let’s introduce the cage.
You can observe that the presence of the cage significantly alters the outcome, despite it not fully separating the two regions.
Now, let's explore an example from fluid dynamics.
Example 3: Fluid Flow in a Pipe
In this scenario, the ends of a 1-meter-wide L-shaped pipe, measuring 3 meters along its sides, are held at velocities of 200m²/s and -200m²/s, with the walls set at 0m²/s. We aim to determine the flow within the pipe. Since V=∇ϕ, the velocity aligns with the increasing ϕ, designating the bottom left boundary at -200m²/s as the pipe’s inlet and the top left boundary at 200m²/s as the outlet.
The relaxation code for this scenario can be found here. Below is the graph of the results:
An intriguing aspect of this potential is its flatness near the corner, which persists even with significantly high potential values at the inlet and outlet. This phenomenon relates to the skew-symmetry of the boundary conditions along the line y=300-x. The graph indicates that a test particle entering the pipe from the middle of the bottom left corner (e.g., at x=0, y=50) initially accelerates, slows down upon reaching the corner, and then speeds up again toward the outlet at the top right.
A streamline plot illustrates some of the trajectories, although it doesn’t capture the speed of the particle at each trajectory point.
An additional fascinating behavior is observed: in the horizontal section of the pipe, test particles are directed toward the wall, while in the vertical section, they are pushed away. Although the streamline plot does not demonstrate this, particles do not cease movement upon reaching the wall; instead, they continue along a path that closely parallels the wall. Similarly, particles do not originate from the vertical walls; the streamlines directed away from the walls represent particles that were flowing along the walls and are subsequently pushed away.
This solution represents the least accurate of the three presented. Specifically, the configuration of the streamlines is more sensitive to grid spacing than the mesh plot of the potential. Unfortunately, I am currently limited to my 2011 Macbook Air until I can acquire a new computer, which restricts my ability to reduce the spacing as much as I would prefer. However, if you have access to a more powerful computer, I highly recommend adjusting the MATLAB code settings, particularly experimenting with a grid spacing of 0.1 centimeters instead of the default 5 centimeters.
Concluding Remarks
While it would be ideal to analyze physical systems with analytical methods to achieve perfectly accurate solutions, the reality is that this is not always feasible. This underscores the significance of numerical techniques for physicists and engineers. Any aspiring professional in these fields is expected to possess at least a foundational understanding of numerical methods.
For those interested in further exploration, the article that inspired Example 3 may provide additional insights.
Copyright and Disclaimers
All images and code referenced in this article are my original creations, and I welcome their sharing or usage, provided proper attribution is given. I encourage downloading, running, and modifying my code. Although MATLAB and Scilab are generally safe and compatible with most modern systems, I cannot guarantee that my code will function flawlessly for everyone.
This essay is part of a series of articles on mathematics-related topics published in Cantor's Paradise, a weekly Medium publication. Thank you for reading!
The first video, titled "Reduce Stress through Progressive Muscle Relaxation (3 of 3) - YouTube," discusses techniques to alleviate stress through progressive muscle relaxation.
The second video, "Box breathing relaxation technique: how to calm feelings of stress or anxiety - YouTube," explores the box breathing method to help reduce stress and anxiety.