# Research on the wood processing method of helium-assisted laser process – Journal of Wood Science

#### ByChunmei Yang, Xinchi Tian, Bo Xue, Qingwei Liu, Jiawei Zhang, Jiuqing Liu and Wenji Yu

Aug 27, 2022

When the laser cooperates with the inert gas emitted from the coaxial nozzle, the initial velocity of the gas is about 10 m/s. The emitted gas absorbs energy in the space and forms turbulence. Meanwhile, the gas is also diffused by the gas flow and the concentration decreases along the direction perpendicular to the axis. As the cutting quality of the machined part is caused by high-temperature ablation, it is vital to employ the mathematical model to solve the ablation area on the machined surface. The gas jet flow during laser cutting is shown in Fig. 1.

### Establishment of model

The three elements of combustion are combustibles, combustion-supporting gas and ignition point. When the laser is processed under the helium package, the scattering degree of the laser would further affect the generation and transfer of heat. In the meantime, the helium would diffuse at a certain initial velocity after being sprayed from the nozzle, so the gas concentration would also be reduced in the plane perpendicular to the nozzle.

In the mixture of helium and oxygen, when the oxygen content is higher than 12.5%, wood can burn, and when it is lower than 12.5%, it is not combustible. Therefore, it is crucial to calculate the proportion of oxygen at any points. The oxygen content is also affected by the diffusion of helium. Therefore, it is critical to calculate the helium concentration at each point. Meanwhile, whether the processed surface burns also depends on the amount of heat transmitted through the wood. The ignition point of the wood is generally between 200 and 300 ℃. The ignition point of the original wood studied in this paper is approximately 280 ℃. Comprehensively, whether the machined surface burns or not is determined by the oxygen concentration and surface temperature.

#### Basic assumptions

In order to simplify the mathematical model and solution process, the following assumptions are put forward:

1. (1)

The material is anisotropic in the heat diffusion direction, namely, the heat diffusion rate at each point remains constant but is randomly distributed among each other.

2. (2)

The mass of inert gas is conserved in the diffusion process. The gas concentration in the diffusion process obeys the normal distribution in the y and z directions, and the flow rate does not exceed 68 m/s, which is regarded as an incompressible fluid.

3. (3)

In the process of laser spraying, each point on the laser beam is supposed to be isothermal, regardless of the temperature drop caused by the emission distance.

4. (4)

The indoor gas flow is stable and the helium concentration does not change with time.

5. (5)

As the laser source is only 1 mm relative to the processing distance, the blocking effect of wood thickness processed by steam plasma on laser is ignored.

6. (6)

Ignore the combustion effects of small sawdust produced by processing.

#### Gas diffusion model

Taking the injection point as the origin and the injection direction as the z-axis direction, a spatial rectangular coordinate is established. The effective point source is located at the coordinate origin o, and the average velocity direction is parallel to and in the same direction as the z-axis. The Reynolds number of helium jet flow is higher than the critical value and flows as turbulent flow. When the jet stream ejects the inert gas from the nozzle at the velocity v0, the surrounding fluid is continuously pumped in the flow, the width of the jet is expanded, and the velocity of the main body of the gas is constantly reduced. The jet diffusion angle is θ. According to the momentum theorem, the momentum along the jet velocity direction remains unchanged at any interface of the jet, so it can be obtained:

$$smallint rho v^{2} dA = rho_{0} v_{0}^{2} A_{0} = rho_{0} pi R_{0}^{2} v_{0}^{2} = 2pi smallint_{0}^{R} rho v^{2} rdr.$$

(1)

In the above formula, A0 is the cross-sectional area at the nozzle, R0 is the nozzle radius and the helium density at the nozzle is (rho_{0}). Since the jet density is the same as that of the surrounding air, that is (rho = rho_{0}), formula (1) can be rewritten as:

$$2intlimits_{0}^{{frac{R}{{R_{0} }}}} {left( {frac{v}{{v_{0} }}} right)^{2} frac{r}{{R_{0} }}dleft( {frac{r}{{R_{0} }}} right)} = 1.$$

(2)

Also (frac{r}{{R_{0} }} = frac{r}{R} cdot frac{R}{{R_{0} }}), (frac{R}{{R_{0} }}) determined by the distance from any cross section to the jet pole. Therefore, (frac{v}{{v_{0} }} = frac{v}{{v_{m} }} cdot frac{{v_{m} }}{{v_{0} }}). So Eq. (2) can be exhibited as:

$$2left( {frac{{v_{m} }}{{v_{0} }}} right)^{2} left( {frac{R}{{R_{0} }}} right)^{2} intlimits_{0}^{1} {left( {frac{v}{{v_{m} }}} right)^{2} frac{r}{R}dleft( frac{r}{R} right)} = 1.$$

(3)

Herein, vm is the velocity on the axis of inert gas injection. In terms of the above formula, the inert gas flow on any section is:

$$q = 2pi intlimits_{0}^{R} {vrdr} = 2pi R_{0}^{2} v_{0} left( {frac{R}{{R_{0} }}} right)^{2} frac{{v_{m} }}{{v_{0} }}intlimits_{0}^{1} {frac{v}{{v_{m} }}frac{r}{R}dleft( frac{r}{R} right)} = 2q_{0} left( {frac{R}{{R_{0} }}} right)^{2} frac{{v_{m} }}{{v_{0} }}intlimits_{0}^{1} {left[ {1 – left( frac{r}{R} right)^{frac{3}{2}} } right]^{2} frac{r}{R}dleft( frac{r}{R} right)}.$$

(4)

Among them, q0 the inert gas flow initially ejected from the nozzle can be calculated at any point according to Eq. (4).

The concentration change of inert gas in the atmosphere is jointly determined by jet and diffusion. As for the gas diffusion, we adopt the steady Gaussian diffusion model. The point source is located at the coordinate origin, so its diffusion can be regarded as gas diffusion at height 0. The diffusion in air is a two-dimensional normal distribution with two coordinate directions of y and z. When the random variables in the two coordinate directions are independent, the distribution density is the product of the one-dimensional normal distribution density function in each coordinate direction. According to the assumption of Gaussian gas diffusion model, taking (mu = 0),the gas concentration distribution function at any points in the velocity direction of point source can be acquired as follows:

$$Cleft( {x,y,z} right) = A(x)exp left[ { – frac{1}{2}left( {frac{{y^{2} }}{{sigma_{y}^{2} }} + frac{{z^{2} }}{{sigma_{z}^{2} }}} right)} right].$$

(5)

In the above formula, C—concentration of pollutants at space point (x, y, z); A(x)—undetermined function; (sigma_{y} ,sigma_{z})—standard deviations in the horizontal and vertical directions, i.e., the diffusion parameters in the y and z directions, mm.

According to the law of conservation of mass and the continuity theorem, on any cross section of air flow perpendicular to the z-axis:

$$q = intlimits_{ – infty }^{ + infty } {intlimits_{ – infty }^{ + infty } {Cud{text{y}}dz} } .$$

(6)

Equation (1) is substituted into Eq. (2) by the velocity stability condition, a is independent of x and y, and (int_{ – infty }^{ + infty } {exp left( {{{ – t^{2} } mathord{left/ {vphantom {{ – t^{2} } 2}} right. kern-nulldelimiterspace} 2}} right)} {kern 1pt} ,dt = sqrt {2pi }). Undetermined coefficients can be obtained by integration:

$$A(x) = frac{q}{{2pi sigma_{y} sigma_{z} }}.$$

(7)

Substituting (3) into (1):

$$C(x,y,z) = frac{q}{{2pi sigma_{{text{y}}} sigma_{z} }}exp left[ { – frac{1}{2}left( {frac{{y^{2} }}{{sigma_{y}^{2} }} + frac{{z^{2} }}{{sigma_{z}^{2} }}} right)} right].$$

The diffusion coefficients (sigma_{y}) and (sigma_{z}) are related to the air stability and the vertical distance z, which increases with the increase of z. Therefore, the concentration of inert gas at each point in the coordinate range can be derived by Eq. (3). Since 1 m3 air mass is about V = 1286 g [14], the proportion of helium at any points can be obtained by substituting (frac{C}{V} times 100{text{% }}), and then utilizing the oxygen content in the air, the proportion of oxygen at any points is (A{text{% }} = left( {1 – frac{C}{V}} right) times 20{text{% }}), so the oxygen concentration and proportion at any points can be calculated to explore whether it could satisfy the combustion-supporting conditions.

#### Laser heat transfer model

Laser is an electromagnetic radiation wave. The laser used for cutting is generally output in the mode of fundamental Gaussian mode. The cross-section spot is circular and the light intensity distribution basically satisfies the Gauss distribution. According to Professor Siegman’s theory [15], heat flux density of laser heat source can be expressed as:

$$q(r) = frac{{3q_{max } }}{{pi w({text{z}})^{2} }}e^{{( – frac{{3r^{2} }}{w(z)})}} .$$

(8)

In this formula, q(r)—heat flux at radius R, W/mm2; (q_{max })—the maximum value of heat flux density of laser beam, viz. the value at the center of beam waist radius, can be replaced by laser power value; W(z)—effective radius of beam at z coordinate, mm.

The beam passing through the lens will have a certain scattering, that is Rayleigh scattering. The divergence of the beam will cause a great impact on the heat transmission. The expression of the effective radius of the beam at the z coordinate is:

$$w(z) = w_{0} sqrt {1 + (frac{{z – z_{0} }}{{z_{R} }})^{2} }$$

(9)

where w(z) is the Gaussian beam radius at the z coordinate, (w_{0}) is the radius at the beam waist, z0 is the ordinate at the beam waist, and zR is the Rayleigh constant. Its value is determined by the scattering degree of the light source.

As there are no combustibles in the air, combustion would not occur even at high temperature, but the surface of processed wood is ablated by laser at high temperature. When the surface temperature is higher than the ignition point of wood and satisfied combustion-supporting conditions, combustion would occur on the surface of wood. Therefore, the boundary between combustion zone and non-combustion zone should be located on the processing surface of wood and perpendicular to the laser centerline. The thickness of wood processed by laser is only 2 mm, so the temperature drop caused by processing is ignored. The heat conduction rate between wood micro elements is roughly the same. Therefore, the heat transmitted at any two points perpendicular to the laser processing axis and equal to the vertical distance of the laser beam in the wood is the same. Meanwhile, the laser is a beam with concentrated energy and small scattering. Therefore, in the wood processed by laser, the heat conduction mode is primarily heat conduction, with small heat convection and ignoring the role of heat radiation. Together, along the direction perpendicular to the laser axis, the heat conduction can be described by the thermal conduction differential equation with constant physical properties, two-dimensional, unsteady state and internal heat source:

$$frac{partial T}{{partial t}} = aleft( {frac{{partial^{2} T}}{{partial x^{2} }} + frac{{partial^{2} T}}{{partial y^{2} }}} right) + frac{{dot{phi }}}{rho c}$$

(10)

In this formula, T—temperature, ℃; t—time, s; a—thermal diffusivity, (a = frac{lambda }{rho c}), m2/s, taken as 1.7 × 10–7; (lambda)—thermal conductivity, W/(m·K), represented as 0.2; (rho)—material density, kg/m3, Fraxinus mandshurica signified as 686 kg/m3; (dot{varphi })—generation heat of heat source per unit volume, expressed as 6.67 × 10–10 J; c—specific heat of material, Fraxinus mandshurica indicated as 1.72 × 103 J/(kg·℃).

The transmission rate of laser heat in any direction of XOY plane is consistent. Therefore, a cylindrical coordinate system can be established in the heat diffusion direction. The two-dimensional heat transfer becomes a one-dimensional heat conduction problem along the radius direction. If the thermal conductivity of wood is a fixed value, the thermal conductivity differential equation can be expressed as:

$$frac{d}{dr}left( {rfrac{dT}{{dr}}} right) = 0$$

(11)

The indoor initial temperature is 18 ℃. Since the laser heating temperature is 1200 ℃ and the ignition point of wood is 280 ℃, when the temperature of mixed gas containing more than 12.5%, viz. above 280 ℃, wood would be ignited. Therefore, the initial conditions of this model are:

$$t = 0,forall {kern 1pt} ,r ne 0,;t = {text{const = 18}} ^circ {text{C,}}$$

$$t = 0,;r = 0,;T = 1200 ;^circ {text{C}}{.}$$

When the oxygen concentration is higher than 12.5%, the mixed gas would burn. The oxygen content ratio of the mixed gas can be described by a certain boundary. In the XOY plane, its function expression is (y = f_{1} (x)). When the temperature is lower than 280 ℃, the gas would not burn on the premise of sufficient oxygen content, and its boundary function expression is (y = f_{2} (x)). Therefore, the boundary conditions of this model are exhibited as follows:

$$A% , = ,{12}.{5}% ,;y, = ,f_{{1}} (x),$$

$$T, = ,{356};^circ {text{C}},;y = f_{2} (x),$$

$$y = f_{1} (x) wedge f_{2} (x),$$

$$y = y_{infty } ,T = const,$$

$$h_{0} (T_{infty } – T(r,t)) = left. { – lambda frac{partial T(r,t)}{{partial r}}} right|_{{r = r_{0} }} ,$$

where h0 is the convective heat transfer coefficient of the laser beam surface, and r0 is the laser heating heat flow radius.

For the heat transfer process of laser processing, the comprehensive model is established as follows:

$$left{ {begin{array}{*{20}c} {{text{Jet equation: }}q = 2q_{0} left( {frac{R}{{R_{0} }}} right)^{2} frac{{v_{m} }}{{v_{0} }}int_{0}^{1} {left[ {1 – left( {frac{r}{R}} right)^{{frac{3}{2}}} } right]^{2} frac{r}{R}dleft( {frac{r}{R}} right)} } \ {{text{Gas diffusion equation: }}C(x,y,z) = frac{q}{{2pi sigma _{y} sigma _{z} }}exp left[ { – frac{1}{2}left( {frac{{y^{2} }}{{sigma _{y}^{2} }} + frac{{z^{2} }}{{sigma _{z}^{2} }}} right)} right];A% = left( {1 – frac{C}{V}} right) times 20% ,} \ {{text{Heat source equation: }}q(r) = frac{{3q_{{max }} }}{{pi w(z)^{2} }}e^{{left( { – frac{{3r^{2} }}{{w(z)}}} right)}} ,} \ {{text{Governing equation: }}frac{{partial T}}{{partial t}} = aleft( {frac{{partial ^{2} T}}{{partial x^{2} }} + frac{{partial ^{2} T}}{{partial y^{2} }}} right) + frac{{dot{phi }}}{{rho c}},} \ {{text{Heat transfer equation: }}frac{d}{{dr}}left( {rfrac{{dT}}{{dr}}} right) = 0,} \ {{text{Boundary condition: }}A% = 12.5% ,y = f_{1} (x);T = 356,} \ {;quad quad y = f_{1} (x) wedge f_{2} (x).} \ end{array} } right.$$

### Solution of model

#### Finite difference method

The finite difference method is to express the differential equation as a difference equation defined on discrete lattice points (i.e., the numerical relationship between lattice points reflected by the differential equation), and iteratively calculate the numerical value on the unknown boundary through the difference relationship between similar lattice points according to the given boundary conditions. The main difference between explicit and implicit difference methods is that the discretization approximate representation of differential terms in the difference equation instead of calculating differential equations is different, so that the explicit difference method can directly deduce the values on the target boundary from the known boundary conditions, while the implicit difference method constantly needs to solve the equations for recursive calculation. However, it elucidates that the results calculated by the difference method are not necessarily stable, and the small error may be amplified, resulting in the error of the calculation results.

To solve the partial differential equation problem by the finite difference method, the continuous problem must be discretized, namely, the solution area is meshed, the continuous solution area is replaced by a finite number of discrete grid nodes, and then the differential quotient of the partial differential equation is substituted by the difference quotient to derive the discretized difference equation system. Finally, the difference equations are solved by direct method or iterative method to achieve the numerical approximate solution of the differential equation.

When established the difference equations of heat conduction and gas diffusion, the solution area needs to be meshed by the internal node method. As shown in Fig. 2, the node temperature and gas concentration can approximately represent the temperature and gas concentration of the whole small area, which is convenient to deal with the uneven problem. Meanwhile, the data of boundary conditions are employed in the boundary points.

Assuming that in the grid, the time step (tau), the space step along the x-axis direction h, the space step along the y-axis direction k, and the step along the radius direction n, the explicit difference scheme is used for discretization in this model to solve the temperature and gas concentration in the jth grid. The discretization scheme in the model is:

$$left{ {begin{array}{*{20}l} {{text{Jet equation: }}q_{j} = 2q_{0} left( {frac{{R_{j} }}{{R_{0} }}} right)^{2} frac{{v_{mj} }}{{v_{0} }}int_{{r_{j} – frac{n}{2}}}^{{r_{j} + frac{n}{2}}} {left[ {1 – left( {frac{{r_{j} }}{R}} right)^{frac{3}{2}} } right]frac{{r_{j} }}{R}dleft( frac{r}{R} right),} } hfill \ {{text{Diffusion equation: }}C_{j}^{n} = frac{{q_{j}^{n} }}{{2pi sigma_{yj}^{n} sigma_{zj}^{n} }}exp left[ { – frac{1}{2}left( {frac{{(y_{j}^{n} )^{2} }}{{(sigma_{yj}^{n} )^{2} }} + frac{{(z_{j}^{n} )^{2} }}{{(sigma_{zj}^{n} )^{2} }}} right)} right],} hfill \ {{text{Governing equation: }}frac{{T_{j}^{n + 1} – T_{j}^{n} }}{tau } = a_{j} left( {frac{{T_{j + 1}^{n} – 2T_{j}^{n} + T_{j – 1}^{n} }}{{h^{2} }} + frac{{T_{j + 1}^{n} – 2T_{j}^{n} + T_{j – 1}^{n} }}{{k^{2} }}} right) + frac{{dot{phi }_{j} }}{{rho_{j} c_{j} }},} hfill \ {{text{Heat transfer equation: }}frac{{T_{j}^{n + 1} – T_{j}^{n} }}{n} + r_{j}^{n} frac{{T_{j + 1}^{n} – 2T_{j}^{n} + T_{j – 1}^{n} }}{{n^{2} }} = 0.} hfill \ end{array} } right.$$

For the explicit difference scheme, the discrete solution of unsteady heat transfer process needs to be considered in the stability condition. The above explicit difference formula displays that the temperature of the time node n + 1 on the space node I is affected by the adjacent points on the left and right sides, and the stability limit condition (Fourier grid number limit) needs to be satisfied, otherwise there would be unreasonable oscillatory solutions[16].

$$left{ {begin{array}{*{20}l} {F_{ODelta } = frac{lambda Delta t}{{rho c(Delta x)^{2} }}({text{grid number}}),} hfill \ {F_{ODelta } le frac{1}{2}({text{inner node constraints}}),} hfill \ {F_{ODelta } le frac{1}{{2left( {1 + frac{hDelta x}{lambda }} right)}}({text{boundary constraints}}).} hfill \ end{array} } right.$$

After discretizing the unsteady process, the gas diffusion and heat transfer equations can be solved by the control equations and boundary conditions.

#### Solution of gas jet diffusion model

Under the idea of finite difference method, MathWorks Inc’s MATLAB2015a software is employed to set the flow rate at the nozzle to 10 m/s and the height to 1 mm. The gas will expand due to energy absorption during injection. Meanwhile, the density of helium is less than that of air. After leaving the nozzle, helium will be subjected to the air buoyancy, resulting in an upward trend, which would further lead to the expansion of the outer diameter of the gas. The distribution of gas trajectory after helium is injected out of the nozzle is shown in Fig. 3.

After the gas is injected through the nozzle, the instability of the gas itself will diffuse and the concentration of helium will decrease along the X and Y directions. The greater the helium concentration is, the smaller the air concentration is. Therefore, along a fixed direction, the gas will diffuse with the flow, and the concentration will gradually decrease during diffusion, resulting in the gradual increase of oxygen content. Therefore, the proportion of oxygen can be achieved by helium concentration. The distribution of helium diffusion concentration is shown in Fig. 4a:

The proportion of oxygen in the air is about 20% and combustion support can be realized when the oxygen content is higher than 12.5%. When the oxygen content is 12.5%, the proportion of air is 62.5%, and the helium concentration is 37.5%. Therefore, when the helium concentration is 37.5%, other gases can realize combustion support. The concentration distribution line of gas diffusion is shown in Fig. 4b, in which the red dotted line is the isoconcentration line of helium concentration of 37.5%. When the oxygen concentration is higher than 12.5%, the combustion-supporting effect can be realized.

#### Solution of laser heat transfer model

According to the assumptions, the laser heating temperature does not decrease with the processing depth, and the wood is a poor conductor with low thermal conductivity. Meanwhile, the wood is an anisotropic material with significant differences in porosity, moisture content and other physical properties. Combustion would occur when the wood surface temperature is higher than its ignition point of 280 ℃, but due to the anisotropy of the wood material, the heat conduction efficiency in each direction is different. When the moisture content in one direction is high, the heat transfer efficiency of wood will also enhance. Set the wood performance parameter as the original wood parameter, and take the processing time t = 2 s to achieve the wood surface temperature distribution, as shown in Fig. 5a.

As the thermal conductivity of wood is relatively poor, the temperature decreases rapidly during heat transfer. In addition, Fraxinus mandshurica is composed of fiber structure, and the transverse and longitudinal thermal conductivity are quite different. Due to the influence of wood fiber structure, the carbon content in the along-grain cutting mode is significantly less than that of the cross-grain cutting mode, whereas the thermal conductivity of longitudinal fiber and transverse wood cells are also very small. Therefore, the temperature drop on the wood surface changes rapidly with coordinates. Chosen various temperature and laser source as the origin, a plane rectangular coordinate is established on the processing surface. The factors affecting heat conduction such as porosity and moisture content on the wood surface are randomly distributed, and the isotherm is shown in Fig. 5b.

The feasible areas of gas concentration and temperature are shown in Fig. 6. Selected 37.5% to mark the critical line of combustion-supporting gas in the concentration change diagram. The oxygen content in the area outside the critical line (marked as area 1) can achieve the combustion-supporting effect. The temperature of 280 ℃ is selected to mark the isotherm on the wood surface. The wood may burn in the inner area of 280 ℃ (marked as area 2), and the intersection area of areas 1 and 2 will produce combustion and forming a combustion area.

It can be signified from the above figures that the oxygen content of the machining surface would be reduced after adding helium, which destroys the combustion conditions. As such, when using inert gas-assisted laser processing, the carbide generated by combustion will be reduced, the slit width declined correspondingly, and the machining accuracy as well as cutting quality will be further improved.