Modeling and simulation of different and representative engineering problems using Network Simulation Method

Mathematical models simulating different and representative engineering problem, atomic dry friction, the moving front problems and elastic and solid mechanics are presented in the form of a set of non-linear, coupled or not coupled differential equations. For different parameters values that influence the solution, the problem is numerically solved by the network method, which provides all the variables of the problems. Although the model is extremely sensitive to the above parameters, no assumptions are considered as regards the linearization of the variables. The design of the models, which are run on standard electrical circuit simulation software, is explained in detail. The network model results are compared with common numerical methods or experimental data, published in the scientific literature, to show the reliability of the model.


Background
In the text of González-Fernández [1] can be found a detailed description of the fundamentals of the method and the first applications in different fields of science and engineering: electrochemical processes, transport through membranes, heat transfer, etc. After this date, the Network Simulation Method hereinafter (NSM), has been employed in new problems developing models not covered by a specific new text, so that the person concerned should refer to the specific scientific publications or doctoral theses in the research group 'Simulation Networks', from Universidad Politécnica de Cartagena (UPCT), or research groups working with this method at the Universities of Granada, Jaen and Murcia: [2][3][4][5][6][7]. In this regard it should be mentioned applications in the fields of fluid flow to transport (of mass or heat), inverse problem in heat transmission, magneto-hydrodynamic flow, mechanical vibrations, tribology, dry friction, membrane transport, elasticity, development of specific numerical calculation programs, etc.
On the other hand, there have been various programs that use the NSM as a tool for numerical calculation: PRODASIM for designing simple fins [28], PROCCA-09 for design and optimization of thermal problems [29], FATSIM-A for simulating flow fluids with solute transport problems [30], FAHET for simulation of flow fluids with heat transfer problems [31], EPSNET_10 for simulation problems elasticity [32] and OXIPSIS_12 for simulating corrosion problems [33].
The equivalence between the physical problem and the network model is that both are governed by the same differential equations in finite difference in space, covering both the elementary cell volume or as to the boundary conditions. However, time remains continuous variable in the design model.
The formal approach, which is the basis for the development of the problems, is the 'Network Theory' of Peusner [34], in which his 'thermodynamics of networks' is supported. Network models are for Peusner an accurate representation of the mathematic characteristic of the processes described. Thus, the variables which characterized of the problem must satisfy Kirchhoff's laws and their relationships determine the corresponding circuit elements. However, in each individual process and once conjugate variables chosen, the information on which circuit elements involved in the network model and how they connect with each and other, is obtained from the mathematical model and not on considerations of physical type on the role playing these variables.
In summary, in network theory, the feasibility of a network model involves: • The existence of an independent network of time • The existence of a conserved quantity called flow associated with each branch, connecting nodes and obeys the Kirchhoff Current Law (KCL) • The existence of a magnitude that satisfies criteria of uniqueness associated with each node, and obeys the Kirchhoff Voltage Law (KVL)

Dry friction between microscope tip and smooth surface at atomic-scale
This study will analyze the atomic scale dry friction phenomenon using models implemented by the network simulation method, namely the Frictional Force Microscope hereinafter (FFM), evolution of the Scanning Force Microscope hereinafter (SFM), increases research capacity of phenomena at the atomic scale, key issue for understanding the origin and nature of the fundamental laws of friction. The analysis of the behavior of these point contacts is done by a two-dimensional model of friction at the atomic scale. The response of the tip is the 'stick-slip' type [35]. After a detailed review of the literature a number of techniques and materials for analysis have been selected. The cases are: • NaF (001): it is a surface with translational symmetry. In this case the FFM is used [36] • Highly Oriented Pyrolytic Graphite hereinafter (HOPG) is well-defined structures sheets, often chemically stable. In this case an SFM and Lateral Friction Microscope hereinafter (LFM) are used [37] • Graphite: the Atomic Friction Microscope hereinafter (AFM) is used [38]

Problems of moving front
In problems with phase transitions in matter, the Stefan problem is a particular kind of boundary problem, in which this phase boundary can move with time. The classical Stefan problem aims to describe the temperature distribution in a phase change of matter. These problems are known like 'problems of moving front' [39].
There are numerous problems of moving front, which mainly include the melting and solidification processes as well as the processes of oxidation at high temperature.

The elastic problem
Theory of Elasticity describes the response of a solid object to external forces [40]. The solution to this problem involves knowing the stresses, strains and displacements generated in the solid. The analysis results allow the engineer to make decisions on the validity or suitability of structures and/or machines. On the other hand, the elasticity theory equations are widely used in other fields of science and technology, such as the study of seismic wave propagation or biomechanical behavior.

Dry friction between atomic force microscope tip and smooth surface at atomic-scale
The AFM, SFM, FFM tip physical scheme on a surface model, Fig 1, considers one rigid sliding body connected by three springs (one for each spatial direction) to the mobile support. Each spring reflects the elastic interaction between the two surfaces during the contact [36,38]. The sliding body involves the interaction with the surface and the inertial effect.
Several forces act on the tip in x-direction: inertial force represented by 'mÁ(d 2 x t /dt 2 )', where x t is the absolute displacement of the tip of mass m; elastic force from the spring, represented by the term 'k x Á(x M -x t )'; and damping force represented by the term 'c x Á(dx t /dt)'. The coefficients of these forces are constants. In addition, the surface-tip interaction force is built from a surface potential [36][37][38]. Thus, the NaF surface with FFM tip potential, Hölscher et al. [36] is: where the potential amplitude, V 0 , is 1 eV, while a x and a y , are 4.62 Å for both parameters.
Hölscher et al. [37] use a slight modification of Eq (1) to study HOPG surface and a SFM tip: where V 0 is 0.5 eV and the lattice constant, a, is 2.46 Å. Sasaki et al. [38] used the Lennard-Jones potential to study the graphite and an AFM tip: where the strength of the potential, ε, is 0.87381Á10 −2 eV, and the value of the distance between the tip and the surface atom at which the potential is zero, σ, is 2.4945 Å. The tip and the i-th atom of the surface distance, r 0i , is calculated from the i-th atom position inside the lattice, which is defined by its constant, 2.46 Å, and the nearest distance between 2 carbon atoms, 1.42 Å. Hölscher et al. employed the same potential with the same surface and microscope [41], but σ is equal to 3.4 Å. As a consequence of applying the balance of the tip forces, Fig 1, the equations could be written: The microscope movement is a set of trips in the x-direction.

Moving front problems
In moving front problems, there is increase, or decrease, of temperature, or a species, such as oxygen, increases its concentration. This depends on the process studied. The result is a new state or phase with different properties. Therefore, a new interface is generated. This interface is moved inside the structure producing the change of phase or state. Its position is represented as a distance function, Fig 2. The problem considered has a semi-infinite structure, whose free surface, x = 0, is exposed to a constant concentration of oxygen or temperature, C 0 , [42,43]. These processes are almost instantaneous. Therefore, as soon as the temperature or the oxygen concentration reaches a determinate value, c cr , at one point, solidification or melting, oxidation interface will move to that point. The process described is a moving boundary problem, whose governing equations are: where D 1 and D 2 are the diffusivities of each phase. The initial and boundary conditions are: The difference between flow values at the boundary between the two phases is responsible for the progress of the interface. Applying the conservation equation to concentrations [44,45], we have: where Δc is the oxygen concentration or temperature jump at the interface, and subscripts 1 and 2 refer to the oxide and metal phases or solid (or liquid and solid phase), respectively. Eq (11) is called the Stefan condition. This equation can be written as:

Governing equation for the elastic problem
Navier's equation [40], as the governing equation, represents the equilibrium forces in terms of displacements u for each point of the elastic body, Fig 3:λ where f is the volume force vector and λ, μ the material properties of the elastic body. For 2D problems in Cartesian coordinates, Navier's equations are reduced to two coupled differential equations: The required boundary conditions can be applied in displacements u b i or tractions terms t b i , Fig 3. The first condition is directly imposed, but the second required a complex relation by coupling of the first partial derivatives of the main unknowns and the outer normal vector of the boundary n: Solved the elastic problem defined by Eq (14), the stress components can be evaluated from the displacement using the Lamé's equations: 4. Network models

Atomic force microscope simulation network models
The initial conditions related to displacements, x M -x t = y M -y t = z M -z t = 0, and velocities, dx t / dt = dy t /dt = dz t /dt = 0, are inserted in the specifications of the capacitors initial conditions (initial voltage value) and of the coils (initial current value), respectively. The whole network model is now run employing PSpice [46,47]. The reliable network model design needs equivalence between the model and the process equations, including the initial conditions. Since Eq (4), the basic network model is corresponded with this equation, to which initial conditions must be added. Modeling basic rules are explained in González-Fernández [1]. The first is to select the equivalence between physical and electrical variables (different choices give rise to different networks). The following equivalence is established for the problem: x t (tip displacement in the x-direction) q (electric charge in the network), or dx t /dt (tip velocity in the x-direction)dq/dt = i (electric current at the network). In the other directions, similar equivalences are applied. Now, each addend of the Eq (4) is represented by an electric potential difference in a circuit element whose constitutive equation is analogous to the addend. Hence, the second Kirchhoff's law application on a network equivalent to the equation allows us to solve the equation. The addend of Eq (4) '(d 2 x t /dt 2 )' is correlated with an inductance whose constitutive equation is V L = L(di/dt) = L(d 2 q/dt 2 ), with L = 1H. Equally, the same electrical element is used for the remainder directions. The non-linear terms cannot be implemented directly, and must be defined by controlled sources or additional circuits. In the event of voltage or current source, their output is defined by a routine as an arbitrary formula of the dependent variables. The most text books do not consider this classical electrical analogy in this way, becoming the NSM more interesting and efficient in the field of numerical computation.
After the integration of 'dx t /dt', equivalent to the current in the main circuit, the term 'x t ' is available. The integration is implemented using a secondary circuit in the models for NaF surface analyzed by FFM, Fig 4(b), for graphite surface analyzed by AFM, Fig 5(c). The terms 'y t 'and 'z t ' are implemented in similar network models. The secondary circuit is implemented by a controlled source, F 1x , which generates the current, 'dx t /dt', obtained from the voltage in an ammeter connected in series in the main circuit, V x . A capacitor with C x1 = C y1 = C x = 1F is used to integrate the current of F 1x , and the voltage in this element, V Cx1 = V Cy1 = V Cx = Cx -1 Á R (dx t /dt)dt, is simply the variable 'x t '. A resistor with a very high value, R INF , improves the algorithm stability.
Once the variable 'x t ' has been determined, the addend 'k x Á(x M -x t )' in the model for NaF surface analyzed by a FFM, is implemented by a voltage from a controlled source, E 2 . For graphite surface analyzed by an AFM, another controlled source, E Tx1 , provides the necessary voltage. The secondary circuit shown in Fig 4(c) provides the variable 'x M ', the product of the constant velocity of the microscope support and time.
A pair of controlled sources, E 1 and E 1a , generates the voltage correspondent to '@V/@x t ', for NaF surface analyzed by a FFM. For graphite surface analyzed by an AFM, the controlled source, E Tx2 , generates the voltage correspondent to potential gradient. The Lennard-Jones' potential, represented in the secondary circuit of Fig 5(b), is selected. The addends associated to @V/@y t and @V/@z t are implemented in a similar way. The last one of the first equation in Eq (4) is equivalent to a resistor, characterized by 'c x ', in the models for NaF surface analyzed by a FFM. For the remainder equations in Eq (4), a similar form is implemented.
The initial conditions corresponding to displacements, x M -x t = y M -y t = z M -z t = 0, and velocities, dx t /dt = dy t /dt = dz t /dt = 0, are implemented by the initial conditions in capacitors and coils, respectively. PSpice [46,47] allows us to simulate the network model.

Moving front simulation network models
The details of the general rules of the NSM, with applications to different types of problems, are explained in González-Fernández and Alhama [1]. However, for a better interpretation the following steps are described for the design of the network model. Firstly, the equivalence between the temperature or chemical concentration (oxygen concentration) and electrical variables: c V e (electric voltage) must be established.
Secondly, it is considered that the process advances parallel to the surface so it is only necessary to consider a dimension. Therefore, the material is discretized in one-dimensional cells, dividing the spatial variable x into n volume elements, Fig 6. Thus, the terms in second derivatives of the above equations can be expressed in finite differences by @c 2 ðx; tÞ @t ¼ c 2;xþDx=2 À c 2;x Dx 2 =2D 2 À c 2;x À c 2;xÀ Dx=2 Dx 2 =2D 2 ; XðtÞ < x < 1 ð18Þ Thirdly, the right terms of the Eqs (8) and (9), dc j /dt are the currents that cross each capacitor, C j . The voltage across each capacitor, V C,j = C j -1 Á R (dθ j /dt)dt, is simply the variable c j when C j = 1F.
Fourthly, the Eq (10) gives the following differential finite difference equation: The time derivative is discretized directly by the Spice Code Software tools [46,47], which automatically adjust the time increment to reach convergence more quickly. As mentioned with spatial discretization, the derivatives are associated with balances on the cell and not on each point of the mesh.
Fifthly, we proceed to the design of the network model. Fig 7 represents a cell inside the domain. In this circuit, each term, in Eqs (17) and (18), equal an electric current that is balanced with the current of another term in a common node.
The term on the left of the above equations, dc j /dt, is the current across the capacitor C j . The voltage in each capacitor, V C,j = C j -1 Á R (dθ j /dt)dt o V C,y = C y -1 Á R (dy/dt)dt is the variable c j when C j = 1F. Since there is a change of state in the cell, from metallic to oxidized, the capacitor used must have two different initial values, according to the state in which it is. In addition, with this software it is not possible to change the initial value of a capacitor, so it has been chosen to use two capacitors whose connection to the circuit is done through two switches controlled by voltage as indicated in Fig 7. The last two linear terms of Eqs (17) and (18), I j,in and I j,out , are represented as simple resistances, R j,in and R j,out , respectively, since the equation of this electrical component is i R = V R /R. The resistance is R j,in = R j,out = (Δx) 2 /2D j . As mentioned, a change of state occurs in the cell, so the resistors used must have two different values. Since it is not possible to change these values with this software, it has been chosen to introduce, again, two resistors in parallel for each addition. An open switch disables the resistance in parallel, leaving only the resistance representing the initial state. When the switch is closed, two parallel resistors whose combined effect represents the final state are enabled , Fig 7. Finally, each term in Eq (19) is considered a voltage equation and therefore a sequence of stresses that are introduced into the circuit by means of successive switches, allowing integrating the current position of the border , Fig 8 [ 44,45].

Network model for the elastic problem
According to the NSM rules [1], the resulting cell is derived from the spatial discretization of the governing equation, Eq (14), and the adequate electro-mechanical analogy. One circuit is resultant from each differential equation. Thus, the resulting cell involves two circuits related to the unknowns u x and u y , respectively. Defining C 1 = (λ + 2μ), C 2 = μ and C 3 = (λ + μ), governing equation can be written in a more compact form: Using the mesh and notation shown in Fig 9, each term in Eq (20) can be expressed in these finite differences: ' u kt;2 À u kt;4 Dx À u kb;2 À u kb;4 Dx 2Dy ¼ u kt;2 À u kt;4 À u kb;2 þ u kb;4 2DxDy The resulting differences equations are: ux k;0 À ux k;2 Dx 2 2 1 C 1 þ ux k;0 À ux k;4 Dx 2 2 1 C 1 þ ux k;0 À ux k;3 Dy 2 2 1 C 2 þ ux k;0 À ux k;1 Dy 2 2 1 C 2 À C 3 uy kt;2 À uy kt;4 À uy kb;2 þ uy kb;4 2DxDy þ f x ! ¼ 0 uy k;0 À uy k;2 Dx 2 2 1 C 2 þ uy k;0 À uy k;4 Dx 2 2 1 C 2 þ uy k;0 À uy k;3 Dy 2 2 1 C 1 þ uy k;0 À uy k;1 Dy 2 2 1 C 1 À C 3 ux kt;2 À ux kt;4 À ux kb;2 þ ux kb;4 2DxDy þ f y Establishing the analogy between the mechanical displacement and electrical voltage, each equation shows a current balance. The resulting network model emerges from the electrical connection, in the whole domain N x ×N y , among the cells shown in Fig 10. The model presents two coupled circuits: circuit 'ux' in corresponding with the rectangular u x displacement and one more for the u y component. Each cell contains the following components for circuits: four resistors corresponding with the primary four addends in both equations; values of the resistors are expressed in the denominator. In other hand, coupled terms between equations and volume forces (last addends in brackets), are implemented by using a voltage-controlled voltage-source. With regard to boundary conditions, while displacements are directly implemented by a constant voltage source, tractions, Eq (15), need the use of voltage-controlled voltage-sources.

Atomic-scale dry friction simulation: Atomic force microscopes network models
Network Simulation Method application of friction processes contained in the introduction is presented. For these processes, it is available partial results obtained by other authors with whom it is possible to verify the method. Furthermore, the effect of the parameters on the behavior of the proposed systems is discussed.
The NSM has been applied to studied surfaces and tips for which data are available: NaF and FFM tip [36], HOPG and SFM tip [37], and graphite and AFM tip [41]. 5.1.1 Friction force microscope on a sodium fluoride surface simulation. The tip characteristics, the surface size and the test conditions for the tests selected are: • The interaction between the NaF surface and the FFM tip is defined m x = m y = 10 −8 kg, c x = c y = 10 −3 NÁs/m and k x = k y = 10 N/m. The tip support velocity is 400 Å/s, the scan range of 20x20 Å and the scan angle of 0˚are used in order to show the results [36]. The potential defined in Eq (1) is used in the model represented by Eq (4).
• The interaction between the HOPG surface and the SFM tip is defined by m x = m y = 10 −8 kg, c x = c y = 10 −3 NÁs/m and k x = k y = 25 N/m. The tip support velocity is 400 Å/s, the scan range of 20x20 Å and the scan angle of 7˚are used in order to show the results [36]. The potential defined in Eq (1) is used in the model represented by Eq (4). The studied surface, with a very high stiffness [37], is integrated by 271 hexagons, defined by six carbon atoms (600 altogether).
• The interaction between graphite surface and the AFM tip is defined by c x = c y = 0 NÁs/m, k x = k y = 0.25-2.5 N/m, k z = 0.25 N/m, and a punching of -6 Å. The scan range of 9.8x8.5 Å and the scan angle of 15˚are used in order to show the results [38]. Besides, m x = m y = 10 −6 kg and a tip support velocity of 10 Å/s, have been assumed in order to compare the results with Sasaki's ones [38].
The potential functions defined in the following equations are used in Eq (4): • Eq (1)  In the last case, for graphite surface, the assumed model parameters are obtained using short computational time. Fig 11(b) indicates the tip elastic force components, F x and F y , in each NaF surface position, whose crystal net is presented in Fig 11(a). The reciprocation between both represented images allows us to handle this code as an image interpretation tool from the microscope. Thus, this code is capable to handle a lot of crystal nets, even with crystallographic defects, providing the appropriate elastic forces, which could be contrasted with those from the microscope. The simulation results correspond to those accomplished by Hölscher et al. [36]. 5.1.2 Atomic force microscope on a highly oriented pyrolytic graphite sample and on a graphite sample simulation. Fig 12(b) indicates the tip elastic force components, Fx and Fy, in each HOPG surface position, whose crystal net is presented in Fig 12(a). The simulation results correspond to those provided by Hölscher et al. [37]. Besides, Fig 13(b) shows the elastic force x-direction, Fx, in each position of the graphite surface, Fig 13(a). In this case, a 3D plot is employed to depict the details better, specifically the stick-slip. The lines are plotted at different y-coordinates, separated by a fixed step. When one of these intersects the surface at a critical position, the instabilities appear, Fig 13(b). The simulation results correspond to those obtained by Sasaki et al. [38].

Moving corrosion fronts in titanium matrix composites
Lagoudas et al. [43] calculated the oxygen distribution throughout the oxide layer into the metal, assuming a rate position of the boundary between metal and oxide phases proportional to the square root of time. We solve this problem without making assumptions. In their analysis for SiC/Ti-15-3, used the values of diffusivity in the oxide, D 1 , and in the metal, D 2 , from the parameters of the Arrhenius equation, for temperatures between 600˚C and 700˚C, [39,48]. These parameters are shown in Table 1 [39,[48][49]. Other authors, as Ariel et al. [49], estimate the diffusivity of oxygen in titanium 99.9% of purity, D 2 , using the parameters of the Arrhenius equation shown in the Table 1. We have used these parameters for modeling the oxidation of 99.9% pure titanium. Parameters of the oxide layer are the same as those used by Lagoudas et al [42,43]. The coefficients to calculate the diffusivity are defined in Table 1.   If Fig 14 is compared with previous work [45], it can be seen that the matrix SiC/Ti 15-3 is oxidized at a faster rate than titanium 99.9%. The explanation is the higher oxygen diffusivity of the first material. As shown, the second material takes nearly 20 hours to achieve an oxide thickness of 2 μm. Fig 15 shows as the progress of the oxidation frontier for a given temperature is stopped when it reaches a certain position, a determinate oxide thickness. After a time interval, the frontier advance is restarted. Furthermore, as expected, an increase in temperature reduces this interval. This figure shows a radically different behavior of matrix SiC/Ti-15-3 [45]. The explanation for this different behavior is due to a titanium diffusivity value lower in three orders of magnitude as the matrix SiC/Ti-15-3.  As expected with increasing temperature is an increase of the oxidation rate. These results are the same as those established for titanium alloys by ASTM B861-10 [50] and 'Asociación de Investigación de la Rama Metalmecánica de Valencia (AIMME)' [51]. Their studies establish that the best alloys are those containing palladium and commercially pure titanium, 99.9%.

Cantilever beam subjected to shear loading
To show the applicability of the NSM in the elasticity field bending of a cantilever is solved for a rectangular cross section, Fig 16, Morales et al. [52]. This plane stress problem has been studied by Timoshenko & Goodier [53] under ideally boundary condition. The resultant shear force at the free end F = 20Á10 3 N is applied following a parabolic distribution according with the elementary strength material theory (p = 30 MPa). For the right end the boundary condition zero-displacement is imposed. Other geometric and material parameters are L = 200 mm, H = 100 mm, thickness = 10 mm, λ = 121 GPa and μ = 80 GPa. Fig 17 presents the results of the network simulation in term of the von Mises stress over the deformed domain draw with a scale factor of 50. Next figures (Figs 18 and 19) exhibit a comparison between NSM and FEM [54] numerical simulations in terms of u x -displacement and σ xx -stress. Fig 20 shows the vertical displacements at y = 0 for the analytical solutions studied by Timoshenko and Goodier (two ideal cases, Eqs n and r [53]), and the numerical solutions obtained by means of NSM and FEM for the boundary condition depicted in Fig 16. The numerical solutions of the deflection curves are very closed and bounded by the two analytical and ideal cases. Fig 21 shows the σ xx -stress solutions for analytical and numerical methods for the upper edge line, y = H/2. These results are according to Saint-Venant's principle, where only significant differences can be appreciated near to the right section. For a better comparison, Table 2 details numerical results in points A, in the origin, and B, at 3/4 of L. (Fig 16).

Conclusions
Some numerical models based on the NSM have been designed and used-with negligible computing times in a suitable circuit simulation code-to successfully simulate different and representative engineering problems, such as dry friction in atomic scale, the moving front problems and elastic and solid mechanics, with no assumptions concerning the linearization  of the governing equations. The influence of the problems' main parameters is studied within a practical range of them. Comparison between the network model results and standard numerical methods or experimental data, which have been published in the scientific literature, confirms the reliability of the model.