Damage level calculation in SMC composite structures

(1)

M A S T E R’S T H E S I S

2006:050 CIV

MARTIN FISK

Damage Level Calculation in SMC Composite Structures

MASTER OF SCIENCE PROGRAMME Mechanical Engineering

Luleå University of Technology

Department of Applied Physics and Mechanical Engineering Division of Polymer Engineering

2006:050 CIV • ISSN: 1402 - 1617 • ISRN: LTU - EX - - 06/50 - - SE

(2)

Abstract

Due to design and weight properties the automotive industry increasingly substitutes ordinary sheet metal with polymer products where one of these new materials is sheet molding compound (SMC). Therefore is an interest from the industry to do Finite Element Analysis (FEA) on these products.

The goal was to create a program in Matlab that calculates the damage level in a tailgate during an open-close case. One solution to that problem is to use a linear material model in the FEA program, write a program that transfers the results from linear to nonlinear within a certain margin, calculate the damage level, and then postprocess the results.

(3)

Acknowledgments

This work has been carried out at Volvo Car Corporation, Gothenburg, Sweden and is an industrial project for the degree of M. Sc. at Luleå University of Technology.

I will express my sincere gratitude to my supervisors PhD. Magnus Oldenbo and M. Sc. Rickard Petersson for their invaluable guidance and support.

I also will like to acknowledge the assistance from the staff at the departments of

“Exterior Engineering” and “CAE Durability and Load Data”.

Finally a great thank you to Professor Janis Varna at the department “Polymer Engineering”, Luleå University of Technology for examining this thesis.

Gothenburg, Sweden January 2006

Martin Fisk

(4)

1 Introduction ...1

1.1 Sheet Molding Compound – Manufacturing ...1

1.2 SMC Characteristics ...2

1.3 BMC ...3

1.4 FEA...3

1.5 Technical information...3

2 Theory...4

2.1 Damage Level and Material Data ...4

2.2 Strain Rate ...5

2.3 Neuber Correction ...7

3 Matlab...7

3.1 Program Structure...9

3.2 Punch Files...9

3.3 Sorted Data ...11

3.4 Patran Format...12

3.5 The Matlab Program ...13

3.5.1 Main...13

3.5.2 prepare ...13

3.5.3 Punch functions ...13

3.5.4 read ...14

3.5.5 hub ...14

3.5.6 principal ...14

3.5.7 neuber ...15

3.5.8 damage...15

3.5.9 Material model functions...15

3.5.10 writeToFile ...15

3.5.11 materialProperties ...16

4 Example and Validation ...17

4.1 Geometry ...17

4.2 FEA Methods...17

(5)

4.3 Results of Neuber Correction ...18

4.4 Results of Damage Level Simulation ...18

5 Discussion and conclusions ...21

References ...22

I Appendix I ...i

I.I Obtaining Tensile Stress-Strain data for BMC...i

II Appendix II...iii

II.II Perl script ... iii

II.III Nastran...v

II.IV Matlab ...vi

II.IV.I startup.m or main.m...vi

II.IV.II Prepare.m ...vii

II.IV.III splitPunch.m ...ix

II.IV.IV splitSort.m... xiii

II.IV.V read.m ...xiv

II.IV.VI hub.m ...xvi

II.IV.VII principal.m ...xvi

II.IV.VIII neuber.m ...xvii

II.IV.IX damage.m...xix

II.IV.X materialModelTime.m ...xxi

II.IV.XI topValues.m ...xxi

II.IV.XII materialModel.m...xxii

II.IV.XIII writeToFile ... xxiii

(6)

1

1 Introduction

During the past decades, the use of polymers in the industry has continuously increased. Ten years ago, most cars had doors, fenders, tailgates, etc.,

manufactured using sheet metal, but today, because of weight and design properties there is an increasing trend for polymers to substitute the traditional material. The industry has shown an interest in simulating the properties of polymer products in finite element analysis (FEA) and is therefore of current value.

A polymer material generally applied to tailgates is SMC. Using today’s

technology, a tailgate can reach approximately 20 % lower weight as compared to steel, and with a proper FEA the concept can offer a very competitive solution.

However, now the problem is how to treat the material properties during the FEA.

1.1 Sheet Molding Compound – Manufacturing

Sheet molding compound is manufactured via a compression molding process, where the “raw material” is SMC prepreg. Using compression molding one can attain products with high surface quality as well as high volume production.

SMC prepreg contains fibers, resin, and additives and is therefore a fully formulated system. The fiberglass is normally chopped E glass fiber

reinforcement with a length and diameter of approximately 25 mm and 14 μm respectively. The fibers are distributed into bundles which contain 180-400 fibers per bundle. A normal SMC panel is approximately 2-4 mm thick and in Table 1 the typical composition of SMC can be seen.

Component Weight Percentage Function

Polystyrene 11 – 13 Reactive monomers – cross

linking function Unsaturated polyester 10 – 16 Reactive polymer

Fiberglass ~ 30 Reinforcement

Calcium carbonate ~ 40 Filler in the matrix

Magnesium oxide 0.5 - 0.7 Increases viscosity during the maturation process

Zinc stearate 1.0 - 2.4 Lubricant

Table 1: A typical example of the composition of sheet molding compounds.

(7)

1.2 SMC Characteristics

As previously mentioned, a normal SMC panel is 2-4 mm thick and contains glass fibers with a length of approximately 25 mm. A sound hypothesis is therefore that all fibers are aligned with the plate and an orthotropic material model can be used.

We further assume that the bundles are randomly distributed in the plane, which gives a transversally isotropic material model, and then allow the 1-2 plane to become the plane of the panel. Current literature states that G

13

≈ G

23,

this is the reason why an in-plane isotropic material model can be applied, see Oldenbo et al.

[1].

To obtain a macroscopic understanding of the behavior of the materials, a stress-strain graph is a good place to start, see Figure 1.

Figure 1: A typical stress-strain graph for SMC. The nonlinearity depends on the micro cracks that develop mainly in the transverse direction in relation to its loading direction. The strain rate is 2x10^-4s^-1. Tensile strain data collected from Oldenbo et al. [2] and pressure data can be estimated as being twice as great as tension. (According to other studies.)

In SMC, the dominant deteriorating mechanism in tension is that cracks develop

during strain. Using its rather large reduction as in Young's modulus, see Figure 1,

we can conclude that cracks mainly develop transversely to its loading direction,

but not parallel to it. Analysis carried out in direction 2 confirms this. Considering

this, damaged composites will be anisotropic after loading, even though they were

isotropic in their original state.

(8)

3

1.3 BMC

Another common composite used to produce structures is bulk molding compounds (BMC). BMC has a fiber length of approximately 5 mm, and therefore responds earlier to crack intensification than SMC, i.e., less macro ductility. The component weight percent composition is similar to SMC.

On the other hand, the shorter fibers result in higher surface quality and styling opportunities. Components such as, headlight houses, spoilers and other extra high form or surface quality products are therefore manufactured using this composite.

No stress-strain curves could be found for BMC whereas specimens where token from a tailgate. Seven specimens formed a stress strain curve for BMC material and the result can bee seen in Appendix I.

The theory is developed for SMC composite and the only changes to be done for BMC is the material property, i.e. Young's modulus, threshold etc.

1.4 FEA

To capture transient waves that arrive in a structure during an open-close case, a dynamic FE analysis is required. If now treating the material as non – linear, every time step must pass some kind of non-linear calculation. With this situation, a realistic answer can bee in practice if magnitude of time or a fast computer can bee used. Perhaps both of it.

One solution of that problem is to use a linear material model, write a program that transfer the result from linear to nonlinear within a certain margin, calculate the damage level, and then postprocessing the results. In this project this is done by Matlab.

Another major benefit is that unnecessary data can be filtered. That gives elements with strain exceeding the damage initiation threshold the opportunity to be

highlighted.

1.5 Technical information

The following programs will be mentioned through the report.

Program Version Program Version

Matlab 7.0.0.19901 (R14) Abaqus viewer 6.3 – 1

Animator 3 - V0.6.4a1 testXpert 6.01

Nastran 2004

Ansa 11.3.6

Abaqus

(9)

2 Theory

2.1 Damage Level and Material Data

In its original state, SMC can be treated as an in-plane isotropic material, but after damage has occurred it becomes in-plane anisotropic. With two symmetry planes (orthotropic) and practicing plane stress, the constitutive relationship for the anisotropic material will be:

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢

⎣

⎡

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢

⎣

⎡

−

− +

−

− − −

−

=

⎥⎥

⎥

⎦

⎤

⎢⎢

⎢

⎣

⎡

12 2 1

2 1 2

2 2

1

2 1 2

1

0 12

2 1

) 1 )(

1 (

) 1 ( 0 2

0

) 0 1 (

1 )

1 )(

1 (

) 0 1 )(

1 ( ) 1 (

1

σ σ σ ν ν

ν

γ ε ε

D D D

D D

D D D

E

(1)

where D

1

and D

2

is the damage level in each direction respectively.

If loading is uniaxial and Hooke's law is applied, (1) becomes

2 1 0 1 1 1 1 2 1 0

1 ~ ~ (1 )

) 1

( D E E E D

E − =

σ

=

ε

⇒ = −

ε (2)

(in the 1 direction) and

0 1 1

~

1 E

D = − E (3)

Using the strain energy expression for the damaged composite in one direction (direction 1) and derivate it with respect to D

1,

the thermodynamic force (sometimes denoted "damaged energy release rate"), Y

i

, will be found [1].

3 1 2 1

1

E ( 1 D )

Y = σ −

;

Y₂ =0

(4)

where

)) ( ( ) ˆ (

1

τ

f Y

D_i = _i

; i = 1, 2 (5)

and

) 0

| ) ˆ ( max(

)

(t D t

D_i = _i

τ

≤

τ

≤

; i = 1, 2 (6)

Oldenbo et al. [2] gives that

0. D

) Y - a(Y ) (Y f D

2

0 1 1

1 1

=

= Y

0

≤ 0.5 (7)

(10)

5

With knowledge of that, use equation (4)

) ) 1

( ( ₃ ₀

1 0

2 1

1 Y

D a E

D −

= σ−

(8)

solve for σ

1

0 2 1 0

1 3

1) ( )

1

( E

a a Y D

D ₊ ₌

σ

−

(9)

with equation (1) and only σ

1

≠ 0 (9) becomes

2 1 0 1 0

1 1

1 ⎟⎟⎠ =

ε

⎜⎜ ⎞

⎝

⎛

− +

aE D

a Y

D

(10)

using Matlab and equation (3) thus Figure 2 is created. Here a = 0.205 MPa

^-1

and Y

0

= 0.0654 MPa, see Oldenbo et al. [2]. In addition, a trend line was estimated from the curve, see equation (11).

052 . 1 1656 . 0 1714 . 0

~

2 0

+

−

= ε ε

E

(11)

In Figure 1 one can see that the material behaves in a linear pattern under a strain of 0.25 % tension and 0.5 % compression. One conclusion arising from this is that non-important damage has occurred before the strain reaches 0.25 % tension and 0.5 % compression. A threshold for the values mentioned above is therefore a sound approximation and are further adopted.

No values of a and Y

0

were found for BMC why a damage curve was estimated from the stress-strain curve. See Appendix I for further information.

2.2 Strain Rate

Compared to quasi-static loading the stress-strain curve differs when the strain rate increases. Jendli et al. [4] have carried out experimental investigations on SMC – R26 with strain rates between 2x10

^-4

to 2x10

²

s

^-1

. The threshold strain for damage initiation between these strain rates increased from 0.25 – 0.77 %, and the ultimate strain increased from 1.4 – 2.1 %. The ultimate strain increased linearly whereas the threshold increased from 0.25 % to 0.3 % in a range from 2x10

^-4

to 1 s

^-1

and then from 0.3 % to 0.77 % in the remaining range.

A graph showing strain rate as a function of time for an element in a dynamic analysis shows that the strain rate is far removed from levels where changes in threshold and ultimate strain changes are. See Figure 3. Therefore any further strain rate effects are ignored.

No strain rate theory for BMC could be found.

(11)

Figure 2: Young's modulus decreases during increasing strain. The curve is important when estimating SMC damage level.

Figure 3: Strain and strain rate as compared to time history from an element in a dynamic structure analysis.

(12)

7

2.3 Neuber Correction

Hooke's law states the elastic relationship between stress and strain as long as the deformations are small. But in cases where the strain is in the over-elastic area the proportionality between stress and strain is lost. For stresses and strains in the elastic area a stress concentrations factor α

k

can be introduced

nom teo

k

σ

α

=

σ (12)

where σ

nom

is the strain that is not affected by any notch and σ

teo

the strain that is affected by the notch, see Figure 4.

In the case where the strain is in the over-elastic area the notch concentrations factor α

k

become invalid and instead we need different concentrations factors for stresses and strains. With the knowledge that real strain is greater than Hookian and real stress is smaller then Hookian, we introduce the stress-concentration factor,

σ

nom

α

_σ =

σ

^max

(13)

where σ

max

is the maximum stress in the notch and σ

nom

is the nominal stress.

Analogical the strain-concentrations factor is also introduced.

ε

nom

α

_ε

= ε

^max

(14)

In 1961, Neuber showed on a shear-strained prismatic body, if the notch angle is neglected, that the relationships between the stress and strain concentration factors are [3],

2

α

k

α

_σ _ε

= ₍₁₅₎

Using the equations (12), (13) and (14) and applying them to equation (15)

2 max

max

⎟⎟ ⎠

⎜⎜ ⎞

⎝

= ⎛

nom teo nom

nom

σ

σ ε

ε σ

σ (16)

further use σ

nom

= E

0

ε

nom

nom nom teo

E ε ε ε σ

σ

0 2 max

max

= (17)

and that σ

teo

= E

0

ε

teo

(13)

2 0 0

2 0 max max

) (

teo

teo E

E

ε ε

ε

σ

= =

(18)

Because α

k

is constant it follows that (18) becomes a hyperbola, as illustrated in Figure 4.

Figure 4: A graph illustrating the use of Neuber's hyperbola to estimate real stress and strain from Hookian stress and strain.

Suppose that theoretical strain (from Hooke's law) and the real stress-strain curve are known, see Figure 4. In order to find where the ε

max

point is for a certain strain, divide (18) by ε

max,

and compare σ

max

with the real stress-strain curve.

When they are equal, the true stress and strain is known. By doing this for all elements in a structure, one can obtain a nonlinear approximation of the strains.

Hookian stress

Real stress

Neuber’s hyperbola

(14)

9

3 Matlab

3.1 Program Structure

The program reads the strain data from the FEA program and transfers it by Neuber correction into non-linear strain after that the principal strains have been calculated. The damage level is then calculated by equation (3) and (11). Figure 5 shows a flow chart of the program structure and the numbers inform the way the program is used. Furthermore a short description of every program will be carried out. For the whole code see Appendix II

Figure 5: A flow chart illustrating the program structure. The numbers inform the way the functions are used.

3.2 Punch Files

The FE program Nastran is used for calculation of stresses and strains in the structure in this project, but of course any FE program could be used. By using Nastran it makes it possible to export so called punch files. The punch files can be printed for stresses or strains and since the damage level is calculated from the strain level, this is exported.

Depending on which analysis method is used dynamic or static, there are mainly two different layouts of the punch files. See Figure 6 for static layout and Figure 7 for the time dependent layout.

In principal the two different layouts are treated in the same way, except that there is no time vector in the static one. In both the figures in the blue, yellow and red fields represent element numbers, ε

xx

, ε

yy

and γ

xy

for Z

1

and ε

xx

, ε

yy

and γ

xy

for Z

2

where Z

n

is the thickness of the shell element. (In this case -1.500000E+00 for Z

1

and 1.500000E+00 for Z

2.

The element is 3 mm thick). The green field in Figure 7

represents the time.

(15)

Figure 6: The numbers highlighted represent data that is selected from the punch files and this data does not contain any time vectors. (Static calculation)

Figure 7: The numbers highlighted represent data that is selected from the punch files and this data contains a time vector. (Incremental solution)

(16)

11

Other differences in the two file types are further used in the Matlab file to separate them from each other. The file representing static calculations only has one new header each time a new element type is found, for example QUAD4 or TRIA3, whereas the file containing time has one new header for each new element. The time vector neither changes its length between each header nor results in frequency repeated headers.

3.3 Sorted Data

In Matlab, the function fopen is used to open the file for reading and writing.

Since fopen has the standard C library implemented, it guarantees only an addressing of 2

³¹

-1 bytes in the memory.

In order to carry out a dynamic analysis with reasonable resolution many time steps are necessary, this is the reason why a punch file in the format as briefly described presumably has a size over 2

³¹

-1 bytes. One way to override this problem is to read the important data using another program, and then sort it in a certain way. This could be done for many program languages, for example, Perl.

A rather simple and short code makes it possible to sort the data from Figure 7 to data as sorted in Figure 8. The eight columns contain: element number, time, ε

xx

, ε

yy

and γ

xy

for Z

1

and ε

xx

, ε

yy

and γ

xy

for Z

2

respectively. After the code has been sorted it is approximately 4.5 times smaller. See Appendix II for the Perl code.

Figure 8: Important data from the punch file. Here sorted to make it smaller.

(17)

3.4 Patran Format

To read the output data in a postprocessor, there must be some known structure.

One simple structure is the Patran format, which has "fef" as a suffix. The postprocessor Animator, for example, can then read the result file. See Figure 9.

The six different columns contain:

1. Element number, damage level for Z

1

in tension and value 0, 1, 2 or 3 depending if the damage value is none, in Z

1

, in Z

2

or both.

2. Damage level for Z

1

in pressure.

3. Damage level for Z

2

in tension.

4. Damage level for Z

2

in pressure

5. Maximum damage level (tension or pressure) for Z

1

. 6. Maximum damage level (tension or pressure) for Z

2

Figure 9: Output file written in Patran format.

(18)

13

3.5 The Matlab Program 3.5.1 Main

This is the file you start from (run file), and contains four important things:

**1. inName: The file or files Matlab is going to read, i.e. *.pch file/files in any** format as described above. Use the path name from root (or write a ~ for home section). Must be in cell format.

inName = {'~/plate/plate_SMC.pch'; 'file no. two'; 'file no. three'; etc …};

2. outName: The file Matlab writes with information about the damage level (the Patran file) This can only contain one name.

outName = '~/plate/plate_SMC.fef'

3. rows: Specify how many rows are contained in the entire file(s). The reason for this is that Matlab is not able to scan text from an entire file at once. If the inName cell contains more then one file, "rows" can be specified either by one value (all the files contain approximately the same no. of rows) or as a vector (the files have different number of rows).

rows = [100000]; or rows = [500000 40000 700000];

4. material: Specify 'SMC' or 'BMC' in cell format. If the inName cell contains more then one file, see no. 3.

material = {'SMC'}; or material = {'SMC'; 'BMC'; 'SMC'};

3.5.2 prepare

1. Enters the header in the outName file and "open" the inName file. (Make a pointer to it).

2. Splits up the rows in a scanLength vector if the file contains too many rows (more than ~100000).

3. Separates ordinary punch files (static and dynamic) and sorted files.

4. Retrieves the splitSort function or splitPunch function.

3.5.3 Punch functions splitPunch

This file has mainly three commissions:

1. Reads the file (for example 70000 rows before it stops) and saves it to a 5*readLength cell.

2. Examines how many headers there are in the file, and in which row. It also

inspects the length of the headers.

(19)

3. Splits the file and investigates where to start and stop reading. It also separates dynamic and static files.

if length(label) >= 4 & label(2)-label(1) == label(4) - label(3) type = 'dynamic'; % If the inName file contains time.

else

type = 'static'; % If the inName file does NOT contain time.

end

Note that the inName file must therefore contain (or read more than) at least four headers (elements) to be able to work.

It then passes on the start and stop information to the read function.

sortPunch

Is quite similar to the splitPunch function but reads the file into an 8*readLength cell. Then it passes on the start and stop information to the read function.

3.5.4 read

Reads the important information (see Figure 6, Figure 7 and Figure 8) and saves it to a cell, Z. See Figure 10 for the structure of Z.

time = 1 if the problem is static.

The information in the file repeats itself every sixth time, for the

non-sorted files, whereas the reading vector must be: start:6:stop. The value of the reading vector is

start:1:stop for the sorted data.

Figure 10: Structure of the Z - cell.

3.5.5 hub

Just passes on information between different functions. In this file is it easy to turn functions on and off, such as Neuber. The functions connected to the hub are written in a way that the input data is Z1 or Z2, you need to retrieve the functions twice.

3.5.6 principal

This transforms the strains to principal strains. The transformation matrix is used.

After that it returns to the hub function and calls on the Neuber function.

(20)

15

3.5.7 neuber

Transforms the linear strain to nonlinear strain by using the Neuber correction.

The neuber function only updates strains that are greater than the threshold value, for example, 0.25 % for SMC tension. Then back to the hub function and further to the damage function.

3.5.8 damage

Investigates if any values are grater than the threshold value and zeros all others.

If the punch file is from a dynamic simulation, it passes on the file to materialModelTime or to materialModel.

3.5.9 Material model functions materialModelTime

The first thing this function does is to pass the information on to the topValues function.

topValues

This function calculates:

• All values that are over the threshold value in an element. These values are sorted separately for each peak.

• The total time a strain peak exceeds the threshold. (If the element contains more then one peak, it comes as separate information for each peak).

• Maximum value for each peak.

• Mean value for each peak.

materialModelTime

From the topValues function, it simply takes the greatest strain value and passes it on to the materialModel function.

materialModel

To calculate the damage level, the materialProperties function needs to receive information about equation (11) or how Young's modulus decreases as a function of strain. The damage level is then calculated by equation (3).

3.5.10 writeToFile

Simply creates a matrix of the damage level information and sends it to the

outName file. This function also calculates the value 0, 1, 2 or 3 depending if the

damage value is zero, in Z

1

, in Z

2

or both.

(21)

3.5.11 materialProperties

Contains information about the material, i.e. Young's modulus, stress-strain curve, threshold value. See Figure 11.

Figure 11: A materialProperty function.

(22)

17

4 Example and Validation

4.1 Geometry

A plate with a hole was designed by Ansa to evaluate the Neuber correction and damage level calculation. The plate was then compared with Abaqus and an incremental material model that Magnus Oldenbo developed in his doctoral thesis.

See Figure 12 for the geometry.

Figure 12: Geometry of the plate on which numerical tests were performed.

The plate contains 774 CQUAD4 linear shell elements and 4 CTRIA3 linear shell elements with a thickness of 3 mm. The plate is then fixed on one short side and pulled with a force of 3000 N from the opposite side.

4.2 FEA Methods

In this project there are procedures to calculate the damage level:

1. Strain exports from Nastran into a punch file.

2. Matlab reads the punch file.

3. Transforms the strains to principal strains.

4. Transforms linear strains to nonlinear using Neuber's hyperbola.

5. The strains are inserted into equation (11) and the then into equation (3)

for damage calculation.

(23)

This is an alternative way of estimating the damage level that is delivered within certain margins. An improved way of calculating the damage level would rather be using an incremental method as Oldenbo et al. [5] did, which is the basis when validating this new routine.

One issue that the incremental method does take into account is the damaged area.

If there is any damage to some elements, Young's modulus would decrease and result in increased strain. In this method the elements are coupled to each other and neighboring elements can reach increased strain, which results in increased damage area.

4.3 Results of Neuber Correction

To evaluate if the Neuber correction gives a correct result, Abaqus is firstly solved using an ordinary linear solver and then Oldenbo's user defined model is applied.

The linear results are then used as input value in equation (19) and compared with the real stress-strain curve. When the values are the same, a nonlinear result is reached, and can be compared with Abaqus' incremental method, see Table 2. As one can see carried the Neuber correction out well and is therefore further

adopted.

Element

number Linear results Nonlinear results Neuber's results

36 5.66*10

^-3

5.95*10

^-3

6.05*10

^-3

1098 4.71*10

^-3

4.95*10

^-3

4.95*10

^-3

162 3.80*10

^-3

3.97*10

^-3

3.93*10

^-3

159 3.04*10

^-3

3.14*10

^-3

3.10*10

^-3

Table 2: Linear and nonlinear results from Abaqus compared with Neuber correction.

4.4 Results of Damage Level Simulation

Nastran is used for the analysis and during the calculation a punch file is exported

with strains as shown in Figure 6. See Appendix II for Nastran's setup file. Matlab

then reads the punch file, makes some calculations and writes a text file in Patran

format. This procedure is repeated twice, once with the Neuber correction and

once without the Neuber correction. Animator is then used as postprocessor. See

Figure 13 for the results without the Neuber correction and Figure 14 when using

the Neuber correction.

(24)

19

Figure 13: Results from a damage calculation without the Neuber correction.

Figure 14: Results from a damage calculation with the Neuber correction.

(25)

As one can see in Figure 13 and Figure 14 is the difference between them only in the damage level in the elements. Some may think it is unnecessary to do the Neuber correction, but in fact it gives a more realistic picture over the damage since the gradient between the elements is bigger. In the reality there are no elements, and therefore no gradients, whereas the damaged area can be easier to visualize.

With the aid of Oldenbo's user defined model, Abaqus solves the damage level and postprocesses in the Abaqus viewer. See Figure 15 for results which have good agreement with figure 14.

Figure 15: Results from a damage calculation using Abaqus and Oldenbo’s user defined method.

Some elements are picked to give the exact value and the result can be seen in Table 3. In some cases is the Neuber correction over estimated compared to Abaqus and the used defined method, and in some cases is it underestimated. The differences are however almost smaller compared to the calculation without Neuber correction.

Element number

User defined method

With Neuber correction

Without Neuber correction

36 0.051 0.052 0.045

37 0.042 0.046 0.040

1098 0.033 0.032 0.028

164 0.015 0.014 0.012

162 0.018 0.015 0.014

160 0.015 0.013 0.012

Table 3: Some damage levels from the plate.

(26)

21

5 Conclusions

In principle anyone can plot the principal strains in a postprocessor and evaluate if a strain is greater than the threshold level. The disadvantage with this solution is that the operator has to evaluate at least four principal strains, two directions for every element face. Some postprocessors do not even have this function. Instead an interpolation between the element faces is used and the strain level can be underestimated. Under these conditions the damaged area on complex geometry can be difficult to discover. This is why an application tool makes life simpler.

The method that has been developed makes it possible to translate strains over threshold levels to a damage level D, which gives the operator a more

comprehensive view of the condition of the structures. There are several benefits with this method, for example, specifically highlighting the damage level. If the model in Nastan is of dynamic character, using this method will select the most severe strain level from the strain history and translates it into a damage level.

Some may think that it is unnecessary to calculate the damage level and it is better to give an integer value of these elements that has a strain level exceeding the threshold instead. This is not a bad idea, and in fact it is used, but only as a complement to the damage level. But if the operator only has this function, any idea of where the cracks begin would vanish.

A static load case is used as an example to validate the method, as seen in chapter

4.

(27)

References

1.

Oldenbo M., Mattson D., Varna J., Berglund L. A. (2004). Global Stiffness of a Panel Considering Process Induced Fiber Orientation. Journal of Reinforced Plastic Composites, Vol. 23 No. 1.

2.

Oldenbo M., J. Varna. (2003). A constitutive model for non-linear behavior of SMC accounting for linear viscoelasticity and micro-damage. Annual Technical Conference of the society of Plastics Engineers (ANTEC´03), Paper no 345.

3.

Neuber H. (1961). Theory of Stress Concentration for Shear-Strained Prismatic Bodies with Arbitrary Nonlinear Stress-Strain Law. ASME Journal of Applied Mechanics 28.

4.

Jendli Z., Meraghni F., Fitoussi J., Baptiste D. (2004). Micromechanical analysis of strain rate effect on damage evolution in sheet molding compound composites.

Composites: Part A 35 (779-785).

5.

Oldenbo M., Varna J. (2005).An incremental 2D constitutive model accounting for linear viscoelasticity and damage development in short fibre composites. International Journal for Numerical Methods in Engineering. Volume 64, Issue 11 (1509-1528)

6.

Jendli Z., Fitoussi J., Meraghni F., Baptiste D. (2005). Anisotropic strain rate effects on the fibre-matrix interface decohension in sheet moulding compound composites.

Composites Science and Technology 65 (387-393).

7.

Thilo Trautwein F. Numerical Validation and Application of the Neuber-Formula in FEA- Analysis. ACES GmbH, Filderstadt, Germany.

8.

Urm H. S., Yoo I. S., Heo J. H., Kim S. C., Lotsberg I. (2004). Low Cycle Fatigue Strength Assessment for Ship Structures. 9th Symposium on Practical Design of Ships an Other Floating Structures Luebech-Travemuende, Germany.

(28)

i

I Appendix I

I.I Obtaining Tensile Stress-Strain data for BMC

Several specimens with dimensions 250 x 23 x 3.6 mm were produced from the tailgates having BMC 1100 material delivered from MCR (Mixt Composites Recyclables). On each sample two circular reflection pieces about one millimeter in diameter was positioned approximately 25 mm from the middle of the sample in each direction. The reflection pieces where then illuminated by two optical extensometers for extension measurement.

The data measured was collected by the software testXpert and then exported to Excel for calculations of mean value and trend line. See Figure i.

Figure i: A stress-strain curve for BMC 1100 material. A trend line was estimated from the mean value curve. The initializing Young's modulus is 14500 MPa. A strain rate of 2x10^-4s^-1 was used.

From the tension trend line, a pressure line was estimated in the same way as for SMC material, i.e., pressure data is twice the value as tension data.

From the stress strain curve several tangential Young's modulus were calculated and represent damaged Young's modulus. This is incorrect, because of the viscoelasticity. However, the viscoelastic effect during a rapid test is believed to be small.

Mean value

(29)

Figure ii: The trend line from Figure i. Pressure data is estimated as twice as much as tension data.

Figure iii: The damage curves are calculated from tangential Young’s modulus in tensions from Figure i. The damage in pressure is twice as in tension. By doing this the viscoelasticity will not be observed and the damage curve is overestimated.

(30)

iii

II Appendix II

II.II Perl script

The Perl script compresses the incremental punch file and writes a startup.m file.

The startup.m file makes it possible to start Matlab, which runs the startup.m file automatically, and then terminates Matlab.

#!/opt/perl/bin/perl

print "Which punch-file do you want to run? ";

chomp($punch=<STDIN>);

open(PCH,"$punch")|| die "Can't open \"$punch\" ! \n";

print "What do you want to call your damage file (In Patran format)? ";

chomp($fefFile=<STDIN>);

$flag="ANY";

while ($flag ne "MAT") {

print "Which material do you have (SMC or BMC)? ";

chomp($mat=<STDIN>);

if ($mat eq "SMC") { $flag="MAT";

} elsif ($mat eq "BMC") { $flag="MAT";

} }

$flagType="AnyType";

while ($flagType ne "TYPE") {

print "Type '(i)ncremental' or '(s)tatic': ";

chomp($type=<STDIN>);

if ($type eq "i") { $flagType="TYPE";

} elsif ($type eq "s") { $flagType="TYPE"

} }

chomp($pwd=`pwd`); #Check in which path you are.

if ($type eq "i") {

$compr="/lfs/temp/compressed.pch";

open (COMP,">$compr")|| die "Can't open \"$compr\" ! \n";

$[ = 1; # set array base to 1

print "\nA file called 'compressed.pch' is now created in

$compr\n";

$rows=0;

(31)

while (<PCH>) {

chomp; # strip record separator @Fld = split(' ', $_);

# header

if ($_ =~ /^\$TITLE/) { $_ = Getline0();

$_ = Getline0();

$elmid = $Fld[4];

printf COMP "%8s\n", 'ele';

$_ = Getline0();

$rows++;

}

$Time = $Fld[1];

$e1xx = $Fld[3];

$e1yy = $Fld[4];

$_ = Getline0();

$e1xy = $Fld[2];

$_ = Getline0();

$e2xx = $Fld[2];

$e2yy = $Fld[3];

$e2xy = $Fld[4];

$_ = Getline0();

$rows++;

printf COMP "%8d %13.6E %13.6E %13.6E %13.6E %13.6E %13.6E

%13.6E\n", $elmid,

$Time, $e1xx, $e1yy, $e1xy, $e2xx, $e2yy, $e2xy;

}

close(PCH);

close(COMP);

} elsif ($type eq "s") {

print "The number of newline characters in the input file calculates\n";

$r=`wc -l $punch`;

@row=split(' ', $r);

$rows=$row[1]

}

print "A 'startup.m' file is now created in $pwd\n";

$matlfile="startup.m";

open (MATL,">$matlfile")|| die "Can't open \"$matlfile\" ! \n";

printf MATL "clc\n";

(32)

v

printf MATL "clear\n\n\n";

if ($type eq "i") {

printf MATL "inName = {'$compr'};\n";

} elsif ($type eq "s") {

printf MATL "inName = {'$punch'};\n";

}

if ($fefFile =~ /^\//i) {

printf MATL "outName = '$fefFile';\n";

} else {

printf MATL "outName = '$pwd/$fefFile';\n";

}

printf MATL "rows = [$rows];\n";

printf MATL "material = {'$mat'};\n\n\n";

printf MATL "path(path,'/vcc/homes/pccf638/Matlab/forslag3');\n";

printf MATL "[Z, eNr ,time, E_Neuber] = prepare(inName, outName, rows, material);\n\n\n";

printf MATL "exit\n\n";

close(MATL);

print "Initializing Matlab\n";

exec("/vcc/ans/sim/Matlab/matlab7.0/bin/matlab");

sub Getline0 {

if ($getline_ok = (($_ = <PCH>) ne '')) { chomp; # strip record separator

@Fld = split(' ', $_);

} $_;

}

II.III Nastran

NASTRAN BUFFSIZE=65537

TIME 999999 $Max analysis CPU run-time SOL 101 $Static solution sequence CEND

TITLE = Plate_hole SUBTITLE = Subtitle

ECHO = NONE $supress input deck echo in .f06-file

$ Output: PLOT=>op2, PRINT=>f06, PUNCH=>pch(ascii-file)

$<---Max 72 characters per row!--->

$SET 99 = 1, 2 THRU 4 $node or element set definition for output DISPLACEMENT(PLOT) = ALL $displacements

STRESS(PLOT) = ALL $vonMises stress STRAIN(PLOT) = ALL $vonMises strain SPC = 1 $constraints SUBCASE 1 $define loadcases

(33)

LABEL = Subcase1 LOAD = 1

REPCASE 3 $use if a result type should be written $to > 1 outfile

STRAIN(PUNCH,FIBER) = ALL $Writes the punch file BEGIN BULK

PARAM,K6ROT,1. $Necessary for models with 1D elements $connected to SHELLs

$PARAM,SNORM,0 $Default=20. See v2001 Release Guide for $discussion

$Stiffness matrix

PARAM,PRGPST,NO $Supress AUTOSPC removed singularities $printout

$PARAM,MAXRATIO,1.E+7 $Use with caution! Default=1.E+7

$Large MAXRATIO indicates ill-conditioned PARAM,GRDPNT,0 $0=Print model mass and inertia data.

Default=-1

$PARAM,COUPMASS,-1 $-1=lumped mass matrix, 1=coupled (full) $mass matrix

$Output

PARAM,POST,-2 $Output format supported by Animator and $HyperMesh

$ Concatenation of bulkdata file.

INCLUDE 'plate.nas'

II.IV Matlab

II.IV.I startup.m or main.m

clear clc

inName = {'/vcc/homes/pccf638/plate/mjukt/plate.pch'};

outName = '/vcc/homes/pccf638/Matlab/plate.fef';

rows = [13258];

material = {'SMC'};

path(path,'"Any path"')

[Z, eNr, time, E_Neuber] = prepare(inName, outName, rows,...

material);

exit

(34)

vii

II.IV.II Prepare.m

function [Z, eNr ,time, E_Neuber] = prepare(inName, outName, ...

rows, material)

%Prepares to read the inName file. It also prepares the outName file.

%---Runs the different inName files after each other--- for i = 1:length(inName)

if length(rows) > 1

scanLength = countLength(rows(i)); %Sub function else

scanLength = countLength(rows) end

if length(material) > 1 %The material can bee one or ar a cell %containing one material for each

%file.

n = i;

else n = 1;

end

if i == 1 %The first time makes a new file action = {material{n}; 'Write'};

else %the other times are the data added to the file.

action = {material{n}; 'Add'};

end

Name = inName{i}

if strcmp(action{2},'Write') == 1;%Makes the outfile and %writes the header.

fidOut = fopen(outName,'w');

fprintf(fidOut,'Crack Initiation\n 7\nC1 = T in Z1, C2 = Pr in Z1, C3 = T in Z2, C4 = Pr in Z2, C5 = max T\n Pr in Z1, C6 = max T or Pr in Z2, C7 = 1,2,3 for damage in Z1 Z2 or Z1 and Z2\n');

fclose(fidOut);

%Send some messages

elseif strcmp(action{2},'Add') == 1;

disp('You are now adding the results to the file') end

%info contains material, Add or Write, and outname.

info = action;

info{3} = outName;

temp{1,1} = 0; %Save relevant data for the next time.

type = splitType(Name); %Sub function

(35)

[fidIn message]= fopen(Name,'r'); %Open the file.

for n = 1:length(scanLength) readLength = scanLength(n);

percent = round((n-1)/length(scanLength)*100);

disp([num2str(percent), ...

' percent of the calculation is finished'])

if fidIn == -1 %send a message if it fails to open the %file.

message end

if strncmp(type, 'sort',4) == 1

[Z,time,temp,eNr,E_Neuber] = splitSort(readLength,...

fidIn, n, temp, info);

elseif strncmp(type, 'punch',4) == 1

[Z,time,temp,eNr,E_Neuber] = splitPunch(readLength,...

end

temp = temp;

end

fclose(fidIn);

end

%---Splitt up the rows in a scanLength vector if the file

contains too many rows--- function scanLength = countLength(rows)

read = 70000; %To read ca 70000 rows is ok.

if rows <= read scanLength = -1;

else

n = rows/read;

integer = round(n);

scanLength = ones(1,integer).*read;

scanLength(length(scanLength)) = -1; %It must end with -1. It %then reads to the end %of the file.

end

function type = splitType(Name)

[fidIn2 message] = fopen(Name); %Open the file for reading if fidIn2 == -1

message

(36)

ix

end

C = textscan(fidIn2, '%s %s %s %s %s',20); %Only 20 rows.

for i = 1:length(C{1});

if strncmp(C{1}(i),'$TITLE',4) == 1;

type = 'punch'; %If the file has NOT been sorted elseif strncmp(C{1}(i),'ele',3) == 1;

type = 'sort'; %If the file has been sorted end

end

fclose(fidIn2);

II.IV.III splitPunch.m

function [Z, time, temp, eNr, E_Neuber] = split(readLength,...

%Split the file and pass it on to the read function and the %hub function.

%---Read all information from the file--- if n == 1 %Only the first time.

%Scans the file as "strings" and make a 5*readLength cell C = textscan(fidIn, '%s %s %s %s %s',readLength);

else

D = textscan(fidIn, '%s %s %s %s %s',readLength);

C = temp; %Use the temp in the beginning

for rows = 1:5

%Add the new information

C{1,rows}([1:length(D{1})]+length(temp{1})) = D{rows};

end end

%---Find where all the headers are--- label = 0;

j = 0;

labelTest = strncmp(C{1}(i),'$TITLE',4); %Find all times the %file contains $TITLE if labelTest == 1;

j = j+1;

label(j) = i;

end end

%---Calculate the length of the header--- if label(1) ~= 0

flag = 0;

(37)

a = 0;

for nn = label(1):label(1)+20 a = a + 1;

test = strncmp(C{1}(nn),'-CONT-',4); %The first time %-CONT- turns up if test == 1 & flag == 0

s(a) = 0;

flag = 1;

end end

stepsInHeader = length(s)-1;

end

%---When the file contains time, the header repeat itself--- if length(label) >= 4 & label(2)-label(1) == label(4) - label(3) type = 'dynamic' %If the inName file contains time.

else

type = 'static' %If the inName file does NOT contain time.

end

if strncmp(type, 'static', 4) == 1

%---Read and write relevant data from C--- if n == 1 %Only the first time.

for a = 1:length(label)

if a <= length(label)-1 %Not the last element if a == 1; %The first element

%There are "stepsInHeader" in relevant data in %the cell C before the first number.

start = stepsInHeader;

stop = label(2)-1; %$LABEL is the first data %The reading function.

[eNr, Z, time] = read(start,stop,type,C);

%The hub function

[Epsilon, E_Neuber] = hub(Z,eNr,info,time);

else

start = label(a) + stepsInHeader - 1;

stop = label(a+1)-1;

end

else %The last header

start = label(length(label)) + stepsInHeader - 1;

%Subfuction

m = steps(C{1}(length(C{1})-5:length(C{1})));

if readLength == -1 %Reads to the end of the file.

stop = length(C{1});%Scans to the end of the %file.

(38)

xi

else

%m is the number of steps there are to a new %element. The element repeats itself every7:th %time so m-7 means that you are at the

%beginning.

stop = length(C{1})+(m-7);

for k = 1:5

%Save the information as stated temp{k} = C{k}(stop:length(C{k}));

end end

else %When n ~= 1

if label == 0 %When no header can be found.

start = 1;

m = steps(C{1}(length(C{1})-5:length(C{1})));

if readLength == -1 stop = length(C{1});

else

for k = 1:5

temp{k} = C{k}(stop:length(C{k}));

end end

[eNr,Z, time] = read(start,stop,type,C);

else %When at least one head can be found.

for b = 1:length(label) + 1 if b <= length(label)-1

if b == 1

start = 1;

stop = label(1)-1;

[eNr,Z,time] = read(start,stop,type,C);

[Epsilon, E_Neuber] =hub(Z,eNr,info,time);

else

start = label(b) + stepsInHeader - 1;

stop = label(b+1)-1;

[Epsilon, E_Neuber] =hub(Z,eNr,info,time);

end

else %Pass the last head.

if b == 1 start = 1;

stop = label(length(label)) - 1;

else

start = label(length(label))...

+ stepsInHeader - 1;

(39)

m = steps(C{1}(length(C{1})- 5:length(C{1})));

if readLength == -1 stop = length(C{1});

else

for k=1:5

temp{k} = C{k}(stop:length(C{k}));

end end end

end end end end

elseif strncmp(type, 'dynamic',4) == 1

%---Read and write relevant data from C--- for a = 1:length(label)

start = stepsInHeader;

stop = label(2)-1; the header.

else start = label(a) + stepsInHeader - 1;

stop = label(a+1)-1;

end

else %The last header

if readLength == -1 %Scans to the end of the file.

start = label(a) + stepsInHeader - 1;

stop = length(C{1});

else

for k = 1:5

start = label(length(label));

temp{k} = C{k}(start:length(C{k}));

end

end end

%---Find how many steps there are to the new element---

(40)

xiii

function m = steps(C) flag = 1;

m = 1;

list = cellstr(['$TIT';'$SUB';'$LAB';'$ELE';'$REA';

'$SUB';'$ELE';'FIBE';'-CON']);

while flag == 1 for i = 1:9

test(i) = strncmp(C(m),list(i),4);

end

if sum(test) == 0 flag = 0;

end

m = m + 1;

end

II.IV.IV splitSort.m

function [Z, time, temp, eNr, E_Neuber] = split(readLength,...

%Slit the file and passes it further to the read function and the

%hub function.

type = 'sort';

%---Read all information from the file--- if n == 1 %Only the first time

%Makes a 8*readLength cell

C = textscan(fidIn, '%s %s %s %s %s %s %s %s',readLength);

else

D = textscan(fidIn, '%s %s %s %s %s %s %s %s',readLength);

C = temp; %Use the temp in the beginning for rows = 1:8

%Add the new information

C{1,rows}([1:length(D{1})]+length(temp{1})) = D{rows};

end end

%---Find where all the heads are--- label = 0;

j = 0;

labelTest = strncmp(C{1}(i),'ele',3); %Find all times the file %contains $TITLE

if labelTest == 1;

j = j+1;

label(j) = i;

end end

%---Read and write relevant data from C--- for a = 1:length(label)

(41)

start = label(1) + 1; %Start to read at the second %row.

stop = label(2) - 1; %Stop one row above label "ele".

%The read function.

%The hub function.

else

start = label(a) + 1;

stop = label(a+1) - 1;

end

else %The last header if readLength == -1

start = label(length(label)) + 1;

stop = length(C{1});

else

for k = 1:8 %Save the data that is between element n %and n+1.

start = label(length(label));

temp{k} = C{k}(start:length(C{k}));

end end end end

II.IV.V read.m

function [eNr, Z, time] = read(start,stop,type,C)

%Reading relevant data from a 5xn cell. Outputs are elementNr, Z

%and time.

Z contains epsilon_x, epsilon_y, epsilon_xy for Z1 and %Z2 respective.

if strncmp(type, 'dynamic' ,3) == 1

%---Makes vectors containing time and element number--- time = str2num(char(C{1}([start:6:stop])));

eNr = str2num(char(C{4}(start-2)));

%---Reads relevant data from the file--- vector1 = [start:6:stop];

vector2 = [start+1:6:stop];

(42)

xv

Z1(:,1) = str2num(char(C{3}([vector])));

Z1(:,3) = str2num(char(C{2}([vector2])));

%---Z contains epsilon_x, epsilon_y, epsilon_xy--- Z{1,1} = Z1;

Z{1,2} = Z2;

elseif strncmp(type, 'sort' ,3) == 1

%---Makes vectors containing time and element number--- eNr = str2num(char(C{1}(start)));

time = str2num(char(C{2}([start:stop])));

%---Reads relevant data from the file--- vector = [start:stop];

Z{1,2} = Z2;

elseif strncmp(type, 'static' ,3) == 1 time = 1;

eNr = str2num(char(C{1}([start:6:stop-5])));

%---The relevant information are 0, 1 and 3 rows below start.-- vector1 = [start:6:stop-1];

Z1(:,1) = str2num(char(C{3}(vector1)));

Z{1,2} = Z2;

end

(43)

II.IV.VI hub.m

function [Epsilon, E_Neuber] = hub(Z, eNr, info, time)

%Pass information between different files.

material = info{1};

outName = info{3};

%Make the strains to principal-strains

[Epsilon{1,1}, theta{1,1}] = principal(Z{1,1}); %Z1 [Epsilon{1,2}, theta{1,2}] = principal(Z{1,2}); %Z2

%Update strains to non-linear with Neuber.

[E_Neuber{1,1}] = neuber(Epsilon{1,1}, material); %Z1 [E_Neuber{1,2}] = neuber(Epsilon{1,2}, material); %Z2

%Calculate the damage level

[D_TensionZ1, D_PressureZ1, D_TotZ1] = damage(E_Neuber{1,1}, ...

eNr, time, material);

[D_TensionZ2, D_PressureZ2, D_TotZ2] = damage(E_Neuber{1,2}, ...

eNr, time, material);

%Write relevant data to the outName file.

[outData] = writeToFile(outName, eNr, D_TensionZ1,...

D_PressureZ1, D_TensionZ2, D_PressureZ2, D_TotZ1, D_TotZ2);

II.IV.VII principal.m

function [Epsilon, theta] = principal(Z)

%Transforming strains to principal strains.

epsilon = Z;

%---Calculate the radii of Mohr's circle--- R = (((epsilon(:,1)-

epsilon(:,2))./2).^2+(epsilon(:,3)./2).^2).^(1/2);

%---Find values where e1 < e2--- sign = ones(length(R),1);

flag = zeros(length(R),1);

sign(find(epsilon(:,1)-epsilon(:,2) < 0),1) = -1;

flag(find(epsilon(:,1)-epsilon(:,2) < 0),1) = 1;

%---Calculate the angle--- notZero = find(R ~= 0); %We cannot divide by zero.

theta = zeros(length(R),1);

(44)

xvii

theta(notZero) = ...

1/2*asin(epsilon(notZero,3)./(2*R(notZero))).*sign(notZero)...

+ pi/2*flag(notZero);

for i = 1:length(theta) c = cos(theta(i,:));

s = sin(theta(i,:));

%---Transformation matrix--- T{i,1} = [c^2 s^2 s*c;

s^2 c^2 -s*c;

-2*s*c 2*s*c c^2-s^2];

%---Principal strains--- Epsilon(i,:) = T{i,1}*epsilon(i,:)';

end

II.IV.VIII neuber.m

function [E_Neuber] = neuber(e1, material)

%Update linear values to nonlinear by Neubers hyperbola.

%---Collect data from the materialProperties function---

%ut = {E0; epsilon; sigma; epsilonPr; sigmaPr; E_tilde; threshold;

epsilonBreak};

in = materialProperties(material); %SMC or BMC material E0 = in{1};

epsilon = in{2};

sigma = in{3};

sigmaPr = in{5};

epsilonPr = in{4};

threshold = in{7}(1);

thresholdPr = in{7}(2);

epsilonBreak = in{8}(1);

epsilonBreakPr = in{8}(2);

%Create a polynomial from the stress-strain curve.

p = polyfit(epsilon,sigma,2);

pPr = polyfit(epsilonPr,sigmaPr,2);

%---Neuber strains--- e_max = threshold:0.000001:epsilonBreak*2;

e_maxPr = [abs(thresholdPr):0.000001:abs(epsilonBreakPr)].*-1;

%---Find all values over threshold for e1 and e2--- for i = 1:2

thres{1,i}(:,1) = find(e1(:,i) >= threshold);

thresPr{1,i}(:,1) = find(e1(:,i) <= thresholdPr);

end

%---Update values over threshold--- for j = 1:2 %e1 and e2

(45)

%Tension

if length(thres{1,j}) >= 1 %If there are any e1/e2 values that

%is over threshold in tension.

for k = 1:length(thres{1,j}) %e1(thres{1,j}(k,1),j)

if e1(thres{1,j}(k,1),j) <= (epsilonBreak-0.001) %Neubers hyperbola

S_euler = e1(thres{1,j}(k,1),j).^2*E0./e_max;

%The stress-strain curve stresses S_nonLin = polyval(p,e_max);

%If the difference between them are <= 0.01

N = e_max(find(abs(S_euler-S_nonLin) <= 0.01));

if length(N) >= 1

%Sometimes can it be more then one value.

Neuber{1,j}(thres{1,j}(k,1),1) = mean(N);

else

end else

disp('You have a value that overrides epsilonBreak')

end end end

%Pressure

if length(thresPr{1,j}) >= 1 for k = 1:length(thresPr{1,j})

if e1(thresPr{1,j}(k,1),j) >= (epsilonBreakPr)

S_euler = e1(thresPr{1,j}(k,1),j).^2*E0./e_maxPr;

S_nonLin = polyval(pPr,e_maxPr);

NPr = e_maxPr(find(abs(S_euler-S_nonLin)...

<= 0.01));

if length(NPr) >= 1

Neuber{1,j}(thresPr{1,j}(k,1),1) = mean(NPr);

else

end else

disp('You have a value that overrides epsilonBreak')

end end end end