Discrete and Continuous Models and Applied Computational ScienceDiscrete and Continuous Models and Applied Computational Science2658-46702658-7149Peoples' Friendship University of Russia named after Patrice Lumumba (RUDN University)2421910.22363/2658-4670-2020-28-2-141-153Research ArticleKinematic support modeling in SageKroytorOleg K.Postgraduate of Department of Applied Probability and Informaticskroytor_ok@pfur.ruMalykhMikhail D.Doctor of Physical and Mathematical Sciences, assistant professor of Department of Applied Probability and Informaticsmalykh_md@pfur.ruKarnilovichSergei P.assistant professor, Ph.d., assistant professor of Institute of Physical Research and Technologykarnilovich_sp@pfur.ruPeoples’ Friendship University of Russia (RUDN University)1512202028214115320072020Copyright © 2020, Kroytor O.K., Malykh M.D., Karnilovich S.P.2020The article discusses the kinematic support, which allows reducing the horizontal dynamic effects on the building during earthquakes. The model of a seismic isolation support is considered from the point of view of classical mechanics, that is, we assume that the support is absolutely solid, oscillating in a vertical plane above a fixed horizontal solid plate. This approach allows a more adequate description of the interaction of the support with the soil and the base plate of the building. The paper describes the procedure for reducing the complete system of equations of motion of a massive rigid body on a fixed horizontal perfectly smooth plane to a form suitable for applying the finite difference method and its implementation in the Sage computer algebra system. The numerical calculations by the Euler method for grids with different number of elements are carried out and a mathematical model of the support as a perfectly rigid body in the Sage computer algebra system is implemented. The article presents the intermediate results of numerical experiments performed in Sage and gives a brief analysis (description) of the results.kinematic supportseismic isolation supportmathematical modelfinite difference methodcomputer algebra systemSagenumerical calculationsSageкинематическая опорасейсмоизолирующая опораматематическая модельМКРсистема компьютерной алгебрычисленные расчёты1. Introduction As one of the types of seismic protection of buildings in seismically active regions of the Earth, among others, auxiliary “seismic suppression” supports are used. A large number of standards of such supports are known with a high level of reliability of operation and a high ability to damp seismic waves of high magnitude. They are high-tech in manufacturing, are sold at high prices and, thanks to efficiency and reliability, are in high demand among large construction companies and in areas with a high standard of living. At the same time, in poor areas of developing countries, citizens and municipalities cannot afford to use them. However, it is precisely in such regions that residents most often suffer from the devastating effects of earthquakes. The Soviet, and then the Russian school of seismic protection architecture was able to offer inexpensive and effective solutions in the form of so-called kinematic supports. In our country and abroad, a large number of active seismic protection systems for buildings have been proposed, developed and applied. Among them are those proposed by A. Kurzanov, S.Yu. Semenov [1], [2], Yu.P. Cherepinsky [3], V.V. Nazina, etc. Some of these systems were practically implemented in real buildings [1], [2], [4], [5], making it possible to assess their workability for building industry. Vibration tests were carried out at many facilities [6], [7], which provided experimental data on the behavior of these systems under dynamic impacts. However, essentially all developed systems need additional analysis under full-scale conditions. Therefore, many aspects of the real behavior of seismic protection systems are difficult to study theoretically or on models due to the very large number of factors affecting the behavior of the structure during an intense earthquake. Kinematic supports are vertically placed cylinders on which the building rests. Neither the place of entry of the support into the ground or concrete slab, nor the contact with the horizontal slab of the building placed on such supports, are fixed rigidly. Supports can be made in the form of short concrete pillars with an outer cage of steel pipe or a reinforcing cage of carbon composite or basalt composite nets. It is promising to use in the construction of concrete racks with dispersed reinforcement basalt fiber, since such concretes have increased resistance to cracking and tensile strength during bending. The essence of the kinematic support is that when the base is displaced by a certain design value, the building slightly rises, receiving some additional kinetic energy. In this case, a returning torque arises, bringing the ”base-building” system to its original state (position before the earthquake). Residential buildings constructed using seismic isolating supports have fullscale confirmation of the reliability of the structure and have proven themselves successful in experimental studies [1], [2]. The aim of our work is to create an adequate mathematical model of the support and its interface with the building, which will help to design kinematic supports taking into account the operational requirements of customers. The solution to this problem can be approached from two sides. First, it is possible to create a model of elastic support in the Ansys, system [8]. With this approach, the main difficulty is the selection of adequate boundary conditions for the place where the support is in contact with the soil and base plate of the building. The second approach proposes to consider the support from the point of view of classical mechanics, assuming it to be a perfectly rigid body, oscillating in a vertical plane above a fixed horizontal absolutely smooth plate, but more adequately describe the interaction of the support with the ground and the building plate. From the point of view of analytical mechanics, the motion of kinematic supports is the motion of a complex system of bodies with non-holding bonds [9]. The mathematical basis of the dynamics of such systems was developed by outstanding mathematicians J. Dalamber, S. Poisson, Yu.A. Arkhangelsky, V.V. Kozlov, A.P. Markeev and others. In our work, we will follow the formalism developed by A.P. Markeev. Models of rigid body mechanics are systems of a large number of ordinary differential equations (ODEs) that are not resolved with respect to derivatives [10], [11]. They, as a rule, do not allow an analytical solution. Therefore, we are going to use the finite difference method. For its successful application, it is important to solve the equations for derivatives by increasing the order. Due to the large number of equations, this procedure turns out to be rather complicated. Therefore, it seems natural to execute it in a computer algebra system (CAS). Computer-aided study of such models is carried out in two stages: 1) symbolic transformations reducing the system to normal form, used in the standard formulation of the Cauchy theorem; 2) numerical solution of this system using the finite difference method (rk4). We chose the Sage system [12], because it can execute both the first and second stages. To test the described approach, we took the simplest model that describes the complete system of equations of motion of a massive rigid body along a fixed horizontal perfectly smooth plane [10], [11]. 2. Description of the mathematical model Let the seismic isolation support be a rigid body and have the shape of a cylinder, in which one of the bases has a spherical shape. The support is installed between the foundation of the building (a horizontal rigid plate) and the building itself. The support touches the foundation of the building always with a spherical end [13]. We assume that when oscillating or when horizontal forces act on the foundation of the structure, the point of contact of the support and the base plate always lies in the plane ????????????. To describe the vibration of the support, we will use the model of motion of a rigid body on the surface [10], [11] and adapt this model to our task. The motion of the body will be considered relative to the fixed laboratory coordinate frame ???????????? with the origin ???? at a certain point of the plane. The axis ???????? is directed vertically, ???? is the unit vector of inner normal to the body surface at point of the axis ????. Let ???????????????? denote the moving coordinate frame rigidly coupled to the body with the origin at its center of gravity ???? and the axes directed along the principal axes of inertia 1. The orientation of the body with respect to the fixed laboratory frame is specified by the Euler angles ????,????,???? or by the matrix of direction cosines ????????????. The unit vector of the ????-axis in the frame ???????????????? is specified by the components ????31, ????32, ????33: ????31 = sin????sin????, ????32 = sin????cos????, ????33 = cos????. Assume ????, ????, ???? to be the principal axes of inertia with respect to the gravity center. Let ???? be the point of contact between the horizontal plane ???????????? and the support (see Figure 1). Its coordinates ????, ????, ???? in the frame ???????????????? will be functions of angles ????,????, determined from the form of equation ????(????,????,????) = 0 that specifies the shape of the body surface [10], [11]. The sign of function ???? is chosen such that ???? = -. Then the quantities ???????????? for the axis ???? are expressed in terms of the Euler angles as 1 ???? ???? = - , |∇????| ∇???? { 1 ???? ⎨????32 = sin????cos???? = -|∇????| ∇????, (1) { { 1 ???? {⎩????33 = cos???? = -|∇????| ∇????. Figure 1. Motion of a perfectly rigid body above the perfectly smooth horizontal plane Following the studies presented in Refs. [10], [11], [14], [15], let us consider the full system of equations of motion for a massive rigid body on a fixed horizontal perfectly smooth plane and introduce the following unknown functions of time ????: center of gravity coordinates ????, ????, ???? of the body in the laboratory frame; Euler angles ????,????,????; components ????, ????, ???? of radius vector ???? of point ???? of contact of the support and the plane (base plate) relative to the gravity center, and magnitude ???? of the normal reaction of the plane. To determine the unknowns listed above the following equations and relations will be used: a) The equations that represent the theorem of momentum variation. Theexternal forces are the reaction of the plane ???? = ???? ⋅ ???? directed vertically (???? ⩾ 0) and the gravity force. The equations are: ⎧{????????̈ = 0, ????????̈ = 0, (2) ⎨{⎩????????̈= -???????? + ????, where ???? is the body mass, ???? is the free fall acceleration. b) The equations that represent the theorem of angular momentum variation: ⎧{???????????????????? + (???? - ????)???????? = ????(????????33 - ????????32), { { ???????? ⎨{???????????? + (???? - ????)???????? = ????(????????31 - ????????33), (3) { ???????? {⎩???????????? + (???? - ????)???????? = ????(????????32 - ????????31), where ????, ????, ???? are the projections of the angular velocity vector ???? on the axes of the coordinate system ????????, ????????, ???????? rigidly bound to the body; ????, ????, ???? are the principal moments of inertia with respect to these axes. c) The relations represented by the Euler kinematic equations: ⎧???? = ????̇ sin ????sin ???? + ????̇ cos ????, { ⎨???? = ???? ????????????????̇ cos ???? - ????̇ sin ????, (4) {⎩???? = ????̇ cos ???? + ????.̇ d) The Poisson equations: ⎧{????̇31 = ????32???? - ????33????, ⎨{⎩????̇????̇3233 == ????????3331???? - ???????? - ????3231????????, (5) indicating the fact that vector ???? defines an invariable direction in the fixed coordinate frame ????????????????; e) The equation of the body surface in the coordinate frame rigidly bound to the body and having the origin at its center of gravity: ????(????,????,????) = 0. (6) In our case, the body surface equation is the equation of a sphere (the support has spherical shape at the point of contact with the plate): ????2 + ????2 + (???? + ????)2 = ????2. (7) f) The constraint equation ???? = -???? ⋅ ???? = -(????31???? + ????32???? + ????33????) (8) means that the support moves contacting with the plate all the time. Equations and relations (1)-(8) determine a complete system of equations of motion of a massive rigid body on a fixed horizontal perfectly smooth plane that can be written in the following form: ⎧????????̈ = 0, {{????????̈ = 0, {{????????̈= -???????? + ????, { ???????? {{???????????? + (???? - ????)???????? = ????(????????33 - ????????32), {{{{???????????????????? + (???? - ????)???????? = ????(????????31 - ????????33), { ???????? {???????????? + (???? - ????)???????? = ????(????????32 - ????????31), { {{????̇31 = ????32???? - ????33????, {{????̇32 = ????33???? - ????31????, {????̇33 = ????31???? - ????32????, (9) ???? + ????32???? + ????33????), { 1 ???????? {{????31 = sin ???? sin ???? = -|∇????| ???????? , {{{????32 = sin ???? cos ???? = -|∇????|1 ????????????????, { {{{????33 = cos ???? = -|∇????|1 ???????????????? , { {{???? = ????̇ sin ????sin ???? + ????̇ cos ????, {{???? = ????̇ sin ????cos ???? - ????̇ sin ????, {⎩???? = ????̇ cos ???? + ????.̇ From relation (1) and equation (7), we get explicit expressions of ???????????? for the spherical base of a kinematic support: ⎧{{????31 = -|∇????|1 ???????????????? = -√4????2 + 4????21+ 4(???? + ????)2 ⋅ 2???? = -????1 ????, { { 1 ???????? 1 1 ⎨{????32 = -|∇????| ???????? = -√4????2 + 4????2 + 4(???? + ????)2 ⋅ 2???? = -????????, (10) 1 ???????? 1 1 ⎩????33 = -|∇????| ???????? = -√4????2 + 4????2 + 4(???? + ????)2 ⋅ 2???? = -????(???? + ????). Using the formulae (10), is it easy to express the components ????, ????, ????, therefore, we can get rid of the Poisson equations in the system (9), since the values of ???????????? and components ????, ????, ???? are already known. According to the two first equations of (9), the center of gravity moves so that its projection on the base horizontal plane moves rectilinearly. Hence, it is obvious that the appropriate equations can be disregarded, too. After a number of executed transformations and assumptions, the system (9) takes the form ⎧???? = ????????̈+ ????????, {{{{???????????????????? + (???? - ????)???????? = ????(????????33 - ????????32), { ???????? {{???????????? + (???? - ????)???????? = ????(????????31 - ????????33), {{{???????????????????? + (???? - ????)???????? = ????(????????32 - ????????31), + (???? + ????)2 = ????2, {{???? = -???? ⋅ ???? = -(????31???? + ????32???? + ????33????), 1 ???????? 1 (11) ⎨????31 = sin ???? sin ???? = -|∇????| ???????? = -????????, { { { 1 ???????? 1 {{????32 = sin ???? cos ???? = -|∇????| ???????? = -????????, { { 1 ???????? 1 {{????33 = cos ???? = -|∇????| ???????? = -????(???? + ????), { {???? = ????̇ sin ????sin ???? + ????̇ cos ????, { {{???? = ????̇ sin ????cos ???? - ????̇ sin ????, {⎩???? = ????̇ cos ???? + ????.̇ Equations (11) are reduced to the form convenient for using the computer algebra system Sage, so that further simplifications and the solution of these equations will be carried out in this CAS. In the next section we investigate the solubility of this system of equations in Sage. 3. Resolving the system with respect to derivatives The system (11) includes both differential equations and algebraic ones (relations). Let us use Sage to reduce it to a simper form. It is obvious that from the system (11) via the values of the direction cosines ????31, ????32, ????33, we can explicitly express variables ????, ????, ???? as follows: ⎧{???? = -???? ⋅ ????31 = -???? ⋅ sin ???? sin ????, ???? = -???? ⋅ ????32 = -???? ⋅ sin ???? cos ????, (12) ???? - ????. Provided that ????, ????, ???? (12) and ????31, ????32, ????33, are known, it is possible to express the values of the function ????. Let us write in Sage the equation ???? = -???? ⋅ ???? = -(????31???? + ????32???? + ????33????) from system (11) and substitute into it the expressions for ????, ????, ???? and ????31, ????32, ????33, then function ???? will be represented by the following symbolic expression: ???? = ???? ⋅ cos????2 ⋅ sin????2 + ???? ⋅ sin????2 ⋅ sin????2 + (???? ⋅ cos???? + ????) ⋅ cos????. We substitute the obtained expression of ???? into equation ???? = ????????̈+ ???????? and calculate the normal reaction of the plane: ???? = (2⋅????⋅cos????2 ⋅cos????2 ⋅????2̇ +2⋅????⋅cos????2 ⋅sin????2 ⋅????2̇ -2⋅????⋅cos????2 ⋅sin????2⋅ ⋅ ????2 - 2 ⋅ ???? ⋅ sin????2 ⋅ sin????2 ⋅ ????2̇ + 2 ⋅ ???? ⋅ cos????2 ⋅ cos???? ⋅ sin???? ⋅ ???? + 2 ⋅ ???? ⋅̈ cos????⋅ ⋅ sin????2 ⋅ sin???? ⋅ ???? - ???? ⋅̈ cos????2 ⋅ ????2̇ + 2 ⋅ ???? ⋅ sin????2 ⋅ ????2̇ - (???? ⋅ cos???? + ????) ⋅ cos????⋅ ⋅ ????2 - ???? ⋅ cos???? ⋅ sin???? ⋅ ???? - (???? ⋅̈ cos???? + ????) ⋅ sin???? ⋅ ????) ⋅ ???? + ???? ⋅ ????.̈ As soon as all quantities ????, ????, ????, ????31, ????32, ????33, ????, and ????, ????, ???? are explicitly expresses, we substitute them into equations 2-4 of the system (11) and of the form ⎧Φ(????,????,????,????,…̇ ????) = 0,̈ { ⎨Ψ(????,????,????,????,…̇ ????) = 0,̈ (13) {⎩Θ(????,????,????,????,…̇ ????) = 0.̈ Explicit expressions for Φ,Ψ,Θ were be found in Sage [12]. The system (13) is linear with respect to ????,̈ ????,̈ ????.̈ Now we resolve the obtained system of differential equations with respect to higher derivatives using the function solve (). Ultimately, we arrive at the system of differential equations of the 6-th order resolved with respect to higher derivatives. The system incorporates three equations of the second order with respect to Euler angles. For further solution of the problem we have to decrease the order of the differential equations to the first one. To reduce the order of the system of differential equations, we perform the changes ???? = ????̇ , ???? = ????̇ , ???? = ????̇ . As a result, we get a system of six first-order differential equations. 4. Numerical experiments in SAGE Let us implement explicit Euler method, in order to confirm the absence of errors related to transformation of data types. The calculations by the Euler method will be performed using the grids with ???? = 400,800,1600 under the following initial conditions: ????0 = 0.1, ????0 = 0.1, ????0 = ????+0.5, ????0 = 1.5, ????0 = 1.5, ????0 = 0, ℎ = 2????/????, ???? = 0, ???? = 304, ???? = 304, ???? = 400, ???? = 1.5, ???? = 0.2, ???? = 500, ???? = 9.8 where ???? is the number of steps; ℎ is the step; ????0, ????0, ????0 are the Euler angles at the initial moment of time; ????0, ????0, ????0 are the initial velocities; ????, ????, ???? are the moments of inertia. The results obtained by the Euler method are presented in Figures 2-7. For ???? = 400 the results are plotted by dash lines, for ???? = 800 by dash-dot lines, and for ???? = 1600 by dot lines. The accuracy of the obtained results is estimated using the Richardson method (see Figures 9, 10). Figure 8 shows graphs. The distance between the center of mass and the plane ground oscillates as it shown on Figure 8. Figure 2. The angle ???? as a function of time Figure 3. The angle ???? as a function of time for the grids with ???? = 400,800,1600 for the grids with ???? = 400,800,1600 Figure 4. The angle ???? as a function of time Figure 5. Rate of angle ???? change with time for the grids with ???? = 400,800,1600 for the grids with ???? = 400,800,1600 Figure 6. Rate of angle ???? change with time Figure 7. Rate of angle ???? change with time for the grids with ???? = 400,800,1600 for the grids with ???? = 400,800,1600 Figure 8. The center of mass height above the plane for the grids with ???? = 400,800,1600 Figure 9. Accuracy of Euler solution Figure 10. Accuracy of Euler solution estimated by Richardson method estimated by Richardson method for the grids with ???? = 400,800,1600 for the grids with ???? = 400,800,1600 The results obtained by the Euler method for the assumed boundary conditions show that the angles ????, ???? linearly increase with time (see Figures 2 and 3), while the angle ???? experiences harmonic-like oscillations with growing amplitude (see Figure 4). Figures 5-7 show that the rates of change of the Euler angles vary according to a harmonic-like law, the amplitude of oscillations increasing with time. 5. Conclusions A crude mathematical model of a kinematic support is constructed, in which the support is considered as a perfectly rigid body oscillating in a vertical plane above a fixed horizontal perfectly smooth plate. The approximate model is rigid and does not take friction into account. A procedure for reducing the system of differential equations of the model to a form suitable for applying the finite difference method is described and implemented in the Sage computer algebra system. An explicit Euler method is implemented for grids with the number of partitions 400, 800, 1600 and the accuracy of the solution is estimated by the Richardson method. For the initial and boundary conditions specified by us, the time dependence of Euler angles and their rates is determined. The results of numerical experiments are consistent with the general idea that small deviations lead to small oscillations of the support. The problems and lines of study that will be addressed at further stages of the research are identified as follows. First, to solve the problem, we have to to determine the correct additional conditions under which the construction is stable or loses stability, i.e., when at strong ground vibrations the center of gravity horizontally shifts beyond the limits of return movement of the support and the support begins to tip over [16]. After finding the correct initial and boundary conditions, to create an adequate mathematical model of the support, we will take friction into account, specifying the reaction of the plane (base plate), and the effect of earthquakes on the movement of the supports.[A. Kurzanov and N. Skladnev, “Seismo-isolating supports for dwellings,” in Earthquake Engineering, Tenth World Conference, Balkema, Rotterdam, 1992, pp. 1951-1954.][A. M. Kurzanov and S. Y. Semenov. (2013). “Pipe-concrete seismic isolating support [Trubobetonnaya seysmoizoliruyuschaya opora],” [Online]. Available: https://patents.google.com/patent/RU2477353C1/enl.][Y. Cherepinsky, “Seismic stability of buildings on kinematic supports,” Soil Mech Found Eng, vol. 9, pp. 164-168, 1972. DOI: 10.1007/ BF01709306.][Y. Cherepinsky, “The seismic isolation of residential buildings,” Res. Seism. Stab. Build. Struct., vol. 23, pp. 438-462, 2015.][T. Belash, U. Begaliev, S. Orunbaev, and M. Abdybaliev, “On the Efficiency of Use of Seismic Isolation in Antiseismic Construction,” American Journal of Environmental Science and Engineering, vol. 3, no. 4, pp. 66-74, 2019. DOI: 10.11648/j.ajese.20190304.11.][V. Smirnov, J. Eisenberg, and A. Vasil’eva, “Seismic isolation of buildings and historical monuments. Recent developments in Russia,” in 13th World Conference on Earthquake Engineering. Vancouver, B.C., Canada, August 1-6, paper no. 966, 2004.][A. Uzdin, F. Doronin, G. Davydova, et al., “Performance analysis of seismic-insulating kinematic foundations on support elements with negative stiffness,” Soil Mechanics and Foundation Engineering, vol. 46, pp. 99-107, 2009. DOI: 10.1007/s11204-009-9052-1.][O. Kroytor, “The penetration modeling of flat obstacles in Ansys Autodyn program,” IOP Conference Series: Materials Science and Engineering, vol. 675, p. 65, 2019. DOI: 10.1088/1757-899X/675/1/012027.][J. Awrejcewicz, “Kinematics of a Rigid Body and Composite Motion of a Point,” in Classical Mechanics. Advances in Mechanics and Mathematics. May 2012, pp. 263-399. DOI: 10.1007/978-1-4614-3791-8_5.][A. Markeev, “On the theory of motion of a rigid body with a vibrating suspension,” Doklady Physics, vol. 54, pp. 392-396, 2009. DOI: 10.1134/ S1028335809080114.][A. Markeev, “On the motion of a heavy dynamically symmetric rigid body with vibrating suspension point.,” Mechanics of Solids, vol. 47, pp. 373-379, 2012. DOI: 10.3103/S0025654412040012.][(2013). “SageMath, the Sage Mathematics Software System (Version 9.0.1), ReleaseDate:2020-02-20,” [Online]. Available: https://www. sagemath.org.][K. Bissembayev, A. Jomartov, A. Tuleshov, and T. Dikambay, “Analysis of the Oscillating Motion of a Solid Body on Vibrating Bearers,” Machines, vol. 7, no. 3, p. 58, 2019. DOI: 10.3390/machines7030058.][A. Markeev, “On the dynamics of a rigid body carrying a material point,” Regular and Chaotic Dynamics, vol. 17, pp. 234-242, 2012. DOI: 10.1134/S1560354712030021.][A. P. Markeyev, “The equations of the approximate theory of the motion of a rigid body with a vibrating suspension point,” Journal of Applied Mathematics and Mechanics, vol. 75 (2), pp. 132-139, 2011. DOI: 10.1016/j.jappmathmech.2011.05.002.][S. Karnilovich, K. Lovetskiy, L. Sevastianov, and E. Shchesnyak, “Seismic stability of oscillating building on kinematic supports,” Discrete and Continuous Models and Applied Computational Science, vol. 27, no. 2, pp. 124-132, 2019. DOI: 10.22363/2658-4670-2019-27-2-124-132.]