Matrix Description of Non-Linear Properties of Materials or Structural Components—Idea and Application Examples

Numerical methods are widely used in structural analysis problems. In the cases of the most complex and practical problems, they are often the only way to obtain solutions, as analytical methods prove ineffective. The motivation for this paper was the desire to extend the scope of numerical methods to cover the problems of creating constitutive models of structural materials. The aim of this research was to develop a matrix or numerical discrete constitutive model of materials. It presents the general assumptions of the developed method for modeling the physical properties of materials. The matrix model is only useful with an appropriate numerical algorithm. Such an algorithm was created and described in this paper. Based on its findings, computer software was developed to perform numerical simulations. Presented calculation examples confirmed the effectiveness of the developed method to create constitutive matrix models of various typical materials, such as steel, but also, e.g., hyper-elastic materials. It also presents the usefulness of constitutive matrix models for simulations of simple stress states and analyses of structural elements such as reinforced concrete. All presented examples involved the physical nonlinearity of the materials. It is proved that the developed matrix constitutive model of materials is efficient and quite versatile. In complex analyses of structures made of nonlinear materials, it can be used as an effective alternative to classical constitutive or analytical models based on elementary mathematical functions.


Introduction
Structural materials used in the construction industry feature a number of parameters that are important primarily in terms of material strength and structure mechanics, but also in terms of other parameters related to, for example, building physics. One of the basic parameters is the material stiffness that is necessary to determine displacements and forces in structural systems.
Basic construction materials (steel and concrete in particular) behave in a non-linear way through the whole operating range. This means that their stiffness changes along with the increase of strain and stress. In the process of designing steel structures, the amount of stress or strain is limited to a certain acceptable level. This level is related to material specification (e.g., a clear yield point) or it is a defined arbitrary limit. This approach means that steel structures almost always operate in constant stiffness range and it is justified to assume the constant material stiffness in static calculations. The situation is completely different in case of reinforced concrete structures. Here, cracks occur very quickly (at a relatively low level of internal forces) in the tensile zones of concrete crosssections. As a result, tensile stresses in the cross-sections of bent reinforced concrete structural elements are mainly transferred with reinforcement bars, and compressive stresses are mainly transferred with concrete. When internal forces increase, concrete or steel becomes plasticized, i.e., it behaves in a non-linear way. Yet, most static analyses of reinforced concrete structures are carried out in a simplified way, assuming that the geometric characteristics of cross-sections correspond to an uncracked concrete crosssection (without any reinforcement). This primarily applies to engineering calculations. It The discussed studies were mainly related to two issues: non-linear analyses of construction elements (mainly reinforced concrete) and non-linear models of hyper-elastic constitutive materials and their utilization. In some cases, the cyclic load of the system that usually generates a hysteretic response is taken under consideration ( [6,10,11,21,30,31]). The studies feature a common theme that is creating models of constitutive materials or elements-using algebraic or differential mathematical relations describing the stressstrain or moment-curvature relation. Whether models in analyses are correct was most often verified by comparison with experiment results or other independent numerical calculations of the same issue. The desired situation was the high compatibility of the obtained results with the experiment, which could not always be achieved. Moving towards the increased compatibility of the constitutive model and actual material behavior during the experiment, Sussman and Bathe [32] proposed an alternative approach to describe the physical properties of the material. They presented the concept of constructing a constitutive relation in the form of a spline function, consisting of segments of a cubic spline. This enables us to achieve a very high consistency with the results of experimental research. The Sussman and Bathe concept was used in other authors' studies [33,34].
This paper presents the concept of the matrix description of material properties, e.g., constitutive relations. This concept is similar to that of Sussman and Bathe in one pointinterpolation approach in both cases is used to ensure high consistency of constitutive relation with the experiment. The difference is that Sussman and Bathe used the interpolation of a one-variable function and the study uses the approximation of a two-variable function. The utilization of two-variable function is the most important difference between the presented solution and those used so far in other studies. The stress-strain (or momentcurvature for bending elements) relation is usually the basis for determining the material stiffness. The stiffness determined this way is a function of one independent variable, i.e., stress or strain. This paper proposes to treat stiffness as a function of two variables: stress and strain. This description of stiffness allows us to include more information about the material, e.g., partial history of previous strain or strain states.
The first part of this paper presents general assumptions for the description of stiffness using matrices and discusses the numerical algorithm to use stiffness matrices for analyses of the displacement and strain of a structure under load. The second part contains calculation examples and some possibilities of describing stiffness using a matrix are presented. The summary of the paper outlines the possible development of the presented concept and potential fields of its use.
The motivation for this paper was an interest in numerical methods and a desire to broaden their application scope in order to create alternative constitutive models of building materials. The main objective of this paper was to present the developed method for modeling the physical characteristics of materials. It is a discrete, matrix-described model. An appropriate numerical algorithm was developed to demonstrate the model effectiveness and some application examples. It was used in the presented calculation examples. The developed method of creating numerical constitutive models is alternative to classical, analytical approaches. No similar approach presented in this paper has been found in the scientific literature. Matrix modeling of constitutive features of materials has both advantages and limitations, compared to classical approaches. It has a potential for researchers involved in nonlinear structural analysis, e.g., made of hyper-elastic materials.

The Matrix Description of Material Stiffness Concept
In order to improve the commonly used approach to describe the properties of building materials, the concept of describing the coefficient of elasticity E as a function of two variables was developed, i.e., strain ε and stress σ. Different values of elasticity coefficients for load and unload may occur in a given state of the body defined by values (σ, ε), therefore the description of material stiffness requires two functions, which are indicated as E p and E n . E p (σ, ε) is a function of Young's modulus at primary load cycle, E n (σ, ε) Materials 2021, 14, 5837 4 of 25 is a function for the unload and secondary load cycle. These functions are related to tangential modules.
Function E(σ, ε) is presented in the form of a discrete array of values that can be conveniently used in computer calculations. This value table is called the stiffness matrix in this paper. The graphical representation of the idea of such a matrix structure is presented in Figure 1.
variables was developed, i.e., strain ε and stress σ. Different values of elasticity cients for load and unload may occur in a given state of the body defined by values therefore the description of material stiffness requires two functions, which are ind as E and E . E (σ, ε) is a function of Young's modulus at primary load cycle, E is a function for the unload and secondary load cycle. These functions are related gential modules.
Function E(σ, ε) is presented in the form of a discrete array of values that conveniently used in computer calculations. This value table is called the stiffness in this paper. The graphical representation of the idea of such a matrix structure sented in Figure 1. In Figure 1, both the X axis (ε values) and Y axis (σ values) were digitized, b division does not have to be uniform. X-axis values are collected in the X vector: X = ε , ε , … , ε , … , ε and Y-axis values-in the Y vector: The graphical interpretation of Young's tangential moduli for the primary loa is shown with solid blue arrows pointing upwards, and for the unload and seconda cycles-with red dashed arrows pointing downwards. Arrow directions are tangen possible σ(ε) graphs. The general form of stiffness matrices is the following: Matrix values describe the following relations: In Figure 1, both the X axis (ε values) and Y axis (σ values) were digitized, but the division does not have to be uniform. X-axis values are collected in the X vector: and Y-axis values-in the Y vector: The graphical interpretation of Young's tangential moduli for the primary load cycle is shown with solid blue arrows pointing upwards, and for the unload and secondary load cycles-with red dashed arrows pointing downwards. Arrow directions are tangential to possible σ(ε) graphs. The general form of stiffness matrices is the following: Matrix values describe the following relations: The α p ij and α n ij angles are shown in Figure 2, that is a fragment of the graph presented in Figure 1. The α and α angles are shown in Figure 2, that is a fragment of the gra sented in Figure 1. Stiffness matrices created according to discussed principles are used later in per.

Calculation Algorithm That Utilizes Stiffness Matrices
Stiffness matrices contain tabular values of tangential moduli and are theref able for use primarily in incremental methods of solving non-linear issues. In t two basic calculation problems arise. The first one is to determine the value of th modulus at any point outside the nodes. The second problem is related to the det tion of the modulus value in a given calculation step (including increment).
A node should be understood as a point defined with a pair of values σ occur in vectors Y and X for which the modulus value E is given in the stiffnes ces. In general, it is necessary to set modulus values for any points. For a selected defined with a pair of values (σ , ε ), the determination of the modulus value quires to run the algorithm in two steps. The first step is to establish the value of Y and X between which the point A lies. This step is just to determine such value j that σ ≤ σ ≤ σ and ε ≤ ε ≤ ε . The second step is to determine the m value at the A point. The principle is shown in Figure 3.  Stiffness matrices created according to discussed principles are used later in this paper.

Calculation Algorithm That Utilizes Stiffness Matrices
Stiffness matrices contain tabular values of tangential moduli and are therefore suitable for use primarily in incremental methods of solving non-linear issues. In this case, two basic calculation problems arise. The first one is to determine the value of the Young modulus at any point outside the nodes. The second problem is related to the determination of the modulus value in a given calculation step (including increment).
A node should be understood as a point defined with a pair of values σ i , ε j that occur in vectors Y and X for which the modulus value E ij is given in the stiffness matrices. In general, it is necessary to set modulus values for any points. For a selected A point, defined with a pair of values (σ A , ε A ), the determination of the modulus value E A requires to run the algorithm in two steps. The first step is to establish the value of vectors Y and X between which the point A lies. This step is just to determine such values i and j that σ i ≤ σ A ≤ σ i+1 and ε j ≤ ε A ≤ ε j+1 . The second step is to determine the modulus value at the A point. The principle is shown in Figure 3.
The α and α angles are shown in Figure 2, that is a fragment of t sented in Figure 1. Stiffness matrices created according to discussed principles are used la per.

Calculation Algorithm That Utilizes Stiffness Matrices
Stiffness matrices contain tabular values of tangential moduli and are able for use primarily in incremental methods of solving non-linear issue two basic calculation problems arise. The first one is to determine the value modulus at any point outside the nodes. The second problem is related to t tion of the modulus value in a given calculation step (including increment) A node should be understood as a point defined with a pair of value occur in vectors Y and X for which the modulus value E is given in the s ces. In general, it is necessary to set modulus values for any points. For a sele defined with a pair of values (σ , ε ), the determination of the modulus quires to run the algorithm in two steps. The first step is to establish the va Y and X between which the point A lies. This step is just to determine such j that σ ≤ σ ≤ σ and ε ≤ ε ≤ ε . The second step is to determine value at the A point. The principle is shown in Figure 3. Taking into account the markings shown in Figure 3 and additionally Taking into account the markings shown in Figure 3 and additionally defining: the requested value is determined with the following relation: As presented in Figure 3 and Formula (8), linear interpolation was used to determine the modulus value outside the nodes.
When determining the elasticity modulus in a given calculation step, it is assumed that the value of stress increase ∆σ is known. The determination of its corresponding ∆ε value, and thus the value of Young's modulus in a given step, is done with successive approximations. The idea of this process is presented in Figure 4.
As presented in Figure 3 and formula (8), linear interpolation was used to de the modulus value outside the nodes.
When determining the elasticity modulus in a given calculation step, it is a that the value of stress increase ∆σ is known. The determination of its correspon value, and thus the value of Young's modulus in a given step, is done with su approximations. The idea of this process is presented in Figure 4. Let us assume that stiffness matrix, point A coordinates and increment value given. The value of the stiffness modulus at the A point shall be determined acco formula (8). During the first, initial approximation, angle α is assumed so that t E and the position of the B point is calculated. Then, with a second approximat E modulus value is determined (i.e., value at the B point) and the α calculate is that tg(α ) = 0.5(E + E ). This enables the calculation of the position of the B and with similar actions, as in case of α , determine the α angle. Such calculat performed to achieve assumed α angle accuracy and thus the accuracy of the Y modulus. The value of the modulus can be determined in the same way if the st crease ∆ε in one step is known.
The basic elements of a calculation algorithm to use the stiffness matrix are ve ple, and based on elementary mathematical operations. However, the algorithm s ity should be considered as its advantage. Figure 5 shows a general basic block diagram of the calculation algorithm. It i plified diagram that focuses on activities relevant to the issues at hand. The activit formed in the operation block marked with a double box are shown in Figure 6 refer to the issues described in Section 2.2. Let us assume that stiffness matrix, point A coordinates and increment value ∆σ are given. The value of the stiffness modulus at the A point shall be determined according to Formula (8). During the first, initial approximation, angle α 0 is assumed so that tg(α 0 ) = E A and the position of the B 0 point is calculated. Then, with a second approximation, the E B0 modulus value is determined (i.e., value at the B 0 point) and the α 1 calculated angle is that tg(α 1 ) = 0.5(E A + E B0 ). This enables the calculation of the position of the B 1 point and with similar actions, as in case of α 1 , determine the α 2 angle. Such calculations are performed to achieve assumed α angle accuracy and thus the accuracy of the Young's modulus. The value of the modulus can be determined in the same way if the strain increase ∆ε in one step is known.

Flowchart of the Calculation Algorithm
The basic elements of a calculation algorithm to use the stiffness matrix are very simple, and based on elementary mathematical operations. However, the algorithm simplicity should be considered as its advantage. Figure 5 shows a general basic block diagram of the calculation algorithm. It is a simplified diagram that focuses on activities relevant to the issues at hand. The activities performed in the operation block marked with a double box are shown in Figure 6. These refer to the issues described in Section 2.2.   For each node of the discretized structure: For each node of the discretized structure: Figure 6. A flowchart that presents the stiffness analysis within the main algorithm.

Results
The first part was limited to the presentation of simulation examples of relatively simple sample tests made of an abstract elastic-plastic material. These tests are the axial tensile stress bar test, cyclic tensile and relief stress bar test, and the tensile stress bar test with relief, compression and another relief stage. During the second part, the tension of three bars connected with a rigid beam was simulated and the reinforced concrete beam was calculated.
For the purpose of the presented simulations, a hypothetical material was adopted, the Young's modulus of which is described with the values of dimensionless deformations ε collected in the X vector: and two stiffness matrices, expressed in (MPa): Reading ∆E . Reading system state (σ , ε ) at the beginning of a specific step. This state corresponds to the position of the A point according to Figure 3. Reading the value of E based on the relationship (8).
Reading the value E .
Determining the system state in relation to point B (σ + ∆σ , ε + E ∆σ ). Reading the value E .
no Inserting E value in the following calculations yes Figure 6. A flowchart that presents the stiffness analysis within the main algorithm.

Results
The first part was limited to the presentation of simulation examples of relatively simple sample tests made of an abstract elastic-plastic material. These tests are the axial tensile stress bar test, cyclic tensile and relief stress bar test, and the tensile stress bar test with relief, compression and another relief stage. During the second part, the tension of three bars connected with a rigid beam was simulated and the reinforced concrete beam was calculated.
For the purpose of the presented simulations, a hypothetical material was adopted, the Young's modulus of which is described with the values of dimensionless deformations ε collected in the X vector:

Example A. Axial Bar Tension
The first analyzed issue is relatively simple. Its solution is to be used primarily to study the course of the σ(ε) function. It is clear that in case of analyzing real issues, laboratory tests must first be performed to identify the σ(ε) relation, and then appropriate stiffness matrices E p and E n must be created. However, stiffness matrices themselves and graphs based on them are not directly represented in the σ(ε) relation. Only the laboratory sample simulation, e.g., the axial bar tension, will allow the corresponding σ(ε) graph to be depicted.
It is assumed that a bar that is stretched has a length L = 2.0 m and a cross-section area A = 0.01 m 2 . The bar material is described with matrices (11) and (12). Tension is static and continues until a tensile force of 105kN is reached.
To solve the issue of the tensile bar simulation, the algorithm discussed in the previous chapter was used. The algorithm was the basis for a computer program written in the Scilab environment. The selection of the ∆σ step (see Figure 4) was carried out adaptively so that the transition of the graph through one cell of the stiffness matrix (the field separated with two adjacent values from the X vector and two adjacent values from the Y vector) took place in about four steps.
The graph based on the simulation is shown in Figure 7. Since a continuous stress increase was present during the simulation, only the E p matrix was used in the calculations. (1

Example A. Axial Bar Tension
The first analyzed issue is relatively simple. Its solution is to be used primarily study the course of the σ(ε) function. It is clear that in case of analyzing real issues, l boratory tests must first be performed to identify the σ(ε) relation, and then appropria stiffness matrices E and E must be created. However, stiffness matrices themselve and graphs based on them are not directly represented in the σ(ε) relation. Only the l boratory sample simulation, e.g., the axial bar tension, will allow the corresponding σ( graph to be depicted. It is assumed that a bar that is stretched has a length L = 2.0 m and a cross-sectio area A = 0.01 m 2 . The bar material is described with matrices (11) and (12). Tension static and continues until a tensile force of 105kN is reached.
To solve the issue of the tensile bar simulation, the algorithm discussed in the prev ous chapter was used. The algorithm was the basis for a computer program written in th Scilab environment. The selection of the ∆σ step (see Figure 4) was carried out adaptive so that the transition of the graph through one cell of the stiffness matrix (the field sep rated with two adjacent values from the X vector and two adjacent values from the vector) took place in about four steps.
The graph based on the simulation is shown in Figure 7. Since a continuous stre increase was present during the simulation, only the E matrix was used in the calcul tions. The graph presented in Figure 7 shows the material behavior, described with X ve tors according to (9), Y vectors according to (10) and the E matrix given in (11).
Based on experimental data, it is quite simple to create a suitable stiffness matri This will be presented in the next example. The graph presented in Figure 7 shows the material behavior, described with X vectors according to (9), Y vectors according to (10) and the E p matrix given in (11).
Based on experimental data, it is quite simple to create a suitable stiffness matrix. This will be presented in the next example.

Example B. Stiffness Matrix Calibration with Axial Tensile Bar Test Simulation
The following example shows the results of stiffness matrix calibration. It was created manually by introducing subsequent corrections in the matrix values and observing the obtained effects. Let us assume that the tensile bar test produced the F(∆L) relation graph. Some characteristic points are indicated in this graph and then the σ stresses corresponding to the F forces and the ε strains corresponding to ∆L elongations are calculated for these points.
Values in theε vector are dimensionless and values in theσ vector are expressed in (MPa). Additionally, the graph in the first section (between points 1 and 2) should be linear.
The vectors X, Y and the E p matrix were produced with successive value adjustments to finally obtain the following: Values in the X vector are dimensionless, and values in the Y vector and the E p matrix are expressed in (MPa). For a specific material stiffness determined with this method, the axial tensile test was simulated until a stress value of 0.9 MPa was reached. A computer program created for example A was used. Figure 8 shows an achieved tensile simulation graph for the calibrated stiffness matrix (solid line). Points defined withε andσ vectors are marked with crosses. It can be seen that the tensile simulation graph matches the base points quite well. The shape of the resulting graph corresponds to a material with a distinct yield strength, e.g., low-carbon steel. The following example shows the results of stiffness matrix calibration. It was created manually by introducing subsequent corrections in the matrix values and observing the obtained effects. Let us assume that the tensile bar test produced the F(∆L) relation graph. Some characteristic points are indicated in this graph and then the σ stresses corresponding to the F forces and the ε strains corresponding to ∆L elongations are calculated for these points. Finally, the σ(ε) graph should go through the vector points ε = 0 0.002 0.0022 0.0031 0.004 0.006 0.009 0.014 , Values in the ε vector are dimensionless and values in the σ vector are expressed in (MPa). Additionally, the graph in the first section (between points 1 and 2) should be linear.
The vectors X, Y and the E matrix were produced with successive value adjustments to finally obtain the following: . (17) Values in the X vector are dimensionless, and values in the Y vector and the E matrix are expressed in (MPa). For a specific material stiffness determined with this method, the axial tensile test was simulated until a stress value of 0.9 MPa was reached. A computer program created for example A was used. Figure 8 shows an achieved tensile simulation graph for the calibrated stiffness matrix (solid line). Points defined with ε and σ vectors are marked with crosses. It can be seen that the tensile simulation graph matches the base points quite well. The shape of the resulting graph corresponds to a material with a distinct yield strength, e.g., low-carbon steel.
These are the experimental data read from Figure 6b in the paper [23].
In this case, stiffness is described with the following vectors and the matrix The graph shown in Figure 9 was obtained with simulating the tensile test. A high degree of compatibility with the base data was also achieved. The shape of the graph is similar to the ones corresponding to hyper-elastic materials that offer a stiffness increase at high deformations. Rubber is a good example of such material.
In the second example, the same actions were performed as before, assuming th lowing set of base data: These are the experimental data read from Figure 6b in the paper [23]. In this case, stiffness is described with the following vectors The graph shown in Figure 9 was obtained with simulating the tensile test. A degree of compatibility with the base data was also achieved. The shape of the gr similar to the ones corresponding to hyper-elastic materials that offer a stiffness in at high deformations. Rubber is a good example of such material.  Figure 6b in Ref. [23] were used).
In both examples, material stiffness matrices were developed to reproduce hyp ical laboratory tensile sample tests of different materials with high accuracy (eval visually, based on graphs). It should be noted that there were different incremen tween the values of X and Y vectors. As a consequence, the ∆σ increments used algorithm are adjusted accordingly. This means that the value density in X and Y v Strech Nominal Stress (MPa) Figure 9. Relation graph achieved for the calibrated stiffness matrix in case 2 (the axis designations used in Figure 6b in Ref. [23] were used).
In both examples, material stiffness matrices were developed to reproduce hypothetical laboratory tensile sample tests of different materials with high accuracy (evaluated visually, based on graphs). It should be noted that there were different increments between the values of X and Y vectors. As a consequence, the ∆σ increments used in the algorithm are adjusted accordingly. This means that the value density in X and Y vectors in the zones with more rapid changes in the elasticity modulus automatically caused a decrease in ∆σ step and thus better algorithm accuracy.
In both cases, the values in vectors X and Y are mostly repetitions of the base data set in vectorsε andσ. Such a relationship is fairly obvious, but not a requirement. Selection of values in vectorsε andσ depends on the person that builds the stiffness matrix of a given material. It must be done in such a way that it closely reproduces the tension graph obtained in the laboratory. In zones of large changes of the inclination angle that is tangent to the graph, the value drops in vectorsε andσ have to be small-this is clearly visible in Figure 8. Incorrect selection of base values in vectorsε andσ will lead to the creation of an inadequate constitutive model. A solution to the problem would also be to use small increments in complete vectorsε andσ, but it would lead to a significant size increase of the stiffness matrix.

Example C. Axial Cyclic Tension and Relief Stress Bar Test
This example shows a simulation of the tension and relief stress test of the bar. The subject of the analysis will be exactly the same bar that is made of the same material as in example A (Section 3.1). This time, however, it was assumed that the bar would be cyclically stretched to a certain force value and then relieved. During individual cycles, tensile forces of 95, 100, 103 and 105 kN, respectively, were assumed to be achieved. The obtained tension graph is shown in Figure 10.
in the zones with more rapid changes in the elasticity modulus automatically caus decrease in ∆σ step and thus better algorithm accuracy.
In both cases, the values in vectors X and Y are mostly repetitions of the base set in vectors ε and σ . Such a relationship is fairly obvious, but not a requirement. S tion of values in vectors ε and σ depends on the person that builds the stiffness m of a given material. It must be done in such a way that it closely reproduces the ten graph obtained in the laboratory. In zones of large changes of the inclination angle th tangent to the graph, the value drops in vectors ε and σ have to be small-this is cl visible in Figure 8. Incorrect selection of base values in vectors ε and σ will lead t creation of an inadequate constitutive model. A solution to the problem would also use small increments in complete vectors ε and σ , but it would lead to a significant increase of the stiffness matrix.

Example C. Axial Cyclic Tension and Relief Stress Bar Test
This example shows a simulation of the tension and relief stress test of the bar subject of the analysis will be exactly the same bar that is made of the same material example A (Section 3.1). This time, however, it was assumed that the bar would be c cally stretched to a certain force value and then relieved. During individual cycles, te forces of 95, 100, 103 and 105 kN, respectively, were assumed to be achieved. The obta tension graph is shown in Figure 10. In this case, only the E matrix was used during the first load cycle. During the and another load stage to the previously achieved stress level, the values were taken the E matrix. When applying load above the previously achieved stress level, th matrix was reused.
The graph obtained in Figure 10 corresponds to the model behavior of such mate as steel. Extreme upper parts of the graph from Figure 10 should match the graph sented in Figure 7. The confirmation of this rule was obtained by comparing both gr in Figure 11. In this case, only the E p matrix was used during the first load cycle. During the relief and another load stage to the previously achieved stress level, the values were taken from the E n matrix. When applying load above the previously achieved stress level, the E p matrix was reused.
The graph obtained in Figure 10 corresponds to the model behavior of such materials as steel. Extreme upper parts of the graph from Figure 10 should match the graph presented in Figure 7. The confirmation of this rule was obtained by comparing both graphs in Figure 11.
The graph shown in Figure 10 shows that the lines corresponding to reliefs and consecutive repeated loads in successive cycles are parallel to one another. This is obvious, as only the E n matrix with constant coefficients is used in the simulation during these phases. Appropriate calibration of this matrix may result in a change of the slope in graph sections related to relief and another application of load. Another method is to change E n stiffness matrix coefficients during the calculation. Figure 12 shows an example graph of cyclic tension and relief with such matrix modification. The graph shown in Figure 10 shows that the lines corresponding to reliefs a secutive repeated loads in successive cycles are parallel to one another. This is obv only the E matrix with constant coefficients is used in the simulation durin phases. Appropriate calibration of this matrix may result in a change of the slope i sections related to relief and another application of load. Another method is to ch stiffness matrix coefficients during the calculation. Figure 12 shows an example g cyclic tension and relief with such matrix modification. In the presented example, after the sample material exceeded certain defo limits (0.07 in this case), the matrix was multiplied by E , by a coefficient who was a certain function of achieved deformation. Such a mechanism can be used, fo ple, to build degradation models.  The graph shown in Figure 10 shows that the lines corresponding to reliefs a secutive repeated loads in successive cycles are parallel to one another. This is obv only the E matrix with constant coefficients is used in the simulation durin phases. Appropriate calibration of this matrix may result in a change of the slope i sections related to relief and another application of load. Another method is to cha stiffness matrix coefficients during the calculation. Figure 12 shows an example g cyclic tension and relief with such matrix modification. In the presented example, after the sample material exceeded certain defor limits (0.07 in this case), the matrix was multiplied by E , by a coefficient whos was a certain function of achieved deformation. Such a mechanism can be used, fo ple, to build degradation models.

Example D. Sample Tension and Compression-Hysteretic Material Response
This example shows the simulation of a tensile sample test described in exa (Section 3.1) with a force of 107 kN, followed by a relief, compression with a forc kN and final relief. In order to exclude the possibility of stability loss by the com sample, this time it was assumed that the sample length will be 0.2 m.
The simulation required to extend the stiffness matrix. First, the vectors w panded: X with the relation (9) and Y with the relation (10) that add specific valu metrically to the 0.00 values and then change their signs to negative. In case of In the presented example, after the sample material exceeded certain deformation limits (0.07 in this case), the matrix was multiplied by E n , by a coefficient whose value was a certain function of achieved deformation. Such a mechanism can be used, for example, to build degradation models.

Example D. Sample Tension and Compression-Hysteretic Material Response
This example shows the simulation of a tensile sample test described in example A (Section 3.1) with a force of 107 kN, followed by a relief, compression with a force of 107 kN and final relief. In order to exclude the possibility of stability loss by the compressed sample, this time it was assumed that the sample length will be 0.2 m.
The simulation required to extend the stiffness matrix. First, the vectors were expanded: X with the relation (9) and Y with the relation (10) that add specific values symmetrically to the 0.00 values and then change their signs to negative. In case of the E p matrix given with a Formula (11), symmetrization relative to the first column (containing values of 200) was performed and an appropriate number of repeating rows was also added. The E n matrix has been expanded to appropriate size and filled with 200 values (see relation (12)). For the data prepared in this way, a simulation was performed that produced the σ(ε) relation graph shown in Figure 13. matrix given with a formula (11), symmetrization relative to the first column values of 200) was performed and an appropriate number of repeating row added. The E matrix has been expanded to appropriate size and filled with (see relation (12)). For the data prepared in this way, a simulation was perf produced the σ(ε) relation graph shown in Figure 13. As it can be seen in Figure 13, a hysteresis loop graph was obtained. The sh from the adoption of the symmetry of values in the E matrix. Precise loop clo (0.0) point additionally indicates high accuracy of the algorithm used.

Example E. Tension and Relief of the Three Bar System
The next analyzed issue is related to the tension of the system of three con shown in Figure 14. The bars have the same length L and the distances between them are co are indicated as a. It is assumed that the upper fastening position of each b with numbers as i = 1, 2 and 3) are not displaced in any way. At the bottom, t attached to a beam, which was treated as a rigid body. The tensile force P is t to the system through the beam. Individual bars have different cross-section ar can also be made of materials that feature different properties. Tensile forces σ (MPa) ε Figure 13. The σ(ε) relation obtained during the tension and subsequent compression sample simulation.
As it can be seen in Figure 13, a hysteresis loop graph was obtained. The shape results from the adoption of the symmetry of values in the E p matrix. Precise loop closure at the (0.0) point additionally indicates high accuracy of the algorithm used.

Example E. Tension and Relief of the Three Bar System
The next analyzed issue is related to the tension of the system of three connected bars shown in Figure 14. As it can be seen in Figure 13, a hysteresis loop graph wa from the adoption of the symmetry of values in the E matr (0.0) point additionally indicates high accuracy of the algor

Example E. Tension and Relief of the Three Bar System
The next analyzed issue is related to the tension of the sy shown in Figure 14. The bars have the same length L and the distances bet are indicated as a. It is assumed that the upper fastening p with numbers as i = 1, 2 and 3) are not displaced in any wa attached to a beam, which was treated as a rigid body. The to the system through the beam. Individual bars have differe σ (MPa) ε Figure 14. A three bar system in tension.
The bars have the same length L and the distances between them are constant and are indicated as a. It is assumed that the upper fastening position of each bar (marked with numbers as i = 1, 2 and 3) are not displaced in any way. At the bottom, the bars are attached to a beam, which was treated as a rigid body. The tensile force P is transmitted to the system through the beam. Individual bars have different cross-section areas A i and can also be made of materials that feature different properties. Tensile forces inside the bars are indicated as P i .
The system being solved is statically indeterminable. Three equations have been formulated to solve the problem: two static equations and one displacement-related equation, resulting from the assumption that the beam is a rigid body (no deformations). The sum of forces equation of the analyzed case is as follows: The equation of the sum of force moments related to any point on bar no. 2 results in the following relation: The last relation on the assumed geometry of the system is in the form of the following: where u i indicates the vertical displacement of the no. i bar end, that is elongation of this bar. Introducing the relation describing the classical Hooke's law (P i = E i A i u i /L) into the above equations, a system of three equations with unknown u i displacements is obtained: Assuming that E i = f(σ i , ε i ), and thus E i = f(u i ), the above system is non-linear. As in previous examples, the equation system was solved using the incremental method with a previously described computational algorithm implemented in the proprietary program created in the Scilab environment. The stiffness of the bar material is described with E p i and E n i matrices. For rods No. 1 and 3 the stiffness matrices identical to problem D were used, while in case of rod No. 2 the matrices from problem D were multiplied by 0.75.
A cycle consisting of a load phase until the P force reaches 3.5 MN and a load relief phase were assumed. The following length L = 2.0 m and bar cross-section areas are assumed: A 1 = 0.17 m 2 , A 2 = 0.20 m 2 , A 3 = 0.09 m 2 . It was assumed that the whole system is weightless and that there were no stresses in the bars before the start of applying the load.
Calculation results are shown in Figures 15 and 16. Figure  Graphs in Figure 16 show the relation between σ i stresses [MPa] and ε i strains. It can be seen that the load curves for bars no. 1 no. 3 overlap (partly as more stress is achieved in bar no. 3), while the curve for bar no. 2 is more inclined. This is also the case during the relief stage. It is caused by different stiffness of the materials from which the bars are made.

load.
Calculation results are shown in Figures 15 and 16. Figure 15 shows the relations between the P forces (MN) in specific bars and their elongations u (m). The graphs are color coded: red for bar No. 1, green for bar No. Graphs in Figure 16 show the relation between σ stresses [MPa] and ε strains. It can be seen that the load curves for bars no. 1 no. 3 overlap (partly as more stress is achieved in bar no. 3), while the curve for bar no. 2 is more inclined. This is also the case during the relief stage. It is caused by different stiffness of the materials from which the bars are made. The presented example differs from the previous ones as its solution required solving a system of non-linear equations. Based on the achieved results, it was concluded that:

•
Obtained solutions (in each calculation step) met the equations of the solved equation system; • σ (ε ) and σ (ε ) stress graphs have the same shape as the graph in Figure 7 (stiffness matrices are the same in each case); a similar check was also performed for the σ (ε ) graph, and in this case, the charts were also confirmed to be consistent. Indicated observations prove the correctness of the algorithm and the usefulness of the stiffness matrix for solving non-linear systems that are statically indeterminable.

Example F. Bending of a Reinforced Concrete Beam
The last example is the most complex and, at the same time, features the most reference for practical aspects. It refers to bending of a double-span reinforced concrete beam shown in Figure 17. The beam cross section is rectangular in shape with the following dimensions b × h = 40.0 cm × 60.0 cm. The bottom reinforcement of 3 Ø20 in positive moment zones, and the top reinforcement of 8 Ø20 in negative moment zones were assumed. C30/37 concrete and 34GS reinforcing steel were used, and concrete features were taken from standard [35] Distances between reinforcement centers of gravity and beam edges σ MPa ε Figure 16. σ i (ε i ) graph (i = 1-red; i = 2-green; i = 3-blue).
The presented example differs from the previous ones as its solution required solving a system of non-linear equations. Based on the achieved results, it was concluded that: • Obtained solutions (in each calculation step) met the equations of the solved equation system; • σ 1 (ε 1 ) and σ 3 (ε 3 ) stress graphs have the same shape as the graph in Figure 7 (stiffness matrices are the same in each case); a similar check was also performed for the σ 2 (ε 2 ) graph, and in this case, the charts were also confirmed to be consistent.
Indicated observations prove the correctness of the algorithm and the usefulness of the stiffness matrix for solving non-linear systems that are statically indeterminable.

Example F. Bending of a Reinforced Concrete Beam
The last example is the most complex and, at the same time, features the most reference for practical aspects. It refers to bending of a double-span reinforced concrete beam shown in Figure 17. The beam cross section is rectangular in shape with the following dimensions b × h = 40.0 cm × 60.0 cm. The bottom reinforcement of 3 Ø20 in positive moment zones, and the top reinforcement of 8 Ø20 in negative moment zones were assumed. C30/37 concrete and 34GS reinforcing steel were used, and concrete features were taken from standard [35] Distances between reinforcement centers of gravity and beam edges are 4.0 cm. The stiffness of intermediate elastic support is k = 700 MN/m.

Example F. Bending of a Reinforced Concrete Beam
The last example is the most complex and, at the same time, features the most refer ence for practical aspects. It refers to bending of a double-span reinforced concrete beam shown in Figure 17. The beam cross section is rectangular in shape with the following dimensions b × h = 40.0 cm × 60.0 cm. The bottom reinforcement of 3 Ø20 in positive mo ment zones, and the top reinforcement of 8 Ø20 in negative moment zones were assumed C30/37 concrete and 34GS reinforcing steel were used, and concrete features were taken from standard [35] Distances between reinforcement centers of gravity and beam edge are 4.0cm. The stiffness of intermediate elastic support is k = 700 MN/m. In this case, the stiffness matrices of two considered beam cross-sections (with different reinforcement) were built first. The algorithm for determining the stiffness of reinforced concrete bend cross-sections was used, which was defined in the paper [36]. It uses crosssectional dimensions, surface area and location of main reinforcement, as well as the physical relations for concrete and reinforcing steel as defined in Eurocode 2. Based on achieved results, stiffness matrices were defined. Section stiffness with the reinforcement of 3 Ø20 is described with vectors The E n 2 (14×15) matrix consist entirely of 98,576 values. Vectors X 1 and X 2 are curvature vectors here (κ) expressed in (m −1 ), Y 1 and Y 2 are bending moment vectors (M) expressed in (kNm), and matrices E p 1_ , E n 1 , E p 2_ and E n 2 describe the bending stiffness of sections expressed in (kNm 2 ). The M(κ) relation graphs generated based on the data compiled above are presented in Figures 18 and 19.  The problem was solved with the finite difference method and the incremental algorithm incorporating stiffness matrices that was discussed earlier and used in previous problems. The beam undergone discretization assuming a digitizing range of 0.1 m. One load cycle for the beam and the subsequent relief was assumed. Maximum load value of q = 160 kN/m was chosen so that the plasticization effect of beam sections could be clearly seen during the simulation. Building the load increment vector, the load phase was divided into two ranges. The first range, covering the force range from 0 to 0.95q (i.e., up to the value of 144 kN/m), was divided into 15 subranges with an increment of 9.60 kN/m. The second range, covering the force range from 0.95 q to q (i.e., from 144 to 160 kN/m), was divided into 15 subranges with an increment of 1.067 kN/m. The load phase (from the value q = 160 kN/m to 0) was divided into 15 subranges with an increment of 10.67 kN/m.
The displacement (deflection) graph obtained during specific steps (for successive load values in the load and relief phases) is shown in Figure 20. The x axis describes the coordinates of the beam sections starting from the leftmost support.   The problem was solved with the finite difference method and the incremental algorithm incorporating stiffness matrices that was discussed earlier and used in previous problems. The beam undergone discretization assuming a digitizing range of 0.1 m. One load cycle for the beam and the subsequent relief was assumed. Maximum load value of q = 160 kN/m was chosen so that the plasticization effect of beam sections could be clearly seen during the simulation. Building the load increment vector, the load phase was divided into two ranges. The first range, covering the force range from 0 to 0.95q (i.e., up to the value of 144 kN/m), was divided into 15 subranges with an increment of 9.60 kN/m. The second range, covering the force range from 0.95 q to q (i.e., from 144 to 160 kN/m), was divided into 15 subranges with an increment of 1.067 kN/m. The load phase (from the value q = 160 kN/m to 0) was divided into 15 subranges with an increment of 10.67 kN/m.
The displacement (deflection) graph obtained during specific steps (for successive load values in the load and relief phases) is shown in Figure 20. The x axis describes the coordinates of the beam sections starting from the leftmost support. The problem was solved with the finite difference method and the incremental algorithm incorporating stiffness matrices that was discussed earlier and used in previous problems. The beam undergone discretization assuming a digitizing range of 0.1 m. One load cycle for the beam and the subsequent relief was assumed. Maximum load value of q = 160 kN/m was chosen so that the plasticization effect of beam sections could be clearly seen during the simulation. Building the load increment vector, the load phase was divided into two ranges. The first range, covering the force range from 0 to 0.95q (i.e., up to the value of 144 kN/m), was divided into 15 subranges with an increment of 9.60 kN/m. The second range, covering the force range from 0.95 q to q (i.e., from 144 to 160 kN/m), was divided into 15 subranges with an increment of 1.067 kN/m. The load phase (from the value q = 160 kN/m to 0) was divided into 15 subranges with an increment of 10.67 kN/m.
The displacement (deflection) graph obtained during specific steps (for successive load values in the load and relief phases) is shown in Figure 20. The x axis describes the coordinates of the beam sections starting from the leftmost support.    The other graphs (Figures 21-25) show only the calculation results corresponding to maximum load (q = 160 kN/m)-solid lines, and the total beam relief-dashed lines.   The other graphs (Figures 21-25) show only the calculation results corresponding to maximum load (q = 160 kN/m)-solid lines, and the total beam relief-dashed lines.     From the achieved results, it can be observed that the beam did not return to its original shape after load relief. The system is statically indeterminable and therefore non-zero internal forces in the beam also remained.
In   From the achieved results, it can be observed that the beam did not return to its original shape after load relief. The system is statically indeterminable and therefore non-zero internal forces in the beam also remained.
In   From the achieved results, it can be observed that the beam did not return to its original shape after load relief. The system is statically indeterminable and therefore non-zero internal forces in the beam also remained.
In From the achieved results, it can be observed that the beam did not return to its original shape after load relief. The system is statically indeterminable and therefore non-zero internal forces in the beam also remained.
In the presented example, there were different values of bending moments and curvatures in individual cross-sections. The most significant plastic rotation of cross sections occurred at beam start and end (at terminal supports). The fact that plastic joints were formed only in the most stressed sections is typical for bar structures. The plasticization of the outermost sections is indicated by large value jumps at these locations seen in Figure 23. In Figure 22, the presence of plastic cross sections in these sections is indicated less clearly, but it is due to the limitations of the finite difference method used. Figure 26 shows the M(κ) relation graph in these cross-sections overlapped with the graph from Figure 19.
Materials 2021, 14, 5837 21 of 2 occurred at beam start and end (at terminal supports). The fact that plastic joints were formed only in the most stressed sections is typical for bar structures. The plasticization of the outermost sections is indicated by large value jumps at these locations seen in Fig  ure 23. In Figure 22, the presence of plastic cross sections in these sections is indicated les clearly, but it is due to the limitations of the finite difference method used. Figure 26 show the M(κ) relation graph in these cross-sections overlapped with the graph from Figure  19. High compatibility of both graphs can be seen in Figure 26. Some differences oc curred at the point of severe changes of graph inclination angle. These differences could be reduced, for example, as a result of reduced load increments.
The analyzed beam was previously solved in study [37]. Calculations performed there took into account stiffness functions determined in accordance with the same algo rithm (study [36]), but without using the stiffness matrix description. One load cycle up to q = 50 kN/m was taken into account. At such load, the maximum deflection value o 0.0031073 m was obtained in [37]. Having solved the same problem using the algorithm discussed in this paper, a maximum deflection value of 0.0031061 m was obtained. Thi means that the relative error of the deflection value is −0.04%. This confirms that the algo rithm works correctly as planned.

Discussion
The purpose of this paper was to present the developed matrix constitutive model o structural materials. The main assumptions of this method and the numerical algorithm to allow for its efficient use were presented in Section 2.
In Section 3, the paper presented calculation examples that confirm the effectivenes of the developed method and show some fields of application. Example B presented two constructed stiffness matrices for the assumed baseline data, with the second matrix (Case 2) referring to the data extracted from the experiment. Despite manual calibration of the stiffness matrix, a good compliance between the tension test simulation results and the baseline data was obtained. The versatility of the method is also apparent here. Stiffnes matrices can be constructed for materials with distinct yield strength, for hyper-elastic materials, and material types not considered in this paper.
Examples C and D show simulations of tension, relief, and compression of a sample in axial compressive stress. Obtained simulation results confirm the high versatility of the method as well as the ease of adaptations. In example E, a degradation model was buil by introducing a coefficient that depends on strain history, but a similar effect could be High compatibility of both graphs can be seen in Figure 26. Some differences occurred at the point of severe changes of graph inclination angle. These differences could be reduced, for example, as a result of reduced load increments.
The analyzed beam was previously solved in study [37]. Calculations performed there took into account stiffness functions determined in accordance with the same algorithm (study [36]), but without using the stiffness matrix description. One load cycle up to q = 50 kN/m was taken into account. At such load, the maximum deflection value of 0.0031073 m was obtained in [37]. Having solved the same problem using the algorithm discussed in this paper, a maximum deflection value of 0.0031061 m was obtained. This means that the relative error of the deflection value is −0.04%. This confirms that the algorithm works correctly as planned.

Discussion
The purpose of this paper was to present the developed matrix constitutive model of structural materials. The main assumptions of this method and the numerical algorithm to allow for its efficient use were presented in Section 2.
In Section 3, the paper presented calculation examples that confirm the effectiveness of the developed method and show some fields of application. Example B presented two constructed stiffness matrices for the assumed baseline data, with the second matrix (Case 2) referring to the data extracted from the experiment. Despite manual calibration of the stiffness matrix, a good compliance between the tension test simulation results and the baseline data was obtained. The versatility of the method is also apparent here. Stiffness matrices can be constructed for materials with distinct yield strength, for hyper-elastic materials, and material types not considered in this paper.
Examples C and D show simulations of tension, relief, and compression of a sample in axial compressive stress. Obtained simulation results confirm the high versatility of the method as well as the ease of adaptations. In example E, a degradation model was built by introducing a coefficient that depends on strain history, but a similar effect could be obtained by modifying the stiffness matrix accordingly. In Example D, the stiffness matrix was modified to easily account for compressive stresses. The simulation of the tensile processes of the sample resulted in a graph in the form of an exactly closed hysteresis loop, which confirms that the method is correct and accurate.
Examples E and F deal with the analysis of statically indeterminate systems. In these cases, the calculation algorithms shown in Figures 5 and 6 were fully utilized. In Example E, a system of three parallel connected bars that were tensioned with a force that causes plastic deformation was analyzed. After relief, the system did not return to its original position and non-zero strains and stresses remained in the system. Example F was related to a double-span reinforced concrete beam. The stiffness matrices created in this case referred to the complete reinforced concrete beam cross-sections that take into account the concrete matrix and reinforcing bars (in earlier examples the matrices described material stiffness). This is another confirmation of the versatility and adaptability of the developed method. In this case, the finite difference method was used to create the mathematical model of the task, but stiffness matrices can be used along with other numerical methods such as the finite element method. Simulation of the beam behavior under stress causing high plasticization of outermost sections and then during relief was carried out. After the load relief, non-zero displacements and internal forces remained in the beam. The example in question was previously solved for another article. Result comparison confirmed the correct functioning of the stiffness matrix and calculation algorithm. During the simulation, high compatibility of the obtained M(κ) function with the assumed one was obtained (see Figure 26), which is another confirmation of the proper operation of the algorithm.
An important advantage of the developed matrix description of material properties is the dependence of the stiffness function on two variables, stress σ and strain ε. The difference between this approach and the classical approach is shown in Figure 27. In the classical approach (Figure 27a), the stiffness function depends on one variable, and it can be represented graphically as a line on a plane graph. Method presented in this paper refers to a description of the stiffness, the graphical interpretation of which is the surface (Figure 27b). This way it is much easier to analyze the loading and unloading processes and take the earlier deformation history of the material into account. was modified to easily account for compressive stresses. The simulation of the tensile processes of the sample resulted in a graph in the form of an exactly closed hysteresis loop, which confirms that the method is correct and accurate. Examples E and F deal with the analysis of statically indeterminate systems. In these cases, the calculation algorithms shown in Figures 5 and 6 were fully utilized. In Example E, a system of three parallel connected bars that were tensioned with a force that causes plastic deformation was analyzed. After relief, the system did not return to its original position and non-zero strains and stresses remained in the system. Example F was related to a double-span reinforced concrete beam. The stiffness matrices created in this case referred to the complete reinforced concrete beam cross-sections that take into account the concrete matrix and reinforcing bars (in earlier examples the matrices described material stiffness). This is another confirmation of the versatility and adaptability of the developed method. In this case, the finite difference method was used to create the mathematical model of the task, but stiffness matrices can be used along with other numerical methods such as the finite element method. Simulation of the beam behavior under stress causing high plasticization of outermost sections and then during relief was carried out. After the load relief, non-zero displacements and internal forces remained in the beam. The example in question was previously solved for another article. Result comparison confirmed the correct functioning of the stiffness matrix and calculation algorithm. During the simulation, high compatibility of the obtained M(κ) function with the assumed one was obtained (see Figure 26), which is another confirmation of the proper operation of the algorithm.
An important advantage of the developed matrix description of material properties is the dependence of the stiffness function on two variables, stress σ and strain ε. The difference between this approach and the classical approach is shown in Figure 27. In the classical approach (Figure 27a), the stiffness function depends on one variable, and it can be represented graphically as a line on a plane graph. Method presented in this paper refers to a description of the stiffness, the graphical interpretation of which is the surface (Figure 27b). This way it is much easier to analyze the loading and unloading processes and take the earlier deformation history of the material into account. The matrix description of materials ensures many more possibilities than those presented in this paper. In particular, it is possible to introduce more variables on which the described parameters will depend. This is explained in Figure 28. The matrix description of materials ensures many more possibilities than those presented in this paper. In particular, it is possible to introduce more variables on which the described parameters will depend. This is explained in Figure 28.  Figure 28a shows a situation where a material parameter depends on two variables. This corresponds to stiffness matrices E(σ, ε) included in this paper. Figure 28b shows the case of a parameter that depends on three variables. This third variable may be, for example, the distortion rate ε , which will enable time-dependent analysis. Figure 28c develops the concept with adding another variable. For practical reasons, the cases shown in Figure 28b,c can be prepared by using matrix structures with more than two dimensions. The use of such multidimensional matrices would require appropriate modification of the calculation algorithm discussed in this paper.
The developed method of describing material properties is alternative to classical constitutive models of constitutive materials. It has both advantages (described above) and limitations related to classical models. These limitations include: • Time-consuming calibration of the stiffness matrix that requires supervision by experienced and dedicated researchers, • Application limitation to nonlinear issues (using this method in linear problems is possible but inefficient), • The need to adapt the incremental method parameters of given problem solution to the distribution of characteristic points in vectors X and Y. Areas of further planned research and activities to improve the method include, among others, the algorithm development for automated calibration of the stiffness matrix. This will enable a quick creation of accurate stiffness matrices based on basic data, e.g., laboratory test results.
Potential areas of application of the method include the analysis of structures made of hyper-elastic materials and structures connected with the ground. The ground behaves in a non-linear and complex way, which is difficult to describe using mathematical expressions, so writing down its properties in the form of a stiffness matrix should give good results. Due to the possibility of direct referencing the stiffness matrix with the results of material laboratory tests, it is possible to apply the method when analyzing structures made of new materials or for which no verified analytical constitutive models exist.

Conclusions
The concept of matrix description of material properties presented in this paper is an original approach to the creation of nonlinear constitutive models. The distinguishing feature of this concept is the exclusive use of numerical description (in the form of one or more matrices) and taking into account the dependence of the stiffness parameter on at least two variables. The developed method offers the accurate description of the material behavior, at the cost of its workload. The computational examples presented in the paper, incorporating the developed numerical algorithm, confirmed the effectiveness and accuracy of the method. It is primarily developed for the description and analysis of nonlinear materials. In particular, it is suitable for the description of materials where effective analytical constitutive models have not been developed. The proposed approach offers great adaptability for describing different types of materials as well as for the analysis of a va-  Figure 28a shows a situation where a material parameter depends on two variables. This corresponds to stiffness matrices E(σ, ε) included in this paper. Figure 28b shows the case of a parameter that depends on three variables. This third variable may be, for example, the distortion rate . ε, which will enable time-dependent analysis. Figure 28c develops the concept with adding another variable. For practical reasons, the cases shown in Figure 28b,c can be prepared by using matrix structures with more than two dimensions. The use of such multidimensional matrices would require appropriate modification of the calculation algorithm discussed in this paper.
The developed method of describing material properties is alternative to classical constitutive models of constitutive materials. It has both advantages (described above) and limitations related to classical models. These limitations include:

•
Time-consuming calibration of the stiffness matrix that requires supervision by experienced and dedicated researchers, • Application limitation to nonlinear issues (using this method in linear problems is possible but inefficient), • The need to adapt the incremental method parameters of given problem solution to the distribution of characteristic points in vectors X and Y.
Areas of further planned research and activities to improve the method include, among others, the algorithm development for automated calibration of the stiffness matrix. This will enable a quick creation of accurate stiffness matrices based on basic data, e.g., laboratory test results.
Potential areas of application of the method include the analysis of structures made of hyper-elastic materials and structures connected with the ground. The ground behaves in a non-linear and complex way, which is difficult to describe using mathematical expressions, so writing down its properties in the form of a stiffness matrix should give good results. Due to the possibility of direct referencing the stiffness matrix with the results of material laboratory tests, it is possible to apply the method when analyzing structures made of new materials or for which no verified analytical constitutive models exist.

Conclusions
The concept of matrix description of material properties presented in this paper is an original approach to the creation of nonlinear constitutive models. The distinguishing feature of this concept is the exclusive use of numerical description (in the form of one or more matrices) and taking into account the dependence of the stiffness parameter on at least two variables. The developed method offers the accurate description of the material behavior, at the cost of its workload. The computational examples presented in the paper, incorporating the developed numerical algorithm, confirmed the effectiveness and accuracy of the method. It is primarily developed for the description and analysis of nonlinear materials. In particular, it is suitable for the description of materials where effective analytical constitutive models have not been developed. The proposed approach offers great adaptability for describing different types of materials as well as for the analysis of a variety of structural analysis issues. Stiffness matrices can be used in conjunction with numerical structural analysis methods such as the finite difference method or the finite element method.
Funding: This research received no external funding.
Institutional Review Board Statement: Not applicable.