JP6048358B2 - Analysis device - Google Patents

Description

  The present invention relates to an analysis apparatus, and more particularly to an apparatus for analyzing viscoelastic material characteristics by a finite element method.

  A viscoelastic material constitutive law in which an elastic element and a viscoelastic element are arranged in parallel is used as a model for analyzing the characteristics of a viscoelastic material such as a rubber material. For example, Patent Document 1 and Non-Patent Document 1 show a method of analyzing the correlation between stress and strain taking into account rubber elasticity and viscoelasticity by setting the relaxation time indicating the damping characteristic of the viscoelastic element as a constant. ing.

JP 2010-256293 A

J.C.Simo and T.J.R Hughes, `` Computational Inelasticity '' Interdisciplinary Applied Mathematics volume 7, Springer Verlag (1998)

  In general, in a harmonic vibration test for a rubber material, the correlation between stress and strain changes depending on the magnitude of the amplitude to be excited. Therefore, it is desirable to use an analytical model that can reproduce the amplitude dependence.

  The present invention has been made in view of these circumstances, and an object thereof is to provide a viscoelastic material analysis technique with improved reproducibility with respect to amplitude dependency.

  In order to solve the above problems, an analysis apparatus according to an aspect of the present invention analyzes stress-strain characteristics of a viscoelastic material based on a nonlinear viscoelastic material constitutive law in which an elastic element and a viscoelastic element are arranged in parallel. A first calculation unit configured to calculate a displacement amount of a node by setting a condition in a viscoelastic material model divided into a finite number of elements having a node as a boundary, and a strain at the node using the displacement amount Using the second calculation unit for calculating the velocity, the third calculation unit for calculating the value proportional to the power of the strain rate as the base, as the relaxation time of the viscoelastic element, and the relaxation time calculated from the strain rate And a fourth calculator for calculating the stress at the nodes.

  According to this aspect, the analysis apparatus can analyze the stress-strain characteristic of the viscoelastic material using the relaxation time expressed as a function of the strain rate. Thereby, the reproducibility of the stress-strain characteristic according to the amplitude amount vibrated by the viscoelastic material can be enhanced.

  ADVANTAGE OF THE INVENTION According to this invention, the analysis technique which improved the reproducibility with respect to an amplitude dependence can be provided by using the value of the power which makes a strain rate the base as relaxation time of a viscoelastic element.

It is a figure which shows typically the viscoelastic material constitutive law which concerns on embodiment of this invention. It is a figure which shows the test conditions in a harmonic vibration test. It is a graph which shows the calculation method of the number of powers typically. It is a figure which shows typically the stress relaxation curve in a stress relaxation test. It is a figure which shows the stress relaxation curve at the time of changing input strain. It is a figure which shows the static characteristic curve of a distortion and stress. It is a figure which shows the relationship curve and static characteristic curve of the distortion | strain and stress in a fixed strain rate test. It is a block diagram which shows the function structure of the analyzer which concerns on embodiment. It is a graph which shows an example of the stress-strain curve obtained by an analyzer. It is a graph which shows an example of the stress relaxation curve obtained by an analyzer. It is a graph which shows an example of the frequency characteristic obtained by an analysis device. It is a flowchart which shows operation | movement of an analyzer.

  Embodiments of the present invention will be described below with reference to the drawings.

  FIG. 1 is a diagram schematically showing a viscoelastic material constitutive law according to an embodiment of the present invention. In the viscoelastic material model shown in the figure, an elastic element having an elastic modulus G0, an elastic element having an elastic modulus Gi (i = 1 to N, N is a natural number) and a viscous element having a viscosity coefficient ηi are connected in series. A viscoelastic element is connected in parallel. This model combining an elastic element and a viscoelastic element is used as a model for representing the dynamic characteristics of rubber parts such as tires and rubber bushes.

  Here, when the rigidity ratio in the elastic element and the viscoelastic element connected in parallel is γi (i = 0 to N) and the relaxation time in the viscoelastic element is τi = ηi / Gi, at the time t in this viscoelastic material model The stress S is expressed by the following formulas (1) and (2).

  In the formula (1), S indicates the second Piola-Kirchhoff stress, and S ○ having a superscript sign indicates that the stress is only the elastic component from which the viscous force component is removed. J represents the volume change rate of the viscoelastic material. The volume change rate J is expressed by J = det [F] using a determinant (det) of the deformation gradient tensor F indicating the relationship between the linear transformation of the position before and after deformation at a certain material point. The operator DEV is represented by the following equation (3) using the right Cauchy-Green tensor C = FT · F. In Expression (3), [•] represents a variable to be operated by the operator DEV.

  Qi indicates the viscous force in each viscoelastic element, and Qi in the equation (1) is expressed by the development equation shown in the equation (2). In equation (2),

Represents a deviation component of strain potential energy in a superelastic body. Also,

Is a modified right Cauchy-Green tensor from which the volume component is removed, and is represented by equation (4).

  The second Piola-Kirchoff stress S expressed by the above equations (1) and (2) is expressed by the following equation (5) by integrating the convolution integral form using the integration factor exp (t / τi). It can be expressed in integral form.

  In addition, in Formula (5)

Is a volume component of strain potential energy. Further, g (t) is a relaxation function and is represented by the above equation (6).

  Here, the second Piola-Kirchoff stress Sn + 1 at the time tn + 1 can be obtained by transforming the expression (5) expressed as a function of the time t into the following expression (7). The function at the calculation step n is expressed as (•) n, and the function at the calculation step n + 1 is expressed as (•) n + 1.

  In equation (7)

Is an intermediate function obtained as an approximate solution obtained by integrating at a time interval of [tn, tn + 1] using the midpoint theorem, and is expressed by the following equation (8).

  here,

Is defined by the following formulas (9) and (10).

  Equation (9) is Kirchhoff elastic stress.

With the following formula (11)

Can be expressed by the following formula (12).

Kirchoff stress tensor

Can be expressed by the following equation (13) using the second Piola-Kirchoff stress Sn + 1.

Therefore, the Kirchoff stress tensor can be expressed by the following equation (14) from the equations (7) and (13).

Note that the operator dev in the equation (14) is defined by the following equation (15). In formula (15), (•) represents a variable to be operated by the operator DEV.

Further, the relaxation function g ^* in the equation (14) is defined by the following equation (16).

  From the above, as an intermediate function

Is maintained for each calculation step, the Kirchoff stress tensor value can be obtained by numerical analysis. The details of the above-described formula modification are described in Non-Patent Document 1, for example.

  In the present embodiment, in the viscoelastic model shown in FIG. 1, the relaxation time τi representing the damping characteristic of the viscoelastic element is used as a variable depending on the strain rate, thereby improving the reproducibility of the amplitude dependence. The analysis method is shown. On the other hand, in the viscoelastic model (Simo model) shown in Non-Patent Document 1 and the like, various dynamic characteristics in the rubber material can be expressed by calculating the relaxation time τi as a constant. There is a problem that the prediction accuracy of the amplitude dependency is not high. The present inventor has found that it is effective to use the relaxation time τi as a function of the strain rate as a method for improving the accuracy of amplitude dependency. In the present embodiment, an analysis method is proposed in which the reproducibility of the amplitude dependence is improved by using the relaxation time τi that satisfies the relational expression shown in the following expression (17).

  In equation (17),

Is a Green-Language strain tensor with the deviation component removed, and a modified right Cauchy-Green tensor.

Is represented by the following formula (18).

Also,

Is the strain rate and represents the time derivative of the strain tensor. Also,

Represents the magnitude of the strain rate, and is expressed by the following equation (19) when a three-dimensional strain tensor is used.

  The power factor mi and the proportionality constant 1 / Ai are values introduced in the present embodiment, and are derived from the relational expression described later using the test result of the harmonic vibration test.

  Next, a description will be given of a modification of the expression necessary for numerical analysis when introducing the relaxation time τi shown in Expression (17). In order to calculate Kirchhoff stress using the above equation (14), it is necessary to modify the equation including the relaxation time τi. Specifically, the equation (8) may be transformed into the following equation (18) and the equation (9) may be transformed into the following equation (19).

Here, the function of the relaxation time τi is defined by the following equations (20) and (21).

Note that the equation (21) is used as the relaxation term τi in the second term on the right side of the equation (18) and the equation (19) because the approximate solution obtained by integrating at the time interval of [tn, tn + 1] using the midpoint theorem is used. This is because it was used.

  Next, a method for identifying a material constant in the viscoelastic material model shown in this embodiment will be described. First, after identifying the identification method of the power factor mi and the proportionality constant 1 / Ai, the identification method of the stiffness ratio γ0 and the elasticity modulus G0 of the elastic elements arranged in parallel is shown, and finally the stiffness ratio γi of the viscoelastic element is identified. The method is shown.

  FIG. 2 is a diagram showing test conditions in the harmonic vibration test. This figure shows a graph of strain E = Epre + εsin (ωt) applied to a viscoelastic material specimen when the horizontal axis is time t and the prestrain Epre, the amplitude ε, and the frequency ω. By measuring the dynamic elastic modulus G of the specimen when such strain E is vibrated, the power mi can be derived using the relationship of the following formula (22).

  Equation (22) is obtained by solving the following equation (23), which is an evolution equation obtained by substituting equation (17) into equation (2) for the input of the harmonic vibration test shown in FIG.

  FIG. 3 is a graph schematically showing a method for calculating the power mi. This figure shows the logarithmic value of the dynamic elastic modulus G, which is the test result value when the prestrain Epre and the frequency ω are changed with the logarithm log (ε) of the amplitude ε to be excited in the harmonic vibration test as the horizontal axis. log (G) is plotted. The approximate line L1 on the graph is an approximate line having a slope l1 corresponding to the plot value when the amplitude values ε1 to ε3 are changed with respect to the pre-strain Epre1 and the frequency ω1. Similarly, the approximate lines L2 and L3 are approximate lines having inclinations l2 and l3 corresponding to plot values when the amplitude values ε1 to ε3 are changed with respect to the pre-strains Epre2 and Epre3 and the frequencies ω2 and ω3. . Since the following equation (24) can be obtained by taking logarithms on both sides of the equation (22), the power mi can be obtained from the slope li.

  From equation (24), the power can be obtained from the relational expression mi = −1−li. In the equation (24), α represents an intercept of the approximate line shown in FIG. Thereby, the power mi corresponding to the excitation frequency ωi in the harmonic vibration test can be obtained. Note that the value of the amplitude ε shown in FIG. 3 is an example, and a value of the amplitude ε of 4 or more may be used to obtain the power mi, or the value of the amplitude ε may be 2.

  When the value of the power mi does not change so much with the excitation frequency ωi, instead of obtaining the mi corresponding to each frequency, the common power m obtained as the average value of each mi is used. Also good. By using the common power m, the calculation load can be reduced as compared with the case where a different power mi is used for each frequency.

  Further, the proportionality constant 1 / Ai can be obtained from the relationship of the following equation (25) using the power mi obtained by the above method. Similar to the derivation method for the power mi, the frequency ωi in the equation (25) is the excitation frequency in the harmonic vibration test, and by changing this frequency, the frequency component ωi corresponds to each different viscoelastic element. The proportionality constant 1 / Ai can be obtained.

  Next, a method for calculating the stiffness ratio of the elastic element γ0 will be described. FIG. 4 is a diagram schematically showing a stress relaxation curve in the stress relaxation test. This figure has shown the test result which measured the time change of the stress Q when input strain Epre is made constant. From the instantaneous stress Q0 at time t = 0 shown in this test result and the relaxation stress Q∞ at time t = ∞, the stiffness ratio of the elastic element can be obtained by γ0 = Q∞ / Q0.

  Next, a method for calculating the elastic modulus G0 of the elastic element will be described. In the present embodiment, the elastic modulus of the elastic element is determined by determining the superelastic coefficients C10, C20, and C30 defined in the Yeoh material model, which is a superelastic material model. These superelastic coefficients can be derived by obtaining a relaxation stress Q∞ for a plurality of input strains Epre in the stress relaxation test shown in FIG.

  FIG. 5 is a diagram showing a stress relaxation curve when the input strain Ek is changed. FIG. 6 is a diagram showing a static characteristic curve QR of strain E and stress Q. As shown in FIG. 5, a stress relaxation test is performed on a plurality of load strains E1 to EN having different values to obtain relaxation stresses QR (E1) to QR (EN) corresponding to each load strain Ek. By graphing this relational expression as a strain-stress curve, a static characteristic curve QR shown in FIG. 6 is obtained. The elastic modulus of the elastic element can be determined by determining the superelastic coefficients C10, C20, and C30 that can approximate the static characteristic curve using a known technique.

  Finally, a method for calculating the rigidity ratio γi of the viscoelastic element will be described. FIG. 7 is a diagram showing a relation curve Qi and a static characteristic curve QR between the strain E and the stress Q in the constant strain rate test. The static characteristic curve QR is the same as the static characteristic curve shown in FIG. Relationship curves Q1 to QN show the test results of strain E and stress Q corresponding to the case where the strain rate Vi is changed to V1 to VN in the constant strain rate test. Here, the i-th relationship curve Qi corresponds to the dynamic characteristic of the i-th viscoelastic element, and corresponds to the viscoelastic element having a proportional constant 1 / Ai calculated from the excitation frequency ωi. In order to obtain such a correspondence relationship, a value satisfying the relationship of strain rate Vi = ωiε may be determined using the excitation frequency ωi and the amplitude ε in the harmonic vibration test.

  At this time, the rigidity ratio γi of each viscoelastic element can be obtained from the area Zi surrounded by the relationship curve Qi between the strain E and the stress Q shown in FIG. If the area surrounded by the static characteristic curve QR is Z0 and the increase in area due to the increase in strain rate from i-1 to i-th is Zi, the rigidity ratio is identified by γi = γ0 × Zi / Z0. .

  By the above method, the material constants mi, Ai, γ0, G0, γi in the viscoelastic material model of the present embodiment can be identified. That is, in the present embodiment, the material constant can be uniquely determined from the result of the material test.

  Other material constant identification methods include an identification method based on optimization calculation, but the optimization method does not converge and the solution cannot be obtained. The obtained optimum solution is a local solution. For this reason, there are problems such as a decrease in reproduction accuracy. On the other hand, in this embodiment, since optimization calculation is not required, the problem that a solution cannot be obtained can be solved. Moreover, since the material constant is identified from the result of the material test, the reproducibility of the analysis can be improved.

  Next, a finite element analysis method using the evolution equation shown in the above equation (23) will be described. FIG. 8 is a block diagram illustrating a functional configuration of the analysis apparatus 100 according to the embodiment. The analysis device 100 includes a construction unit 10, a calculation unit 20, a storage unit 30, and an output unit 40.

  Each block shown in the block diagram of the present specification can be realized in terms of hardware by an element such as a CPU of a computer or a mechanical device, and in terms of software, it can be realized by a computer program or the like. The functional block realized by those cooperation is drawn. Therefore, those skilled in the art will understand that these functional blocks can be realized in various forms by a combination of hardware and software.

  The construction unit 10 constructs a viscoelastic material model to be calculated by the calculation unit 20 by replacing the viscoelastic material to be analyzed with a finite number of small elements. Each element constituting the viscoelastic material model is defined such that numerical analysis is possible. Specifically, for each element, a node coordinate value, an element shape, a material characteristic, and the like in the coordinate system are defined. When a two-dimensional model is targeted, a three-node element having a triangular shape or a four-node element having a quadrangular shape may be used as each element. When a three-dimensional model is targeted, a 4-node element having a tetrahedral shape, a 6-node element having a hexahedral shape, or the like may be used. The construction unit 10 can be realized, for example, by using general-purpose software based on a known technique called a preprocessor.

  The storage unit 30 stores input values necessary for processing by the calculation unit 20 and stores output values at each calculation step calculated by the calculation unit 20. The storage unit 30 stores at least the power factor mi, the proportionality constant 1 / Ai, the stiffness ratios γ0 and γi, and the superelastic coefficients C10, C20, and C30 of the elastic elements as input values. These material constants are values set by the above-described material constant identification method according to the characteristics of the viscoelastic material to be analyzed. For example, a plurality of types of material constants are stored in the storage unit 30 in order to enable analysis of a plurality of types of viscoelastic materials, and corresponding material constants are determined according to the type of viscoelastic material to be analyzed. Is read by the calculation unit 20.

  The calculation unit 20 uses the viscoelastic material model constructed by the construction unit 10 to calculate values such as displacement, strain, and stress at the nodes of each element for each calculation step. The calculation unit 20 repeatedly executes the calculation process until a predetermined calculation end condition is satisfied, and causes the storage unit 30 to store the obtained value. The calculation unit 20 includes a first calculation unit 21, a second calculation unit 22, a third calculation unit 23, and a fourth calculation unit 24.

  The first calculation unit 21 sets boundary conditions for the viscoelastic material model divided into a finite number of elements, and calculates the displacement amount U of the node in each element from the input conditions for each calculation step. The 1st calculation part 21 produces the whole structure matrix showing the whole structure of a viscoelastic material model, after producing the rigidity matrix for solving the rigidity equation in each element. A displacement amount U of an unknown node is calculated by performing an analysis process by introducing a displacement amount and a node force of a known node, which are input conditions, to the entire stiffness matrix. In addition, the 1st calculation part 21 is realizable by using the general purpose software based on the well-known technique called a solver, for example.

  The second calculator 22 calculates the strain rate at the node of each element using the displacement amount U obtained by the first calculator 21 as an input value. In addition, the 2nd calculation part 22 calculates the strain rate in step n + 1 which is the next calculation step by using the distortion amount obtained in calculation step n.

  The second calculation unit 22 calculates the total displacement amount at the calculation step n + 1 from the total displacement amount ψn of the node at the calculation step n and the displacement amount U obtained by the first calculation unit 21 by ψn + 1 = ψn + U. Using this total displacement amount ψn + 1, the deformation gradient tensor Fn + 1, volume change rate (Jacobiane) Jn + 1, right Cauchy-Green tensor Cn + 1, and deviation component in the calculation step n + 1 are represented by the relationship of the following equations (26)-(31). Modified deformation gradient tensor removed

And modified right Cauchy-Green tensor

Calculate

Note that D in Equation (27) is a differential operator for obtaining the deformation gradient tensor F.

  From the corrected right Cauchy-Green tensor obtained by equation (31) and the strain amount at calculation step n, the strain amount and strain rate at calculation step n + 1 can be obtained by the following equations (32) and (33). Note that Δtn in Equation (18) is a time step discretized corresponding to each calculation step.

  The third calculator 23 calculates the relaxation time τi necessary for the stress calculation from the strain rate obtained by the second calculator 22. The third calculating unit 23 reads the power mi and the proportionality constant 1 / Ai stored in the storage unit 30, and calculates the relaxation time τi using the above equations (20) and (21).

  The fourth calculator 24 obtains Kirchhoff stress at the node using the relaxation time τi calculated from the strain rate. First, the fourth calculation unit 24 obtains Kirchoff elastic stress in the calculation step n + 1 using the equation (11). Next, an intermediate variable is obtained using equations (12) and (18). Finally, Kirchhoff stress is obtained using equations (14) and (19). The fourth calculation unit 24 stores the intermediate variables and stress values obtained in the respective calculation steps in the storage unit 30.

  The output unit 40 outputs the calculation result obtained by the calculation unit 20. The output unit 40 outputs numerical data for obtaining a stress-strain curve or a stress relaxation curve in the viscoelastic material model as a calculation result.

  FIG. 9 is a graph showing an example of a stress-strain curve obtained by the analysis apparatus 100, and corresponds to the test result of the constant strain rate test. In addition, in this figure, the "test value" obtained by the corresponding constant strain rate test is shown with the "calculated value" which is an output result from the output part 40. FIG. Calculation values C1 to C4 are calculation results corresponding to different strain rates, and test values T1 to T4 are test results when a constant strain rate test is performed at the same strain rate as the calculation conditions. By using the analysis apparatus of the present embodiment, as shown in the figure, it is possible to obtain a calculation result that accurately reproduces the result of the constant strain rate test.

  FIG. 10 is a graph illustrating an example of a stress relaxation curve obtained by the analysis apparatus 100. Also in this figure, the “test value” obtained in the stress relaxation test is indicated by a broken line together with the solid line indicating the “calculated value” that is the output result from the output unit 40. Each stress relaxation curve shows a calculated value and a test value when different instantaneous stresses are applied. By using the analysis apparatus of the present embodiment, as shown in the figure, it is possible to obtain a calculation result that accurately reproduces the result of the stress relaxation test.

  FIG. 11 is a graph illustrating an example of frequency characteristics obtained by the analysis apparatus 100. This figure shows the result of obtaining the value of the absolute elastic modulus when the excitation frequency ω and the amplitude ε are changed in the harmonic vibration test. This figure shows the calculation results when the magnitude ε in the harmonic vibration test is 0.2% with respect to the pre-strain and 2.0%. As shown in this figure, it was possible to represent a change in absolute elastic modulus due to a change in excitation frequency and a decrease in absolute elastic modulus accompanying an increase in amplitude. Generally, in viscoelastic materials such as rubber materials, it is known that the absolute elastic modulus decreases as the amplitude increases. Therefore, according to the analysis apparatus 100 in the present embodiment, the amplitude dependence in the viscoelastic material is reduced. Can be reproduced.

The operation of the analysis apparatus 100 configured as above will be described.
FIG. 12 is a flowchart showing the operation of the monitoring device 2. First, a viscoelastic material model to be analyzed is constructed (S10), and material constants corresponding to the viscoelastic material model to be analyzed are set (S12). Next, the nodal displacement of each element constituting the viscoelastic material model is calculated for each calculation step based on the input conditions (S14). The strain rate is calculated from the obtained set displacement amount (S16), the relaxation time is calculated from the strain rate (S18), and the stress is calculated using the relaxation time obtained from the strain rate (S20). If the calculation end condition is not satisfied, the steps S14 to S20 are repeated to calculate the stress for each calculation step (N in S22). If the end condition is satisfied (Y in S22), the obtained calculation result is output (S24).

  In the above, this invention was demonstrated based on embodiment. It will be understood by those skilled in the art that the present invention is not limited to the above-described embodiment, and various design changes are possible, various modifications are possible, and such modifications are within the scope of the present invention. By the way.

  In the above-described embodiment, the method for identifying the proportionality constant 1 / Ai using the relational expression (25) from the result of the harmonic vibration test has been described. In the modification, the proportionality constant 1 / Ai may be identified from the result of the constant strain rate test. Hereinafter, a method for identifying the proportionality constant 1 / Ai in the modification will be described.

  The following formula (34) is obtained by modifying the left side of the above formula (25).

  Here, since ωiε of the denominator is the maximum strain rate in the harmonic vibration test, it can be associated with the strain rate Vnom in the constant strain rate test. Moreover, since the amplitude ε of the molecule is the maximum strain amount in the harmonic vibration test, it can be associated with the measured strain E * in the constant strain rate test. Therefore, using the strain rate Vnom in the constant strain rate test and the measured strain E *, the equation (34) can be rewritten as the following equation (35).

  Therefore, in the modification, the proportionality constant 1 / Ai can be identified from the test result in the constant strain rate test using the equation (35).

  DESCRIPTION OF SYMBOLS 10 ... Construction part, 20 ... Calculation part, 21 ... 1st calculation part, 22 ... 2nd calculation part, 23 ... 3rd calculation part, 24 ... 4th calculation part, 30 ... Memory | storage part, 100 ... Analysis apparatus.

Claims (4)

An apparatus for analyzing stress-strain characteristics of a viscoelastic material based on a nonlinear viscoelastic material constitutive law in which an elastic element and a viscoelastic element are arranged in parallel,
A first calculation unit configured to calculate a displacement amount of the node by setting a condition in a viscoelastic material model divided into a finite number of elements having a node as a boundary;
A second calculator for calculating a strain rate at the node using the displacement;
A third calculation unit that calculates a value proportional to a power value based on the strain rate as a relaxation time of the viscoelastic element;
A fourth calculation unit for calculating stress at the node using a relaxation time calculated from the strain rate;
An analysis device comprising: As a constant determined according to the viscoelastic material to be analyzed in order to obtain the relaxation time, a power factor mi having the strain rate as a base, and a proportionality constant 1 / Ai multiplied by the power value, Further comprising a storage unit for storing
The analysis apparatus according to claim 1, wherein the third calculation unit reads the power mi and the proportionality constant 1 / Ai stored in the storage unit and calculates the relaxation time. The storage unit
Store the power mi obtained by the following equation (1) with the measurement results of the absolute elastic modulus G and amplitude ε in the harmonic vibration test using the viscoelastic material to be analyzed as input,
3. The proportionality constant 1 / Ai obtained by the following equation (2) is stored with the measurement result of the excitation frequency ω and amplitude ε in the harmonic vibration test and the value of the power mi as input. The analysis device described in 1.

A program for obtaining stress-strain characteristics of a viscoelastic material based on a nonlinear viscoelastic material constitutive law in which an elastic element and a viscoelastic element are arranged in parallel,
On the computer,
A function for calculating a displacement amount of the node by setting a condition in a viscoelastic material model divided into a finite number of elements having a node as a boundary;
A function of calculating a strain rate at the node using the displacement;
A function of calculating a value proportional to a power value with the strain rate as a base as a relaxation time of the viscoelastic element;
A function of calculating stress at the node using a relaxation time calculated from the strain rate;
A program to realize

Images