Numerical Simulation of Two-Dimensional Natural Convection in a Confined Environment with the Vorticity-Stream Function Formulation ()

Mame Khady Kane^{}, Cheikh Mbow^{}, Mamadou Lamine Sow^{}, Joseph Sarr^{}

Department of Physics, Cheikh Anta Diop University of Dakar, Dakar, Senegal.

**DOI: **10.4236/ojfd.2018.81004
PDF
HTML XML
1,106
Downloads
2,539
Views
Citations

Department of Physics, Cheikh Anta Diop University of Dakar, Dakar, Senegal.

A numerical study is presented on the problem of 2D natural convection in a differentially heated cavity. The equations governing this unsteady flow phenomenon were solved using the vorticity-stream function formulation of the Navier-Stokes equations and heat. The results obtained are compared with the results of the literature and make it possible to validate this approach. In this work, we studied the heat transfer in a cavity and we determined the variation of the local Nusselt number which allows obtaining the rate of thermal transfer by convection in an enclosure. We analyzed thermal fields for different Rayleigh numbers by selecting two points to visualize temperature fluctuations over time. Thus, the creation of the ascending and descending movements of the fluid inside the cavity was analyzed. We have also established temperature histograms for the graphical presentation of the temperature distribution. The modeling of the two-dimensional problem was established using a “Fortran 90” calculation code. The results also show the different vorticity contour maps in laminar flow regime. We have presented our results of numerical simulations using a visualization tool. The Rayleigh number varies in the range of 10^{3} to 10^{6} for a Prandtl number equal to 0.72 corresponding to air.

Keywords

Numerical Study, Differentially Heated Cavity, Vorticity, Modeling

Share and Cite:

Kane, M. , Mbow, C. , Sow, M. and Sarr, J. (2018) Numerical Simulation of Two-Dimensional Natural Convection in a Confined Environment with the Vorticity-Stream Function Formulation. *Open Journal of Fluid Dynamics*, **8**, 44-58. doi: 10.4236/ojfd.2018.81004.

1. Introduction

In general, numerical resolutions of heat transfer problems [1] mainly related to convection [2] and defined in various geometries, are performed by finite element methods [3] or finite volume [4] . Some studies have focused on natural convection and we have referred to some important works that can serve as a backdrop for this work. The study of the thermal and vorticity fields has been presented. We must know that the vorticity being connected to the circulation; a line of vorticity is defined analogously to a stream function. It is at a given instant a curve which is tangent at all points to the local vorticity. Differently from the square cavity [5] , natural convection has also been studied in a rectangular enclosure with partially active walls for nine different heating locations [6] . For this, the hot and cold regions have been analyzed (above, middle and bottom), to identify, where the rate of heat transfer is maximum and minimum. Natural convection has important applications in solar energy collectors, nuclear reactors, cooling of electronic equipments, oil recovery, ventilation of buildings during fires, double glazing for thermal insulation, aerospace, etc.

We are interested here in the numerical 2D simulation of natural air convection (Pr = 0.72) confined in a square cavity heated on the central halves of its vertical walls. It was interesting to conduct research in the study of the behavior of the flow for different numbers of Rayleigh and showing the effect of time on the properties of the fluid. In this study, we used the finite volume method with quadrilateral control volumes. To solve the theoretical model, it is necessary to go through a numerical processing; for this, the equations were rendered first in conservative form. We began our discretization by adopting the notion of the staggered mesh [7] , which means that the pressure and the temperature are indicated in the center of the mesh while the velocities are indicated on the centers of the four faces of the staggered mesh and we have U_ {e}, U_ {w} on the abscissa and V_ {n}, V_ {s} on the ordinate. For temporal discretization, we used Adams Bashforth’s method for time integration. We recall that we have essentially given the results of calculations for air, because it is one of the most used fluids [8] [9] . However, the assumptions on which our mathematical model is based make it applicable to other Newtonian fluids.

2. Materials and Methods

2.1. Basic Equations

We consider the two-dimensional natural convection phenomenon in a confined environment differentially heated over the centers of its vertical walls as shown in Figure 1.

The governing equations, the fluid flow are written in dimensionless form taking into account the hypotheses following simplifiers:

-The fluid is Newtonian, viscous and transparent.

-The regime is laminar.

-The flow is unsteady and the fluid is incompressible.

-The two-dimensional flow is in the plane (Ox, Oy).

-The physical properties of the fluid are constant, except the density in the gravitational term of the equations of motion, in which it varies linearly with

Figure 1. Cavity studied with conditions to the limits: boundary conditions (dynamic and thermal) of the system under consideration.

temperature according to Boussinesq’s hypotheses [10] :

$\rho ={\rho}_{0}\left[1-\beta \left(T-Tr\right)\right]$ (1)

We have dimensionless quantities:

$\begin{array}{l}X=\frac{x}{A},\text{\hspace{0.17em}}Y=\frac{y}{A},\text{\hspace{0.17em}}\tau =\frac{t}{\left(\frac{{A}^{2}}{\alpha}\right)},\text{\hspace{0.17em}}\theta =\frac{T-Tr}{Tg-Td},\text{\hspace{0.17em}}U=\frac{u}{\left(\frac{\alpha}{A}\right)},\\ V=\frac{v}{\left(\frac{\alpha}{A}\right)},\text{\hspace{0.17em}}\Psi =\frac{\psi}{\alpha},\text{\hspace{0.17em}}\zeta =\frac{\omega}{\left(\frac{\alpha}{{A}^{2}}\right)}\end{array}$ (2)

In the vorticity-stream function formulation, the dimensionless Navier-stokes equations are written:

$\frac{\partial \zeta}{\partial \tau}-\frac{\partial \Psi}{\partial X}\frac{\partial \zeta}{\partial Y}+\frac{\partial \Psi}{\partial Y}\frac{\partial \zeta}{\partial X}=Pr\cdot di{v}^{2}\zeta +RaPr\frac{\partial \theta}{\partial X}$ (3)

$\frac{\partial \theta}{\partial \tau}-\frac{\partial \Psi}{\partial X}\frac{\partial \theta}{\partial Y}+\frac{\partial \Psi}{\partial Y}\frac{\partial \theta}{\partial X}=di{v}^{2}\theta $ (4)

$\zeta =-di{v}^{2}\Psi $ (5)

In this study, we used the finite volume method with quadrilateral control volumes and a staggered mesh. The latter is the subdivision of the field of study into longitudinal and transverse grids whose intersection represents a node, where the variables P and θ are located while the components U and V of the velocity vector lie in the middle of the segments connecting two nodes adjacent. After discretization of the differential transport equations we obtain a system of nonlinear algebraic equations, these equations describe the discrete properties of the fluid at the nodes in the solution domain.

2.2. Validation of the Code

To examine the reliability of our computer code in Fortran 90, using the finite volume method, we compared our results with those of the literature. The results are graphically displayed and show the effect of internal heat generation. We first used a classical test of G. De Vahl Davis Benchmark [11] and P. M. Gresho [12] , we studied the fluid flow in a square cavity consisting of two adiabatic horizontal walls and vertical walls (left and right) subjected to constant hot and cold temperatures respectively at their centers.

In addition, numerical results in terms of Nusselt numbers were compared with those in the literature. We found a good agreement for a 81 × 81 mesh as shown in Table 1.

2.3. Numerical Procedure

The governing Equations (3)-(5) are discretized using the finite volume formulation with the central difference scheme. The region of interest is covered with m vertical and n horizontal uniformly spaced grid lines. The numerical solution is truly transient. An iterative process is employed to find the temperature and vorticity fields. The process is repeated until the following convergence criterion:

$\mathrm{max}\left\{\frac{{\displaystyle {\sum}_{i,j}\left|{\theta}_{i,j}^{n+1}-{\theta}_{i,j}^{n}\right|}}{{\displaystyle {\sum}_{i,j}\left|{\theta}_{i,j}^{n+1}\right|}},\frac{{\displaystyle {\sum}_{i,j}\left|{\Psi}_{i,j}^{n+1}-{\Psi}_{i,j}^{n}\right|}}{{\displaystyle {\sum}_{i,j}\left|{\Psi}_{i,j}^{n+1}\right|}}\right\}\le {10}^{-5}$ (6)

In the above expression n is any time level: in the square cavity, the flow depends on the time. Unsteady laminar natural convection in a cavity with partially thermally active side walls is studied numerically (code Fortran 90). The calculation time was of the order of:

-for the mesh 41 * 41

$\{\begin{array}{l}12\text{\hspace{0.17em}}\mathrm{minutes}\left(Ra=1000\right)\\ 11\text{\hspace{0.17em}}\mathrm{minutes}\left(Ra=10000\right)\end{array}$ (7)

Table 1. Local Nusselt number for Ra = 100,000.

-and for the mesh 81 * 81

$\{\begin{array}{l}57\text{\hspace{0.17em}}\mathrm{minutes}\left(Ra=100000\right)\\ 01\text{\hspace{0.17em}}\text{h}\text{\hspace{0.17em}}05\text{\hspace{0.17em}}\mathrm{minutes}\left(Ra=1000000\right)\end{array}$ (8)

3. Results and Discussion

We performed our simulations by varying the Rayleigh number.

3.1. Nusselt Number

Thermal exchanges during the flow of the fluid are characterized by the number of Nusselt. It is a dimensionless number which represents the ratio between the heat flux exchanged by convection to that by conduction. The rate of heat transfer by convection in a cavity is obtained from the calculation of the Nusselt number. We are interested in heat transfer at the heated part. The Nusselt number measures the effectiveness of convection. The convective heat transfers (local and average) on the hot active part are expressed respectively by:

$Nu={\frac{\partial \theta}{\partial X}|}_{X=0}$ (9)

$\{\begin{array}{l}\stackrel{\xaf}{Nu}={\displaystyle {\int}_{s}Nu\text{d}Y}\\ \text{Withs}=A/2\text{\hspace{0.17em}}\text{isthelengthofthethermallyactivepart}\end{array}$ (10)

We note in Figure 2 that the convective heat transfer increases with the number of Ra. In the beginning, the growth is weak (conduction regime), then it becomes important when the number of Ra exceeds 1000. The differences between the local numbers of Nusselt become accentuated as the number of Ra increases. We notice that this number reaches its maximum when the first contact between the fluid and the hot wall is made (heat exchange becomes important).

Figure 2. Variation of the local Nusselt number along the active part for different Rayleigh numbers: middle active locations.

The local Nusselt numbers for various grid sizes are presented to develop an understanding of the grid fineness that is necessary for accurate numerical simulation. Hence considering the accuracy of the results required and computational time involved, an 41 × 41 grid size is chosen for (Ra = 1000, Ra = 10,000) and an 81 × 81 grid size is chosen for (Ra = 100,000, Ra = 1,000,000). There is a considerable change in the local Nusselt number from Ra = 1000 to Ra = 1,000,000. The Nusselt number profiles gradually decrease until reaching a minimum value. This is explained by the fact that the temperature difference between the fluid and the hot wall starts to drop, which gives a low heat flow. If Nu = 1, there is no convection and if the Nusselt increases, convection cooling is more and more effective. The relative variation of the local Nusselt number along the active left portion is illustrated in Figure 2 for different Rayleigh numbers.

3.2. Influence of the Variation of Time

Natural convection is a movement whose origin is a thermal imbalance. It disappears when the temperature gradients are zero. In the cavity, the flow is unsteady and therefore depends on the time. The established flow is a spatial notion while the notion of unsteady (non-permanent) flow is a temporal notion. It would be interesting to see how temperature changes over time at points on the cavity. To study the temporal behavior of the flow, we studied the temperature fluctuations at two points (A and B) located respectively on the left and right walls. These latter points are discussed below. We have chosen these points, although many more can be established, because they constantly come up in our study and, as we see it, represent the essential issues relating to heat transfer. The central parts of the vertical walls are at imposed temperatures (horizontal gradient).

Our point A is marked with a pink color (Figure 3). Figure 4 shows the temperature variations of a point A at the top of the left vertical wall as a function of time for different Rayleigh numbers. We have temperature profiles as a function of time that increases before stabilizing. The curves are upward before becoming constant, which is explained by the rise of hot air. In Figure 4(a), for a given time (t = 50 s), the temperature profiles are almost constant whereas for Figure 4(b), they are almost constant at t = 25 s. In Figure 4(c), for a given time (t = 20 s), the temperature profiles are almost constant whereas for Figure 4(d), they are almost constant at t = 15 s.

Our point B is marked with a pink color (Figure 5). Figure 6 shows the temperature variations of a point B at the bottom of the right vertical wall as a function of time for different Rayleigh numbers. We have temperature profiles as a function of time that decreases before stabilizing. The curves are descending before becoming constant, which is explained by the descent of cold air. The movement of the air in the form of a clockwise vortex is explained by the fact that we neglect the density variations in all the terms except in the Archimedes

(a) (b) (c) (d)

Figure 3. Temperature field and point A for different Rayleigh numbers: (a) Ra = 1000; (b) Ra = 10,000; (c) Ra = 100,000; (d) Ra = 1,000,000 at the end of simulation.

Principle. In Figure 6(a), for a given time (t = 50 s), the temperature profiles are almost constant whereas for Figure 6(b), they are almost constant at t = 25 s. In Figure 6(c), for a given time (t = 20 s), the temperature profiles are almost constant whereas for Figure 6(d), they are almost constant at t = 15 s. We can also point out that the temperature profiles on the upper left and lower right walls vary in opposite ways, reflecting incessant movements of cold and hot particles from the walls to the inside of the cavity. This agitation is much more prevalent if the Rayleigh number increases.

3.3. Histograms of Temperature

The histogram is a graphical presentation of the distribution of a variable. In this part we have analyzed temperature graphs over time for different Rayleigh numbers. In statistics, a histogram is a type of column chart that shows the distribution of the data. They are frequently used because their standardized format makes them easily understandable and promotes communication between users, even those unfamiliar with statistical methods. Thus, we clearly see the temperature that dominates every moment in the cavity thanks to the

(a) (b) (c) (d)

Figure 4. Temperature fluctuations at point A for different Rayleigh numbers: (a) Ra = 1000; (b) Ra = 10,000; (c) Ra = 100,000; (d) Ra = 1,000,000 over time.

(a) (b)(c) (d)

Figure 5. Temperature field and point B for different Rayleigh numbers: (a) Ra = 1000; (b) Ra = 10,000; (c) Ra = 100,000; (d) Ra = 1,000,000 at the end of simulation.

height of the peaks.

Histograms are useful for illustrating changes in data over a period or illustrating comparisons between elements. Each histogram is composed of adjoining columns of varying heights. The ordinate (y-vertical axis) receives the values and the abscissa (x-horizontal axis) determines categories. Columns are drawn from the x-axis to the value of the y-axis. Categories are grouped into classes or groups. As the cavity is heated differentially, temperature histograms are established for different Rayleigh numbers as shown in Figure 7.

-For Ra = 1000 and for Ra = 10,000: the histograms 7(a) and 7(b) show that the temperature measures have their peaks around 295 K.

-For Ra = 100,000: the histogram 7(c) shows that the temperature measures have two peaks, one around 294 K, and another higher peak around 296 K.

-For Ra = 1,000,000: the histogram 7(d) shows that the temperature measures have two peaks, one around 280 K, and another higher peak around 310 K.

3.4. The Isovorticity Contours

Some works of natural convection in a cavity are also studying vorticity [13]

(a)(b)(c)(d)

Figure 6. Temperature fluctuations at point B for different Rayleigh numbers: (a) Ra = 1000; (b) Ra = 10,000; (c) Ra = 100,000; (d) Ra = 1,000,000 over time.

(a) (b)(c) (d)

Figure 7. Temperature histograms: (a) Ra = 1000; (b) Ra = 10,000; (c) Ra = 100,000; (d) Ra = 1,000,000 at the end of simulation.

[14] . A vorticity represents a vector quantity that locally measures the rotation of the fluid (whatever its state). In mathematical terms, the vorticity is equal to the rotational velocity vector of the wind.

We have varied the temperature variation of the heated parts until a Rayleigh number belonging to the range of 1000 to 1,000,000 is obtained in a laminar regime. We solved the problem from the finite volume method by solving the system for Prandtl Pr = 0.72 (air). Let us note that the computation times for a finer mesh 81 * 81 (Ra = 100,000 and Ra = 1,000,000) are here equal to approximately 5 times those in case the mesh size is 41 * 41 (Ra = 1000 and Ra = 10,000). The computation times become increasingly longer with the increase of number of cells of the mesh since the solver has more values to be determined. The simulation is based on a modeling of the problem which required assumptions that simplify. Since the flow is unsteady, numerical simulations have been carried out for each Rayleigh number at different times. In Figure 8, we show the vorticity fields corresponding to the beginning of the simulation for each Rayleigh number.

In Figure 9, we show the vorticity fields corresponding to the end of the

(a) (b)(c) (d)

Figure 8. Contour maps of vorticity ω at the beginning of simulation: (a) Ra = 1000; (b) Ra = 10,000; (c) Ra = 100,000; (d) Ra = 1,000,000 inside the cavity.

(a) (b)(c) (d)

Figure 9. Contour maps of vorticity ω at the end of simulation: (a) Ra = 1000; (b) Ra = 10,000; (c) Ra = 100,000; (d) Ra = 1,000,000 inside the cavity.

simulation (final time). We can observe the existence of a stagnant zone in the center of the cavity, which means that the heat exchange takes place in an intense way at the corners of the cavity (Figure 9(a) and Figure 9(b)). In Figure 9(c) and Figure 9(d), for higher Rayleigh numbers, we observed the existence of two recirculation cells inside the cavity. In the course of time, we observe that the hot air is ascending for the left vertical wall and the cold air descends for the right one (Boussinesq approximation). The formation of the two cells shows that the solution has a symmetrical flow. The plot of the vorticity shows the existence of a two-dimensional centro-symmetry in the shape of the simulations. This property is described in the literature, it concerns the three fields U, V and θ and it is stated as follows:

$\left(U,V,\theta \right)\left(X,Y\right)=-\left(U,V,\theta \right)\left(1-X,1-Y\right)$ (11)

4. Conclusions

In the present study, we studied laminar natural convection in a two-dimensional confined environment differentially heated by numerical simulation and filled with a fluid. The geometric configuration of the physical model is a square cavity, subjected to a horizontal temperature gradient. Temperatures have been imposed on the central halves of the vertical walls, knowing that the left one is warmer than the one on the right. The inactive parts being adiabatic like the horizontal walls, fluidic movements have been created inside the cavity. Based on the Boussinesq approximations, a mathematical model describing the problem was developed. It was approached using a numerical approach, based on the finite volume method, using a code Fortran and a visualization tool ParaView.

The model developed makes it possible to determine the number of Nusselt, the fields, profiles and histograms of temperature, as well as the isovorticity contours in all the field studied. The analysis consisted in studying the influence of Rayleigh number variation and time too on temperature and vorticity fields. As the Rayleigh number increases, the heat transfer rate also increases. The results show the existence of a two-dimensional centro-symmetry. In this study, we studied the influence of certain parameters such as Rayleigh number on (Nusselt number and isovorticity contours) and also the effect of time on temperature profiles. From the set of numerical results we can conclude that: the Rayleigh number has a great influence on the dominant heat transfer mode in the cavity especially for the number Ra = 10,000 where the convective motion begins and for Ra = 1,000,000 where convection is dominant. So the heat exchange in the cavity increases with the increase of Ra. The determination of the flow field allowed us to analyze the behavior of the fluid inside the cavity. For a small number of Rayleigh, of the order of 1000, we have noticed the dominance of the mode of heat transfer by conduction. Beyond this value, convection dominates and appears more clearly for Ra = 1,000,000. We validated for each case our code based on the results of other authors. These results illustrate the importance of convective flow in many natural phenomena and industrial applications such as the development of materials, energy production and storage, solar collectors, double glazing for thermal insulation, etc.

Acknowledgements

We thank all the contributing members of SOLMATS group for everything.

Competing Interests

We author(s) declare that we have no competing interests.

Nomenclature

A, Side of the cavity, m

max, Maximum

Nu, Nusselt number

Pr, Prandtl number

Ra, Rayleigh number

s, Length of active parts, m

t, Dimensional time, s

T, Dimensional temperature, K

Td, Dimensional temperature of right active part, K

Tg, Dimensional temperature of left active part, K

Tr, Reference temperature, K

u, Dimensional horizontal velocity, m/s

U, Dimensionless horizontal velocity

v, Dimensional vertical velocity, m/s

V, Dimensionless vertical velocity

x, Dimensional abscissa, m

X, Dimensionless abscissa

y, Dimensional ordinate, m

Y, Dimensionless ordinate

Greek Symbols

α, Thermal diffusivity, m^{2}/s

β, Coefficient of thermal expansion, 1/K

ζ, Dimensionless vorticity

θ, Dimensionless temperature

ρ, Density of fluid, kg/m^{3} τ, Dimensionless time Ψ, Dimensionless stream function ψ, Dimensional stream function, m^{2}/s ω, Dimensional vorticity, s⁻¹

Subscript

i, j, Indices

local, Local value

max, Maximum value

min, Minimum value

0, Reference.

Superscript

−, Average value

n, Any time level.

Conflicts of Interest

The authors declare no conflicts of interest.

[1] | Patankar, S.V. (1980) Numerical Heat Transfer and Fluid Flow. Hemisphere Publishing Corporation, Taylor and Francis Group, New York. |

[2] |
Elder, J.W. (1965) Laminar Free Convection in a Vertical Slot. Journal of Fluid Mechanics, 23, 77-98. https://doi.org/10.1017/S0022112065001246 |

[3] |
Ismail, K.A.R. and Scalon V.L. (2000) A Finite Element Free Convection Model for the Side Wall Heated Cavity. International Journal of Heat and Mass Transfer, 43, 1373-1389. https://doi.org/10.1016/S0017-9310(99)00225-2 |

[4] |
Kandaswamy, P., Sivasankaran, S. and Nithyadevi, N. (2007) Buoyancy-Driven Convection of Water near Its Density Maximum with Partially Active Vertical Walls. International Journal of Heat and Mass Transfer, 50, 942-948.
https://doi.org/10.1016/j.ijheatmasstransfer.2006.08.013 |

[5] |
Nithyadevi, N., Kandaswamy, P. and Sivasankaran, S. (2006) Natural Convection in a Square Cavity with Partially Active Vertical Walls: Time-Periodic Boundary Condition. Mathematical Problems in Engineering, 2006, 1-16.
https://doi.org/10.1155/MPE/2006/23425 |

[6] |
Kandaswamy, P., Nithyadevi, N. and Ng, C.O. (2008) Natural Convection in Enclosures with Partially Thermally Active Side Walls Containing Internal Heat Sources. Physics of Fluids, 20, 1063-1834. https://doi.org/10.1063/1.2981834 |

[7] |
Harlow, F.H. and Welch, J.E. (1965) Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface. Physics of Fluids, 8, 2182-2189. https://doi.org/10.1063/1.1761178 |

[8] |
Lord Rayleigh, O.M.F.R.S. (1916) On Convection Currents in a Horizontal Layer of Fluid, When the Higher Temperature Is on the Under Side. Philosophical Magazine Series 6, 32, 529-546. https://doi.org/10.1080/14786441608635602 |

[9] | Winters, K.H. (1980) A Numerical Study of Natural Convection in a Square Cavity. United Kingdom Atomic Energy Authority, AERE-R9747, Harwell. |

[10] |
Bae, J.H. and Hyun, J.M. (2004) Time-Dependent Buoyant Convection in an Enclosure with Discrete Heat Sources. International Journal of Thermal Sciences, 43, 3-11. https://doi.org/10.1016/S1290-0729(03)00102-9 |

[11] |
De Vahl Davis, G. (1983) Natural Convection of Air in a Square Cavity: A Bench Mark Numerical Solution. International Journal for Numerical Methods in Fluids, 3, 249-264. https://doi.org/10.1002/fld.1650030305 |

[12] | Gresho, P.M., Upson, C.D. and Lee, R.L. (1980) Finite Element Simulations of Thermally Induced Convection in an Enclosed Cavity. Lawrence Livermore Laboratory Report UCID-18602, Livermore. |

[13] |
Lo, D.C., Young, D.L., Murugesan, K., Tsai, C.C. and Gou, M.H. (2007) Velocity-Vorticity Formulation for 3D Natural Convection in an Inclined Cavity by DQ Method. International Journal of Heat and Mass Transfer, 50, 479-491.
https://doi.org/10.1016/j.ijheatmasstransfer.2006.07.025 |

[14] |
Lo, D.C., Young, D.L. and Murugesan, K. (2005) GDQ Method for Natural Convection in a Cubic Cavity Using Velocity-Vorticity Formulation. Numerical Heat Transfer, Part B, 48, 363-386. https://doi.org/10.1080/10407790591009251 |

Journals Menu

Contact us

+1 323-425-8868 | |

customer@scirp.org | |

+86 18163351462(WhatsApp) | |

1655362766 | |

Paper Publishing WeChat |

Copyright © 2024 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.