UPTEC F 19023
Examensarbete 30 hp
10 Juni 2019
A SIMULATION WITH FINITE
ELEMENTS TO MODEL STEEL SHEET
SLITTING
A Master Thesis in Engineering Physics
Adam Ahlgren Peters
Teknisk- naturvetenskaplig fakultet UTH-enheten Besöksadress: Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0 Postadress: Box 536 751 21 Uppsala Telefon: 018 – 471 30 03 Telefax: 018 – 471 30 00 Hemsida: http://www.teknat.uu.se/student
Abstract
A SIMULATION WITH FINITE ELEMENTS TO
MODEL STEEL SHEET SLITTING
Adam Ahlgren Peters
A steel slitting process is simulated using FEM (Finite Element Method) in order to see potential defects along the edge in a steel sheet after it has been cut. The model's results were compared to microscope images of the steel sheet in order to verify accuracy. The purpose is conceptual and to find a model that successfully simulates a steel cutting process and (hopefully) how the edge depends on different parameters. The model developed seems to achieve this task, and a more thorough calibration of the model could result in (more) optimal parameters for the machine to use.
Popul¨
arvetenskaplig sammanfattning
Det ¨okande kravet p˚a allt mer h¨ogh˚allfasta st˚al (st˚al som t˚al att deformeras mycket) har lett till att st˚alen utvecklas till att bli allt mer h˚allbara och t˚aliga. Detta kan t¨ara p˚a utrustningen vid klippning, men vad som framf¨orallt unders¨oks i denna rapport ¨ar huruvida en modell kan utvecklas f¨or en specifik klipp-process (spaltning), d¨ar fokus ligger p˚a hur klippkanten ser ut n¨ar st˚alet klipps. Kantens skada p˚averkar st˚alets
h˚allfasthet mycket. F¨or att inse detta kan man j¨amf¨ora ett papper med ett litet hack i kanten med ett papper som har hel kant; papperet med skadan kommer rivas is¨ar mycket enklare n¨ar man b¨orjar dra i papperet l¨angs med kanten. Det ¨ar d¨arf¨or av intresse att veta mer om kantskador.
De simulerade resultaten har j¨amf¨orts med mikroskopbilder p˚a klippkanten av st˚alet hur det ser ut n¨ar det klippts. En komplett kalibrerad modell skulle kunna anv¨andas f¨or att optimera parametrar f¨or klippningen. F¨orhopnningsvis skulle d˚a klippskadan minimeras, vilket skulle inneb¨ara ¨
okad h˚allfasthet f¨or st˚alet eftersom sprickor l¨attare propagerar d¨arifr˚an det redan finns sm˚a sprickor eller oj¨amnheter i materialet.
Acknowledgements
I would like to thank my official supervisor Dr Saed Mousavi for all the help with material-related topics, my unofficial supervisor Dr Lars Olovsson (IMPETUS’s CEO) for all the help with simulation-related
questions. Also a thanks to Prof. Per Isaksson for the valuable feedback on my report.
Table of contents
1 Introduction 5
1.1 Background . . . 5
1.2 Purpose . . . 5
2 Theory 6 2.1 Parameter convention and definition . . . 6
2.2 Slitting. . . 6
2.3 Stress-strain . . . 8
2.4 Material model . . . 9
2.5 Damage model . . . 9
2.6 Residual stresses . . . 10
2.7 Finite Element Method . . . 12
2.7.1 Derivation . . . 13
2.7.2 Finite difference scheme . . . 15
2.8 FEM software . . . 15
2.8.1 Time stepping . . . 16
2.8.2 Mass scaling . . . 17
3 Method 18 3.1 Material . . . 18
3.2 Measurements and parameters . . . 18
3.2.1 Stress-strain . . . 18
3.2.2 Slitting machine . . . 18
3.2.3 Residual stress . . . 19
3.3 Material model adaptation . . . 20
3.3.1 Parameter calibration. . . 20
3.3.2 Choice of damage model . . . 20
3.4 Software setup. . . 20
3.4.1 Element type . . . 21
3.4.2 Element order . . . 21
3.4.3 Element density . . . 21
3.4.4 Time stepping and mass scaling . . . 23
3.4.5 Boundary conditions . . . 23
3.4.6 Initial conditions . . . 23
4 Results 24 4.1 Stress-Strain curve . . . 24
4.2 Scanning Electron Microscope . . . 24
4.2.2 Shear zone (bright zone) . . . 25
4.2.3 Secondary shear zone (fracture zone) . . . 26
4.2.4 Burr zone . . . 27
4.2.5 Plastic deformation zone . . . 27
4.3 Two-dimensional model . . . 29 4.3.1 Energy balance (2D) . . . 30 4.4 Three-dimensional model . . . 31 4.4.1 Energy balance (3D) . . . 32 4.4.2 Different sharpness . . . 33 4.4.3 Different gap . . . 34
4.4.4 Different damage tolerance . . . 35
5 Discussion 36 5.1 Software . . . 36
5.2 Material model . . . 37
5.3 Stress-Strain curve . . . 37
5.4 Scanning Electron Microscope . . . 37
5.5 Two-dimensional model . . . 38
5.6 Three-dimensional model . . . 38
5.7 Future improvements . . . 39
1
Introduction
Metallurgy has until recently depended solely upon trial and error. The recent development in computer science -and computational power
specifically- has led to a new complement to the more traditional approach. By successfully implementing a Finite Element Model (FEM) [1], one can narrow down the necessary experiments needed for development and hence speed up the development process. Ultimately, the strength along with other material properties of the steel alloy is what truly matters. If these can be precisely estimated numerically, the trial and error process would become far more effective.
1.1
Background
Due to the increased demand in high-strength steels (mainly in the automotive industry), the need to precisely simulate the manufacturing process has increased as well. Steel plates are among the most common types of steel elements produced in the world. During the last two decades the yield strength of the hardest steel grades has had an extraordinary improvement. Materials with yield strength of 1.5 GPa and higher, are now commonly used, especially in automotive industry. As a result, the
equipment used to cut such steel has seen a significant increase in wear.
1.2
Purpose
In short, the purpose is to construct a model based on observation of an industrial process that involves steel slitting. The simulations are to be compared to images of the steel’s edge in order to verify accuracy. The goal is to give a hint of the optimal parameters for the machine involved (e.g slitting gap and slitting overlap).
The purpose of this thesis is to lay a foundation for further material studies and process improvements. The study should be seen as
conceptual; the method matters more than the precision of the calibration that depends upon experimental data. Ultimately, a well-developed
method would improve the pace at which improvements can be made, since it narrows down the experiments considered necessary.
As a first initial step towards minimising the damage on the produced steel, the goal is to simulate the slitting process as precisely as possible. Once such a model exists, it should provide an easier way to minimise the waste (and wear) through simulation experiments, in union with real-life experimental controls.
2
Theory
Since this project has roots in both the material and the computational sciences, the relevant theory will be presented for both of these.
2.1
Parameter convention and definition
For simplicity the variable, abbreviation and parameter names are defined as in Table 2.1 unless specified otherwise. SI-units are used. Text in bold symbolises a matrix of any size. Some variables may be reused from section to section and might not refer to the same value.
Symbol Meaning Unit ∆ Difference Operational
∆l/L N umeric
ρ Density kg/m3
ν Poisson’s ratio N umeric
σ Stress P a
ω Ang. frequency rad/s a Acceleration m/s2 c Wave velocity m/s E Young’s modulus P a F Force N L Length m r Radius m m Mass kg t Time s v Velocity m/s
2.2
Slitting
Slitting process (knowledge) is of high importance, as it is an essential part of the manufacturing process. Since customers’ need are in focus and different costumers have different needs, the material that is sold might have to be of different sizes. The slitting machine ensures the width is correct and according to the customers’ specifications. The slitting process is shown in Figure 1 schematically and in Figure2 in reality, where one can see how a steel sheet is split into several slimmer sheets.
Figure 1: Slitting process shown schematically. The metal sheet (or mate-rial to be slitted) is the grey, horizontal line and exits out of the plane if one considers this paper to define a 2D plane. The red, vertical lines repre-sent the circular cutters.
Figure 2: Slitting process. The material is slitted into slimmer chunks by the circular cutters.
The process is important to understand, since the process in fact inflicts some damage to the edges of the material. Edge damage is kept at a minimum (for obvious reasons), but still has a high impact on the
durability of the material. For instance, the burr as seen in Figure 3 should be kept inward if the material is to be folded. Otherwise the small burr may fracture and severely impact the durability of the material, since a larger fracture then more easily may propagate deeper into the material.
Figure 3: Characteristic edge of a slitted material.
2.3
Stress-strain
One of the most, if not the most essential parts of material characterisation is the stress-strain test. This test shows the tensile strength, yield strength, elasticity, plasticity, and failure for the material tested. A typical
stress-strain curve and setup can be seen in Figure 4.
Figure 4: Typical engineering stress (green) plotted together with the true stress (blue). As can be seen, they are identical up until the necking occurs. At 1 we have the yield strength and at 2 the maximum tensile strength. Red cross symbolises fracture. E is the slope of the elastic zone. Figure based on stress-strain figures in [2] and should be seen as purely schematic.
The material is attached in opposing ends. A linearly increasing force is applied to drag the material apart. Initially, this will result in an elastic response of the material (much like a spring) and will follow Hooke’s Law [3]. The maximum value when the material is still elastic is called yield
strength. Typically, high performance steel (such as SSAB’s Docol series [4]) have a very high yield strength. Needless to say, this property is very desirable in products that needs to withstand high degrees of deformation without failure or plastic deformation. At some point however, the material will no longer be elastic but will instead become plastic. This means the material will start to (non-reversibly) deform as more force is applied, until it ruptures. An ideally plastic material will have a plastic part of the stress-strain curve that is parallel to the x-axis. Steel however, continues to increase in stress until it reaches a certain maximum (ultimate tensile strength). At failure the material ruptures.
For calibration of damage models, true stress-strain curve has to be used. The most common stress-strain curve however is the engineering
stress-strain curve [2] [3].
2.4
Material model
Over the years different material models have been developed with one purpose; describe the material(s) properties mathematically. In this project we consider a model with parameters E, ρ, ν, ξ, tresca, and 0. Parameters
not considered are parameters that refers to temperature dependencies, such as thermal softening, thermal volymetric expansion, and adiabatic heat development. The material model used is a linear combination of Hooke’s law and Ludwik’s empirical formula [5]. It can be expressed as
σ() = ( E, if ≤ p σ(p) + C0( − p)C1, if > p (1) E, p, C0, C1 ∈ IR and constant(s).
where p is the value of where the elastic zone ends. p and C0 are
calculated based on empirical data, C1 is an arbitrary constant that is
chosen by trial and error, given the empirical data. Other parameters and variables as defined in section 2.1.
2.5
Damage model
In this project it was decided to work with a so called uncoupled damage model. Damage growth is driven by plastic straining of the material (damage grows faster in tension than in compression). There is no gradual material softening/degradation due to a increasing damage. However, as damage reaches 100% (D=1) the material is suddenly failing and the
element housing the failed integration point will be removed (eroded) from the model.
The Cockcroft- Latham model was used as damage model because of its high suitability for ductile solids [6]. The model can be expressed1 as
D = 1 Wc Z pf f 0 max(0, σ1)dpf f+ 1 ts Z t 0 (σ¯1 σs )αsdt (2)
where D is the damage level of the element, σ1 is maximum principal
stress, and in this case σ1 = ¯σ1, p
f f the effective plastic strain, d p f f the
effective plastic strain increment, αs an exponent controlling the time for
spall fracture, ts is time to develop spall fracture at threshold stress, and σs
is the spall strength. The element is eroded (deleted) when D = 1. The most important parameter to calibrate here is the normalisation constant Wc, since this determines how much the element may deform before it is
eroded.
2.6
Residual stresses
Residual stresses2 are stresses that remain in a solid material. Such stresses
are created in almost every manufacturing process; from casting to welding. They can occur without any external loads and have high impact on the strength of the material. Overlooked and neglected these stresses may cause major failures in constructions (for instance, the famous crack in the Liberty Bell as seen in Figure 5 is the result of neglected and escalated residual stresses). Accounted for, these stresses may be used to increase the durability of the material, which makes the subject highly relevant to study for the manufacturers.
1This notation might upset mathematicians since it looks like we integrate from zero to the same variable as the function depends on (exR0tf (t)dt), however it is used for convenience. The more correct way to express it is to change variables (so that Rt
0f (t)dt becomes R0tf (s)ds).
Figure 5: As a famous example of residual stresses’ potential harm, the crack in the Liberty Bell occurred as a result of residual stresses.
The stresses are self-equilibrating (hence the name residual). This basically means that local areas of tensile and compressive stresses sum to create a force and moment component that sums to zero, much like a conservational mechanical system3. This leads to that the tensile stresses in the central region of the material is cancelled/countered by equal (after summing) compressive stresses at the surfaces. If the stresses are large enough, the material will deform elastically in order to preserve any plastic
deformation’s effect on the residual stresses. During the lifetime of the material external loads such as bending, temperature changes,
transformation stresses, surface corrosion, and intergranular stresses all play a role in how the residual stress might change over time and might cause a failure after some time instead of directly after manufacturing, albeit some residual stresses are so large initially that the cause the material to fail almost immediately.
All force components must always sum to zero when no external loads are applied. Because of this, the stresses might not be very apparent and hence overlooked or neglected. In material design, these stresses may be used to enhance the performance of the material. In brittle materials for example
3A classic example of this is planets that orbit stars thanks to the conservation of energy and angular momentum; forces sum to zero.
(glass), a residual stress that compresses the surface at the expense of tensile stresses inside it will increase the performance since brittle materials usually can withstand high compressive forces and are sensitive to cracks. By hardening (compression hardens the material) the brittle material’s surface, the resistance to any cracks will increase at the surface. A practical example of this is a car’s door glass. The glass is tempered to such a high degree that it will withstand a high degree of external violence without any cracks to be seen. However, once a crack do form, the entire glass panel will shatter. In a car this i very desirable, since it basically eliminates the risk of fatal bleeding injuries as a result of cutting on the shattered glass, since all the shattered pieces are so small.
The residual stress components can be expressed on matrix form, S = σxx σxy σxz σyx σyy σyz σzx σzy σzz (3)
where each component is some function of space and time that symbolises the distribution in that corresponding direction. Since the stresses are equal and opposite, we have in (3) that ST = S, i.e S is symmetric
(σij = σji, where i, j = {x, y, z}).
2.7
Finite Element Method
As a foundation of the Finite Element Software (see section 2.8) lies the finite element method (FEM) [1]. Basically, this method uses the fact that any geometry can be approximated with a linear combination of a finite number of elements that are composed of so-called basis functions, as is shown in Figure 6.
Figure 6: The blue line is the exact solution to some diffusion-problem on [0,1] in 1D. The red line consists of two linear elements (three nodes) and approximates the solution. The green line consists of three elements, and the yellow line four. As can be seen, more elements is better when it comes to accuracy. Note that the elements have for the sake of simplicity been chosen to have the same interval lengths.
2.7.1 Derivation
A minor derivation of the FEM the software uses is necessary in order to understand the advantages and drawbacks of different element types. Given that one wants to approximate the function u(x) with the linear combination of some set of basis functions
u(x) =
n
X
i=1
ciφi, u ∈ [0, 1] (4)
where i is number of grid points located at xi and ci some constant defined
for node i. Transforming this into a local coordinate system so that the elements are defined for x = [0, 1] yields
ξ = x − xi xi+1− xi
, (5)
hence u is locally approximated by ξ, i.e for the 1D linear case
u(ξ) = c1+ c2ξ, ξ1,2 = [0, 1] (6)
which in turn leads to
c1 = u1c2 = −u1+ u2 (7)
that can be expressed on matrix form as
A = 1 0 −1 1 .
The approximate function with node values then becomes
u(ξ) = u1+ (−u1+ u2)ξ = u1(1 − ξ) + u2ξ = u1N1(ξ) + N2(ξ)ξ (9)
where N1, N2 are linear basis functions for 1D elements. This derivation
can be performed for higher order elements as well, in higher dimensions. For second order 1D one makes the quadratic ansatz
u(ξ) = c1 + c2ξ + c3ξ2, ξ1,2,3 = [0, 1] (10)
and a cubic ansatz for third order elements. For 2D and 3D the derivation becomes a bit more tricky, but the same method still applies (coordinate transformation, followed by an ansatz, which leads to some basis functions N1, N2, ... ,Ni. Some basis functions for the 2D second order is shown in
Figure 8, and basis functions for the 1D first, second, and third order is shown in Figure 7.
Figure 7: The first, second, and third order basis functions for 1D.
Figure 8: Two of the quadratic basis functions for rectangular elements. N1 and N2 are derived in section 2.7.1.
2.7.2 Finite difference scheme
In order to solve time-dependent problems, the equations need to be solved for each timestep. This is derived in [8]. IMPETUS Afea’s Solver uses a second order explicit scheme. For example, given the differential equation (harmonic oscillator) in [8]
∂2x
∂2t + ω
2x = 0 (11)
with ω ∈ IR, x0 = 1 and ∂x∂t0 = 0, the integration scheme can be written on
matrix form as xn+1= Axn (12) where xn+1 = xn+1 ∂xn+1 ∂t , xn= xn ∂xn ∂t , and A = 1 −12∆t2ω2 ∆t −∆tω2+1 2γ∆t 3ω4 1 − γ∆t2ω2 . The integration scheme (11) is stable if the eigenvalues of A λ1, λ2 satisfy
|λ1|, |λ2| ≤ 1. (13)
Solving (11) with (13) as constraints for ∆t, one gets that the critical timestep ∆tcr is ∆tcr = 1 ω r 2 γ (14)
If the time step would be larger than this value, the solution would start to oscillate more and more until it diverges completely. This may lead to long simulation times, if the grid is very fine (i.e a lot of elements requires smaller time step).
2.8
FEM software
The FEM software used in this project is IMPETUS Afea’s solver. The IMPETUS Afea Solver is a non-linear explicit finite element solver,
specialised to model material deformations under load. The the solver uses a J2 (von Mises) yield criterion and associated plastic flow rule [9] [10], and has slightly modified (optimised4) quadratic and cubic functions than what
is presented in Section 2.7.1. It is designed for a high level of robustness (double precision), capable of simulating advanced material distortion,
4As the exact source code is not available to the public, these will and can not be presented. The derivation is there to create a basic understanding of the project.
hydrodynamics, welding, explosives and collisions [9]. IMPETUS Afea Solver works with Green-Naghdi stress rate, where material rotations are defined through R in the relation (see also [8] [9])
F = RU (15)
where F is the deformation gradient, U is a symmetric stretch tensor (describes deformation), and R is a rotation tensor (describes rigid body rotation). For contact between objects, a penalty based contact algorithm and a Coulomb model of friction is used. A penalty based contact
algorithm allows for small geometric errors (penetrations between the bodies in contact). The applied contact pressure5 is a direct and linear function of the penetration error.
The solver also includes an advanced post-processing visualisation tool and relies entirely on Lagrangian equations of motion. The IMPETUS AFEA Solver is GPU accelerated and is reliant on NVIDIA’s CUDA architecture. Thanks to most parts of the problem being parallelisable, the benefit of the GPU’s many cores is huge and speeds up calculations a whole lot. The solver was developed to simulate (amongst other things) large mechanical deformations and suits this project very well [9].
2.8.1 Time stepping
As stated in section 2.7.2, the time stepping scheme is of second order accuracy. However, since it is an explicit scheme, this is only true if the time step size is less than a certain threshold. The time step size in an explicit model is dependent on the size of the element and the wave speed. Using too large a time step means the speed of sound is large enough to propagate through an entire element in a single time step, which will result in a unstable solution since this produces shockwaves.
The critical time step is limited by the maximum eigenfrequency of the system6. Given that γ = 2 in IMPETUS’s software, (14) becomes
∆tcr =
2 ωmax
. (16)
5The contact stiffness in this work was set to 5.0 · 1015 Pa/m. That is, a penetration error of 1 µm will generate a contact pressure of 5.0 · 1015Pa/m · 1.0 · 10−6 m = 5 · 109 Pa = 5 GPa.
6This in turn can be derived from the central difference scheme,
https://ftp.lstc. com/anonymous/outgoing/jday/olovsson2005_selective_mass_scaling.pdfin turn referencing ”The Finite Element Method: Linear Static and Dynamic Finite Element Analysis”.
However, an exact calculation of ωmax is expensive. The Impetus Afea
Solver instead estimates the critical time step ∆tc as
∆tcr = min(∆te, ∆tcont, ∆tblast, ∆tSP H) (17)
where ∆te is the critical element time step, ∆tcont is critical contact time
step, ∆tblast is the critical particle blast time step size, ∆tSP H is the critical
SPH (Smoothed-particle hydrodynamics) time step size.
For a linear model ∆tcr is approximately equal to the Courant equation
C = ∆t n X i=1 uxi ∆xi ≤ Cmax, (18)
where uxi is the magnitude of velocity in dimension xi, ∆t is the time step
size, and ∆xi is the step size in dimension xi. This will yield a rough
estimate for more complicated elements where it is unclear how the element length should be measured. If one lets c be the dilatational wave speed through the material, Lc be the smallest element size, then one gets
that ∆t ≈ Lc c (19) c = s E ρ (20)
where E is Young’s modulus and ρ is the mass density. This implies that a material with a higher density and equal modulus E will have a lower wave speed, leading to a larger stability limit [11] [12].
∆tstable = Lc
r ρ
E (21)
2.8.2 Mass scaling
In order to avoid simulation times that are excruciatingly long, the
software exploits the fact that the critical time step is related to the speed of sound (c) in the material [9] [8]. Equation (21) basically states that if one adds artificial mass to the system, increasing the density ρ whilst leaving E unchanged, ∆tstable will increase. This effectively reduces c and
hence increases the critical timestep size. The approximation leads to shorter simulation times, however this may come at the cost of increased energy in the system that derives from the added mass. If the added mass is small enough this is negligible, however if the energies from the added mass becomes large some results might not be very accurate.
3
Method
Overall, the method presented should be seen as conceptual. This means that the precise choice of material and material parameters should matter less and the choice of method and proposed experiments matter more, as the method should be applicable to some arbitrary, similar process. First, the process was observed. The process was then implemented in IMPETUS Afea’s IDE, along with all measurements. After the simulation was complete, the results were objectively compared with material samples of the chosen material. The edge of the slitted material (both from the simulation and the samples) was also compared to well-established characterisation images to further measure the accuracy of the model.
3.1
Material
High strength steel for the sheet was SSAB’s Docol 930S [4]. The material was observed in the slitting machine, and test pieces were obtained for microscopic study.
3.2
Measurements and parameters
3.2.1 Stress-strain
Stress-strain data was achieved by a classic stress-strain test of the
material. However the cross-section area was not measured and hence the achieved curves could not be transformed into true stress-strain curves. Instead, the material model (1) was fitted to the part of the plastic zone before any real difference had occurred between the engineering
stress-strain and the true stress-strain. That is, the curve was optimised to fit experimental data until onset of diffuse necking.
3.2.2 Slitting machine
The slitting machine considered is depicted in Figure 9, where it is drawn schematically with each part scaled down by an equal amount. In Figure
Figure 9: The slitting machine at SSAB shown schematically from the side.
Figure 10: The slitting process at SSAB shown schematically from the side. 3.2.3 Residual stress
Residual stresses were applied arbitrary to a model of a large sheet to see what effect a certain distribution had on the material. The functions applied and the size of the residual stresses where chosen by studying [2] and [7].
3.3
Material model adaptation
Since this is a conceptual study parameter calibration is not of the essence. However, some calibrations where necessary to perform in order to be able to be verify the accuracy of the model.
3.3.1 Parameter calibration
By obtaining stress-strain curves for the steel, parameters C0, C1, and p
could be calculated7 in MATLAB [13] by first transforming the non-linear
equation (1) into a linear regression problem on matrix form
Ax = b (22) where A = C1 0 1 C1 1 1 C1 2 1 .. . ... C1 N −1 1 , x = C0 σ(p) , and b = b0 b1 b2 .. . bN −1 .
Solving for x (minimisation problem for each guess of C1), one arrives at
an optimum rather quickly, given that the size of b is not too large8.
3.3.2 Choice of damage model
The Cockcroft- Latham (CL) damage model was chosen for this study. The CL damage model is not very precise, however it is generally the go-to option when there is limited amounts of material data. CL damage model is appropriate for ductile materials and suits the needs for this project well [9].
3.4
Software setup
IMPETUS uses a command structure very similar to LS Dyna [14] to setup the problem. The source code for this may be found in the Appendix. The parameters used are stated in Table 1 unless otherwise specified. The material model used is the function resulting from the calibration. See Section 2.4 and Figure 13for these.
7See appendix for source code.
8In practise this typically means that N < 104, although this is more a rule of thumb and may vary a lot.
Parameter Value Unit E 210 · 109 P a gap 5 · 10−5 m height 5 · 10−4 m overlap 3 · 10−3 m width 2 · 10−2 m R 0.1 m r 0.2 · 10−3 m v 3.33 m/s ρ 7800 kg/m3 ν 0.3 N umeric Wc 1 · 109 P a ∆tcr 4 · 10−8 s
Table 1: Gap corresponds to the slit space between cutters (x-direction), height the sheet thickness (y-direction), overlap the maximum overlap of the cutters in 3D (y-direction), width is cutters’ contact side width, R is cutter radius, r is radius of the edge that cuts the sheet, v is the velocity at which the sheet enters the cutters. Other variables defined as before. All parameters based on real-life measurements.
3.4.1 Element type
In higher dimensions, such as 2D or 3D, the elements may be based on rectangles, triangles, tetrahedrons, hexahedrons, and (in more rare cases) pentahedrons (basically prisms in their appearance). In Figure 8, examples of second order (quadratic) rectangular elements are shown.
3.4.2 Element order
Most common by far is the first order, which basically consists of piece-wise linear polynomial basis, i.e the functions described in section
2.7.1 of the first order. These are less computationally expensive than second or third order elements, but does not approximate non-linear shapes or plastic deformation in general as well as the higher order types do. 3.4.3 Element density
The density of the elements is ultimately what limits the ability to see details in the simulated material. The more elements per volume, the higher the detail (see Figure 6for an example of this). However, the number of elements used is limited by the time it takes to simulate them.
By putting more elements where the fracture will occur and not as many where no fractures will occur, the resources are used more efficiently. This distribution can be seen for the 2D-case in 11and in Figure 12 for the 3D-case. For the 2D model at most 28590 elements were used. For the 3D model the element count was 907698 at most.
Figure 11: 2D model of the process. Note the denser distribution of ele-ments around the centre.
Figure 12: 3D model of the process. Note the denser distribution of ele-ments around the centre.
3.4.4 Time stepping and mass scaling
The time step size was set to be no smaller than 40 ns. If the critical timestep was smaller than this, mass-scaling was to be used. This led to simulation times of 2-7 days for the 3D model and a few hours for the 2D model. If no mass-scaling had been used, the 3D-model with the most elements would have taken more than one month to complete, maybe more.
3.4.5 Boundary conditions
As boundary condition for the 3D model, the edges of the sheet (the edges with highest |x|-coordinate values) were chosen to only be able to move along the z-axis. The cutters had boundary conditions that prevented any translational motion and were set to rotate around their centre with a constant angular velocity |w|, such that the relative velocity with respect to the sheet was 0 at the point of contact (|w × r| = 0), i.e the cutters rolled on(to) the sheet.
3.4.6 Initial conditions
Due to the effects of mass-scaling, it is better for accuracy and simulation time to initialise an object that has very high density in the mid-section
(such as the simulated sheet) at a certain speed, instead of applying an acceleration. The sheet was set to an initial velocity v = vz = 3.3 m/s.
4
Results
4.1
Stress-Strain curve
The material model was fitted to the stress-strain curve for Docol 930S successfully. The resulting function can be seen in Figure 13.
Figure 13: Fitted curve with least squared error to the data.
4.2
Scanning Electron Microscope
In these SEM images, the direction ”up” is where the shearing tool first impacted the steel, and ”bottom” where the tool last was in contact with the steel. This corresponds to the y-axis in the model if one considers the part of the sheet that is on the negative side of the x-axis. The different zones depicted here corresponds to what is shown in Figure 3, where the cutted zone is referred to as the bright zone. Note that some SEM images are rotated for convenience. Figure 14 depicts the slitted edge from the side.
Figure 14: Overview of the slitted edge. The bright zone is found in the up-per part, with a distinct line separating it from the fracture zone.
4.2.1 Rollover zone
The rollover zone is characterised by its more or less rounded/compressed edge at the top and a vertically straight, smeared-out zone beneath this. An effect of the initial compression of the material. See the upper part of Figure 18for a clear view of this. The rollover zone measured 173 µm in width and 57 µm in height.
4.2.2 Shear zone (bright zone)
The primary shear zone9 is the result of a compression in such a way that
the bright zone becomes smooth and very even. The bright zone is seen zoomed in in Figure 15, and Figure 14 depicts the slitted edge from the side, where the bright zone is found in the upper part.
Figure 15: A zoomed-in image of the bright zone. Note the smoothness. 4.2.3 Secondary shear zone (fracture zone)
Secondary shear zone is located beneath the primary shear zone, see Figure
14. Here, the material has suffered (almost) exclusively tensile forces and has hence been ripped apart. The zone is characterised by its moon-like landscape with mountain-like pikes and high variation in height compared to the primary shear zone. A zoomed-in image of the fracture zone is seen in Figure 16
Figure 16: Fracture zone. Note how the material has been ripped apart and created the very non-smooth surface.
4.2.4 Burr zone
At the very bottom the burr zone is located. This is the result of tensile forces that has dragged the material out to such an extent that it has plastically deformed and formed a sharp edge. The burr zone may be seen in Figure 18 at the very bottom left, as well as in Figure17. The burr zone measured 40 µm in height.
Figure 17: Overview of the slitted edge a bit from the side. The burr zone is the entire, outermost part of the edge.
4.2.5 Plastic deformation zone
The plastic deformation zone is defined as the zone that contains parts that has been plastically deformed as a result of the slitting process. It is typically a semi-circle seen from the positive direction of the z-axis that extends some distance into the negative x-direction of the sheet. In Figure
18 the left sample has been slitted with a gap of 5 · 10−5 m between the two cutters, whereas the right sample depicts the exact same material, cut with more gap ( 5 · 10−4 m) between the cutters.
Figure 18: Profile of the slitted material, achieved with EBSD. Note the high impact of the slitting gap; the right one has suffered a much higher deformation as a result of this.
Figure 19: Profile of the slitted material, achieved with EBSD with plastic deformation as overlay. The plastic deformation colour depicts the degree of deformation (i.e damage), where red is highest and blue lowest. Note that the overlay is normalised to its picture; that is, it not necessarily true that the red zone to the left has the same plastic deformation degree as the the red zone on the right.
Figure 20: Profile of the slitted material, achieved with EBSD. The plastic deformation is depicted and lets one see the direction of the deformation. Figure 19depicts the samples from Figure 18, but with an overlay that shows the degree of plastic deformation. Both Figure 18 and 19were achieved with EBSD (Electron BackScatter Diffraction) [16].
4.3
Two-dimensional model
The 2D model was performed in another project, see [17] for a complete project description. Residual stresses were not considered in this model. The result for the 2D model may be seen in Figure 21. The 2D model acted as an uncalibrated simplification in order to see if the simulation was possible or not.
Figure 21: The damage zone is shown graphically. Note that the fracture now follows the path of maximum damage thanks to the altered mesh. To see how the mesh was altered, please see [17].
4.3.1 Energy balance (2D)
Energy levels for the added mass in 2D are very small, as may be seen in Figure 22.
Figure 22: Here some different types of energies can be seen as a function of time in the simulation. For accuracy, it is important (when mass scal-ing is used) that the kinetic energy created by added mass (Ekinetic−added)
does not increase much. It is especially important if one wants to ne-glect the added mass completely that Ekinetic−added << Einternal−plastic and
Ekinetic−added << Eexternal−mechanical.
4.4
Three-dimensional model
In order to improve the model the first 3D models were simulated with values from the slitting machine (Wc chosen arbitrarily) and then
compared to different simulations with slightly tweaked parameters. The first simulation used as reference point may be found in Figure 23. Movies of the process may be found here10 and here11.
10The zoomed-out version: https://drive.google.com/open?id=
129UPT9rz7DUrF6Xx_s7BhNkbAiw-Hxgp
11The zoomed-in version: https://drive.google.com/open?id=
Figure 23: Simulation with Wc = 1 GPa, r = 0.2 mm, overlap = 3 mm,
and gap = 0.1 mm.
4.4.1 Energy balance (3D)
Energy levels for the added mass in 3D are large, as may be seen in Figure
24. This means that the added mass may not be negligible in the same way as for 2D, and measurements made in the model needs to take into account the added mass.
Figure 24: Here some different types of energies can be seen as a function of time in the simulation. For accuracy, it is important (when mass scal-ing is used) that the kinetic energy created by added mass (Ekinetic−added)
does not increase much. It is especially important if one wants to ne-glect the added mass completely that Ekinetic−added << Einternal−plastic and
Ekinetic−added << Eexternal−mechanical. This is not true here, since the energy
from the added mass is very large. 4.4.2 Different sharpness
Simulations with more sharp and rounded edges of the tools are depicted in Figures 25 and 26, respectively. The more rounded tools appear to deform the steel more.
Figure 25: The edges of the tools have a sharpness radius r = 0.1 mm.
Figure 26: The edges of the tools have a sharpness radius r = 0.35 mm. 4.4.3 Different gap
A simulation with a smaller (infinitesimal) gap between the cutters is found in Figure 27. Worth to note is that the edge appears to be more damaged than when a small gap is used as space.
Figure 27: The edges of the tools have a sharpness radius r = 0.2 mm. The gap between the two cutter tools is infinitesimal.
4.4.4 Different damage tolerance
Figure 28depicts the case when the damage normalisation parameter Wc
(see section 2.4 for definition) was set to half the value of what was
initially simulated (Wc= 0.5 · 109 Pa). Figure 29when the value was set to
double the initial value (Wc= 2.0 · 109 Pa). As can be seen, the elements
are allowed to deform to a much higher extent with higher values of Wc.
Figure 28: Wc = 0.5 GPa. The edges of the tools have a sharpness radius
Figure 29: Wc = 2.0 GPa. The edges of the tools have a sharpness radius
r = 0.1 mm.
5
Discussion
Achieved results seem promising for future development. The 3D model appears to simulate the slitting process very well, and might be able to accurately predict new optimal parameters with the right calibration.
5.1
Software
IMPETUS Afea has proven itself to be extremely accurate and fast to work with. The software’s robustness opens up for future, larger simulations on clusters. The mass-scaling’s effect on the 3D model are a bit unclear, the only noticeable difference appears to be an increased pressure on the cutters, although this remains to be investigated. The time-step is quite hard to decrease due to the increase in simulation time, and by comparing the SEM and 2D results with the 3D model one may assume that the model is accurate despite the sinister energy plot seen in Figure 24. The largest problem by far with this study has been the fact that there are almost no research or material models that handle fracture very well. Typically these (FEM) studies instead simulates the physics up until fracture point. As follows from Section 2.7.1, one may come to realise that the higher order elements simulate plastic deformation much better than the linear elements do [8]. However, since we are more interested in how the fracture propagates and not so much in how the plastic deformation
depicts itself, the trade-off was made to increase the element density at the cost of the elements’ polynomial order.
5.2
Material model
As of now, the only material model considered has been the one specified in Section 2.4. A development of a new material model12 (along with
damage model) would probably enhance the model further, however the current material model must be seen as a very well-working compromise -especially since it may be calibrated further. The current material and damage models are, after all, based on excessive experimental results and may be expected to be very accurate.
5.3
Stress-Strain curve
Calibration of the material model may be seen in Figure 13. The fit was made due to the lack of true stress-strain data (which is hard to come by) and is an approximation of the true curve. The approximation relies on the fact that the true stress-strain relationship continues to grow for ductile materials as seen in Figure 4.
5.4
Scanning Electron Microscope
SEM images were beyond expectations in accuracy and showed very high-detailed images of the test pieces. This opened up for a more precise and objective verification process of the simulated results, especially due to the EBSD images seen in Figures 18, 19, and 20. The EBSD images shows that the plastic deformation and damage simulated in the model resembles what actually occurs in the real-life process. Figure 14 shows that the rollover zone and bright zone together fluctuates around 130 to 200 µm in depth, whereas the rollover zone itself measures around 50-60 µm. This verifies the results from the profile image in Figure 18, where the rollover zone was measured to be 57 µm. The variation in bright zone depth may be expected due to the fact that the steel in real life is a bit heterogeneous, as well as the fact that the cutters may vary a bit as well.
There are major differences between the bright zone and the fracture zone, as can be seen in Figures 15and 16. The difference is due to the fact that the material (on a microscopic case) has been compressed in the former
12A new material model is not an easy task to develop, and might even exceed a PhD in span.
case, and ripped apart in the latter. It is a bit unclear how much this affects the propagation of micro-cracks in the material, but the hypothesis is that the fracture zone is weaker than the bright zone due to the fact that force applied is more likely to be focused into one of the already existing cavities of the fracture zone.
5.5
Two-dimensional model
The two-dimensional model was developed in the project specified here [17]. During this time there was no material data available to calibrate the model and the constants C0, C1 and p were chosen arbitrarily. As a result
of this, the produced model results are not to be considered trustworthy, but rather a first, simple, proof of concept that the simulations indeed can be made. The mesh in Figure 21is a result of first simulating the process with straight mesh, and then, the timestep before any elements are eroded, the most damaged elements were chosen and their coordinates extracted into another, new simulation. This resulted in a very accurate prediction of where the fracture will propagate, and is an effective way to reduce the mesh bias when the element density is high. See [17] for a complete description of this process. The 2D model became the basis for the 3D model.
5.6
Three-dimensional model
Compared to the 2D version, the 3D model is far more complex since it includes depth dimension along with rotation of cutters and other parameters that correspond to those used in the real process. Albeit the calibration made was kept to a minimum, the results does indeed seem very promising for further studies. The influence of residual stresses seems to be negligible compared to the models that doesn’t take residual stresses into account. Also, since residual stresses typically appear in the
magnitude of 200 MPa at most (which is on the surface) [2] [7], the sheet would have to be larger in surface with equal thickness for stresses to have any noticable effect. All tests performed in Section 4.4 were simulated both with and without residual stresses with no noticeable difference.
Hypothetically, the model would have to include a larger portion of the sheet for these effects to appear. Simulations of such large models were skipped mainly due to lack of memory on GPU, but also because the results would be irrelevant due to the uncalibrated Wc.
5.7
Future improvements
The two most necessary experiments needed for an increase in accuracy are a true stress strain curve for the desired material, along with a calibration of Wc. A true stress strain curve would make one able to skip the data
fitting step in Section 3.3, since the raw data then can be imported directly into IMPETUS’s IDE. For better model verification the pressure may be measured in real life on the cutters and then compared to those in the model. Slitting is a continuous crack propagation process, while
Cockcroft-Latham is a criterion for ductile crack initiation. Hence, it might not be the ideal model for this process. However, shifting to something more accurate is not easy without access to results from a material
characterisation program. A future model might also consider the wear on the cutters (currently they are considered to be rigid bodies) and be able to predict when they should be replaced for better performance with regard to cost. A better calibrated model might also be able to predict optimal parameters (cutter overlap, cutter gap etc.) for damage
minimisation. CL damage model may cause mesh bias and should therefore be analysed in future development. An alternative might be to use an energy based crack propagation criteria (energy per unit crack). This would likely reduce mesh dependency greatly while also ensure convergence (see parameter GI in [9]).
A fully calibrated model should be able to estimate the spread of micro-cracks in the material by looking at the damage zone and the
plastically deformed zone in the model, as these have proven themselves to be accurate (compare, for instance, Figure 28with Figures 18and 19). Given that the micro-cracks may be estimated by the model, they could in theory be minimised based on the model and the model could then act as a basis for real-life experiments. The rapid development of GPU technology gives hope of an even more fine mesh, even on desktops.
5.8
Conclusions
Overall the model appears to successfully simulate the real-world process. Further calibrations may enhance the accuracy and perhaps make the simulation able to predict a new, better optimal set of parameters for the process. The most accurate parameters, based on the results, appear to be with Wc= 0.5 GPa and r = 0.1 mm. This is objectively achieved by
comparing the EBSD results in Figure 18 with the model results in Figure
28. It appears to be (based on the results) that a less sharp cutting tool deforms the steel more (Figure 26versus 25), too small a gap damages the
steel more (Figure 27versus 23), and it also seems like Wc should lie in the
neighbourhood of 0.5 GPa in size. The results should be verified to ensure this, although based on the EBSD images (Figures 18, 19, and 20) the rollover zone does indeed worsen with a more round tool and much larger gap.
One may say, based on these results, that the simulation as of now needs some more calibration before it can be seen as a complete success, but a success nevertheless. Due to the uncertainty in the simulation and the need for more essential parameters to be calibrated first, the effect of residual stress cannot be trusted as the model stands now.
The 2D model may be seen as a successful simplification of the 3D model, as their results deviate very little from one another. This basically means that for better crack propagation and prediction, a finer mesh may be used in the 2D model and used for comparison with the 3D model.
Furthermore, the method the model was developed with may be applied to other industrial processes. The applications are uncountable, and albeit improvements may not be guaranteed with this method, they are certainly far more likely to occur.
References
[1] Mathematics Handbook
By Lennart R˚ade and Bertil Westergren. Published 2004 by Studentlitteratur AB. [2] Mechanical Behaviour of MATERIALS
Second edition
Marc Meyers and Krishan Chawla Cambridge University Press
University Printing House, Cambridge CB2 8BS, UK
[3] Teknisk h˚allfasthetsl¨ara by Tore Dahlberg. Published 2001 by Studentlitteratur AB.
[4] SSAB’s dual-phase high-performance steel. Fetched 12/5 2019.
https://www.ssab.se/produkter/varumarken/docol
[5] Ludwik’s empirical relation. Fetched 12/2 2019
[6] Cockcroft-Latham damage model.
https://www.dynalook.com/conferences/
12th-international-ls-dyna-conference/keynotep-a.pdf
[7] Practical Residual Stress Measurement Methods By Gary S. Schajer
Published 2013 by WILEY
[8] On the Arbitrary Lagrangian- Eulerian Finite Element Method By Lars Olovsson
Division of Solid Mechanics
Departement of Mechanical Engineering
Link¨opings universitet, SE-58183 Link¨oping, Sweden [9] IMPETUS Afea’s homepage.
https://www.impetus-afea.com/
IMPETUS’s manual, Cockcroft Latham:
https:
//www.impetus-afea.com/support/manual/?command=PROP_DAMAGE_CL
IMPETUS’s manual, Material properties:
https:
//www.impetus-afea.com/support/manual/?command=MAT_METAL
[10] Von Mises yield criterion.
https://www.engineersedge.com/material_science/von_mises.htm
[11] On solid mechanics.
https://classes.engineering.wustl.edu/2009/spring/mase5513/ abaqus/docs/v6.6/books/gsa/default.htm?startat=ch09s03.html
[12] Physics Handbook
By Carl Nordling and Jonny ¨Osterman. Published 2004 by Studentlitteratur AB. [13] MATLAB’s homepage.
https://se.mathworks.com/products/matlab.html
[14] LS Dyna’s homepage. Fetched 13/4 2019
http://www.lstc.com/products/ls-dyna
[15] Example images of how sheared steel usually look. Fetched on 9/1 2019.
https://www.thefabricator.com/article/stamping/ blanking-questions-have-you-on-the-edger
[16] The EBSD technique. Fetched 12/5 2019.
https:
//en.wikipedia.org/wiki/Electron_backscatter_diffraction
[17] Project report of the previous study. Finished 19/2 2019
https:
APPENDIX
3D model script. 1 ∗UNIT SYSTEM 2 S I 3 ∗PARAMETER 4 L = 0 . 1 5 D = 0 . 0 8 6 h = 0 . 5 e−3 7 o f f = 0 . 0 0 9 5 8 gap = 0.2∗%h 9 R = 0 . 1 10 width = 0 . 0 2 11 o v e r l a p = 0 . 3 e−3 12 v = 3 . 3 3 13 w = %v/%R 14 t e n d = 0.5∗(%D+%o f f )/%v 15 ∗TIME 16 [% t e n d ] , 0 , 4 . 0 e−8 17 # 18 # −−− MESH −−− 19 # 20 ∗COMPONENT BOX 21 ” S h e e t − c e n t e r f i n e ” 22 3 , 3 , 4 0 , 3 0 , 750 23 [−%gap −1.5∗%h ] , [−%h / 2 ] , [−%D/4−% o f f ] , [% gap +1.5∗%h ] , [%h / 2 ] , [−% o f f ] 24 ∗COMPONENT BOX 25 ” S h e e t − c e n t e r c o a r s e ” 26 4 , 4 , 1 , 1 , 1827 [−%gap −1.5∗%h ] , [−%h / 2 ] , [−%D−%o f f ] , [% gap +1.5∗%h ] , [%h / 2 ] , [−% D/4−% o f f ] 28 ∗COMPONENT BOX 29 ” S h e e t − l e f t ” 30 5 , 4 , 1 0 , 1 , 24 31 [−%L / 2 ] , [−%h / 2 ] , [−%D−%o f f ] , [−%gap −1.5∗%h ] , [%h / 2 ] , [−% o f f ] 32 ∗COMPONENT BOX 33 ” S h e e t − r i g h t ” 34 6 , 4 , 1 0 , 1 , 24 35 [% gap +1.5∗%h ] , [−%h / 2 ] , [−%D−%o f f ] , [%L / 2 ] , [%h / 2 ] , [−% o f f ] 36 ∗MERGE DUPLICATED NODES
37 P , 4 , P , 4 , [%h / 1 0 ] 38 ∗INCLUDE 39 ” C u t t e r 1 ” 40 . . / c u t t e r r 0 . 2 . k 41 1 , 1 , 1 42 0 , 0 , 0 , [% gap / 2 ] , [%R−%o v e r l a p / 2 ] , 0 43 ∗INCLUDE 44 ” C u t t e r 2 ”
45 . . / c u t t e r r 0 . 2 . k 46 1 , 1 , 1 , 1 0 0 0 0 , 1 0 0 0 0 , 1 47 0 , 0 , 0 , [−%gap / 2 ] , [−%R+%o v e r l a p / 2 ] , 0 48 −1, 0 , 0 49 ∗CHANGE P−ORDER 50 PS , 1 2 4 , 2 51 ∗SMOOTH MESH 52 PS , 1 2 4 , 6 0 . 0 53 ∗SET PART 54 124 55 1 , 2 , 4 56 ∗MERGE 57 P , 3 , P , 4 , [%h / 1 0 0 ]
58 ∗REDISTRIBUTE MESH CARTESIAN 59 1
60 P , 3 , X, 0 , 1 61 0 , 0 , [−% o f f ]
62 ∗TRANSFORM MESH CARTESIAN 63 1 , P , 3 , 0 , 4444 64 ∗FUNCTION 65 4444 66 0 . 4 ∗ y ∗ ( 1 − abs ( x ) /(%gap +1.5∗%h ) ) 67 ∗GEOMETRY BOX 68 123 69 [−%L / 4 ] , [−%h / 2 ] , [−%D−%o f f ] , [%L / 4 ] , [%h / 2 ] , [−% o f f ] 70 # 71 # −−− MATERIALS −−− 72 # 73 ∗MAT METAL 74 1 , 7 8 0 0 . 0 , 2 1 0 . 0 e9 , 0 . 3 , 1 75 1 76 ∗FUNCTION 77 1 78 6 4 7 . 8 1 e6 + 1 4 3 8 . 5 e6 ∗ e p s p ˆ 0 . 3 79 ∗PROP DAMAGE CL 80 1 , 1 81 1 . 0 e9 82 ∗MAT RIGID 83 ” C u t t e r ” 84 2 , 7 8 0 0 . 0 85 # 86 # −−− PARTS −−− 87 # 88 ∗PART 89 ” C u t t e r 1 ” 90 1 , 2 91 ” C u t t e r 2 ” 92 2 , 2 93 ” S h e e t − c e n t e r ”
94 3 , 1
95 ” S h e e t − e d g e s ” 96 4 , 1
97 #
98 # −−− BOUNDARY CONDITION and CONTACT −−− 99 # 100 ∗CONTACT 101 ” g e n e r a l c o n t a c t ” 102 1 , 1 103 PS , 3 4 , PS , 1 2 , 0 . 1 , 5 . 0 e 1 5 104 0 , 0 , 333 105 8 3 , 8 3 , 0 , 0 , 1 106 ∗FUNCTION 107 83 108 p∗ vtang 109 ∗GEOMETRY BOX 110 333
111 [−%width ] , [−%R ] , [−%D−%o f f ] , [% width ] , [%R ] , 0 112 ∗SET PART 113 34 114 3 , 4 115 ∗SET PART 116 12 117 1 , 2 118 ∗BC MOTION 119 ” C u t t e r 1 ” 120 1 121 P , 1 , XYZ, YZ 122 V, RX, 1001 123 ∗BC MOTION 124 ” C u t t e r 2 ” 125 2 126 P , 2 , XYZ, YZ 127 V, RX, 1 0 0 1 , −1.0 128 ∗FUNCTION 129 1001 130 −%w 131 ∗INITIAL VELOCITY 132 PS , 3 4 , 0 , 0 , [%v ] 133 ∗BC MOTION 134 ” S h e e t l e f t ” 135 3 136 G, 1 , Y 137 V, Z , 1002 138 ∗BC MOTION 139 ” S h e e t l e f t ” 140 4 141 G, 2 , Y 142 V, Z , 1002
143 ∗FUNCTION 144 1002 145 %v
146 ∗GEOMETRY SEED COORDINATE 147 ” S h e e t l e f t ”
148 1
149 [−%L / 2 ] , 0 , [−%D/ 2 ]
150 ∗GEOMETRY SEED COORDINATE 151 ” S h e e t r i g h t ”
152 2
153 [%L / 2 ] , 0 , [−%D/ 2 ] 154 ∗END
MATLAB code for curvefit:
1 c l o s e a l l 2 c l e a r a l l 3 d a t a = x l s r e a d ( ’ d o c o l 9 3 0 s . x l s x ’ ) ; 4 d a t a 2 = x l s r e a d ( ’ docol1400m . x l s x ’ ) ; 5 e p s i l o n 2 = d a t a 2 ( 1 : end , 1 ) ; 6 sigma2 = d a t a 2 ( 1 : end , 2 ) ∗1 e +6; 7 8 9 e p s i l o n = d a t a ( 1 : end , 1 ) ; 10 sigma = d a t a ( 1 : end , 2 ) ∗1 e +6; 11 y = sigma−sigma ( 1 ) ; %s e t x ( 0 ) =0 12 13 p = [ e p s i l o n . ˆ ( 0 . 3 ) o n e s ( s i z e ( y ) ) ] \ sigma ; % E s t i m a t e P a r a m e t e r s 14 % i f y = k∗x . ˆ 0 . 5 +C, t h e n p ( 1 )= k and p ( 2 )=C 15 z = e p s i l o n . ˆ 0 . 3 ∗ p ( 1 )+p ( 2 ) ; 16 o r i g = e p s i l o n . ˆ 0 . 3 ∗ 0 . 1 ∗ 1 0 ˆ 9 + 1 . 5 ∗ 1 0 ˆ 9 ; 17 18 d i s p ( ’ k = ’ ) 19 d i s p ( p ( 1 ) ) 20 21 d i s p ( ’C = ’ ) 22 d i s p ( p ( 2 ) ) 23 24 f i g u r e 25 h o l d on 26 p l o t ( e p s i l o n , sigma ) 27 % p l o t ( e p s i l o n 2 , sigma2 ) 28 p l o t ( e p s i l o n , z ) 29 % p l o t ( e p s i l o n , o r i g )
30 % l e g e n d ( ’ Docol 930 S ’ , ’ Docol 1400M’ , ’ F i t t e d c u r v e y=k∗ e p s ˆ{0.3}+ C ’ , . . .
31 % ’ O r i g i n a l m a t e r i a l model ’ )
32 l e g e n d ( ’ Docol 930 S ’ , ’ F i t t e d c u r v e y=k∗ e p s ˆ{0.3}+C ’ ) 33 x l a b e l ( ’ e p s i l o n ’ )
34 y l a b e l ( ’ sigma ’ )
MATLAB code for mesh alteration:
1 c l e a r a l l 2 format long 3 d a t a=l o a d ( ’ d a m a g e l i s t 4 1 . t x t ’ ) ; 4 x = d a t a ( : , 1 ) ; 5 y = d a t a ( : , 2 ) ; 6 z = d a t a ( : , 3 ) ; 7 dmg = d a t a ( : , 4 ) ; 8 9 %k = l e n g t h ( y ) ; 10 11 % f o r k = 1 : l e n g t h ( y ) 12 % i f z ( k ) ˜=0 && y˜=0 13 % y ( k ) =100; 14 % x ( k ) =100; 15 % end 16 % 17 % end 18 19 i =1; 20 j =1;
21 % Take away e l e m e n t s where z !=0 22 while i <l e n g t h ( y ) 23 i f ( z ( i ) ˜=0) 24 X( j )=x ( i ) ; 25 Y( j )=y ( i ) ; 26 DMG( j )=dmg( i ) ; 27 j=j +1; 28 end 29 i=i +1; 30 31 end 32 33 h=max(Y) ; 34 TOL = 2 . e −4; 35 Ycord = u n i q u e (Y) ; 36 maxD = z e r o s ( 1 , l e n g t h ( Ycord ) ) ; 37 Xcord = z e r o s ( 1 , l e n g t h ( Ycord ) ) ; 38 i =1; 39 % C c o n t a i n s t h e u n i q u e v a l u e s o f Y. For e a c h v a l u e o f Ycord t h a t e q u a l s
40 % some Y, c h e c k a l l t h o s e damages and s a v e t h e one w i t h h i g h e s t dmg .
41 % A l s o s a v e t h e x−c o o r d f o r t h i s v a l u e . 42 f o r k =1: l e n g t h ( Ycord )
43 f o r i =1: l e n g t h (Y)
45 i f (maxD( k ) == 0 ) 46 maxD( k ) = DMG( i ) ; 47 end 48 %e l s e 49 i f k==1 50 i f (maxD( k )<DMG( i ) ) 51 maxD( k )=DMG( i ) ; 52 Xcord ( k )=X( i ) ; 53 end 54 end 55 i f ( k>1)
56 i f (maxD( k )<DMG( i ) && ( abs (X( i )−Xcord ( k−1) )<TOL) ) 57 maxD( k )=DMG( i ) ; 58 Xcord ( k )=X( i ) ; 59 end 60 end 61 %end 62 end 63 end 64 end 65 66 % For e a c h h e i g h t v a l u e , p l o t t h e x−c o o r d i n a t e w i t h h i g h e s t dmg 67 f i g u r e 68 p l o t ( Xcord , Ycord ) 69 a x i s ( [ − 1 . e−3 1 . e−3 0 h ] ) 70 x l a b e l ( ’ H o r i s o n t a l p o s i t i o n [m] ’ ) 71 y l a b e l ( ’ V e r t i c a l p o s i t i o n [m] ’ ) 72 t i t l e ( ’ P r e d i c t e d c r a c k r e s u l t , 33743 c u b i c e l e m e n t s ’ ) 73 74 X2 = Xcord ’ ; 75 Y2 = Ycord ’ ; 76 SavedCord = [ Y2 X2 ] ; 77 %d l m w r i t e ( ’ SavedCord . t x t ’ , SavedCord , ’ d e l i m i t e r ’ , ’ \ t ’ )
![Figure 6: The blue line is the exact solution to some diffusion-problem on [0,1] in 1D](https://thumb-eu.123doks.com/thumbv2/5dokorg/5387162.137396/15.892.304.590.196.378/figure-blue-line-exact-solution-diffusion-problem-d.webp)
