Finite Element Investigation of Load Acting on the Hotspot Detector Located inside the Silo Caused by Material Discharge

Abstract

Spontaneous ignition caused by material discharge inside a silo causes considerable economic damage. To prevent this, we developed a silo hotspot detector that can be installed inside the silo to monitor the temperature according to the depth of the silo. However, if the silo hotspot detector located inside the silo is destroyed because of the pressure and load generated during material discharge, it could lead to a larger accident. Therefore, the structural safety of the silo hotspot detector should be evaluated based on material discharge; currently, there is no particular method to achieve this. Therefore, in this study, the theoretical formula is obtained through Eurocode, and the pressure and tensile force acting on the silo hotspot detector are predicted through the finite element method (FEM) using the Coupled Eulerian–Lagrangian(CEL) method. These result were verified by comparing the load measurement data acting on the silo hotspot detector when the silo material was discharged. It was confirmed that simulation using the CEL method can sufficiently simulate the behavior of the silo according to material discharge. Additionally, we confirmed that the structural safety of the silo hotspot detector inside the silo can be evaluated through FEM.

1. Introduction

In plants using coal or biomass fuels, there is always a risk of fire due to the spontaneous combustion of silo [1]. As a countermeasure against this, IR cameras and gas detectors are used to prepare for fire in order to detect the natural ignition of the silo in advance [1]. However, since the start of spontaneous ignition starts with the increase in the temperature of the center of the silo, it is difficult to detect the starting point of spontaneous ignition by measuring the temperature of the silo surface through an IR camera or measuring gas leakage through a gas detector. The silo hotspot detector, which can measure the temperature of each level, developed for this purpose, is installed inside the silo and monitors the internal temperature in real time to prevent fires caused by spontaneous ignition. However, the silo internal material is frequently charged and discharged, and the pressure and tensile force generated therefrom also frequently changed. If the silo hotspot detector is destroyed by this, it will cause a big accident, so it is necessary to check the possibility of the silo hotspot detector destruction due to the discharge of material inside the silo.

The pressure and load on an internal structure caused by material discharge is governed by the flow of material. The flow of material inside silos is a phenomenon that occurs in many industrial applications; this phenomenon has unique flow and stress characteristics, and it must be considered when designing a device that meets operational and safety requirements. The study of the flow of material has been an active topic of study over the past century [2,3,4,5].

Jansen [6] and Jenike [7,8] described the general flow of material inside a silo. However, the method of accurately evaluating the stress field and material discharge rate occurring inside the silo remains a popular research topic. This is because the theory described by Jansen and Jenike is limited to silos that have a shape where the silo outlet and the center line of the hopper coincide. However, silos that do not actually meet this condition are also used in many industrial fields [9]. Even in silos with specific conditions, accurate evaluation may not be possible based on the scenario. Ding et al. [10] suggested that, as assumed by the pressure theory [6,11,12], the theory is applicable only to steep hoppers, as the wall friction effect of these hoppers develops completely during outflow. For shallow hoppers in which the wall friction develops only partially, significant modifications to the theory are required.

Discrete and continuum numerical methods are commonly used to model and simulate the flow of materials. The discrete element method (DEM) achieves particle-scale simulation and simulates the motion of each individual particle under certain conditions. Furthermore, it can provide micromechanical information such as trajectories and transients acting on individual particles. Therefore, the DEM has been increasingly used to study phenomena and behaviors occurring in hoppers such as mass flow rate [13], particle flow pattern [14], and internal bulk stress [15,16]. However, DEM simulation is computationally demanding and can handle only small systems. Methods such as multidomain continuum and other discrete methods [17,18,19] have been proposed to solve these problems; however, there remains difficulties in using DEM simulation with large systems.

In the continuum type method, the dynamics of the material are governed by the momentum and energy conservation equations, and the flow is modeled on a macroscopic or global scale, with a construction model describing the unique properties of the material. The finite element method (FEM)—a continuum approach—has been used in hopper studies. Thus, hopper internal stress and wall pressure have been predicted satisfactorily. However, to simulate the flow conditions of a material, specific techniques or treatments such as remeshing and rezoning plans [20] and adaptive meshing techniques [21] are required.

In a recent study [22], the dynamic behavior of materials including the interaction between the material flow and a rigid barrier was investigated through the coupled Eulerian–Lagrangian (CEL) technique. This analysis method provides reliable information about material flow, and it can provide reliable information for estimating the impact force caused by the material. Furthermore, the method demonstrated the feasibility of CEL technology in actual engineering applications.

This study aims to evaluate the structural safety of silo hotspot detectors under material discharge in silos. To this end, the discharge process was simulated using the CEL technique, and the results were compared with theoretical equations and experimental measurements. By comparing the results, we confirmed that the CEL method is suitable for the discharge process. Subsequently, the structural effect of the silo hotspot detector due to the release of material inside the silo was analyzed through the CEL method, and it is confirmed that FEM is effective as a method for evaluating the structural safety of the silo hotspot detector.

The rest of this paper is organized as follows. First, the theoretical value is obtained in Section 2; the calculation method and the constitutive relationship and boundary conditions used in the simulation are explained in Section 3; the simulation results and experimental results are presented in Section 4; finally, the conclusions are discussed in Section 5.

2. Load on Silo Hotspot Detector

The pressure acting on the silo hotspot detector because of the material discharge was assessed to be equal to the horizontal pressure acting on the silo wall via a force equilibrium relationship, and a prediction equation was obtained. For the process of obtaining the prediction equation, refer to Eurocode 1 [23]. The shape, cross-sectional shape, pressure, and coordinate system of the target silo are shown in Figure 1. The capacity of the silo was 900 tons and the type of material was coal.

Table 1. Material properties.

Material	Unit Weight $(\mathit{\gamma }, \mathit{k}\mathit{N}/{\mathit{m}}^3)$		Internal Friction Degree $\left({\mathit{\varphi }}_{\mathit{i}}, \mathit{d}\mathit{e}\mathit{g}\mathit{r}\mathit{e}\mathit{e}\right)$		Lateral Pressure Ratio $\left(\mathit{K}\right)$		Wall Friction $\left(\mathit{\mu }\right)$
Coal	${\gamma }_l$	${\gamma }_u$	${\varphi }_{im}$	$a_{\varphi }$	$K_m$	$a_K$	$D2$	$a_{\mu }$
	Lower	Upper	Mean	Factor	Mean	Factor	Mean	Factor
	7.0	10.0	31	1.16	0.52	1.15	0.49	1.12

summarizes the values of the physical properties of the material and friction coefficients according to the wall classification of the target silos, and

Table 2. Characteristics of target silo and silo hotspot detector.

Silo	Vertical Wall Height $\left({\mathit{h}}_{\mathit{c}}, \mathit{m}\right)$	Diameter $\left({\mathit{d}}_{\mathit{c}}, \mathit{m}\right)$	Hopper Half Apex $\left(\mathit{\beta }, \mathit{d}\mathit{e}\mathit{g}\mathit{r}\mathit{e}\mathit{e}\right)$	Hopper Height $\left({\mathit{h}}_{\mathit{h}}, \mathit{m}\right)$	Capacity $\left(\mathit{t}\right)$	Classification of Silo
	13.015	8.7	20.649	10.5	900	AAC 2
Detector	Diameter $\left(m\right)$			Tensile strength $\left(kN\right)$
	0.035			173

summarizes the classification values according to the dimensions and capacity of the target silos and developed products. The property values in Table 1 and Table 2 were set with reference to Eurocode [23]. In addition, the tensile strength of the silo hotspot detector shown in Table 2 cited the value provided by the manufacturer.

2.1. Normal Pressures on Vertical Wall

To find the tensile force applied to the product in the vertical wall section, pressure $P_{hf}$ applied to the vertical wall (Figure 1) must be obtained. $P_{hf}$ can be obtained as

$$P_{hf}\left(z\right)=P_{ho}{Y}_J\left(z\right)$$

(1)

where $z$ is the depth from the load equivalent surface, $Y_J\left(z\right)$ is the Janssen function, and $P_{ho}$ is the pressure when $z=z_o$. Furthermore, $z_o$ indicates the specific depth at which the pressure becomes stable. To determine the maximum pressure generated, it is necessary to assume that the silo is full and extreme. In addition, in order to consider maximum normal pressures on the vertical wall with reference to Eurocode, a friction coefficient $\mu$ is the minimum value (lower limit value), a lateral pressure ratio $K$ is the maximum value (upper value), and an internal friction angle angle_${\varphi }_i$ with a minimum value (lower limit value).

• Lower characteristic value of $\mu =\frac{{\mu }_m}{a_{\mu }}$
• Upper characteristic value of $K=a_K{K}_m$
• Lower characteristic value of ${\varphi }_i=\frac{{\varphi }_{im}}{a_{\varphi }}$

Through the above values, the pressure $P_{hf}\left(z\right)$ applied normally to the vertical wall of the target silo according to depth $z$ can be obtained. However, this study aimed to find the tensile force generated in the product by the outflow, and therefore, the outflow coefficient $C_h$ for the horizontal pressure was considered. The discharge coefficient $C_h$ is a factor that considers many phenomena that occur when the silo load is discharged, and it is determined by classification based on the capacity of the silo. It can be confirmed that the target silo is classified as AAC 2 [23] and has a value of a discharge factor for a horizontal pressure $C_h=1.15$. Therefore, the pressure acting perpendicular to the silo wall by the outflow is given by

$$P_{hf}\left(z\right)\times C_h=1.15P_{hf}\left(z\right)$$

(2)

This pressure acts equally on the product. Therefore, to determine the tensile force applied to the detector, $1.15 P_{hf}\left(z\right)$ is integrated by the outflow against depth $z$ and the coefficient of friction is multiplied by the circumference of the product. The integral value for $z$ is

$${\int }_0^{z}1.15 P_{hf}\left(z\right)dz=1.15P_{ho}\left[z-z_o{Y}_J\left(z\right)\right].$$

(3)

Multiplying this value by the coefficient of friction $\mu$ and the circumference of the detector gives the tensile force acting on the detector by the silo vertical wall section.

2.2. Normal Pressures on Hopper Wall

Eurocode divides the slope of the hopper into “flat bottoms”, “steep hoppers”, and “shallow hoppers” to calculate the pressure acting on the hopper wall [23]. The target silo satisfies, and therefore, the slope of the hopper is classified as “steep”.

$$\therefore \tan \beta <\frac{1-K}{2{\mu }_h}$$

(4)

To find the tensile force applied to the product in the hopper section, pressure $P_{ne}$ (Figure 1) applied perpendicularly to the hopper wall must be obtained. The formula for $P_{ne}$ is

$$P_{ne}=F_e{P}_v$$

(5)

where $F_e$ is the hopper pressure ratio according to the discharge, and $P_v$ is the vertical stress acting on the load. As in the vertical wall section, to determine the maximum pressure acting on the hopper wall, it is necessary to assume that the silo is full and extreme. In addition, a minimum value (lower limit value) for the friction coefficient $\mu$ and a maximum value (upper value) for the internal friction angle ${\varphi }_i$ were applied to consider maximum hopper pressures on discharge with reference to Eurocode.

• Lower characteristic value of $\mu =\frac{{\mu }_m}{a_{\mu }}$
• Upper characteristic value of ${\varphi }_i=a_{\varphi }{\varphi }_{im}$

To determine $P_{ne}$ applied to the hopper wall vertically, $F_e$ and $P_v$ must be obtained. $P_v$ is given by

$$P_v=\left(\frac{\gamma h_h}{n-1}\right)\left\{\left(\frac{x}{h_h}\right)-{\left(\frac{x}{h_h}\right)}^n\right\}+P_{vft}{\left(\frac{x}{h_h}\right)}^n$$

(6)

where x is the height from the bottom of the hopper. The vertical pressure $P_{vft}$ (Figure 1) required to calculate $P_v$ is given by

$$P_{vft}=C_b{P}_{vf}$$

(7)

where $C_b$ is a load expansion coefficient that describes a phenomenon wherein the load in the vertical wall section is transferred to the bottom of the hopper or silo. The target silo is classified as AAC 2, and therefore, $C_b=1.0$. $P_{vf}$ is the vertical pressure obtained by the load in the vertical wall section, and it can be obtained by

$$P_{vf}\left(z\right)=\frac{P_{ho}}{K}{Y}_J\left(z\right)$$

(8)

where $P_{vf}$ is the pressure generated in the vertical wall section. In order to obtain the maximum value of $P_{vf}$, the maximum vertical load on hopper situation was considered with reference to Eurocode. These conditions were: friction coefficient $\mu$ is the minimum value; lateral pressure ratio $K$ is the minimum value; internal friction angle ${\varphi }_i$ is the maximum value. These values were only used to obtain the value of $P_{vf}$.

• Lower characteristic value of $\mu =\frac{{\mu }_m}{a_{\mu }}$
• Lower characteristic value of $K=\frac{K_m}{a_K}$
• Upper characteristic value of ${\varphi }_i=a_{\varphi }{\varphi }_{im}$

Next, we needed to find the exponent $n$ value required to obtain the $P_v$ value. This can be obtained using

$$n=S\left(F{\mu }_{heff}\cot \beta +F\right)-2$$

(9)

where $S$ is the hopper shape factor, ${\mu }_{heff}$ is the effective friction coefficient, and F is the hopper pressure ratio. The value of $S$ varies depending on the shape of the hopper. Because the hopper shape of the target structure is conical, it can be confirmed that $S=2$. The expressions for ${\mu }_{heff}$ and F are different based on the slope classification of the hopper. Furthermore, because the hopper slope of the target silo is classified as “steep,” it can be confirmed that ${\mu }_{heff}={\mu }_h$ and $F=F_e$, where ${\mu }_h$ denotes the hopper wall friction coefficient minimum value (Lower).

To determine $F_e$, the friction angle ${\varphi }_{wh}$ of the hopper wall must first be obtained. This can be obtained using

$${\varphi }_{wh} = {\tan }^{-1}{\mu }_h$$

(10)

Substituting the values, $F_e$ can be obtained as

$$F=F_e=\frac{1+\sin {\varphi }_i\cos {\epsilon }_h}{1-\sin {\varphi }_i\cos \left(2\beta +{\epsilon }_h\right)}$$

(11)

In addition, the ${\epsilon }_h$ value required to obtain the $F_e$ value can be obtained using

$${\epsilon }_h = {\varphi }_{wh}+{\sin }^{-1}\left\{\frac{\sin {\varphi }_{wh}}{\sin {\varphi }_i}\right\}$$

(12)

Thus, the index n value can be obtained, and then, it can be substituted into Equation (6) to obtain the $P_v$ value. However, because $P_{ne}$ obtained above is not the vertical pressure applied to the product, $P_{neh}$, which is the horizontal component of $P_{ne}$, must be obtained as shown in Figure 2.

$P_{neh}$ is given as

$$P_{neh}=\cos \beta P_{ne}.$$

(13)

This value is the pressure $P_{neh}\left(x\right)$ applied vertically to the product according to $x$. To obtain the tensile force applied to the product in the silo hopper section, this pressure must be integrated over $\mathrm{x}$ and multiplied by the coefficient of friction and circumference of the detector.

$${{\displaystyle \int }}_0^x{P}_{neh}dx\times \mu \times 2\pi r$$

(14)

The pressure and load acting on the entire silo section and detector are shown in Figure 3.

2.3. Experimental Studies

An experiment was performed to measure the tensile load acting on the silo hotspot detector. The installation location of the detector was in the center of the inside of the silo as shown in Figure 4a. The specification of the target silo is as shown as Figure 1. The silo hotspot detector was installed through the flange located at the top of the silo to place it in the center of the inside of the silo; the silo hotspot detector and tensile load measuring device were installed using a chain block. The self-weight of the silo hotspot detector was neglected to measure only the tensile load caused by the load spill. To achieve this, the scale was set to zero after the installation was completed. The state in which the experimental setup was completed is shown in Figure 4b; the load was measured using a balance meter during material discharge after the setup. As a load measuring device, a suspension-type scale with a maximum capacity of 5 tons with a model name of 5THB of CAS scale Korea was used.

The experiment was conducted six times, and the results are shown in Figure 5. The third experiment was excluded because the emission per hour was small, unlike the other cases. The tensile load in the case of a storage rate of 55–60% was inferred using the measured value, and the results were compared with results obtained in Section 2.1 and Section 2.2, as shown in Figure 6.

Figure 6 indicates that the two results are quite similar. The experimental value at the storage ratio of 83% was different from the theoretical value because the discharge started at this value, and it was not the state at which the discharge was released.

3. Discharge Simulation

3.1. Constitutive Model

The accuracy of the analysis depends on the constitutive model used to describe the behavior of the target material in the continuum approach. Therefore, selecting an appropriate constitutive model is the first step in interpreting the continuum approach. Several models for determining the composition of materials have been proposed: classical plastic theory [24,25,26], pole elastoplastic theory [27], and others. However, this analysis is considered to be suitable for expressing the properties of granular materials in the classical elastic–plastic model because micromechanics such as shear banding are not as important as bulk scale [28,29].

Thus, in this study, the mechanical behavior of the discharge of the material stored in the silo was simulated using the traditional Mohr–Coulomb elastoplastic model. This model comprises a linear isotropic elastic law that defines material behavior at small loads, a yield criterion that determines the transition of material to yield, and an unrelated plastic flow potential that determines the direction of flow after material yield. The exact description of the Mohr–Coulomb elastoplastic model can be found in Abaqus 2018 [30], and it is briefly described below.

The elastic behavior is determined by Young’s modulus E and Poisson’s ratio ν in the law of isotropic elasticity as

$${\sigma }_{ij} = D_{ijkl}^{el}{\epsilon }_{kl}^{el}$$

(15)

where ${\sigma }_{ij}$ is the total stress, ${\epsilon }_{kl}^{el}$ is the elastic strain, and $D_{ijkl}^{el}$ is the fourth-order elastic tensor.

The yield condition is

$$R_{mc}q-p\tan \varphi -c = 0$$

(16)

where

$$R_{mc}=\frac{1}{\sqrt{3}\cos \phi }\sin \left(\mathsf{\Theta }+\frac{\pi }{3}\right)+\frac{1}{3}\cos \left(+\frac{\pi }{3}\right)\tan \phi$$

(17)

Furthermore, $p=-\frac{1}{3}trace\left({\sigma }_{ij}\right)$ is the first invariant value of the stress representing the equivalent pressure, $\mathrm{q}=\sqrt{\frac{3}{2}\left(S_{ij}{S}_{ij}\right)}$ is the Mises equivalent stress, and $S_{ij}$ is the deviatoric stress tensor; $\mathsf{\phi }$ and $\mathrm{c}$ are the internal friction angle and the cohesive force of the material, respectively; $\mathsf{\Theta }$ is the deflection angle defined by $\cos \left(3\mathsf{\Theta }\right)={\left(r/q\right)}^3$, where $r$ is the third invariant of the bias stress that is defined as $r={\left(\frac{9}{2}\left(S_{ij}{S}_{jk}{S}_{ki}\right)\right)}^{1/3}$. The flow potential $\mathrm{G}$ is selected based on the hyperbolic function in the meridian stress plane and the smooth elliptic function in the deflection stress plane. The mathematical expression is

$$G=\sqrt{{\left({ϵc|}_0\tan \psi \right)}^2+{\left(R_{mw}q\right)}^2}-p\tan \psi$$

(18)

where

$$R_{mw}=\frac{4\left(1-e^2\right){\cos }^2\mathsf{\Theta }+{\left(2e-1\right)}^2}{2\left(1-e^2\right)\cos \mathsf{\Theta }+\left(2e-1\right)\sqrt{4\left(1-e^2\right){\cos }^2\mathsf{\Theta }+5e^2-4e}}\frac{3-\sin \phi }{6\cos \phi }$$

(19)

where $\psi$ is the dilatancy angle of the material, ${c|}_0$ is the initial cohesive yield stress, and $ϵ$ is a variable that characterizes eccentricity. In the analysis in this study, a fixed value of $\mathsf{\epsilon }=0.1$ was used. The $e$ value was determined by the internal friction angle as $\mathrm{e}=\left(3-\sin \phi \right)/\left(3+\sin \phi \right)$.

The volumetric expansion of the material can be controlled by the expansion angle $\psi$. Since the commonly used material has small expandability, the $\psi$ value was set to 5° in this study. Furthermore, because it was assumed that there was no cohesive force, the cohesive force $\mathrm{c}=0$. Other material property values are listed in

Table 3. Parameters used in finite element method (FEM) simulation.

Parameters	Basic Values
Silo height, H (m)	23.515
Silo diameter, dc (m)	8.7
Orifice diameter, d0 (m)	0.786
Silo hotspot detector diameter, dt (m)	0.035
Half angle of hopper, β (deg)	20.65
Bulk density, ρ (kg/m3)	1019
Young’s modulus, E (MPa)	16
Poisson ratio, ν	0.3
Internal friction angle, φ (deg)	35.96
Dilatancy angle, ψ (deg)	2.5
Cohesion, c (Pa)	0
Wall friction, μ	0.49

.

3.2. Numerical Analysis of Discharge Process Using Coupled Eulerian–Lagrangian (CEL) Method

From the numerical analysis point of view, there are several Lagrangian analysis and Eulerian analysis methods to simulate the behavior of an element with a volume over time. The Lagrangian analysis method defines the behavior of the continuum as a function of the physical coordinates and the time of each element. Furthermore, because the nodes of each mesh move with the properties of the elements, the nodes must be accurately shared on the contact surface between the two objects (Figure 7a).

The Eulerian analysis method defines the behavior of a continuum as a function of spatial coordinates and time. Here, the behavior of an object is defined as the amount of movement of a substance within a specific region (Figure 7b). Therefore, there is an advantage in that the distortion or warping of the mesh does not occur during Eulerian analysis. Therefore, in this study, the silo discharge process was analyzed using the CEL analysis method.

Figure 8 shows the calculation procedure of the CEL analysis method. In the CEL analysis method, elements that generate a large strain are processed by the Eulerian technique, and the remaining elements are processed by the Lagrangian technique. As shown in Figure 8, the volume fraction tool created a discrete field based on the Eulerian part. The Eulerian volume fraction (EVF) was calculated between the Eulerian part and the reference part. The EVF had a value between “0” and “1,” which represented the percentage of the area occupied by the specified material. Materials were thus assigned to predefined fields, and material boundaries were determined based on the EVF of all elements.

Areas that were not defined in the Eulerian part were designated as voids. These voids had no material properties at first; however, as the analysis proceeded further, the material flowed or passed through the void region. During each calculation step of the CEL method, the calculation results were stored in each node and mapped to the computational mesh (step 2, Figure 8b). Subsequently, unknown variables were calculated using the Lagrangian formula; the velocity, acceleration, and displacement of each node were updated using these values (step 3, Figure 8b), and the stress and strain were calculated by the applied structural model. After mesh adjustment, the position of the node was updated (step 4, Figure 8b). To solve the Eulerian step, the deformed mesh was remapped to the original fixed mesh, and the material volume between adjacent elements was calculated by the EVF (Figure 8c). Finally, a Lagrangian step process was performed to explain phenomena such as energy, stress, and strain between adjacent elements, as shown in Figure 8b.

3.3. CEL Modeling

The structure of the silo to be analyzed is shown in Figure 9a. The vertical wall section consisted of a cylindrical bunker with diameter d_c, and the hopper section was a conical hopper with hopper angle β and diameter d₀ orifice consisting of a total height H. The values for these geometric variables are listed in Table 3. The shape of the target silo was axially symmetric, and only 1/4 of the silo was modeled for efficiency in the numerical analysis. Because only 1/4 was modeled, plane symmetry conditions were applied to each side as shown in Figure 9b. In fact, it is noted that the load flow caused by silo material discharge may be asymmetric because of some transient influence, asymmetric arrangement, or operation. In this case, the above symmetric processing was invalid. However, in the case of a general cylindrical silos and conical hoppers such as those considered in this study, the computational time is effectively reduced by modeling only 1/4 of the target silos because a symmetric model is applicable.

The silo discharge process simulation was performed using the CEL analysis method described in Section 3.2 using the commercial package Abaqus 2018. The silo wall and silo hotspot detector were modeled as rigid shells (element type R3D4), assuming that both the thickness and sag were negligible under the lateral pressure of the inner material. Translation and rotation were completely limited throughout the simulation. The material inside the silo was modeled using a 3D Eulerian element (EC3D8R) integrated with the basic hourglass control. The number of mesh was 8798.

The procedures and boundary conditions used in the analytical modeling in this study are shown in Figure 8. To facilitate the identification of results such as stresses and velocities occurring near the silo wall, the computational mesh was constructed in the same shape as the silo and slightly increased in width to surround the silo wall. Figure 9b shows the mesh used in the simulation. For silo discharge analysis (or case) with different specifications, the FEM mesh should be changed accordingly. In the basic environment of the CEL method, the Eulerian element is free of materials. In other words, the element is occupied by “voids” (there are no physical properties in the void, but they are treated as materials in the Eulerian FEM system). Therefore, before starting the simulation, materials must be assigned to these silo elements as initial conditions, which is similar to the actual silo material filling process.

As applied boundary conditions (Figure 9b), two Eulerian boundary conditions were created at the top and bottom of the structure so that the material could flow freely in and out of the computational space. As the analysis proceeded, the material flowed out because the discharge slowly replaced the voids from the top. The substances exiting the outlet were not considered in subsequent calculation steps because they were outside the calculation area.

In summary, the silo discharge process simulating consists of four main steps: (i) creating silo walls, a silo hotspot detector, and Eulerian meshes using desired parameters; (ii) feeding material to the expected levels of the mesh; (iii) applying boundary conditions and gravity as described; finally, (iv) running the simulation.

3.4. Verification

To verify the finite element model, the results of the pressure measurement experiments on the silo walls performed by Couto et al. [32] were compared with the numerical results. The analytical model for verification was modeled by applying the same conditions as those applied to the actual experiment, except for one condition. In the actual experiments, because of the presence of humidity, there was a certain level of cohesive interaction between the material particles. However, because the cohesion value c was not specifically provided in the experiment, it was assumed to be “0” in the FE model. Other relevant parameters used the values as listed in

Table 4. Parameters used in the FEM simulation.

Parameters	Basic Values
Silo height, H (m)	2.48
Silo diameter, dc (m)	1.0
Orifice diameter, d0 (m)	0.35
Half angle of hopper, α (deg)	34.3
Bulk density, ρ (kg/m3)	856
Young’s modulus, E (MPa)	10
Poisson ratio, ν	0.3
Internal friction angle, φ (deg)	30
Dilatancy angle, ψ (deg)	5
Cohesion, c (Pa)	0
Wall friction, μ	0.3

. The specifications and analytical models of the silo to be verified are shown in Figure 10. The number of elements was 15,725.

Figure 11a shows the pattern of internal stress in the silo under material discharge obtained using the numerical analysis. Figure 11b shows the experimental values and numerical analysis values graphically. As a result of comparing the results of the FEM with the experimental values, it can be seen that the results showed quite similar results in the transition section where the maximum pressure occurred. However, a difference occurred in some sections of the vertical wall. Thus, it can be said that the numerical analysis for evaluating the pressure generated by the silo material discharge is reasonable.

4. Result and Discussion

4.1. Silo Wall Pressure

Pressure on the silo wall caused by material discharge is a subject studied intensively using the first algebraic equation proposed by Janssen [6]. Janssen proposed an equation to determine the wall pressure by dividing the cylindrical silo into several elements and by analyzing the equilibrium conditions of each element. The characteristic of this equation is that the frictional resistance of the wall increases with the depth of the silo, and thus, the wall pressure is “saturated” at a certain depth. Walker [11] improved this equation to determine the pressure applied to the hopper by applying it to a conical hopper, and Walters [12] considered the hopper shape factor for the rapid stress change between the silo vertical wall section and the hopper section. Nedderman [3] reviewed the equation for obtaining this hopper stress and calculated the pressure distribution in the silo.

Figure 12a shows the silo internal stress pattern (horizontal component σ_xx) obtained via the silo material discharge simulation. The shape of the internal stress distribution varied depending on this section. In the vertical wall section, the frictional force generated on the wall could be balanced with gravity, and therefore, the internal stress of the silo increased with the depth of the material. In the hopper section, because the frictional force generated on the inclined wall balanced gravity with a certain level, the internal stress was gradually released toward the discharge port.

The maximum stress value occurred in the transition region from the vertical wall section to the hopper section. In the theoretical equation, the transition zone is simply regarded as a plane, which is not as realistic as the “arch” in FE modeling, which caused a pressure difference at a silo depth of 9–17 m, as shown in Figure 12b.

The pressure acting on the silo wall can be observed as the horizontal component of the stress (σ_xx) in the vertical wall section (variable SVAVG11 of Abaqus output). In the hopper section, it must be obtained through the transformation of σ_xx, σ_yy (SVAVG22), and σ_xy (SVAVG12) because of the slope of the hopper wall. Figure 12b is a graph plotting the pressure acting on the silo wall obtained from the FEM simulation when a flow field was formed inside the silo. As shown in Figure 12b, the results of the theoretical equation and the results of the FEM simulation are almost identical; however, some differences were observed in some sections.

There is another difference between the theoretical and FEM simulation results. Unlike in the theoretical approach, the FEM results shown in Figure 12 considered the nonuniformity of the frictional force depending on the wall position. The wall pressure graph shows the vibration phenomenon, as in Figure 12b. In addition, the FEM simulation can be used to check the silo internal velocity and stress results in more detail. The silo wall pressure is affected by internal particle flow patterns and velocities [15,33,34]. In theory, it is difficult to fully address this phenomenon; however, it can be explained more easily with the FEM simulation.

4.2. Pressure Acting on a Silo Hotspot Detector

If a combustible object such as coal is stored inside the silo, an unexpected fire may occur, and this can lead to considerable damage; therefore, it is necessary to monitor the temperature inside the silo with the silo hotspot detector. However, if the silo hotspot detector is destroyed by the pressure and load generated when the material inside the silo discharges, another accident can occur. Therefore, it is necessary to measure the load generated in the silo hotspot detector under material discharge. In this study, the generated load was measured through experiments, obtained through an analytical solution, and confirmed through numerical simulation.

Figure 13 shows the pressure acting on the detector because of the material discharge. The pressure profile in Figure 13 is the result of the analysis at 1, 3, and 5 s, the time after which all particles start flowing. Although it did not coincide with the pressure acting on the wall, the pressure pattern acting according to the section was the same. In the vertical wall section, the pressure increased with depth, and in the hopper section, it gradually decreased toward the outlet. As shown in Figure 13, the difference between the equation result and the FEM remained in the entire interval. In addition, there was no transition section unlike the silo wall, and it can be confirmed that it was different from the pressure applied to the silo wall.

The magnitude of the tensile force acting on the silo hotspot detector is summarized in

Table 5. Tensile force applied to the detector.

No.	Range of Depth (m)	Tensile Force by FEM (kN)	Tensile Force by Predictive Formula (kN)
1	0–23.515	38	37.51
3		37.3
5		37

. Although the types of pressures predicted by the FEM and the theoretical formulas are different, the magnitudes have almost similar values. Figure 14 shows the distribution of the internal stress caused by the outflow of load when the storage rate of the internal load of the silo was 83.8%, 71.7%, 62.2%, or 48.8%. The location and size where the maximum stress occurred varied depending on the storage capacity. The number of elements was 7778 at 83.8%, 6928 at 71.7%, 6333 at 62.2%, and 5398 at 48.8%.

In Figure 15, the tensile force generated in the silo hotspot detector for storage capacities of 83.8%, 71.7%, 62.2%, and 48.8% were compared with the theoretical formula, experimental values, and numerical analysis results. Figure 15 indicates that the numerical analysis and theoretical results are similar. In addition, it can be confirmed that the numerical analysis results and experimental results generally coincide; however, there are sections with differences. The experimental results at the storage capacity of 85% are different from the numerical analysis results, similar to the theoretical results. Since this was the value at the start of the experiment, the results of the static state would have been measured, which caused the above difference.

Based on the above results, it can stated that the numerical analysis method used in this paper is suitable for silo discharge process analysis and can be used to perform structural safety evaluation of silo hotspot detectors.

5. Conclusions

In this study, a method for evaluating the pressure and tensile force generated in the silo hotspot detector by the material discharge inside the silo by numerical analysis was proposed. Before using the FEM, the pressure and tensile force generated through the Eurocode were predicted. The CEL method was applied to simulate the release process, and the material was treated with an elasto plastic solid described by the Mohr–Coulomb model. The proposed FEM was verified by comparison with the results measured through experimental methods. The main results of this study are summarized as follows.

(1) The CEL approach applied in this study can be applied to the overall process modeling of silo discharge. Unlike Lagrangian analysis, simulation is possible without the distortion and warping of the mesh. As such, the CEL approach proposed in this study can be applied to engineering fields related to silo discharge.

(2) The pressure acting on the silo hotspot detector tends to be different from the pressure acting on the silo wall. Because of the flow pattern, the pressure applied to the silo wall tends to be different. Nevertheless, it can be seen that the magnitude of the applied tensile load does not differ significantly.

(3) The actual behavior of silo hotspot detectors can be displaced from the center depending on the flow pattern. Therefore, there is a difference between the experimental value and FEM and theoretical values. More accurate safety measures are expected using the results of the analysis.

(4) Currently, there is no other method than the experimental method proposed herein to confirm the structural safety of the silo hotspot detector installed inside the silo. It is hoped that the proposed method can replace the experimental method.

The limitations in this study are as follows. The silo applied to the analysis and experiment of this study has a geometrically symmetrical shape. However, silos with different geometric properties may differ significantly from the results of this study. In addition, in this study, assuming that the friction traction due to the normal pressure in the regular discharge contributes the most to the silo hotspot detector tensile force, only two cases, from the prediction equation to the analysis result, are considered. However, in reality, phenomena such as shaking, vibration, and voids in the material occur. It is a situation that requires research considering this.

Future study will attempt to compare the stress distribution on the internal stress distribution and the silo wall and silo hotspot detector through an analysis that considers influence factors not considered in this study. In addition, the pressure change acting on the silo hotspot detector in the silo with a geometrically different shape from the silo targeted in this study will be compared.