Journal of the Civil Engineering Forum. January 2026, 12. :1-10 DOI 10. 22146/jcef. Available Online at https://jurnal. id/v3/jcef/issue/archive Cold-Formed Steel Truss Roof Structure Failure Considering Seismic Load and Buckling Analysis Muslikh1* . Miftahul Iman2 1 Department of Civil Engineering. Universitas Muhammadiyah Yogyakarta. Yogyakarta. INDONESIA 2 Department of Civil Engineering Master Program. Universitas Janabadra. Yogyakarta. INDONESIA *Corresponding author: muslikh@umy. SUBMITTED 28 April 2024 REVISED 05 September 2024 ACCEPTED 19 June 2025 ABSTRACT There were many incidents of cold-formed steel roof truss structures in the last 5 years in Indonesia. Various kinds of allegations have been addressed to cold-formed steel material applications especially in the case of seismic resistance. Some of them concern the authenticity of the steel material itself and the selection of cold-formed steel material. On the other hand, recently, people have installed . cold-formed steel trusses without involving a certified cold-formed steel applicator. This research is based on a numerical study that modeled the collapse pattern of cold-formed steel truss roof structures by considering buckling failure and the seismic load capacity. The cold-formed steel roof truss structure was modeled with 3D-truss elements in two model types: the overall structure and a single compression member element in 3D solid idealization. Buckling analysis with eigenvalue and nonlinear static analysis was performed to evaluate the critical load (Pcr ). The buckling mode shape also was also compared with the mode shape of modal analysis. This research also evaluated the effect of seismic load on the overall cold-formed steel truss structure and the slenderness of the compression member. The numerical simulation of cyclic loading on the single compression member was conducted in this research. The numerical analysis results showed that cold-formed steel roof truss structure have high vulnerability to seismic hazard effect. The cold-formed steel material has lower ductility than hot rolled steel material. This causes the lateral displacement that occurs to be lower than the displacement produced by the hysteretic curve of numerical cyclic simulation. This research also evaluated the dynamic properties, such as frequency, periods, and mode shapes, of some typical cold-formed steel for roof truss structure. KEYWORDS Buckling. Cold-formed steel. Earthquake. Hysteretic. Nonlinear. A The Author. This article is distributed under a Creative Commons Attribution-ShareAlike 4. 0 International license. 1 INTRODUCTION Some cold-formed steel truss failure cases have been recorded in the last five years. Several roof constructions that were made of cold-formed steel failed immediately. Various kinds of presumptions have been addressed to cold-formed steel material properties, especially in the case of seismic resistance. Some people even doubt the quality of cold-formed steel materials that are available on the market. On another hand, many people have recently installed . coldformed steel truss without involving a certified coldformed steel applicator, which caused the slenderness ratio of the member to not be considered in the design. The geometrical failure often occurs early before material failure of the structure. Hui et al. conducted a study by carrying out analysis and monitoring of stresses and delays on a 54 m long span steel roof truss. Numerical and theoretical analysis was carried out to illuminate the roof truss points experiencing critical stress. The research results showed that the maximum deflection met the requirements. The maximum stress indicates that the stress that was elastic. Miftahul et al. conducted a nu- merical finite element analysis on the roof steel truss Tyfekci et al. investigated the failure of the roof that was constructed by cold-formed steel by numerical analysis had been conducted. The results showed that the cold-formed steel truss structure was adequate to support the gravity load but had a high vulnerability to lateral load actions, such in seismic load. Johnson et al. performed the experimental test to evaluate the reliability of cold-formed steel roof truss structure. Pranoto and Jepriani . conducted study on cold-formed steel material by comparing several cold-formed steel popular brands that available through distributors. Based on investigation, all popular brands of cold-formed steel were sufficient to support the loads carried by roof truss structure. Pise and Sawai . also conducted study focused on the seismic performance of cold-formed steel and hot rolled steel for truss structure. The research results showed that significant reduction of chord capacity and the deflection due to the incremental increase of ultimate load carrying capacity. Journal of the Civil Engineering Forum Vol. 12 No. 1 (January 2. Miftahul et al. conducted a numerical study by comparing two numerical models of slender tubular compression members respectively: 3D shells and 3D solid idealization. By considering eigenvalue analysis . inear buckling analysi. and P Oe nonlinear analysis the critical load (Pcr ) could be evaluated. Ismail et al. also studied on buckling behaviour of coldformed stub channels under compression load by conducting experimental test and numerical analysis. The results showed that the cold-formed stub channel could be considered as an ideal stub compression member. Sumit et al. also studied on buckling analysis of C-sectional cold-formed steel using numerical analysis. Sani and Muftah . conducted another comparison study on Warren-type trusses from cold-formed steel and hot rolled steel structure. vides the numerical software formula to determine linear buckling load, as shown as in Equation . [K Oe G. ] O = 0 where K is material stiffness matrix, is eigenvalue . umerical buckling loa. is geometric stiffness matrix, and O is eigen vector . uckling mod. Equation . is called linear eigenvalue analysis. Next, the buckling load formula can be evaluated by determining the nonlinear P Oe asymptotic curve. Suhendro . considered the second order of strain displacement function in Equation . , which can be re-written in Equation . [K0 ] Pcr [N1 ] Pcr2 [N2 ] {Q } [OIQ] = 0 Anishma and Anjusha . conducted an experimental study to review the cold-formed steel roof truss structure system. The results showed the cold from steel had problems related to stability, such as buckling Based on the previous studies on cold-formed steel buckling failure, the slender compression member potentially caused buckling. Iman et al. conducted experimental test and numerical study on overall buckling failure on tubular steel pipe structure. This research compared the single compression member and the truss system structure buckling failure. Kang et al. conducted the numerical study on buckling behaviour of lipped and unlipped compression member of cold-formed steel. One of the research results showed that the lipped and unlipped canal profile in back-back connections had the greatest column capacity. Several studies on cyclic test that were carried on coldformed steel experimentally and numerically have been Wang et al. conducted experimental and numerical study on the hysteretic cold-formed steel truss composite floor behaviour. This research evaluated some dynamic parameters such as failure modes type, hysteretic curve, stiffness, and energy dissipation capacities. Arthy and Aruna . conducted an experimental study on cold-formed steel member by conducting monotonic and reversed cyclic test. The classic buckling load formula was derived by Leonhard Euler . 7Ae1. as shown as Equation . Pe = A 2 EI Paz . showed the dynamical properties such in frequency, periods, and mode shape can be determined by eigen value as shown as in Equation . |[K] Oe O 2 [M ]| {I} = 0 where [K] is material stiffness matrix. O is natural rotation frequency, [M ] is mass matrix, and I is mode shape from Equation . the fundamental frequency . ) and fundamental period (T ) can be determined by Equation . OT = 2A T = where A is 3. Chopra . stated that the hysteresis curve shows the initial loading curve is nonlinear at large deformation amplitudes, and the unloading and reloading curves differ from the initial loading branch. where Pe is buckling load. A is 3. E is modulus of elasticity. I is moment of inertia, and le is effective length. Buckling load can be evaluated in numerical analysis by performing computer calculations based on the finite element method. Two methods are used to evaluate buckling load: linear buckling analysis . igenvalue analysi. and nonlinear analysis (P O. CSI . pro2 where [K0 ] is material stiffness matrix. Pcr is critical load . uckling loa. , [N1 ] is geometric stiffness matrix first order, [N2 ] is geometric stiffness matrix second order, {Q0 } is displacement at minimum load, and [OIQ] is eigen vector . uckling mod. 2 METHODS This research is based on a numerical study that modeled the collapse pattern of cold-formed steel truss roof structures by considering elastic buckling failure and the seismic load capacity. The cold-formed steel roof truss structure 3D model is shown in Figure 1. Figure 1. shows the 2D model of cold-formed roof steel truss structure. Based on 2D-truss model, the ef- Vol. 12 No. 1 (January 2. Journal of the Civil Engineering Forum The material properties (Table . obtained from tension tests conducted according to ASTM E-8M. The cold-formed steel properties, such as Fy . Fu , and E were derived from the stress-strain curve using the 0. 2 strain offset method. The seismic zone in this research was selected, and the soil site class was classified as D for determining response spectrum acceleration coefficients, namely: SDS and SD1 . Based on equivalent linear force and modal analysis the lateral displacement, frequency, period, and mode shape were evaluated. 3D-Truss Structure Model . 2D-Truss Structure Model Figure 1 3D-Model Cold-Formed Roof Truss Structure. fects of gravity load and seismic load were compared and then verified by the 3D-truss model. In this case, the dead load was defined as the gravity load, consisting of the roof cover load, purlin, and self-weight. The equivalent linear force (ELF) analysis was conducted on both of models . D and 3D mode. , and lateral displacement and deformation were evaluated. The three-model of cold-formed roof steel truss were constructed using 3D truss elements for the chord The geometric and material properties of cold-formed steel truss structure are presented in Table 1. Table 1. Geometric and Material Properties Data Geometry. C75-0. Material Parameter Number Properties Number Cross section. A . m2 ) Yielding stress. Fy (MP. Moment of Inertia. I . m4 ) Ultimate stress. Fu (MP. Modulus of elasticity. E (MP. Radius of gyration, r . Figure 2 shows two numerical models of single compression member with the simple supports, representing an idealized bolt connection. Figure 2. shows the 3D frame idealization for single compression member of C75-0. On the other hand, the compression member was also modelled in 3D solid as shown as in Figure 2. and in the 3D-overall truss model idealization, as shown as in Figure 2. The three models were analyzed in linear . and nonlinear P Oe buckling analysis. The eigen value analysis was conducted to find out the lowest load that caused the member buckled, referred as critical load . uckling loa. Nonlinear static analysis was performed to evaluate the critical load (Pcr ) by determining the asymptote of nonlinear P Oe curve. The numerical buckling analysis was verified by 2D-Euler buckling theoretical analysis for a single simple supported compression member. This research also evaluated the effect of seismic load on the single C75-0. 75 cold-formed steel the compression member of truss structure. The compression member was identified as the slender element by considering its slenderness ratio. The single compression member C75-0. 75 was evaluated by modal analysis. The dynamical properties such in: frequency, period, and mode shape were evaluated in this research. The buckling mode shape was also compared with the mode shape of modal analysis. Additionally, the numerical simulation of cyclic loading on the single compression member was conducted in this research. Figure 2. shows the 3D lateral cyclic load on single compression member. The time-displacement was calculated for 45 cycles with displacement control of 2. This number was determined by the maximum displacement that was calculated in elastic condition . 75 m. The numerical cyclic load test on single compression member C75-0. 75 produced the hysteretic 3 RESULTS Based on Equation . the Euler theoretical critical load (Pe ) can be calculated (Pe = 28497. 43 N). The first buckling mode occurred in the x-direction because of the weaker axis of C75-0. 75 cold-formed steel section pro3 Journal of the Civil Engineering Forum . 3D-Frame . 3D-Solid Vol. 12 No. 1 (January 2. 3D-Overall Truss . Cyclic Model Figure 2 3D-Single Compression Member Model. The x-axis of the C75-0. 75 had the lowest radius of gyration . x = 12. 35 m. which produced the highest slenderness ratio. That is why buckling occurred in the x-axis of the C75-0. 75 cold-formed steel section profile. Theoretically, the buckling mode shape is a single curvature because of the boundary condition . imply supporte. Self-tapping screws had been utilized for the cold-formed steel truss structure connection systems. Hence, this is appropriate for assuming that the boundary condition of the cold-formed steel truss structure is a simply supported boundary condition idealization. The cold-formed steel truss structure had been modelled in 2D and 3D truss element by utilizing SAP2000. The cold-formed steel truss structure model had been analyzed by the finite element method to find out the free body diagram. In the case of a truss structure, the Figure 3 Check Design Structure. axial load is the most important factor to be evaluated. This research emphasized the compression axial member load of the cold-formed steel truss structure system. By considering the gravity load and based on the structure analysis, it showed that all top chord members are compression member load. The highest compression load is 37341. 42 N, which occurred on the edge top chord member of the cold-formed steel truss structure. On another hands, the cold-formed steel structure model also had been analyzed by considering the AISIASD 96. The gravity load combination had been involved in this case by assuming that the cold-formed steel truss structure only supported dead load . oof cover, purlins, and self weigh. The dead load was carried at every node of the roof truss structure top chord The service load . ctual loa. had been con- Vol. 12 No. 1 (January 2. Journal of the Civil Engineering Forum sidered as the load combination that was carried by the cold-formed steel truss structure. The results showed that some chord members were identified as Not Calculated (N/C), especially in the case of the top chord member and two adjacent vertical chord members as shown as in Figure 3. Figure 3 also showed that the top chord had been defined as an unsafe (N/C) chord member that wes identified in red color. Based on the analysis results, the top chord member is overstressed. That means the compression axial load member (Pa ) is larger than the nominal compression member load capacity (Pn ). In other words, the cross-section area of the C75-0. 75 coldformed steel profile is not adequate to support the service load that it was carrying. The N/C chord member might have been caused by an undefined alpha . lenderness rati. , so the design could not be completed for this section. The critical load (Pe ) was lower than the nominal compression member load capacity (Pa ). This situation had been validated by some numerical solutions for the single top chord member with buckling analysis as shown as in Figure 4. Figure 4 showed the first elastic buckling pattern of a single compression member C75-0. 75 consisting 3 models, namely: 3D-frame model, 3D-solid model, and 3D-overall truss model. The models had been analyzed by eigenvalue analysis and nonlinear P Oe analysis by performing the Riks method. Figure 4 also showed that the first mode of numerical elastic buckling pattern is a single curvature for the deformation. This satisfies the first Euler buckling mode shape . ingle curvatur. The nonlinear P Oe curve of C75-0. 75 that was produced by nonlinear static analysis. The critical load (Pe ) had been evaluated by determining the asymptotic of the nonlinear P Oe curve . 17 N). Figure 4. and Figure 4. showed the first buckling mode shape of the single compression member C750. The numerical model of the single C75-0. 75 had . 3D-Frame been idealized as 3D-frame in Figure 4. and a 3Dsolid idealization in Figure 4. The critical load had been recorded as 28232. 42 N . D-fram. and 26414 N . D-soli. by eigenvalue analysis . inear buckling analysi. Based on nonlinear analysis by utilizing the Riks method, the critical load (Pe ) was 25408 N for the 3Dtruss and 22390 N for the 3D-solid model of single compression member of C75-0. 75 idealization. The overall buckling failure also can be observed by modelling an overall 2D-truss system as shown as in Figure 4. The truss model consisted of two member types: tension member . and compression member . This model represented the coldformed steel truss structure for evaluating the overall buckling failure that occurred in the compression In this case, the observed C75-0. 75 compression member was the horizontal chord. Figure 4. also showed that the first mode buckling pattern is a single curvature with a critical load (Pe ) of 28640. 63 N by eigenvalue analysis. This number is larger than that calculated by EulerAos formula . 49%). The nonlinear analysis results also showed that the critical load (Pe ) is larger than EulerAos value for linear buckling analysis. However, this situation explains that buckling failure can occur in the C75-0. 75 compression element both in single form and in the overall truss structure system. Both numerical results can verify each other and explain that the buckling load is lower than the actual (Pa ) and limit load (Pu ) that occurred on the compression member, in this case, the C75-0. 75 cold-formed steel truss of top chord member. The critical load (Pe ) comparison can be shown in Table 2. Table 2 also showed that buckling stress is lower than yielding stress in the case of the cold-formed steel structure . MP. This means that buckling failure also occurred in the elastic range of the material. That is why the buckling failure is also called geometric failure. 3D-Solid . 3D-Overall Buckling Mode Shape Figure 4 First buckling Mode shape. Journal of the Civil Engineering Forum Vol. 12 No. 1 (January 2. Table 2. Critical Load (Pe ) Comparison C75-0. Critical Load (Pe ). N Euler 3D-Truss 3D-Solid 3D-Overall-Truss EV* NL** EV* * EV is eigen value analysis. ** NL is nonlinear analysis 4 DISCUSSION These results also showed the drift level of the coldformed steel truss structure. The truss structure had been defined in the steel system is not specifically detailed for seismic resistance. By considering the linear equivalent linear force for seismic load combination, the drift level can be evaluated for the 2D case, as shown as in Table 3. Based on Table 2 the 3D-overall truss model had a critical load number is closer to the Euler critical number . 59%). This number had been found by performing nonlinear analysis on a 3D-truss idealization that builds the 3D-overall truss structure model. Table 2 also showed that the 3D solid idealization of the C75-0. 75 single compression model using the nonlinear analysis with the Riks method had the lowest critical load . 90 N). The percentage different is 21. compared to the theoretical Euler critical load . 97 N). Several factors can be evaluated. Mesh discretization is one of the dominant factors that should be considered. Too many nodes need to be calculated which caused too many iterations in the 3D-solid idealization. Hence, mesh discretization should consider the memory capacity of the computer. In this case, the load combination for the equivalent linear seismic force load combination is 0. 9D Ex for the 2D-cold formed steel truss structure model where. D is the dead load and Ex is the equivalent linear force seismic load in the x-direction. The seismic load combination refers to ASCE for the limit state condition. Table 3 also showed the drift level of the edge vertical chord did not satisfied the allowable the drift level of the vertical chord . ,98 m. This means that the edge vertical chord needs a bracing system for the seismic Table 2 also showed that the critical load obtained by performing eigenvalue analysis is close to the exact solution for the elastic buckling load that was calculated by EulerAos formula. The eigenvalue buckling analysis is a linear analysis. Hence, it doesnAot give any information on the post-buckling behaviour of the structure. However, in certain cases, it can calculate with a good approximation the actual buckling load factor of the structure. The percentage different between eigenvalue analysis on Euler buckling formula is, successively: 0. 93% . D truss-single elemen. 50% . truss-overall elemen. 31% . D solid-single elemen. Practically, for conservative reasons, engineers should consider the nonlinear P Oe curve from 3D solid Table 3 also showed that displacement increased on: the 4th level . 66 m. , 3rd level . 55 m. , and 2nd level . 99 m. of the vertical chord by load combination 9D Ex . The situation showed that the equivalent linear force seismic load combination had a significant effect on the cold-formed steel truss structure, especially in the case of end supported failures. This can be understood as the truss structure not having been designed for lateral loads such as seismic loads. In other words, the cold-formed steel truss structure has a high vulnerability to collapse when a seismic load occurs. Furthermore, the slenderness ratio is the most important parameter that causes the structure instability. Table 3. The Drift Level of Cold-Formed Steel Truss Structure . D-Cas. Displacement OI . OIcd /l OIn . Remarks Safe Safe Safe Unsafe Vol. 12 No. 1 (January 2. Journal of the Civil Engineering Forum The single C75-0. 75 compression member also had been evaluated by performing modal analysis. The frequency had been recorded in Table 4. Table 4. Frequency. Period, and Buckling Load Mode Buckling Load (N) Figure 5 Critical Load C75-0. 75 Curve Comparison. The linear buckling analysis . is usually used for preliminary buckling load prediction. Figure 5 showed the comparison of the C75-0. 75 critical load (Pe ) that was plotted as curve. The Euler formula and the eigen value analysis for the critical load had been presented as horizontal line curves, which lie above the non-linear curve. In other words, the Euler formula and eigenvalue analysis provides the highest critical load estimation. Figure 5 also showed the nonlinear load-displacement curve which consisted of 3Dsingle truss, 3D-single solid, and 3D-overall truss idealizations. The critical load is the asymptotic of the nonlinear load-displacement curves that were generated by the Riks method, as shown in Table 2. In this case, the nonlinear critical load of the 3Dsingle solid idealization is the lowest critical load value. On the other hand, the highest nonlinear critical load had been shown by the 3D-overall truss idealization . 50 N). However, this number is lower than the highest compression load . 42 N). This situation explains that the C75-0. 75 compression member experienced buckling failure even though the compression member is part of a steel truss configuration. Figure 6 showed the comparison of the buckling pattern and modal analysis mode shape at the 1st mode of the C75-0. 75 single compression member model. Figure 6. showed the mode shape of the 3D-truss single compression member model . = 92. 65 H. , and Figure 6. showed the mode shape of 3D-solid single compression member model . = 89. 96 H. by modal Both mode shapes had similar mode shapes to the first buckling mode shape, as shown in Figure 4. The mode shape is a single curvature mode shape. However, each had different parameters that were involved in the analysis. Specifically: geometric stiffness matrix for buckling analysis and mass matrix for modal analysis. Figure 6 also showed the buckling load-frequency comparison curve. Frequency (H. 3D-Truss 3D-Solid Table 4 also showed the comparison of dynamical properties . and the buckling load values for 6 modes shape for the single compression member, in 3D-truss and 3D-solid idealization. As the buckling load increased, the frequency also increased. Conversely, the period decreased when the buckling load Table 4 also showed that the frequency of 3D-truss and 3D-solid idealization are similar with an average percentage different of 6. However, in a fact, the buckling load in the first mode occurred early. The buckling modes shape is related to the modal analysis modes shape. This is because the two types of modes shape are eigen vectors that were produced by eigen problems analysis. The eigen problem analysis had been developed to find out the lowest critical load that caused instability as shown in Equation . On the other hand, the dynamical properties such as frequency and period also had been produced by eigen value analysis as shown in Equation . Hence, the buckling mode and modes shape of modal analysis had been determined by nontrivial solution to find out the mode shape vector. However, the key difference is that buckling analysis does not involved the mass matrix in the analysis. There is no requirement for 90% mass absorption in the structure as is the case in modal This research also performed a numerical cyclic test on the 3D-solid C75-0. 75 single compression member By performing lateral displacement control at one boundary condition (Ux = 2. 00 m. the hysteretic curve as shown in Figure 7, can be evaluated. Figure 7. showed the numerical lateral cyclic deformation of the 3D solid C75-0. 75 model as a representative top chord member. Figure 7. also showed the stress range was recorded between 358 MPa Ae 509 MPa. Journal of the Civil Engineering Forum . 3D-Truss Vol. 12 No. 1 (January 2. 3D-Solid . Buckling Load-Frequency Comparison Curve Figure 6 Buckling Mode and Mode Shape Modal Analysis for Single Compression Model. This stress range showed that the stress was in elastic Figure 7. showed the hysteresis curve of the single compression cold-formed truss member (C750. The hysteresis curve had been produced by a numerical cyclic test on the 3D solid C75-0. 75 single compression model. The analysis had involved 45 cyclic repetitions in the lateral direction. Figure 7. also showed that the hysteresis curve had a maximum axial load that was evaluated by the displacement control method. Based on the hysteretic curve the maximum axial compression load is 19159. N with a displacement of 25. 2 cm for the displacement. This number is lower than the minimum buckling load . 90 N) that was evaluated by the nonlinear P Oe curve of the 3D-solid idealization. This means that the seismic failure in the C75-0. 75 single compression member occured earlier than the buckling failures. other words, we can say that the lateral loads, such as seismic loads, can trigger buckling failure immediately. 5 CONCLUSION Based on the numerical studies, some points can be concluded below: The C75-0. 75 cold-formed steel profile had not satisfied the compression member design requirements (Pe < Pn ). The C75-0. 75 had a very high vulnerability to buckling failure hazards, especially at the top chord of the cold-formed steel roof truss structure. Numerical solutions for buckling load had results that are close to the theoretical buckling load calculated with EulerAos formula . 9% - 21%). In the case of the C75-0. 75 cold-formed steel compression member the 3D-Truss model with eigenvalue analysis is more appropriate than other numerical models for buckling . Lateral Cyclic Mode Figure 7 3D-Solid Compression Member Deformation by Lateral Cyclic Test. Modal Analysis Mode Shape Vol. 12 No. 1 (January 2. The single compression member buckling load can be verified by comparing it with the overall buckling load . 43%). The frequency was observed in the first mode as 92. 65 Hz in the 3D-truss model and 96 Hz in the 3D-solid model through modal analysis. The first six buckling modes of the single compression element corresponded to the first six mode shapes of modal analysis. As the frequency increased, the buckling load also increased. Based on the equivalent linear force seismic load combination analysis with the load combination . Ex ), the displacement caused the largest drift level of the vertical chord near the boundary condition. This means the equivalent linear force seismic load combination had a significant effect on the cold-formed steel truss structure, potentially causing collapse when seismic load occurred. By considering the numerical cyclic load analysis, the stress range was recorded to be about 358 MPa Ae 509 MPa . ithin the elastic rang. Based on the hysteresis curve, the maximum compression load 50 N with a displacement of 25. 2 cm. This means that when seismic load occurred, the C75-0. cold-formed steel profile as a compression member can become unstable is potentially buckling failure. This number is lower than the numerical buckling load number . 90 N) that was evaluated by the nonlinearP Oe curve using the Riks method. The situation showed that cold-formed steel is highly vulnerable to seismic In other words, lateral seismic loads can trigger buckling failure immediately, especially in cases involving compression member with high slenderness ratios. DISCLAIMER The authors declare no conflict of interest. ACKNOWLEDGMENTS The authors would like to thank to the Institute of Research and Community Outreach. Universitas Muhammadiyah Yogyakarta, for funding and supporting the research. REFERENCES