The 9/11 Forum

All you need is Love? Same elastic collision problem, different results.


I think I mentioned it on physorg many moons ago, before Gregory Urich appeared on the scene, but at the end of Love's treatment of problem 281, 'Rod fixed at one end and struck longitudinally at the other', he states,

It may be shown also that the greatest compression of the rod occurs at the fixed end, and that, if m < 5, its value is 2(1 + exp(-2/m)V/a, but, if m > 5, its value is approximately equal to (1 + sqrt(m))V/a

V is the impact velocity, taken to be 8.5 m/s in the BZ scenario.

I used Greg's spreadsheet at http://www.cool-places.0catch.com/911/calcMassAndPeWtc1_v2.5.xls
to calculate m, which in the formulas above is the mass ratio between the striking mass (considered rigid) and the rod (considered elastic, with bar velocity a)


total mass 98 - 111
32,815.80 kg


from Greg's posts here:
http://the911forum.freeforums.org/numerical-solution-of-elastic-rods-colliding-t75.html#p1002

"Mass of columns B6-97 = 55,700 metric tons."

which is 55,700,000 kg

so taking m to be the ratio of all the stories above floor 97 to the mass of the columns in floor 97 and below, we get

m = 32,815,800 / 55,700,000 = .589

being less than 5, we use the first formula mentioned above, viz.,


2(1 + exp(-2/m))V/a

substituting, we get

2(1 + exp (-2/ .589 )) V/a
=2(1 + .0335) V / a

= 2.067 (8.5 / 5200)

= .0038


While this is > .002, considering that cross area of the steel columns in the basement level is at least 4x the cross area at floor 97, and spring constants add for springs in parallel, if we ignore the fact that the problem is for an "unloaded" (gravity free) rod, then surpassing the elastic limit at the base is impossible. (Not sure about the rest of the rod.)

For a more realistic "loaded" solution, this has to solved numerically.

Now, as to the peculiar sub-title I have chosen, I am confused as to the difference in stated results between A. E. H. Love's A Treatise on the Mathematical Theory of Elasticity, problem 281, and the results stated by Johnson in Contact Mechanics, Section 11.1

http://books.google.com/books?id=Do6WQlUwbpkC&pg=PA342&lpg=PA342&dq=%22maximum+stress%22+longitudinal+elastic+rod&source=web&ots=gneklsfh9_&sig=wXw0T7Tys1vyxFgIOtDYNaWXw-k&hl=en&sa=X&oi=book_result&resnum=4&ct=result#PPA341,M1


As far as I can see, they are describing the same problem - rigid mass hitting an elastic rod, the rod is fixed at the opposite end. Gravity is ignored.

However, as I mentioned above, Love gives maximum compression at the fixed end, for all mass ratios, while Johnson states that, for low mass ratio, maximum stress is at the end of the rod being struck.

Love doesn't actually show the derivation of his maximum compression claim, though he does state his boundary condition explicitly in his initial statement of the problem in a manner that looks "well posed", and he also actually solves the wave equation, mathematically. The boundary condition at the free, impacted end is : lim t-> +0 ( d_ w / d_ t) = - V
( "d_" is for partial derivatives) . Does this at least imply that the lim t-> +0 stress -> - rho * c0 * v, for small m? (see below).

Johnson, OTOH, doesn't even do a full derivation of the wave equation in 11.1 (though he gives some references). He uses a momentum argument to show that stress = - rho * c0 * v, where v is the (local) velocity of the rod itself, c0 is the velocity of the wave front, and rho is density. In a scenario such as this, where you consider the impactor to be rigid, it does seem like a good first approximation to take v = V, where V is the velocity of the impactor, at the instance of impact.

However, absent a full treatment, I am wondering if Johnson is simply wrong in his claims.

I also wonder if I am simply wrong in assuming that maximum stress is completely determined by maximum compression!

Interesting post. It is a nice way to approach the physics problem.


Here metamars finds 2 bodies of text on the subject we are studying which have different conclusions.

It is very important to see that the hypothetical situation in the BZ paper actually already has detailed solutions in physics literature.


This means we don't have to "reinvent the wheel" using physics 101 equations and "back of the envelope" calculations.

I've never understood why people didn't just refer to the existing physics literature on the subject before now.

Notice how there are much more complex factors involved than are considered or even mentioned in BZ.

What about the striking mass being not totally rigid? This means there would be even more "give" by both objects. The strike would be more distributed over time.

From Wave Motion in Elastic Solids, by Graff, p. 102, the ratio of initial force between the elastic rod impactor (Impedance Z1) and a rigid rod impactor on a collinear elastic rod of impedance Z2 is

1 / 1 + (Z1 / Z2) (Eq. 2.4.25)

Impedance is defined as Z = rho * A * c0

rho - density
A - cross sectional area
c0 - bar velocity ("speed of sound")

in the WTC scenario, I would take Z1 = Z2 as a first approximation ( or guess ). Therefore, the initial force will only be 1/2 of what you get when you consider the impactor rigid, as Bazant and Zhou do.

BTW, it looks like I can easily compute stress values for the OP problem via Eq. 3.92 (sigma), in Goldsmith's Impact. I'll do this sometime in the near future.

in the WTC scenario, I would take Z1 = Z2 as a first approximation ( or guess ). Therefore, the initial force will only be 1/2 of what you get when you consider the impactor rigid, as Bazant and Zhou do.

Interesting. Classical physics has existed for a few hundred years. Many of our basic physics relations have already been worked out. Why not refer to the existing literature when working out the BZ hypothetical problem?????????????

metamars wrote:From Wave Motion in Elastic Solids, by Graff, p. 102, the ratio of initial force between the elastic rod impactor (Impedance Z1) and a rigid rod impactor on a collinear elastic rod of impedance Z2 is

1 / 1 + (Z1 / Z2) (Eq. 2.4.25)

Impedance is defined as Z = rho * A * c0

rho - density
A - cross sectional area
c0 - bar velocity ("speed of sound")

in the WTC scenario, I would take Z1 = Z2 as a first approximation ( or guess ). Therefore, the initial force will only be 1/2 of what you get when you consider the impactor rigid, as Bazant and Zhou do.

BTW, it looks like I can easily compute stress values for the OP problem via Eq. 3.92 (sigma), in Goldsmith's Impact. I'll do this sometime in the near future.

We mention this type of thing in the Missing Jolt paper. Bazant was off by a factor of ten in his total axial stiffness calculation. He calculated 71 GN/m while the total axial stiffness we calculated, using the actual core column cross sections and a good estimate of the perimeter column cross sections, was 7.1 GN/m. The result of this would be that the energy would be absorbed over a longer period of time with a lower deceleration and correspongly lower force amplification. However, this is somewhat moot as there is no deceleration to be observed in the first few seconds of the falls of the upper blocks of the North Tower and WTC 7, which have been looked at for this thus far.

What happens when a steel structure collides with another steel structure has been studied extensively in Japan, e.g.:

The performance of buffer bow structures against collision. 1st report: Collapse strength of the simplified structure models.

Abstract;In order to prevent a disastrous damage induced by collision accident which might cause cargo oil spill from a struck oil tanker, the double hull system in side hull of oil tanker has been recognized as an effective countermeasures. However, when considering that ocean-going vessels are increasing not only in size but also in speed, the threat of disastrous collision accident should be further mitigated also on the responsibility of striking ships. The concepts of buffer bow have been proposed supposing its installation to fast passenger ships in the late 1960s and to large cargo vessels equipped with bulbous bow in the 1990s. RR76 panel of Shipbuilding Research Association of Japan and Ship Research Institute have started, as a joint work, to study the necessary items for the buffer bow structure since 1998. SRI has conducted a series of crush test using three kinds of simplified scale models of the buffer bow. A simple analysis method has been developed to evaluate the crashworthiness of buffer bow structure and its accuracy has been verified by comparing the results with the test and FEM calculation as well. The paper describes the performance of the buffer bow focusing on the collapse mechanism and the P/.DELTA. relation, then demonstrates the effectiveness of the present analysis method. (author abst.)

http://sciencelinks.jp/j-east/article/2 ... 728396.php

I reviewed these reports in European technical magazines back in 2001. One test was full scale - a barge with a special steel bow rammed another steel barge (in a Dutch canal) and produced plenty of structural damages to both structures. And then the same thing was simulated/calculated using FEM - or 1000+ FEM analysises of every step of the collision. It took the PCs at the time several weeks to to the calculations. That's why I always suggest that Bazant & Co can learn from ships' collisions.

BTW - these collision tests/models indicate that my oil tanker design performs better (less oil spills) than the only US approved design of the US Oil Pollution Act 1990. USCG decided to ignore this knowledge.

Well, that's interesting,.... both Metamars and Heiwa have presented discussions of collisions in which gravity is not a factor.

In W. Johnson's book "Impact Strength of Materials" two examples of longitudinal impact are considered:

(i) A mass M1 striking a long uniform bar perfectly fixed at its far end.

The equation of motion is:

M1 dV1/dt = - A2 sigma2

where V1 is the velocity of the interface, A2 the cross sectional area and sigma2 is the stress intensity at the interface.

(ii) A mass M1 falling under gravity and striking a vertical bar perfectly fixed at its far end.

The equation of motion is:

M1 dV1/dt = M1g - A2sigma2

So we have two different equations of motion giving, (surprise, surprise!), two different results!

metamars wrote:
BTW, it looks like I can easily compute stress values for the OP problem via Eq. 3.92 (sigma), in Goldsmith's Impact. I'll do this sometime in the near future.

I've now calculated this, and gotten a maximum stress far beyond the yield strength of A36 steel. In fact, it is so far beyond (3 orders of magnitude), that I am wondering if I have done something wrong.

Below is my commented VB.NET code. Results are:
First ksi interval: -10076310173.3565
Second ksi interval: -76617747629.4889
Third ksi interval, Maximum Stress: -367,536,625,701.757 Pa >> Yield strength of A36 steel, which is 250 MPa

I calculated 10,000 points in each interval, and about 140 exceeded the yield strength of A36 steel in the last interval that had maximum stress. It's a little hard to interpret this since ksi encompasses points in both time and space, but since they seem to be contiguous in ksi-space, I have to assume that there's a spatial range where stress is exceeding yield(A36).

///////////////////////////////////////////////////////////////
Code
///////////////////////////////////////////////////////////////

Private Sub Button7_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button7.Click


' From "Impact", by Werner Goldsmith
' the problem is for a rigid mass hitting an elastic rod of length L, opposite end fixed, gravity free
' Eqn 3.92 gives the stress, sigma, as:
' -E[f'<ksi> + f'<ksi - 2 L>] (note: Goldsmith uses ksi-bar)
' f' is determined iteratively, on intervals [0 , 2L], [2L , 4L], [4L , 6L], [6L , 8L]
' the mass ratio of impactor to struck rod is m = 32,815,800 / 55,700,000 = .589; Goldsmith uses M = 1 / m
' according to table 1 on p.49, contact between impactor and rod ceases in the interval [4L , 6L] since m ~ .5, thus M^-1 ~ 2, so I will not calculate beyond this interval
' ksi is defined to be c0*t -+ x , where c0 is bar velocity, 5200 m/s, t is time, and x is displacement from the impacted end


Dim ksi As Double

Dim c0 As Double = 5200 'the bar velocity
Dim L As Double = 413.0 * (98 / 111) 'length of the rod representing the lower WTC, in meters
Dim V As Double = 8.5 'impact velocity in m/s ; Goldsmith denotes this v2,0
Dim El As Double = 200000000000 '200 GPa


'Dim t_per As Double
'Dim s_per As Double
'Dim t As Double = 1 '1 second
'Dim s As Double = L 'see above
'Dim dt As Double
'Dim ds As Double

Dim M As Double = 1 / 0.589 'the mass ratio of he mass of the columns in the WTC base divided by the total mass of the WTC top >= Floor 98
Dim factor As Long = 10000

' I write s where Goldsmith uses f'; the integer following s indicates the interval it is valid on
's1 is the solution in [0 , 2L]
's2 is the solution in [2L , 4L]
's3 is the solution in [4L , 6L]
's1 = (V / c0) * (Exp((-M * ksi) / L))
's2 = s1 + (V / c0) * (1 - 2 * M * ((ksi / L) - 2)) * Exp(-M * ((ksi / L) - 2))
's3 = s2 + (V / c0) * (1 - 4 * M * ((ksi / L) - 4) + 2 * M * M * ((ksi / L) - 4) * ((ksi / L) - 4)) * Exp(-M * ((ksi / L) - 4))


Dim MaxCompression(3) As Double

'Dim ksiVal As Double
Dim i As Long
'find maximum stress in [0 , 2L], = -E(s1<ksi> + s1<ksi - 2L>)
ksi = 0
For i = 1 To factor
ksi = ksi + i * ((2 * L) / factor)
If MaxCompression(1) < s1(ksi) + s1(ksi - (2 * L)) Then
MaxCompression(1) = s1(ksi) + s1(ksi - (2 * L))
End If
Next
MaxCompression(1) = -MaxCompression(1) * El

'find maximum stress in [2L , 4L], = -E(s2<ksi> + s2<ksi - 2L>)
ksi = 2 * L
For i = 1 To factor
ksi = ksi + i * ((2 * L) / factor)
If MaxCompression(2) < s2(ksi) + s2(ksi - (2 * L)) Then
MaxCompression(2) = s2(ksi) + s2(ksi - (2 * L))
End If
Next
MaxCompression(2) = -MaxCompression(2) * El


Dim count As Long
'find maximum stress in [4L , 6L], = -E(s3<ksi> + s3<ksi - 2L>)
ksi = 4 * L
For i = 1 To factor
ksi = ksi + i * ((2 * L) / factor)

'Debug.WriteLine(s3(ksi) + s3(ksi - (2 * L)))
If s3(ksi) + s3(ksi - (2 * L)) > 250000000 / El Then
count += 1
End If
If MaxCompression(3) < s3(ksi) + s3(ksi - (2 * L)) Then
MaxCompression(3) = s3(ksi) + s3(ksi - (2 * L))
End If

Next
MaxCompression(3) = -MaxCompression(3) * El

For i = 1 To 3
Debug.WriteLine(MaxCompression(i))
Next
Debug.WriteLine("count is " & count)

End Sub
Private Function s1(ByVal ksi As Double) As Double

Dim c0 As Double = 5200 'the bar velocity
Dim L As Double = 413.0 * (98 / 111) 'length of the rod representing the lower WTC, in meters
Dim V As Double = 8.5 'impact velocity in m/s ; Goldsmith denotes this v2,0
Dim M As Double = 1 / 0.589 'the mass ratio of he mass of the columns in the WTC base divided by the total mass of the WTC top >= Floor 98

s1 = (V / c0) * (Exp((-M * ksi) / L))

End Function
Private Function s2(ByVal ksi As Double) As Double

Dim c0 As Double = 5200 'the bar velocity
Dim L As Double = 413.0 * (98 / 111) 'length of the rod representing the lower WTC, in meters
Dim V As Double = 8.5 'impact velocity in m/s ; Goldsmith denotes this v2,0
Dim M As Double = 1 / 0.589 'the mass ratio of he mass of the columns in the WTC base divided by the total mass of the WTC top >= Floor 98

s2 = s1(ksi) + (V / c0) * (1 - 2 * M * ((ksi / L) - 2)) * Exp(-M * ((ksi / L) - 2))

End Function
Private Function s3(ByVal ksi As Double) As Double

Dim c0 As Double = 5200 'the bar velocity
Dim L As Double = 413.0 * (98 / 111) 'length of the rod representing the lower WTC, in meters
Dim V As Double = 8.5 'impact velocity in m/s ; Goldsmith denotes this v2,0
Dim M As Double = 1 / 0.589 'the mass ratio of he mass of the columns in the WTC base divided by the total mass of the WTC top >= Floor 98

s3 = s2(ksi) + (V / c0) * (1 - 4 * M * ((ksi / L) - 4) + 2 * M * M * ((ksi / L) - 4) * ((ksi / L) - 4)) * Exp(-M * ((ksi / L) - 4))
End Function

End Class

The standard (conservative) approach to engineering "rods" which are subject to sudden (i.e., impact) loads is to assume that the shock wave at the fixed end is entirely reflected. Then the load due to an impact, at the impacted end, is twice that of appliying the same loading mass quasi-statically (very slowly).

Engineers designing in steel are required to use the ASCI DCR (demand to capacity ratio) method. There is a bit about this in NCSTAR1--2A. NIST states that the core columns were at 0.5 DCR before the aircraft impacts. For just before the collapse their calculated DCRs are given in figures in NCSTAR1--6D; the DCRs are on average much higher due to reduced capacity.

The load previusly bourne by the exterior walls was suddenly, transfered by the hat truss, imposed upon the core columns. This load, for stories 98 and up, appears to be higher than even the intact core columns from story 105 and up could handle; some or all of the core columns could fail there. However this might have been, it is easy to show that the sudden load, if it reached that far, would impose a load far in excess of the ability of the reduced capacity column members around story 98 could possibly bear, a DCR in excess of 2.

It turns out that DCR analysis is rather conservative. One design manual, for GSA if I recall correctly, requires the steel design engineer to calculate the DCR of all components. If the DCR is greater than 1.1, redesign to lower that is required. If the DCR is less than 0.9, the engineer is required to raise that or document why the extra strength is considered necessary. (This is for multi-story structures only a few stories tall.) The assumption used here, and which I've seen elsewhere, is that some inelastic behavior is thought to occur at a DCR of 1.15.

The other conservative aspect is the assumption of total reflection at the interface (base). For the WTC towers I worked it out on the second part of the PhysOrgForum thread(s) devoted to this topic; it is more like 1/3 than 1. However, there were also all the reflections at all the interfaces caused by changing size or shape of the column members. Each of these approximately 98/3 interfaces would induce substantial reflections, so doing the geometric progession will suggest to you that the conservative factor of two for sudden loads is about right.

The conclusion of this DCR analysis is that the core columns fail, at any place from around story 98 up to the hat truss connections under the roof.

metamars wrote: ' according to table 1 on p.49, contact between impactor and rod ceases in the interval [4L , 6L] since m ~ .5, thus M^-1 ~ 2, so I will not calculate beyond this interval

I goofed. M^-1 = m.

Besides that, the interval in which contact ceases is a time interval, [2L / c0 , 4L / c0], not a ksi interval [2L , 4L]. However, considering a fixed point on the bar such that x ~ 0 , ksi ~ c0t , so

2L < ksi < 4L

implies

2 L / co < c0t / c0 = t < 4L / c0

Since I found peak stress in the intervals that I computed not in [2L , 4L ], but in [4L , 6L], this is either another sign that I have an error in my calculation, or else my understanding of the relationship between "no contact" and "peak stress has already occurred" is wrong.

David B. Benson wrote:The load previusly bourne by the exterior walls was suddenly, transfered by the hat truss, imposed upon the core columns. This load, for stories 98 and up, appears to be higher than even the intact core columns from story 105 and up could handle; some or all of the core columns could fail there. However this might have been, it is easy to show that the sudden load, if it reached that far, would impose a load far in excess of the ability of the reduced capacity column members around story 98 could possibly bear, a DCR in excess of 2.

The hat truss was designed to spread the wind loads from the antenna among all of the core columns and then used the outriggers to the perimeter for buttressing and stiffness.

I wish I had the actual sizes of the hat truss outriggers and their connection to the exterior columns. In looking at photos it seems highly unlikely that the hat truss outriggers could have transmitted 12 stories of the perimeter's gravity load to the core before failing themselves.

T_Szamboti wrote:The hat truss was designed to spread the wind loads from the antenna among all of the core columns and then used the outriggers to the perimeter for buttressing and stiffness.

This is widely misunderstood, I am sorry to say. Note that WTC 2, with no antenna towerr, also had a hat truss.

Read NCSTAR1--2A to understand the function of the hat truss in distributing wind loads from the exterior walls to the core. The antenna tower made little difference.

I wish I had the actual sizes of the hat truss outriggers and their connection to the exterior columns. In looking at photos it seems highly unlikely that the hat truss outriggers could have transmitted 12 stories of the perimeter's gravity load to the core before failing themselves.

There is enough data in NCSTAR1--6D to determine that indeed, such was possible without even reaching the elastic limit of the outrigger welds to the exterior walls. (Note that NIST calculated that one of these outrigger connections failed completely in WTC 2.)

There are even some graphics of the hat truss, showing where it was connected and important load paths.

A standard principle of engineering design is to make the connection as strong as the weakest of the members being connected. So we can be quite sure that the hat truss was at least as strong as the welded connection to the exterior walls.

David B. Benson wrote:

T_Szamboti wrote:The hat truss was designed to spread the wind loads from the antenna among all of the core columns and then used the outriggers to the perimeter for buttressing and stiffness.

This is widely misunderstood, I am sorry to say. Note that WTC 2, with no antenna towerr, also had a hat truss.

Read NCSTAR1--2A to understand the function of the hat truss in distributing wind loads from the exterior walls to the core. The antenna tower made little difference.

I wish I had the actual sizes of the hat truss outriggers and their connection to the exterior columns. In looking at photos it seems highly unlikely that the hat truss outriggers could have transmitted 12 stories of the perimeter's gravity load to the core before failing themselves.

There is enough data in NCSTAR1--6D to determine that indeed, such was possible without even reaching the elastic limit of the outrigger welds to the exterior walls. (Note that NIST calculated that one of these outrigger connections failed completely in WTC 2.)

There are even some graphics of the hat truss, showing where it was connected and important load paths.

A standard principle of engineering design is to make the connection as strong as the weakest of the members being connected. So we can be quite sure that the hat truss was at least as strong as the welded connection to the exterior walls.

While the hat trusses had a couple of purposes, it is incorrect to say their primary purpose was not to deal with the antenna loads. WTC 2 would have gotten a hat truss in case it was used for an antenna mount in the future. Did WTC 7 have anything like a hat truss?

I will look at the NIST information on the hat truss more closely before continuing discussion on this issue.

T_Szamboti wrote:Did WTC 7 have anything like a hat truss?

No. While not exactly conventional design, it was a steel space frame without trussed floors. That, combined with "only" 47 stories, meant it did not need one.

metamars wrote:

metamars wrote:
BTW, it looks like I can easily compute stress values for the OP problem via Eq. 3.92 (sigma), in Goldsmith's Impact. I'll do this sometime in the near future.

I've now calculated this, and gotten a maximum stress far beyond the yield strength of A36 steel. In fact, it is so far beyond (3 orders of magnitude), that I am wondering if I have done something wrong.

I went ahead and calculated the stress values for specific times, to see what the pressure wave would look like. Code for this is below, using the same s (f' in Goldsmith) functions I used earlier. In a subsequent post, I will show a sanity check calculation comparing energy*. At t = .03 seconds, I get an average compression of 8756 Mpa, over just 1.787m of the bar. (The rest of the values were zero) Later outputs show a much more variation of pressure, but occurring over shorter distances.


Sorry, the sanity check did indeed fail, so things are not consistent. The note below is wrong.

( * note: edited to remove " which is clearly inconsistent. Can anybody tell me what I'm doing wrong?" I double checked the sanity calculation just before posting, and when I corrected it, it came out close enough)

///////////////////////////////////////////////////
Code
//////////////////////////
Private Sub Button8_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Handles Button8.Click


' From "Impact", by Werner Goldsmith
' the problem is for a rigid mass hitting an elastic rod of length L, opposite end fixed, gravity free
' Eqn 3.92 gives the stress, sigma, as:
' -E[f'<ksi> + f'<ksi - 2 L>] (note: Goldsmith uses ksi-bar)
' f' is determined iteratively, on intervals [0 , 2L], [2L , 4L], [4L , 6L], [6L , 8L]
' the mass ratio of impactor to struck rod is m = 32,815,800 / 55,700,000 = .589; Goldsmith uses M = 1 / m
' according to table 1 on p.49, contact impactor and rod ceases in the interval [4L , 6L] since m ~ .5, thus M^-1 ~ 2, so I will not calculate beyond this interval
' ksi is defined to be c0*t -+ x , where c0 is bar velocity, 5200 m/s, t is time, and x is displacement from the impacted end


Dim ksi As Double

Dim c0 As Double = 5200 'the bar velocity
Dim L As Double = 413.0 * (98 / 111) 'length of the rod representing the lower WTC, in meters
Dim V As Double = 8.5 'impact velocity in m/s ; Goldsmith denotes this v2,0
Dim El As Double = 200000000000 '200 GPa

Dim M As Double = 1 / 0.589 'the mass ratio of he mass of the columns in the WTC base divided by the total mass of the WTC top >= Floor 98
Dim xfactor As Long = 10000
Dim tfactor As Long = 1000
' I write s where Goldsmith uses f'; the integer following s indicates the interval it is valid on
's1 is the solution in [0 , 2L]
's2 is the solution in [2L , 4L]
's3 is the solution in [4L , 6L]
's1 = (V / c0) * (Exp((-M * ksi) / L))
's2 = s1 + (V / c0) * (1 - 2 * M * ((ksi / L) - 2)) * Exp(-M * ((ksi / L) - 2))
's3 = s2 + (V / c0) * (1 - 4 * M * ((ksi / L) - 4) + 2 * M * M * ((ksi / L) - 4) * ((ksi / L) - 4)) * Exp(-M * ((ksi / L) - 4))


'Dim MaxCompression(3) As Double

' declare arrays which will record compression along the length of the bar, at fixed time intervals
Dim Compression(tfactor, xfactor) As Double 'use ksi = c0t - x
Dim Compression2(tfactor, xfactor) As Double 'use ksi = c0t + x


Dim t As Double = 0


Dim j As Long
Dim i As Long

For j = 1 To tfactor
t = t + (j * 0.6) / tfactor

x = 0
For i = 1 To xfactor
x = x + (i * L) / xfactor


ksi = c0 * t - x
If ksi >= 0 And ksi < 2 * L Then
Compression(j, i) = s1(ksi) + s1(ksi - (2 * L))

ElseIf ksi >= 2 * L And ksi < 4 * L Then
Compression(j, i) = s2(ksi) + s2(ksi - (2 * L))

ElseIf ksi >= 4 * L And ksi < 6 * L Then
Compression(j, i) = s3(ksi) + s3(ksi - (2 * L))
End If
Compression(j, i) = El * Compression(j, i)


ksi = c0 * t + x
If ksi >= 0 And ksi < 2 * L Then
Compression2(j, i) = s1(ksi) + s1(ksi - (2 * L))

ElseIf ksi >= 2 * L And ksi < 4 * L Then
Compression2(j, i) = s2(ksi) + s2(ksi - (2 * L))

ElseIf ksi >= 4 * L And ksi < 6 * L Then
Compression2(j, i) = s3(ksi) + s3(ksi - (2 * L))
End If
Compression2(j, i) = El * Compression(j, i)


Next
Next

'examine some results, dividing by 1,000,000 to make plotting in Excel easier
j = 5
Debug.WriteLine("======== t = " & (j / tfactor) * 0.6 & " seconds =======")
Debug.WriteLine("======== Distance in meters | Compression in MPa =======")
For i = 1 To xfactor 'Step 5
If Compression(j, i) > 0 Then
Debug.WriteLine(Round((i * L) / xfactor, 3) & vbTab & Round(Compression(j, i) / 1000000, 0))
End If
Next


j = 10
Debug.WriteLine("======== t = " & (j / tfactor) * 0.6 & " seconds =======")
Debug.WriteLine("======== Distance in meters | Compression in MPa =======")
For i = 1 To xfactor 'Step 5
If Compression(j, i) > 0 Then
Debug.WriteLine(Round((i * L) / xfactor, 3) & vbTab & Round(Compression(j, i) / 1000000, 0))
End If
Next
End Sub

///////////////////////////////////////////////////
Results
///////////////////////////////////////////////////
======== t = 0.003 seconds =======
======== Distance in meters | Compression in MPa =======
0.036 8107
0.073 8110
0.109 8114
0.146 8120
0.182 8127
0.219 8135
0.255 8145
0.292 8156
0.328 8168
0.365 8182
0.401 8197
0.438 8214
0.474 8232
0.51 8252
0.547 8273
0.583 8295
0.62 8319
0.656 8345
0.693 8372
0.729 8400
0.766 8430
0.802 8462
0.839 8495
0.875 8530
0.912 8566
0.948 8604
0.985 8643
1.021 8685
1.057 8727
1.094 8772
1.13 8818
1.167 8866
1.203 8916
1.24 8968
1.276 9021
1.313 9077
1.349 9134
1.386 9193
1.422 9254
1.459 9317
1.495 9382
1.531 9449
1.568 9518
1.604 9590
1.641 9663
1.677 9739
1.714 9817
1.75 9897
1.787 9980
1.823 10065