![]() |
deal.II version GIT relicensing-2879-g253433159b 2025-03-21 00:10:00+00:00
|
This program was contributed by Hongfeng Ma <hongfeng.mark.ma@gmail.com>.
It comes without any warranty or support by its authors or the authors of deal.II.
This program is part of the deal.II code gallery and consists of the following files (click to inspect):
![]() | ![]() | ![]() | ![]() |
![]() |
In this tutorial, the studied problem is to simulate temperature distributions of a sample under a moving laser. Light penetrates the substrate without loss. The top-covered thin-film is, however, a light absorbing material. For simplicity, the thin-film is assumed to be TiO \(_2\) mixed with silver nanoparticles, which supports light heating of the material by absorbing light energy. For this tutorial, we only consider the isotropic absorption. The Figure below illustrates the sketch of the problem. The absorption coefficient is assumed to be \(10^4 m^{-1}\). The substrate is glass. The thickness of the thin-film is assumed to be 400 nm. The spot size at the top of thin film is \(20 \mu m\) at \(e^{-2}\). The writing speed is assumed to be 10 mm/s. The laser power is 0.4 W. The time step is set as 10 \(\mu s\). The initial position of laser center is \(-50 \mu m\) with 50 \(\mu m\) distance away from the left boundary to avoid boundary effects.
In general, the non-uniform isotropic heat equation is as following
\begin{align*} \rho C_m \frac{\partial T}{\partial t} -\nabla \cdot (k\nabla T) = f(\textbf{x}) \end{align*}
Now, we discretize the equation in time with the theta-scheme as
\begin{align*} \rho C_m \frac{T^{n} - T^{n-1}}{dt} - [(1-\theta) \nabla \cdot (k \nabla T^{n-1}) + \theta \nabla \cdot (k \nabla T^n)] = (1-\theta) f^{n-1}(\textbf{x}) + \theta f^n(\textbf{x}) \end{align*}
where \(\theta\) is a parameter; if \(\theta = 0 (\text{or} = 1)\), it becomes forward (backward) explicit Euler method; the Crank-Nicolson method is when \(\theta = 0.5\). Integrating by the test function \(T^{'}\) and do integration by parts
\begin{align*} \int T^{'} \nabla \cdot (k \nabla T) dV = T^{'} k \nabla T |^a_b - \int k \nabla T \cdot \nabla T^{'} dV \end{align*}
since the test function is time invariant (because the grid is not changed), we have \(T^{'n}\) = \(T^{'n-1}\).
\begin{align*} T &= \sum{_i} {u_i} T^{'}_{i} \\ \int T T^{'}_{i} dV &= u_{i} \end{align*}
Let
\begin{align*} M &= \int \rho C_m T^{'_i} T^{'_j} dV \\ A & = \int k \nabla T^{'_i} \cdot \nabla T^{'_j} dV \\ F^n & = \int T' f^{n}(\textbf{x}) \end{align*}
we have the following term,
\begin{align*} \int T^{'} \rho C_m [T^{n} - T^{n-1}] - dt \int T^{'} [(1-\theta) \nabla \cdot (k \nabla T^{n-1}) + \theta \nabla \cdot (k \nabla T^n)] \\ = dt \int T^{'} (1-\theta) f^{n-1}(\textbf{x}) + dt \int T^{'} \theta f^n(\textbf{x}) \end{align*}
\begin{align*} M U^n - M U^{n-1} - dt \int T^{'} [(1-\theta) \nabla \cdot (k \nabla T^{n-1}) + \theta \nabla \cdot (k \nabla T^n)] \\ = dt \int T^{'} (1-\theta) f^{n-1}(\textbf{x}) + dt \int T^{'} \theta f^n(\textbf{x}) \end{align*}
\begin{align*} M U^n - M U^{n-1} + dt [(1-\theta) A U^{n-1} + \theta A U^n] \\ = dt (1-\theta) F^{n-1} + dt \theta F^{n} \end{align*}
the final equation becomes
\begin{align*} (M + dt \theta A) U^n = [M - dt (1-\theta) A] U^{n-1} + dt (1-\theta) F^{n-1} + dt \theta F^{n} \end{align*}
The initial temperature can be interpolated over each vertex as follows,
\begin{align*} M_0 &= \int T^{'_i} T^{'_j} dV \\ T_0 &= \sum_i u_i T^{'i} = g_0(x) \\ M_0 U &= \int g_0(\textbf{x}) T^{'i} dV \end{align*}
which is robust for general use. In fact, Deal.II provides a function (VectorTools::interpolate) doing the same thing, which is, however, may not necessary work for parallel version.
To simplify the question, the heat equation is solved in two-dimensions (x-y) by assuming that the z-axis is homogeneous. Following is part of the running results in 4-threads:
Solving problem in 2 space dimensions. Number of active cells: 6567 Total number of cells: 6567 Number of degrees of freedom: 11185 9 CG iterations needed to obtain convergence. initial convergence value = nan final convergence value = 2.31623e-20 Time step 1 at t=1e-05 time_step = 1e-05 80 CG iterations needed to obtain convergence. initial convergence value = nan final convergence value = 1.66925e-13 +------------------------------------------+------------+------+ | Total wallclock time elapsed since start | 2.98s | | | | | | | | Section | no. calls | wall time| % | +------------------------------+-----------+------------+------+ | assemble_rhs_T() | 1 | 0.107s | 3.6% | | assemble_system_matrix_init() | 1 | 0.245s | 8.2% | | make_grid() | 1 | 1.11s | 37% | | refine_mesh_at_beginning | 1 | 0.652s | 22% | | setup_system() | 1 | 0.276s | 9.3% | | solve_T() | 2 | 0.426s | 14% | +-------------------------------+-----------+-----------+------+ Time step 2 at t=2e-05 time_step = 1e-05 79 CG iterations needed to obtain convergence. initial convergence value = nan final convergence value = 2.06942e-13 +------------------------------------------+------------+------+ | Total wallclock time elapsed since start | 0.293s | | | | | | | Section | no. calls | wall time | % | +---------------------------------+--------+------------+------+ | assemble_rhs_T() | 1 | 0.0969s | 33% | | solve_T() | 1 | 0.161s | 55% | +------------------------------+-----------+------------+------+ Time step 3 at t=3e-05 time_step = 1e-05 80 CG iterations needed to obtain convergence. initial convergence value = nan final convergence value = 1.71207e-13
The colors stand for different cores.
All the include files have already been discussed in previous tutorials.
grid refine
remember to use the namespace dealii before defining a class that is inheritaged from Function<dim>
In general, it is more clear to separate boundary and initial condictions, as well as the right hand side function from a file holding all the things. To do so, in this work, the globalPara.h file defines physical constants, laser parameters, and heat characteristics of materials involved. Boundary and initial conditions are defined in boundaryInit.h. The rightHandSide.h defines the heat source, which in this work is a moving Gaussian beam.
Now the main class is defined as following
system_matrix
for storing left matrix
for storing right matrix
System_rhs, only locally owned cells
Solutions Old Solutions with ghost cells, for output
Old Solutions only with locally owned cells
New Solutions only with locally owned cells
Dynamic assembling of the righthandside terms
the constructor
the destructor
make grid by importing msh file, and rescale
initialization
we want to output solution, so here should have ghost cells
locally owned cells
copy to system_matrix_T and system_rhs_T for projecting initial values
store M + dt*theta*A as the left_system_matrix
store M - dt*(1-theta)*A as the right_system_matrix
The right hand side is assembled each time during running, which is necessary as the laser source is moving. To separate the heat source and the right hand side assembling, the right hand side function is defined as RightHandside<dim>.
Solving the equation is direct. Recall that we have defined several matrices and vectors, to avoid ambiguous, here, we only use system_matrix_T as the system matrix, system_rhs_T as the system right hand side. The vector completely_distributed_solution is used to store the obtained solution.
Print the number of iterations by hand.
only refine mesh 5um above the TiO2 and glass interface
the length of output numbering
output the overall solution
This is the function which has the top-level control over everything. Apart from one line of additional output, it is the same as for the previous example.
projection of initial conditions by solving. solution stored in new_solution_T;
reinitialization
output initial values; need ghost cells
the dynamic solving part
old_solution_T is used for output, holding ghost cells old_solution_T_cal is used for calculation, holding only locally owned cells.
main
function#-------------------------------------------------------— #
# #-------------------------------------------------------—
#-------------------------------------------------------—
#-------------------------------------------------------—
#-------------------------------------------------------—
#-------------------------------------------------------—
thin film
substrate
#-------------------------------------------------------— #
# #-------------------------------------------------------—
#-------------------------------------------------------—
#-------------------------------------------------------—