Non-destructive test method of wood moisture content based on multiple varying bounds integral numerical method

Assumptions of the model and notation description

Some notations and basic assumptions are given in this paper:

1) $B_{1}$ and $B_{2}$ denote two thin electrode plates of the capacitive sensor.

2) $B$ denotes the measured wood, and $B$, $B_{1}$, $B_{2}$ do not intersect each other.

3) $\varepsilon = \varepsilon(x, y, z)$ denotes the dielectric constant of the wood at the point $(x, y, z)$.

4) $\tilde{\varepsilon}$ denotes the dielectric constant of air.

5) $V = V(x, y, z)$ denotes the potential at the point $(x, y, z)$.

6) $V_{1}$ and $V_{2}$ denote the potentials of $B_{1}$ and $B_{2}$, respectively.

7) $E = E(x, y, z)$ denotes the electric field strength at the point $(x, y, z)$.

8) $C$ denotes the capacitance of the capacitive sensor.

9) $q_{1}$ and $q_{2}$ denote the charge carried by $B_{1}$ and $B_{2}$, respectively, with $q_{1} = q_{2}$. Additionally, $q_{1}$ also satisfies

10) $\partial B_{in}$ and $\partial B_{out}$ denote the interior and exterior of the wood boundary, respectively.

1.   Introduction

Wood moisture content (MC) detection is one essential problem in the modern world. All wood products contain some level of moisture, and too high or too low MC can cause quality issues in wood products. Thus, inaccurate MC detection severely impacts quality control in the wood processing industry. It also impacts the maintenance of wooden historical buildings and the protection of rare tree species ^{[1,2,3,4,5,6,7]}. Once wood products or wooden structures are completed, their shape and material no longer change. The key factors determining their internal quality are primarily the MC and drying stress. When the equilibrium moisture content is reached, cracking and deformation are least likely to occur. In recent years, inaccurate MC detection has resulted in a wood utilization rate of only 50$\%$–60$\%$. In contrast, technologically advanced countries have achieved rates as high as 90$\%$. This has led to frequent significant economic losses in the national economy.

Due to the significant social and economic benefits of MC detection, scholars from various countries have conducted extensive and in-depth research on this topic. With continuous research, methods for MC detection ^[8] have also been continuously developed and improved. Martin et al.^[9] used the resistance method to study the relationship between MC and average resistivity of silver fir. Van Blokland and Adamopoulos ^[10] analyzed the electrical resistance characteristics of Norway spruce. Through simple linear regression, they derived a first-order polynomial function relationship between MC and resistance of the wood. However, the resistance method can obtain only ideal measurement results when the MC is between 6$\%$ and 25$\%$. The measurement error is severely limited by the value of MC. Edwards and Jarvis ^[11], Batranin et al. ^[12], and Penttila et al. ^[13] used the ray technique to detect changes in MC, which are not limited by specific MC values. However, the ray technique has high detection costs, slow speeds, and safety hazards, posing challenges for inspectors. Qin et al. ^[14] studied the relationship between MC and $\varepsilon$ but noted that different tree species require different dielectric models. Compared to these methods, we aim to find a more widely applicable method that can meet more measurement requirements. Capacitive sensors ^[15] are electronic components that can convert the measured physical quantities and their variation laws into capacitance and capacitance change laws. They have the outstanding advantage of performing non-contact, continuous, and high-precision measurements. Therefore, we consider measuring the capacitance between the electrode plates to calculate MC ^{[16,17,18,19]}. On the other hand, we found that most researchers use only experimental methods to MC detection. We propose combining experiments with mathematical methods to build a mathematical model and complete the MC detection.

For the mathematical model, it is also very important to choose the appropriate numerical method, which will have a great influence on numerical accuracy and speed. To some extent, the finite volume method ^[20,21,22] has absorbed some advantages of the finite element method and the finite difference method. The finite volume method is mainly applied to solve fluid flow and heat transfer problems. With this method, the numerical scheme obtained can keep some properties of original differential equations, such as the conservation of mass, momentum, and energy.

The multiple varying bounds integral method ^{[23,24,25,26,27]} is a new numerical method developed based on the finite volume method. First, by means of multiple integrals, all the derivatives in the space direction of the differential equation can be eliminated. In this way, the differential equation can be equitably represented by another new equation that contains only the unknown function than its derivatives. On this basis, we begin to construct the corresponding numerical scheme for this new equation. This avoids the occurrence of large errors greatly, especially where the change rate of the derivative is large. Furthermore, multiple integrals with varying bounds are important. Different Integral bounds can help us get different numerical schemes. Depending on the physical properties of differential equation, we can choose an appropriate scheme from them.

The layout of the paper is as follows. In Section 2, we provide the methodology for MC detection. In Section 3, we discuss the steps for MC detection. In Section 4, by multiple varying bounds integral method, we construct second-order compact discrete schemes in fields where the potential is continuous and discontinuous, respectively, and conduct error estimation and numerical experiments. In Section 5, a conclusion is provided.

2.   Methodology for wood moisture content detection

A core issue for wood moisture content detection is to construct the function relationship between the capacitance C and the dielectric constant $\varepsilon$. The methodology is shown in Figure 1. This is a multi-physics inverse problem with partial differential equations. We choose many suitable values of $\varepsilon$. On this basis, through the mathematics model we constructed, we obtain the values of $C$. The data $C$ and $\varepsilon$ are regressed to fit a function $\varepsilon = \varepsilon(C)$. Then, by capacitive sensor, we measure the value of $C$, and substitute it into the aforementioned fitted function to calculate $\varepsilon$. Ultimately, by the one-to-one correspondence between $\varepsilon$ and MC, we obtain the value of MC.

3.   Steps for wood MC detection

We focus on constructing the function relationship between $\varepsilon$ and $C$. The steps are as follows:

1) The measured wood $B$ and the electrode plates $B_{1}$ and $B_{2}$ are shown in Figure 2.

During practical wood detection, MC is always indirectly measured. This is accomplished by finding a physical quantity correlated with MC, measuring this quantity, and establishing its relationship with MC to determine MC. In this paper, we use the dielectric constant as the physical quantity related to MC.

As shown in Figure 2, two metal plates of identical material and specifications are used as the electrode plates of the capacitive sensor and are fixed at specific positions, labeled $B_{1}$ and $B_{2}$. According to physical law, the capacitive sensor's internal dielectric material affects the amount of charge it can store. Therefore, the capacitance between the two plates is influenced by the surrounding environment. If the surrounding environment remains unchanged, and the measured wood is placed on top of the plates, an electromagnetic field is generated around it. If the dimension of measured wood is known, the value of $C$ will be uniquely determined by $\varepsilon$, establishing a one-to-one correspondence between $C$ and $\varepsilon$. Thus, by measuring $C$, we can obtain $\varepsilon$. Additionally, we have found that there is a one-to-one relationship between $\varepsilon$ and MC. Therefore, by $\varepsilon$ as an intermediary, we first establish the relationship between $C$ and $\varepsilon$, and then leverage the corresponding relationship between $\varepsilon$ and MC to determine MC.

2) According to the mathematical model in ^[28], satisfied by the potential $V(x, y, z)$,

If $V_{1}$ and $\varepsilon$ are given, then $V = V(x, y, z)$ and $V_{2}$ are uniquely determined.

3) If $V(x, y, z)$ is determined, the electric field strength $E(x, y, z)$ can be determined by

It can be seen that $E(x, y, z)$ is related to $\varepsilon(x, y, z)$.

4) By integral,

and

the charge $q_1$ (or $q_2$) of $B_1$ (or $B_{2}$) is obtained ^[29]. Obviously, $q_1$ (or $q_2$) is also related to $\varepsilon(x, y, z)$.

5) By the equation

the capacitance $C$ can be obtained, and it can be seen that $C$ is also related to $\varepsilon(x, y, z)$.

Thus, $V(x, y, z)$, $E(x, y, z)$, and $C$ are related to $\varepsilon(x, y, z)$. Therefore, when $V_{\mathrm{l}}$ and $\varepsilon$ are given, it means that the system is determined. In view of this, when we regard the values of $\varepsilon$ as integers from 1 to 20, a series of corresponding $C$ can be obtained. Through regression analysis, we can construct the function relationship between $\varepsilon$ and $C$.

4.   The multiple varying bounds integral method

To maintain the advantages of the finite volume method and meet the high-precision demands of detection problems, we develop a new method called the multiple varying bounds integral method to calculate the capacitive sensor model.

The two electrode plates occupy the domain

and

The measured wood occupies the domain

This is an unbounded problem. To achieve numerical computation, we set a reasonable artificial boundary $[x_{L}, x_{R}]\times[y_{L}, y_{R}]\times[z_{L}, z_{R}]$. The computational domain is discretized by taking $h_x = \frac{x_R-x_L}{N_x}, h_y = \frac{y_R-y_L}{N_y}$ and $h_z = \frac{z_R-z_L}{N_z}$ as the step size in the direction of the $x$-axis, $y$-axis, and $z$-axis, respectively. Here $N_x, N_y$, and $N_z$ are three given positive integers, so we get $\left(N_x+1\right) \times\left(N_y+1\right) \times\left(N_z+1\right)$ grid nodes. Let each grid node be $(x_{i}, y_{j}, z_{k})$, where $x_i = x_L+i h_x\left(i = 0, 1, \cdots, N_x\right), y_j = y_L+j h_y\left(j = 0, 1, \cdots, N_y\right)$ and $z_k = z_L+k h_z\left(k = 0, 1, \cdots, N_z\right)$. Set $V(x_{i}, y_{j}, z_{k}) = V_{i, j, k}$, which approximately represents the corresponding potential at each node.

4.1.   Establishment of a compact discrete scheme

In fact, the mathematical model to be solved is a multi-physics inverse problem. Therefore, for different fields, we should select appropriate numerical methods for different differential equations. We divide physical fields into two regions: Field 1, where $V_{i, j, k}$ is continuous, and field 2, where $V_{i, j, k}$ is discontinuous.

4.1.1.   Field 1

In the computational domain, the variable $V_{i, j, k}$ is continuous except at the faces, vertices and edges of the wood. Therefore, the multiple varying bounds integral method can be directly applied to discretize the equation in field 1.

When $(x, y, z)\in(R^{3}-B_{1}-B_{2}-\overline{B})$, the potential $V(x, y, z)$ satisfies

and is continuous, so we can directly apply the multiple varying bounds integral method in this field to integrate Eq (4.1). We perform multiple varying bound integrals with respect to $x$, $y$, and $z$ on both sides of Eq (4.1), so we have

By Eq (4.2)

Then, integrating both sides of Eq (4.3) within the control volume, we have

The three terms on the left-hand side of Eq (4.4) have identical structures. Therefore, we compute only the first term on the left-hand side. The other two terms can be handled similarly.

where, $\alpha_{x}, \alpha_{y}$ and $\alpha_{z}$ are undetermined parameters. In the case of following the wood physical properties, we set $\alpha_{x} = \alpha_{y} = \alpha_{z} = \alpha$ and $h_{x} = h_{y} = h_{z} = h$.

Through the integral, the derivative function of an unknown function $V(x, y, z)$ can be eliminated, thus reducing numerical errors, especially where the rate of change of the derivative is large. Next, we interpolate the original function of $V(x, y, z)$. Here, let the node step size in each direction be $h$. We select three interpolation nodes in each of the $x$, $y$, and $z$ directions, generating a total of 27 grid nodes $(x_{i}, y_{j}, z_{k})$, $(x_{i}, y_{j-1}, z_{k})$, $(x_{i}, y_{j+1}, z_{k})$, $(x_{i-1}, y_{j}, z_{k})$, $(x_{i-1}, y_{j-1}, z_{k})$, $(x_{i-1}, y_{j+1}, z_{k})$, $(x_{i+1}, y_{j}, z_{k})$, $(x_{i+1}, y_{j-1}, z_{k})$, $(x_{i+1}, y_{j+1}, z_{k})$, $(x_{i}, y_{j}, z_{k-1})$, $(x_{i}, y_{j-1}, z_{k-1})$, $(x_{i}, y_{j+1}, z_{k-1})$, $(x_{i-1}, y_{j}, z_{k-1})$, $(x_{i-1}, y_{j-1}, z_{k-1})$, $(x_{i-1}, y_{j+1}, z_{k-1})$, $(x_{i+1}, y_{j}, z_{k-1})$, $(x_{i+1}, y_{j-1}, z_{k-1})$, $(x_{i+1}, y_{j+1}, z_{k-1})$, $(x_{i}, y_{j}, z_{k-2})$, $(x_{i}, y_{j-1}, z_{k-2})$, $(x_{i}, y_{j+1}, z_{k-2})$, $(x_{i-1}, y_{j}, z_{k-2})$, $(x_{i-1}, y_{j-1}, z_{k-2})$, $(x_{i-1}, y_{j+1}, z_{k-2})$, $(x_{i+1}, y_{j}, z_{k-2})$, $(x_{i+1}, y_{j-1}, z_{k-2})$, $(x_{i+1}, y_{j+1}, z_{k-2})$, as shown in Figure 3.

Based on these 27 nodes, we construct the interpolating function for the unknown function $V(x, y, z)$, as follows

Substituting Eq (4.6) into Eq (4.5), we have

where ${\lambda_{1}} = \frac{\alpha^2(4h+3\alpha)^2}4,$ ${\lambda_{2}} = -\frac{3\alpha(-32h^3-24h^2\alpha+4h\alpha^2+3\alpha^3)}2,$ ${\lambda_{3}} = \frac{\alpha^2(-16h^2+9\alpha^2)}4, {\lambda_{4}} = 9(-8h^2+\alpha^2)^2,$ ${\lambda_{5}} = -\frac{3\alpha(32h^3-24h^2\alpha-4h\alpha^2+3\alpha^3)}2$ and ${\lambda_{6}} = \frac{(4h-3\alpha)^2\alpha^2}4$. To this end, we have completed the discretization of Eq (4.1) by the multiple varying bounds integral method.

4.1.2.   Field 2

We divide the wood boundary into three parts: Faces, vertices, and edges. We found that $\frac{\partial V}{\partial n}$ is discontinuous in these fields, so direct integral cannot be performed and special treatment is required.

1) Analysis about faces of wood.

The wood has a total of six faces. We use the upper face as an example. The other faces follow in a similar fashion. By the Taylor expansion exterior to the upper face, we have

and

From Eqs (4.8) and (4.9), it follows that

Similarly, within the interior of the upper face of the wood, there are

Since the upper face is horizontal, the $V_{i, j, k}'$ direction is vertical. Thus, $V_{i, j, k}' = \frac{\partial V}{\partial n}$. Consequently, $\varepsilon\:\frac{\partial V}{\partial n}\Big|_{\partial B_{in}} = \tilde{\varepsilon}\:\frac{\partial V}{\partial n}\Big|_{\partial B_{out}}$ can be rewritten as

Similarly, we can obtain the discrete equations for the other 5 faces of the wood.

2) Analysis about vertices of wood.

The wood has a total of eight vertices. We use the upper front right vertex of the wood as an example. By model (3.1), there is

This is assuming that $\Omega$ represents a spherical surface with its center at the upper front right vertex and a radius of $h$. $\Omega_{1}$ and $\Omega_{2}$ denote the parts of $\Omega$ interior and exterior the wood, respectively. There are $\Omega_{1} = \frac{1}{8}\Omega\:$ and $\Omega_{2} = \frac{7}{8}\Omega\:$. $\Omega_{2}$ can also be divided into three parts based on the vertices, edges, and faces into three parts, that is $\Omega_{2} = \Omega_{21}+\Omega_{22}+\Omega_{23}$. Moreover, $\Omega_{21} = \frac{1}{8}\Omega\:, \Omega_{22} = \frac{1}{4}\Omega\:$, and $\Omega_{23} = \frac{1}{2}\Omega\:$, as shown in Figure 4.

From $\frac{\partial V}{\partial n} = \frac{\partial V}{\partial x}\cos\alpha+\frac{\partial V}{\partial y}\cos\beta+\frac{\partial V}{\partial z}\cos\gamma$, we have Eq (4.13), which can be written as

Next, we discuss each of the items in Eq (4.14).

(ⅰ) Consider the interior of the sphere

Take $\iint_{\Omega_{1}}\frac{\partial V}{\partial x}dydz$ as an example. When calculating $\iint_{\Omega_{1}}\frac{\partial V}{\partial x}dydz$, we take $(x_{i}, y_{j}, z_{k})$ as the center and extend two steps in the negative $x$, $y$, and $z$ directions, with a step length of $2h$. We define the upper front right vertex as $(x_{0}, y_{0}, z_{0})$. At this point, 27 points $(x_{0}, y_{0}, z_{0})$, $(x_{0}, y_{0}, z_{-1})$, $(x_{0}, y_{0}, z_{-2})$, $(x_{0}, y_{-1}, z_{0})$, $(x_{0}, y_{-1}, z_{-1})$, $(x_{0}, y_{-1}, z_{-2})$, $(x_{0}, y_{-2}, z_{0})$, $(x_{0}, y_{-2}, z_{-1})$, $(x_{0}, y_{-2}, z_{-2})$, $(x_{-1}, y_{0}, z_{0})$, $(x_{-1}, y_{0}, z_{-1})$, $(x_{-1}, y_{0}, z_{-2})$, $(x_{-1}, y_{-1}, z_{0})$, $(x_{-1}, y_{-1}, z_{-1})$, $(x_{-1}, y_{-1}, z_{-2})$, $(x_{-1}, y_{-2}, z_{0})$, $(x_{-1}, y_{-2}, z_{-1})$, $(x_{-1}, y_{-2}, z_{-2})$, $(x_{-2}, y_{0}, z_{0})$, $(x_{-2}, y_{0}, z_{-1})$, $(x_{-2}, y_{0}, z_{-2})$, $(x_{-2}, y_{-1}, z_{0})$, $(x_{-2}, y_{-1}, z_{-1})$, $(x_{-2}, y_{-1}, z_{-2})$, $(x_{-2}, y_{-2}, z_{0})$, $(x_{-2}, y_{-2}, z_{-1})$, $(x_{-2}, y_{-2}, z_{-2})$ inside the wood are selected for interpolation, as shown in Figure 5.

The interpolation function is

Since $\frac{\partial V}{\partial x}|_{i, j, k}$ is unknown, it is represented by the difference quotient, and we have

In that case,

then there is

The integration domain of $A_{1}$ can be seen according to the $yoz$ plane of Figure 6, so there is

The same reasoning leads to

and

Further simplification and organization gives

where,

Similarly, when calculating $\iint\limits_{\Omega_{1}}\frac{\partial V}{\partial y}\:dxdz$ and $\iint\limits_{\Omega_{1}}\frac{\partial V}{\partial z}\:dxdy$, it is also necessary to construct an interpolation function. The interpolation nodes used for this purpose are illustrated in Figure 5. Thus, there are

and

among them,

and

Therefore,

We have completed the discretization at $\Omega_{1}$.

(ⅱ) Consider the exterior of the sphere

We use the method for calculating $\Omega_{1}$ to compute the integral on $\Omega_{2i}\; (i = 1, 2, 3)$. First, consider field $\Omega_{21}$. We need to interpolate around it by selecting 27 points $(x_{0}, y_{0}, z_{0})$, $(x_{0}, y_{0}, z_{1})$, $(x_{0}, y_{0}, z_{2})$, $(x_{0}, y_{-1}, z_{0})$, $(x_{0}, y_{-1}, z_{1})$, $(x_{0}, y_{-1}, z_{2})$, $(x_{0}, y_{-2}, z_{0})$, $(x_{0}, y_{-2}, z_{1})$, $(x_{0}, y_{-2}, z_{2})$, $(x_{-1}, y_{0}, z_{0})$, $(x_{-1}, y_{0}, z_{1})$, $(x_{-1}, y_{0}, z_{2})$, $(x_{-1}, y_{-1}, z_{0})$, $(x_{-1}, y_{-1}, z_{1})$, $(x_{-1}, y_{-1}, z_{2})$, $(x_{-1}, y_{-2}, z_{0})$, $(x_{-1}, y_{-2}, z_{1})$, $(x_{-1}, y_{-2}, z_{2})$, $(x_{-2}, y_{0}, z_{0})$, $(x_{-2}, y_{0}, z_{1})$, $(x_{-2}, y_{0}, z_{2})$, $(x_{-2}, y_{-1}, z_{0})$, $(x_{-2}, y_{-1}, z_{1})$, $(x_{-2}, y_{-1}, z_{2})$, $(x_{-2}, y_{-2}, z_{0})$, $(x_{-2}, y_{-2}, z_{1})$, $(x_{-2}, y_{-2}, z_{2})$ on the exterior of the wood, as shown in Figure 7.

At this point, the interpolating polynomial constructed are

and

Thus, we can get

where,

and

Next, consider field $\Omega_{22}$, as shown in Figure 8.

Thus, we obtain

where,

and

Finally, consider region $\Omega_{23}$, as shown in Figure 9.

We can get

among it,

and

Consequently, Eqs (4.28), (4.32), (4.36), and (4.41) constitute the discrete scheme of the upper front right vertex Eq (4.13). After rearrangement, we have

Similarly, discrete schemes for the other 7 vertices can be obtained.

3) Analysis about edges of wood.

The wood has twelve edges. We use the upper right edge of the wood as an example. By model (3.1), there is

where $\Omega$ denotes a cylindrical surface with the upper right edge as the axis and $h$ as the radius. $\Omega_{\mathrm{l}}$ and $\Omega_{2}$ denote the parts of $\Omega$ in the interior and exterior of the wood, respectively. There are $\Omega_1 = \frac{1}{4}\Omega$ and $\Omega_2 = \frac{3}{4}\Omega$. $\Omega_{2}$ is also divided into two parts according to the edges and faces, that is $\Omega_{2} = \Omega_{21} + \Omega_{22}$, where $\Omega_{21} = \frac{1}{4}\Omega$ and $\Omega_{22} = \frac{1}{2}\Omega$, as illustrated in Figure 10.

First, consider region $\Omega_1$. We need to interpolate around it by selecting 27 points $(x_{0}, y_{0}, z_{0})$, $(x_{0}, y_{0}, z_{-2})$, $(x_{0}, y_{0}, z_{-1})$, $(x_{0}, y_{-1}, z_{0})$, $(x_{0}, y_{-1}, z_{-1})$, $(x_{0}, y_{-1}, z_{-2})$, $(x_{0}, y_{-2}, z_{0})$, $(x_{0}, y_{-2}, z_{-1})$, $(x_{0}, y_{-2}, z_{-2})$, $(x_{1}, y_{0}, z_{0})$, $(x_{1}, y_{0}, z_{-1})$, $(x_{1}, y_{0}, z_{-2})$, $(x_{1}, y_{-1}, z_{0})$, $(x_{1}, y_{-1}, z_{-1})$, $(x_{1}, y_{-1}, z_{-2})$, $(x_{1}, y_{-2}, z_{0})$, $(x_{1}, y_{-2}, z_{-1})$, $(x_{1}, y_{-2}, z_{-2})$, $(x_{2}, y_{0}, z_{0})$, $(x_{2}, y_{0}, z_{-1})$, $(x_{2}, y_{0}, z_{-2})$, $(x_{2}, y_{-1}, z_{0})$, $(x_{2}, y_{-1}, z_{-1})$, $(x_{2}, y_{-1}, z_{-2})$, $(x_{2}, y_{-2}, z_{0})$, $(x_{2}, y_{-2}, z_{-1})$, $(x_{2}, y_{-2}, z_{-2})$ on the interior of the wood as shown in Figure 11.

The interpolating polynomial constructed are

and

Therefore,

where,

and

Second, we consider $\iint_{\Omega_{21}}\frac{\partial V}{\partial n}ds$. The interpolation nodes of field $\Omega_{21}$ are shown in Figure 12.

Equally,

among it

and

Finally, we consider $\iint_{\Omega_{22}}\frac{\partial V}{\partial n}ds$. The interpolation nodes of field $\Omega_{22}$ are shown in Figure 13.

With the same way,

where,

and

Thus, Eqs (4.52), (4.57), and (4.62) constitute the discrete scheme of the upper right edge Eq (4.48), organized as

Similarly, we can obtain discrete equations for the other 11 edges of the wood.

4.2.   Error estimation of discrete schemes

We take the first equation and third equation in system (3.1) as examples for error estimation, and the remaining equations follow a similar approach.

1) Error estimation of the first equation

Since the first equation $\nabla^2V(x, y, z) = 0$ employs the Lagrange interpolation function, the remainder term of the Lagrange interpolation is given by

where, $\xi \in\left(x_{i-1}, x_{i+1}\right), \eta \in\left(y_{j-1}, y_{j+1}\right), \zeta \in\left(z_{k-1}, z_{k+1}\right)$. If $\left|\frac{\partial^{9} V(\xi, \eta, \zeta)}{\partial x^{3} \partial y^{3} \partial z^{3}}\right| \leq M$, we have

Therefore, the error of the first equation is

Due to the identical structure of the three terms, we perform only the calculation for $\frac{\partial^{2} R}{\partial x^{2}}$, as the other two terms can be computed similarly. First, integrating with respect to $x$, we have

Next, integrating with respect to $y$, we have

At this point, integrating with respect to $y_2$, we have

Similarly, the integral for $y_1$ is

Likewise, the integral with respect to the $z$ direction is also $O(h_z^6)$. Therefore,

However, to ensure accuracy, it is also necessary to divide by the integration factor. The integration factor is

Therefore, the accuracy of Eq (4.76) should be $O(h_x^4+h_y^3+h_z^3)$. Thus, the error of Eq (4.71) is

Hence, the local truncation error of the discrete scheme is $O(h_x^3+h_y^3+h_z^3)$, and the global error is $O(h_x^2+h_y^2+h_z^2)$.

2) Error estimation of the third equation

By the Taylor expansion, we have

and

By $4\times(4.78)-(4.79)$ and $(4.81)-4\times(4.80)$, we eliminate $V''_{i, j, k}$ and obtain

and

Thus,

In practical situations, $\tilde{\varepsilon}\neq{\varepsilon}$. Therefore, the error is $O{\left(h_{z}^2\right)}$. Similarly, it can be concluded that the error in the $x$ and $y$ directions are $O{\left(h_{x}^2\right)}$ and $O{\left(h_{y}^2\right)}$. As a result, the third equation error is $O(h_x^2+h_y^2+h_z^2)$.

4.3.   Numerical analysis

We illustrate the theory in this paper through a numerical example using the multiple varying bounds integral method.

The two electrode plates occupy the domain

and

The measured wood occupies the domain

Moreover, we set the artificial boundaries to $[-10, 10]\times[-10, 10]\times[-10, 10]$. In the above domains, the unit of length is centimeters and the step size is $\alpha_x = \alpha_y = \alpha_z = \alpha = h_x = h_y = h_z = h = 1$. It is clear that the domain is symmetric. Based on the physical principles, $V_2 = -V_1$. Suppose $V_1 = 1$ and $V_2 = a$ (where $a$ is an unknown to be determined). If, through numerical computation, the value of $a$ is close to -1, it indicates that the multiple varying bounds integral method is reasonable.

We perform numerical experiments on the above example and obtain the data in Table 1. These data are consistent with the physical principle that $C$ increases as $\varepsilon$ increases. Based on the data, the corresponding image is shown in Figure 14(a). It can be seen that, compared with other methods ^[30], our method eliminates the initial oscillations and maintains the original trend of $C$ variation. Figure 14(b) shows the fitted function image derived from Reference ^[30]. Then, we perform a regression analysis on these data to establish a regression model:

First, we evaluate the goodness of fit for Eq (4.85). A higher goodness of fit indicates a stronger ability of the model to predict the dependent variable. As calculated from the data presented in Table 2, the determination coefficient $R^{2} = \frac{SSR}{SST} = \frac{603.42}{665} = 0.91$, where SSR represents the sum of squares regression and SST denotes the sum of squares total for Eq (4.85). This indicates that, when $C$ is given, the regression model can more accurately predict the value of $\varepsilon$.

In addition, we consider the confidence degree, which represents the reliability of $\varepsilon$. The higher the confidence degree, the stronger the reliability. Table 2 shows that $F(1, 18) = 176.39$, and according to the F-distribution table, $F(1, 18) = 8.29 \; (\alpha = 0.01)$. Because 176.39 $>$ 8.29, it can be concluded that the confidence degree of the model is greater than 0.99. This means the values of $\varepsilon$ obtained through (4.85) are highly reliable and more concentrated around the true value. The results of this regression analysis indicate that the functional relationship model between $\varepsilon$ and $C$ that we established can effectively predict the value of $\varepsilon$ from $C$, with small errors. It can be seen that the multiple varying bounds integral method is a valid numerical method for this kind of problem.

5.   Conclusions

Wood moisture content detection is an inverse problem in multi-physics fields. Based on the high precision demand for moisture content detection, we propose a new numerical method, that is, the multiple varying bounds integral method. Moreover, in different physics fields, we choose different discrete methods to construct numerical schemes. For the physical field where the unknown function is discontinuous, this field has to be divided into several additional parts. For each smaller part, such as faces, vertices, and edges, we build corresponding interpolation functions and handle the integral to meet the precision requirements for the engineering problem. Moreover, we carry out numerical experiments and perform regression analysis to obtain the function relationship between $\varepsilon$ and $C$. This model obeys the physical principle that $C$ increases with an increase in $\varepsilon$. Moreover, the established regression model $R^2$ is greater than 0.91, indicating a high goodness of fit and effectively reflecting the relationship between $\varepsilon$ and $C$. This demonstrates that the data derived from this method are valid. Additionally, this numerical method is applicable to other engineering problems.

On the other hand, the discrete scheme constructed in this paper has second-order precision. There are 2982 unknowns, and the computation time is 875.57 s. If a higher precision is desired, we can consider adding interpolation nodes to improve precision. Moreover, this would also increase the computational cost.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to thank all the references which arehelpful for this article and the valuable suggestions which are put out by experts andreaders. These have improved this paper greatly and made it perfect. This paper was supported by the National Natural Science Foundation of China (Grant No. NSFC11526064), the Fundamental Research Funds for the Central Universities (Grant No. 3072024XX2402), and Harbin Engineering University (Grant No. KYWZ220240710).

Conflict of interest

The authors declare there are no conflicts of interest.