Plane impact equivalence to force in time

[javimm]

The problem comes up as part of a final master's project in FEM analysis. The project consists in calculating the oscillation response of the World Trade Center towers after the plane impacts. The particular problem I'm having trouble is that I need to transform the plane impact of which mass and speed is known, to data the FEM software can use: force in a period of time. The software allows you to define a table of points of F and t.
The data I have is:
Plane speed: 166.67 m/s
Mass: 170 metric tonnes.
Length= 50m

The plane penetrates 30 m inside the building before it stops.

I've tried two suggested solutions.
- Aceleration is constant. I use the classic equations:
F=m*a
a=a_0 (constant)
v=a_0*t+v_0 (v_0=166.67 m/s)
s=.5*a_0*t^2+166.66*t+s_0 (s_0=0 as I use the very moment of impact as reference for t=0)

Solving with the data:
v_final^2-v_0^2=2*a_0*s; a_0=(0-166.67^2)/2*30=-462.98 m/s^2
From velocity equation: 0=-462.98*t+166.67; Solving for t gives t=0.36 s (this is the impact time)

Equivalent Force:
F=170000*462.98=7870660 N
So the plane impact is equivalent to a Force of 7870660 applied during 0.36s. Using the software, it's equivalent to a "pulse": A horizontal line during .36 seconds.

Is this correct?.
After running the FEM sim, maximum deflection at top of tower gives 0.4 m, that I assume as correct, as the real tower swinged 0.6 m and after all simplifications of my model I think is close enough.

- Second solution: deceleration is linear (not constant). 0 at the moment of impact and maximum when the plane stops.
Our acceleration will be:
a=a_plane*t;(Assuming a=0 at t=0)
Integrating the expression we can obtain velocity and distance:
v=a_plane*t^2/2+v_0
s=a_plane*t^3/6+v_0*t
We know that at s=30→v=0. We also know that v_0=166.67 m⁄s . We obtain a non-linear two equation system. Solving it we obtain than a_plane=-4572.4188 m⁄s^2 and t=0.27 s. The plane stops at 0.27 s from impact.
The force we need to define in is:
F=m∙a=170000∙-4572.4188∙t=-7.773∙10^8∙t N, applied during 0.27s, gives a triangular shape in the software F-t graphic.

This solution gives a deflection of several meters, so it's not correct, but I don't see where I made a mistake.

The area of the graphics for both assumptions should be the same, but it isn't.

Any ideas?

Thank you very much.

 

[voko]

All of the figures after the decimal point in your value for a_plane are wrong. The correct value is more like -4572.7.

 

[javimm]

All of the figures after the decimal point in your value for a_plane are wrong. The correct value is more like -4572.7.

That is caused by rounding error from converting the plane speed from km/h to m/s. It dooesn't change the result by much though.

 

[voko]

I do not see any other errors in your second solution. If the final result that follows from the second solution is too far off, then the assumptions in the second solution are wrong.

 

[javimm]

I do not see any other errors in your second solution. If the final result that follows from the second solution is too far off, then the assumptions in the second solution are wrong.

That's what's driving me crazy. In theory both solutions should lead to very similar results, but they don't and by a large margin. I spotted a mistake on my part using the software, and the deflection is not several meters, but they one is twice the other (see below).
The area under the F-t graphics should be more or less the same, but the area of the first assumption (constant acceleration) is 78706600*.36=28334376 and the area of the linear deceletarion gives 0.5*.27*(7.773∙10^8*.27)=28332585 that is the same.
The difference in deflection for the two cases are double: less than 0.4 m in the first assumption and 0.7 m in the second. One is double of the other.
Still I don't spot the mistake.
Thank you for your time.

 

[voko]

The area under the Ft graph must be equal in magnitude to the original momentum of the plane. Because the mass stays constant, it means that the area under the acceleration-time graph must be equal in magnitude to the original speed of the plane. Now, the speed of the plane was 166.67 m/s. In the constant acceleration case, the area under the a-t graph is 462.98 * 0.36 = 166.6728; in the linear deceleration case, it is 4572.7 * 0.27^2 / 2 = 166.6749. Both match the original speed fairly closely. There must be an error in your multiplying the acc. by mass that gives you such a huge (and round) error.

 

[javimm]

I've found a typing error in my first and fifth posts that I'm going to correct. As you said, both areas are the same. I've checked the numbers I've used in the software and they seem to be correct, but there is a difference between the two deflections. I guess using these two assumptions can lead to that, but I'm not completely sure why the difference in deflection is that big.
EDITED TO ADD: I've checked the areas directly in the software and they are the same, but the generated deflections are different. Could be the linear application of force and the instant "pulse" that is causing this.

 

[voko]

So how big is the difference? You gave two very different magnitudes earlier.

 

[javimm]

Yes, sorry about that. As I said there was an error in the data I entered into the program. Now I've tested the areas and values directly in the software and they do match what they are supposed to be.
- With constant acceleration, maximum deflection is 0.4 m
- With linear acceleration, maximum delection is 0.7 m that is quite the same value as the real deflection of the tower.
I think a slight difference in the estimates of the plane variables can lead to ample results.

 

[voko]

Given the extreme simplicity of both models, I would say the agreement with the observed facts is remarkably good.

 

[AlephZero]

EDITED TO ADD: I've checked the areas directly in the software and they are the same, but the generated deflections are different. Could be the linear application of force and the instant "pulse" that is causing this.

If you want to investigate that, either mathematically or numerically with your FE software, compare your impact times with the natural vibration frequencies of the building.

It might be instructive to look at a simpler model, for example the response of single degree of freedom mass on a spring to a constant amount of impulse (force x time) applied with different profiles of force against time. Then try to relate that to what happens in your large model.

 

[javimm]

Given the extreme simplicity of both models, I would say the agreement with the observed facts is remarkably good.

Yes, indeed.
The tower is modeled as a simple cantilever beam, using the first normal mode of T=10 s in order to obtain the physical properties. No floors or columns are modeled, just the beam with a fictional material that makes the tower has the desired mass and first normal mode. It's a very interesting exercise, as it approaches the real event behaviour with very good accuracy. The model oscilates for 4 minutes, just as the real towers did. The exercise tries to prove that they absorbed the impact of the planes no problem. Project wind made the towers swing more or less the same that the planes did, as far as I know. Fire was the real problem there.
Thank you for your answers.

 

[javimm]

If you want to investigate that, either mathematically or numerically with your FE software, compare your impact times with the natural vibration frequencies of the building.

First and second natural frequency of the model is 0.1 Hz (T=10s). It's one of the imposed data of the exercise (see my previous post). Third mode is 0.571 Hz.

 

[AlephZero]

So the structure takes 5 sec for half a cycle of vibration of mode 1 and 0.87 sec for mode 3, compared with about 0.36 sec of impact time (calculated from average speed / distance).

I would expect the differences in the max amplitude to be mainly caused by amplitude and phase difference in the response of mode 3 relative to mode 1. The impact is probably short enough for the mode I response to be dominated by the total impulse rather than the time-variation of the force, but the force is applied over about 40% of the half-cycle of mode 3 so the response in mode 3 will be more sensitive to the time-variation.

You haven't said what is the main objective of your project, but decomposing the response of your model into modal components would be an interesting way to get more insight into what it is telling you (unless all this is just a detail in terms of your main objective, of course).

 

[javimm]

AlephZero, thanks for the insight.
The object of the problem is modelling a simplified version of the towers, obtaining their approximate physical properties (density, E, etc) knowing their dimensions and the period of the first mode. With that, a model of the tower can be made. More data is provided, (damping is 5% of critical).
Secondly we have the planes, and the exercise asks for the response during the first 20 seconds after impact (maximum stresses at the bottom of the tower), using a transient analysis. It also asks to compare this response to a static wind load. I suppose this is a way of putting in perspective the magnitude that the impact had in comparison of what they were designed for.
The exercise is short, but as it's the project for the FEM master, I'm gathering a lot of info and explaining every detail I can think of.
About decomposing the response of the model into modal components, I'm forcing the model to only include the first mode. There's an option for that in the damping section you can define.
Using a w_min=0.05 Hz and w_max=0.2 Hz (min and max freq of the structure) gives the results I've been posting.
Testing with a w_min=0.05 Hz and w_max=0.6 Hz gives a maximum deflection of 0.8 m instead of 0.7 m and the damping is less apparent as in the first case where it goes from 0.7 m to 0.5 m to 0.3 m to 0.2 m and including the third mode goes from 0.8 m to 0.7 m to almost 0.6 m to 0.5 m.
Is it more accurate to include the third mode too?.
Thank you.

 

[AlephZero]

Is it more accurate to include the third mode too?.

That depends very much on whether the third mode is an accurate model of the structure.

If you think about superimposing two sine-wave responses at different frequencies and amplitudes, the resultant peak amplitude is sensitive to when the two peaks occur (i.e. the frequency ratio of the modes) as well as their relative amplitudes.

I'm not a civil engineer, so I don't know whether or not a uniform cantilever beam is a realistic model of a tower. Tweaking a model to get one frequency correct is easy. Getting several frequencies correct, and the amplitudes of the response in the correct ratio, is a harder problem.

It would certainly be more accurate to include more modes in an accurate structural model, but the key part of that statement is the worlds in bold.

 

[javimm]

That's exactly was my initial toughts when I purposely left the third frequency out. As the model is way simplistic, I tought that the only frequency that I had to consider was the first as it was the one I imposed from the beginning. As you mentioned, I had some doubts if the third one could be included, but seeing the model I guess it's better not to. That I'm going to explain in the final article I'm writing.
I'm a Civil Engineer that happens to be studying this master a lot of years after finishing my main studies and having been working in several civil engineering tasks other than structures and FEM, so I'm more than rusty on my dynamic theory, but a cantilever beam is far from a correct model for a building. It can be used for simple calculations, but natural modes for a multi-story building differs greatly from those of a simple beam. Here's a little graphic with an example of a real multi-story building and its first three modes (taken randomly from google):

That's why I didn't include it at first.
Thanks.

 

[AlephZero]

Those mode shapes make the point rather well. I'm assuming the plane was traveling from left to right. If the impact was at 75% height you would hardly excite mode 2 at all, because the impact is near a nodal point. On the other hand it the impact was at 50% you would excite mode 2 with the initial motion in the opposite direction to the diagram, and you might find the top of the tower started moving in the opposite direction to impact before it swung in the same direction, because mode 2 is a higher frequency than mode 1.

Mode 3 would be a similar story to mode 2, except 50% is now near the nodal position.