NASA/TM-2009-215704

A Computational Approach for Probabilistic 

Analysis of Water Impact Simulations 

Lucas G. Horta, Brian H. Mason, and Karen H. Lyle 

Langley Research Center, Hampton, Virginia 

 

April 2009 
NASA STI Program . . . in Profile

 Since its founding, NASA has been dedicated to 

the advancement of aeronautics and space science. 

The NASA scientific and technical information (STI) 

program plays a key part in helping NASA maintain 

this important role. 

 The NASA STI program operates under the 

auspices of the Agency Chief Information Officer. It 

collects, organizes, provides for archiving, and 

disseminates NASA’s STI. The NASA STI program 

provides access to the NASA Aeronautics and Space 

Database and its public interface, the NASA Technical 

Report Server, thus providing one of the largest 

collections of aeronautical and space science STI in 

the world. Results are published in both non-NASA 

channels and by NASA in the NASA STI Report 

Series, which includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 

completed research or a major significant phase 

of research that present the results of NASA 

programs and include extensive data or 

theoretical analysis. Includes compilations of 

significant scientific and technical data and 

information deemed to be of continuing 

reference value. NASA counterpart of peerreviewed formal professional papers, but having 

less stringent limitations on manuscript length 

and extent of graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 

and technical findings that are preliminary or of 

specialized interest, e.g., quick release reports, 

working papers, and bibliographies that contain 

minimal annotation. Does not contain extensive 

analysis. 

• CONTRACTOR REPORT. Scientific and 

technical findings by NASA-sponsored 

contractors and grantees. 

• CONFERENCE PUBLICATION. Collected 

papers from scientific and technical 

conferences, symposia, seminars, or other 

meetings sponsored or co-sponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 

technical, or historical information from NASA 

programs, projects, and missions, often 

concerned with subjects having substantial 

public interest. 

• TECHNICAL TRANSLATION. Englishlanguage translations of foreign scientific and 

technical material pertinent to NASA’s mission. 

 Specialized services also include creating custom 

thesauri, building customized databases, and 

organizing and publishing research results. 

 For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page at 

http://www.sti.nasa.gov

• E-mail your question via the Internet to 

help@sti.nasa.gov

• Fax your question to the NASA STI Help Desk 

at 443-757-5803 

• Phone the NASA STI Help Desk at 

443-757-5802 

• Write to: 

 NASA STI Help Desk 

 NASA Center for AeroSpace Information 

 7115 Standard Drive 

 Hanover, MD 21076-1320
NASA/TM-2009-215704

A Computational Approach for Probabilistic 

Analysis of Water Impact Simulations 

Lucas G. Horta, Brian H. Mason, and Karen H. Lyle 

Langley Research Center, Hampton, Virginia 

National Aeronautics and 

Space Administration 

Langley Research Center 

Hampton, Virginia 23681-2199 

 

April 2009 
Available from: 

NASA Center for AeroSpace Information 

7115 Standard Drive 

Hanover, MD 21076-1320 

443-757-5802 
A Computational Approach for Probabilistic Analysis of 

Water Impact Simulations 

Lucas G. Horta, Brian H. Mason, and Karen H. Lyle 

Abstract 

NASA’s development of new concepts for the Crew Exploration Vehicle Orion presents 

many similar challenges to those worked in the sixties during the Apollo program. However, 

with improved modeling capabilities, new challenges arise. For example, the use of the 

commercial code LS-DYNA1, although widely used and accepted in the technical 

community, often involves high-dimensional, time consuming, and computationally intensive 

simulations. Because of the computational cost, these tools are often used to evaluate specific 

conditions and rarely used for statistical analysis. The challenge is to capture what is learned 

from a limited number of LS-DYNA simulations to develop models that allow users to 

conduct interpolation of solutions at a fraction of the computational time. For this problem, 

response surface models are used to predict the system time responses to a water landing as a 

function of capsule speed, direction, attitude, water speed, and water direction. Furthermore, 

these models can also be used to ascertain the adequacy of the design in terms of probability 

measures. 

This paper presents a description of the LS-DYNA model, a brief summary of the response 

surface techniques, the analysis of variance approach used in the sensitivity studies, equations 

used to estimate impact parameters, results showing conditions that might cause injuries, and 

concluding remarks. 

1. Introduction 

During the development of the Apollo command module (CM), many studies were 

conducted to understand the behavior of the CM upon returning to earth. Because of the 

limited computational capabilities at the time, engineers made extensive use of experimental 

 

1 Trade names and trademarks are used in this report for identification only. Their usage does not 

constitute an official endorsement, either expressed or implied, by the National Aeronautics and Space 

Administration

1 
data to complement the relatively simple analyses. These studies were performed to 

understand the landing behavior under a variety of conditions. Interestingly, Apollo was 

initially to include a land-landing architecture. Thus, many of the early studies focused on 

assessing attenuation systems for landing on soil, sand, and clay surfaces [1-3]. However, the 

land-landing capability was abandoned for water splashdown, similar to the Mercury and 

Gemini capsules [4-5]. The exceptions to the water landing studies, were assessments of the 

capsule response in the event of a pad abort [6]. 

The Apollo water landing studies built on the extensive experience of the Mercury and 

Gemini programs. A wide range of experimental parameter studies, using subscale test 

articles, was conducted to evaluate the Apollo capsule performance. A primary interest of 

these studies was to understand the effect of landing attitude and longitudinal or horizontal 

velocity, as well as heatshield flexibility, on capsule response [7-11]. In addition, numerous 

full-scale tests were conducted to assess the capsule performance on more flight-like vehicles, 

e.g., References [12-14]. Finally, based on knowledge of probable landing conditions, a 

design envelope was established [15]. The purpose of the present work is to provide a 

probabilistic framework for the Orion Crew Module water landing performance boundary 

evaluation. 

Results from two independent studies are reported here, using data from two sets of LSDYNA simulations, to illustrate the computational approach. One of the simulations used 4 

parameters, whereas the second used 7 parameters to study landing of a rigid capsule on 

water. Although both efforts used the same LS-DYNA model, because the initial conditions 

used in each study are not identical, results cannot be easily combined into a single database. 

Instead results, in the form of time histories, from each study are treated separately and used 

to predict the Brinkley index [16] and to conduct statistical analysis. Fundamental differences 

in the approaches will be pointed out as results are presented. 

For readers not familiar with the Brinkley index, it is a criterion, initially developed for 

emergency escape systems, to assess the probability of injury to a crew member when 

exposed to sudden acceleration pulses, like those observed during landings. Using local 

accelerations at the crew locations in all three directions, the criterion divides the risk of 

injury into three categories: 1) low risk corresponds to a 0.5% probability of injury, 2) 

moderate risk corresponds to a 5% probability of injury and, 3) high risk corresponds to a 

2
50% probability of injury. Although there have been efforts to modify the thresholds set 

initially by Brinkley, the calculation of the index is unchanged and is what is used throughout 

the paper. 

The paper presents a computational approach to assess impact conditions (velocities and 

orientation of capsule and impact plane) that could result in injuries to the crew when landing 

a rigid capsule on water. For this problem, a rigid LS-DYNA model is used and described 

first, followed by a detailed discussion of the formulation, the process to create the LS-DYNA 

database, sampling of the parameter space, sensitivity using conventional and analysis of 

variance, development of response surface models, and predictions of injury boundaries using 

the Brinkley index [16]. 

2. Model Description 

The finite element (FE) model of Orion consists of three LS-DYNA [17] parts: the Orion 

capsule, a volume of air, and a volume of water. The capsule is built from Lagrangian rigid 

shell elements using 10307 elements, with an average element size of 3 in. Six nodes in the 

capsule model representing astronaut locations are used for extraction of acceleration results. 

For the water, a total of 1,282,500 cubic elements (4 in. on a side) are used whereas the air 

uses 967,500 elements. Mass and mass moment of inertia properties for the capsule are given 

in Table 1. The volumes of air and water are composed of Eulerian eight-node brick elements 

with dimensions given in Fig. 1. The air is defined with a vacuum material, and the water is 

defined with a null material with a density of 9.3365e-5 lbf s

2

 / in4 and a viscosity of 1.63e-7. 

The FE model is analyzed using the Arbitrary Lagrangian-Eulerian (ALE) formulation in 

LS-DYNA. The element time step for this model is 22.8 μs, resulting in an analysis time of 

20 hrs per 100 ms of simulation time on a single node of a LINUX-based computer. 

Acceleration results are recorded to the nodout text file. 

Because the impact conditions for the Orion are uncertain, seven critical impact parameters 

have been identified for study, as defined in Table 2, and are assumed to have a uniform 

distribution between the bounds provided in the table. Four of these parameters (horizontal 

velocity, vertical velocity, capsule pitch, and wave slope) are illustrated in Fig. 2. Also noted 

on the capsule are approximate crew locations 1-6. Transformations for both the capsule and 

the Euler grid (water and air) are defined in a LS-DYNA keyword file that is separate from 

3
the rest of the model. Separate keyword files for the Euler grid and the capsule are rotated, 

translated, and linked to the analysis run-stream using the *INCLUDE_TRANSFORM 

option. 

 

2.1 Mesh Convergence Study

In this section, the effects of mesh refinement on the FE model are investigated. FE models 

with a rigid capsule at a 38-degree pitch angle and a vertical velocity of 50 ft/s are analyzed 

using 4 different Eulerian models. In all 4 cases, only the thickness of the brick elements in 

the vertical (global X) direction are varied; the in-plane dimensions for the brick elements are 

fixed at 4 in. wide by 4 in. long. Dimensions for the 4 cases are given in Table 3. Also in 

Table 3 are the number of “interface” elements between air and water which are within 8 in. 

of the boundary surface; “far-field” elements are further away from the air/water boundary. 

Results of the mesh convergence study are shown in Fig 3 in terms of time histories of the 

Brinkley index for cases 1 to 4. Time histories for case 2 are within 0.01% of case 1 (the 

baseline case), but the run time for case 2 is 30% less than for case 1. The maximum Brinkley 

index for cases 3 and 4 are 4.0% and 1.6% higher than the baseline case, respectively. 

Unfortunately, the analysis time for cases 3 and 4 are 4.25 times longer than case 1. Because 

of the computational expense of the more refined analyses (cases 3 and 4), the 4% difference 

in the refined and baseline analyses is considered small, and the solution is considered to be 

converged for case 1. Hence, results to follow will use the parameters in Case 1. 

3. Response Surface Description 

Response surface techniques provide a mathematical framework to synthesize models of 

systems for which users have only input/output data. For our purposes, these techniques will 

be used to provide time history predictions of water landings given a limited number of LSDYNA simulations. It is assumed that in order to solve the problem at hand, the amount of 

time to complete all required LS-DYNA runs is prohibited. Therefore an alternate way to 

obtain solutions, i.e., response surface models, is needed. Two techniques have been 

implemented and used; the first one was developed by Sacks [18] is often referred to as 

Kriging and the second technique is the Extended Radial Basis Functions [19-20]. Both 

4
techniques are briefly discussed to aid the reader in understanding some of the parameter 

settings used in our study. 

 

3.1 Kriging Response Surface

In this approach, estimates of the system response are computed using the following 

expression 

ˆ 1 ˆ( ) ( ) ( 1) T y r RY ˆ ν βν β − =+ − (3.1) 

with 

11 1

2

1

1 2

1 2

ˆ (1 1) 1

( , ) exp[ ( ) ]

( ) [ ( , ) ( , ) ( , )]

[ ]

1, 2,..., 1, 2,...,

T T

m k k

i j i j

k

T

N

T

N

R RY

R

rR R R

Y yy y

i N and j N

β

νν θ ν ν

ν νν νν νν

−− −

=

=

=− −

=

=

= =

∑

...

...

where is k

th component of the ith parameter vector, v k

i v m×1 ∈\ is a vector representing the 

parameters, N is the population size, y and ˆ

β are scalars, N N R × ∈\

N 1

is a correlation matrix for 

the parameters, is a correlation vector, Y 1 ( ) N

r ν × ∈\ × ∈\ is a vector with N response 

values, is a vector whose elements are all ones, ×1 1 N ∈\ θ is a scalar parameter yet to be 

computed, and m is the number of parameters. Although Eq. (3.1) is written for a system with 

a scalar output, the formulation is applied sequentially for cases with multiple outputs. A 

subtle but critical aspect of this formulation is the fact that the unknown parameterθ , used to 

define the correlation matrix R, must be computed using maximum likelihood estimation. 

That is, θ is the value that maximizes the scalar metric 

( ) 1 2 ( ) ln( ) ln , [0, ] 2

LN R θσ θ = −+∈ ∞ (3.2) 

where the standard deviation is computed as 

2 1 1 ˆ ( 1) ( T Y RY

N

ˆ σ β 1) − =− − β 

Solving this scalar optimization problem produces a response surface fitting error with a 

Gaussian distribution that matches the N response values. However, estimation of the 

parameter θ requires one to solve a scalar optimization for each output vector; which is 

5
conducted here using the genetic optimization algorithm in [21]. For problems containing q

outputs sampled N times, to generate a Kriging surface requires qN optimizations of the 

metric in Eq.(3.2). Although computationally intensive, the solutions are very accurate. 

3.2 Extended Radial Basis Functions

The Extended Radial Basis Function (ERBF) method developed by Mullur and Messac in 

[19], approximates the system response as follows 

1 1

() ( ) ( )

N N

i ii

i i

y ν i τψ ν ν φ ν ν

= =

= −+ ∑ ∑ −  (3.3)

where ψ is a radial basis function, i τ is a scalar parameter solved as part of the least-square 

fitting process, and the function φ (known as the extended function) takes on different values 

depending upon the value of ν −ν i , see [19] and [20]. To solve for φ requires solving a least 

square problem with 3mN parameters. Since there are N i τ values, the total number of 

parameters needed to compute an ERBF model is N(3m+1). In addition, to completely 

defineφ , two free parameters must be chosen: 1) the order of a local polynomial (set to 4 in 

this example), and 2) the smoothness parameter γ (set to 0.15). Finally, the radial basis 

function ψ is chosen as an exponentially decaying function with characteristic 

radius set to 0.15. 

2 2 / 2 i − − ( ) ν ν rc e

cr

In both Kriging and ERBF, scaling of the parameter variables ν is extremely important. In 

our implementation the absolute value of the maximum parameter value is used for scaling 

such that when variables are at their maximum value the normalized variable equals 1. 

4. Parameter Sampling 

Parameter sampling is perhaps one of the most important steps when creating response 

surfaces. Techniques such as Latin Hypercube [22], random sampling, Hammersly [23], Doptimal[24], Halton-leap [25-26], and even manual setting of parameter variables are just a 

few of the possible means to prescribe a sample parameter. In our study D-optimal, Haltonleap, and manual parameter selection have been used at different times for different reasons. 

For example, in the beginning of the study the D-Optimal sampling approach was used 

because it would be followed by an analysis of variance study and this technique is ideally 

6
suited for that purpose. On the other hand, for the single variable sensitivity study, being 

conducted in parallel, manual sampling (heuristic approach) made the most sense. 

Nonetheless, for good coverage of the parameter range, Halton-leap deterministic sampling 

approach is also used. 

Later in the paper, results for two independent studies using 4 and 7 independent variables 

are presented. For the 4-variable study, sampling of the parameter variables is conducted 

manually to study specific areas of concern and to facilitate sensitivity studies one parameter 

at a time. In contrast, the 7-parameter study used D-Optimal, Halton-leap, and manual 

sampling. 

5. Analysis of Variance 

Parameter sensitivity in most engineering fields in often associated with derivative 

calculations at specific points in the parameter space. However, for analysis of nonlinear 

systems under uncertainty, sensitivity studies are often conducted using Analysis of Variance 

(ANOVA). In classical ANOVA studies, data is collected from multiple experiments while 

varying all parameters (factors) and also varying one parameter at a time. These results are 

then used to assess the output response variance due to variations of a particular parameter as 

compared to the total output variance when varying all the parameters simultaneously. The 

ratio of these two variance contributions is a direct measure of the parameter importance. 

Sobol in [27] and others [27-31] have studied the problem of global sensitivity analysis 

using variance measures. Sobol developed the "so called" Sobol indices to provide a measure 

of parameter sensitivity. In his formulation the variance is computed using the following 

expressions; 

1

2 2

1

2

1

( )

( )

( )( , )

N

i

i

N

i

i

N

x i ii

i

y

D y

D y yx

μ ν

μ ν

μ ν

=

=

=

≅

+ ≅

+ ≅

∑

∑

∑ z





 

(5.1)

where μ is the mean value, D is the total variance computed using N samples of the function, 

and Dx is the single parameter variance due to parameter x. Before using Eq. (5.1) it is 

important to point out a notation change. As before, ν i refers to the i

th sample of the 

7
parameter vector but xi is a parameter within ν i about which the variance is being evaluated 

and zi are all other parameters that comprise ν i not including xi. To properly evaluate Eq. 

(5.1), one needs at least 2N function (output) evaluations; one where all variables are 

randomly sampled, and a second set where all but xi are re-sampled. With this information 

the 1st index in Sobol is computed as / x S D= x D . Of course, this is only one of the m possible 

factors. Incidentally, Eq. (5.1) is easily extended to study two or more parameter (factor) 

interaction, as described in [27], simply by adjusting the number of parameter that gets resampled. 

A final note on ANOVA using Sobol's approach is regards convergence of the variance 

estimates. Although Ref. 27 discussed the asymptotic behavior of the variance estimates 

using Eq. (5.1), for cases where a surrogate model is used instead of LS-DYNA, the variance 

estimates are only as good as our surrogate model. Nonetheless, this approach provides an 

excellent way to rank variable importance when only a limited number of LS-DYNA runs are 

possible. 

6. Computational Framework 

To conduct this statistical study, it is preferable to automate the generation of LS-DYNA 

simulations with different parameter values. For this work the computational framework 

described in [32] and [33] is used based on the commercially available MATLAB1

 program. 

MATLAB scripts are used to modify the LS-DYNA keyword input file automatically, to 

update the parameter values from a known distribution, to control execution of LS-DYNA, 

and to read LS-DYNA output files. Afterwards, results are stored within the MATLAB 

environment to conduct ANOVA and optimization studies. Although comparable response 

surface and variance tools are also available in the LS-Opt [17,34], processing of the time 

history data with the existing MATLAB toolbox provides more freedom. 

As mentioned earlier, a companion set of LS-DYNA solutions while varying 4 parameters 

was created manually, outside the aforementioned computational approach, using a heuristic 

sampling approach to prescribe parameter values, to execute LS-DYNA, and to create a 

solution set. These results are also loaded into MATLAB for analysis and are presented later 

in the paper. 

8
7. Discussion of Results 

Numerical results are discussed in sections to simplify the presentation and to highlight 

specific issues associated with the various aspects of the problem. First, results for response 

surface (RS) models are presented along with examples of LS-DYNA time history predictions 

using the surrogate model. Following this, results using conventional sensitivity analysis are 

presented. In contrast, results from variance analysis using the approach developed by Sobol 

are discussed using both the 4 and 7 parameter surrogate models. Finally, the RS model is 

used to create contours of constant Brinkley values as a function of impact conditions. 

7.1 Response Surface Modeling Results

Before RS models are developed, a database must be created with solutions computed 

using LS-DYNA. For this, two independent efforts took place that produced two databases. 

One database contains results for the case when 4 parameters (1- vertical velocity Vx, 2- 

horizontal velocity Vz, 3- capsule pitch angle θy, and 4- wave slope ψy) are varied manually 

(i.e. heuristic sampling), and the second approach were 7 parameters (1- vertical velocity Vx, 

2- horizontal velocity Vz, 3- capsule pitch angle θy, 4- capsule yaw angle θz, 5- capsule roll 

angle θx, 6- Wave slope ψy, and 7- Wave direction ψx) were varied using a combination of 

sampling techniques which included Halton-leap, D-Optimal, and manual. Each LS-DYNA 

time history output file is read into MATLAB and stored for post-processing. Because the 

initial conditions in the LS-DYNA input key files are set using different logic, the two 

databases cannot be easily combined into one. 

Figure 4 shows results from 231 LS-DYNA runs while varying 7 parameters. To keep the 

amount of data manageable, only 0.2 seconds of simulation time is created and stored, which 

is sufficient to capture the primary impact pulse. At the top plot of Fig. 4a is a summary of 

the time it takes to reach peak acceleration upon impact. This plot is part of the initial 

screening of the database because those cases where the peak time equals 0.2 seconds 

correspond to cases where the capsule is either skimming over the water or cases where the 

horizontal speed is so high that the capsule does not hit the water within 0.2 seconds. Fig. 4b 

shows the resultant maximum acceleration norm in units of gs, and Fig. 4c shows estimates 

for the Brinkley index, computed from predicted local accelerations at the crew location 5. 

9
Cases where Brinkley values are greater or equal to 1 correspond to landing conditions with 

probabilities of injuries. 

Once a database is constructed, RS models are created to predict the time responses as a 

function of the input parameters. Figure 5 shows an example of predictions using Kriging 

response surface to predict velocities; solid blue lines correspond to two LS-DYNA solutions 

and the solid green lines correspond to response surface predictions. Admittedly, this plot is 

only a qualitative comparison of the ability of the surrogate model to predict the time 

responses but it shows that the response prediction is within the LS-DYNA bounds. 

A more quantitative comparison illustrating the adequacy of the RS model is shown in Figs. 

6a and 6b where two LS-DYNA solutions are removed from the 4-parameter database and 

time histories are predicted with the resulting RS model. The parameter vector removed for 

each case is shown with the figure. Figure 6a and 5b shows LS-DYNA results in solid blue 

and RS predictions in dashed green. In Figure 6a the error difference between LS-DYNA and 

the RS model is 0.002; whereas, in Fig. 6b, after removing a different parameter set, the error 

difference of 0.03. Although, these are excellent results for these two cases, as is often the 

case with nonlinear problems, there might be regions of the solution space where the 

interpolation is not as good. 

A note regarding the use of Kriging and/or ERBF is in order. Although both approaches 

provide excellent and comparable results, in cases like our 4-parameter case where sampling 

of the parameter space is heuristic, Kriging has problems when the sampling produces ill 

conditioned (, ) R ν i j ν in Eq.(3.1). Results with the 4-parameter database are obtained 

exclusively using ERBF, but with the 7-parameters database, both Kriging and ERBF provide 

practically identical results. 

7.2 Conventional Sensitivity Analysis of Brinkley Index

A standard approach in design is to perform a parameter study, in which one parameter at a 

time is varied while the other parameters are fixed. This approach is sometimes useful to 

show the sensitivity of responses to a given parameter. However, it is highly dependent on the 

selection of the nominal values. As shown in Fig. 7, the Brinkley index is a non-linear 

function of the relative pitch angle (capsule pitch angle minus wave slope). In addition, the 

shape of the Brinkley index curve is different at different combinations of tangential and 

10
normal velocity. Clearly, a simple polynomial curve fits will not be able to predict the 

Brinkley index for Orion; more advanced response surface methods must be used. 

7.3 Analysis of Variance 

A more general approach to study the sensitivity of nonlinear system responses to 

parameter changes is the approach developed by Sobol. Normally, Sobol indices are 

computed by direct evaluation of the function values (i.e. peak accelerations, Brinkley values, 

etc) as the parameters are varied. For example, if one needs to evaluate the sensitivity of the 

peak acceleration to variations in pitch angle, LS-DYNA would have to be executed multiple 

times as required by Eq. (5.1). This approach is time consuming with models whose 

simulation times are between 12 to 24 hours. Instead, the RS model is used to provide 

predictions for the Sobol formula. 

To aid in understanding the process, Fig. 8 shows a notional chart where the abscissa and 

ordinate contains Brinkley values for two crew locations 1 and 2. Each dot in the chart 

corresponds to a possible LS-DYNA solution when all the parameters are varied. The outer 

blue circle represents to the total variance, whereas the single parameter contribution is 

depicted using the inner black circle. Finally, the ratio of these two variances provides for a 

direct measure of the parameter importance. 

7.4 Four Parameter Variance Analysis

To demonstrate the approach for the capsule impacting the water, the 4-parameter response 

surface model is used to estimate the parameter variance using Eq. (5.1). Figure 9 shows the 

fractional contribution from each parameter variation to the total variance of the Brinkley 

index. This fractional contribution is stacked as a percentage of the total variance 

(normalized to 1) observed at each of the six crew locations. Colors are used to distinguish 

among the different parameters, as indicated by the color-bar. Across the top of Fig. 9 are the 

computed mean values for the Brinkley index at each crew location. After reviewing Fig. 9, it 

is clear that the vertical velocity Vx is the most important parameter followed by the wave 

slope ψx. 

11
7.5 Seven Parameter Variance Analysis 

In contrast to the variance results using the 4-parameter surrogate, Fig. 10 shows the results 

using a 7-parameter model. The results are significantly different from the 4-parameter case 

and show that in this case both vertical and horizontal velocities play a significant role along 

with the wave slope. After further review, it is observed that yaw θz and pitch θy show an 

increased contribution to the variance, when compared to the 4 parameter results, and 

significantly change the capsule behavior on impact. 

To better understand these results, it is important to highlight the difference between the 4-

parameter and 7-parameter databases. Because these were two independent efforts, the LSDYNA initial condition set up is different. The main difference is that in the 4-parameter 

database, the capsule is always placed such that the closest distance from the capsule to the 

water is the same. Hence, during free fall, the capsule always travels the same distance. 

However, in the 7-parameter database this is not the case. Thus results need to be adjusted 

based on impact attitude and speed as opposed to the initial conditions in the LS-DYNA key 

file. 

7.6 Brinkley Contours

Water impact analysis is aimed at estimating the risk of injury to the crew upon landing. 

For this, the Brinkley index is an industry accepted criteria to assess the risk of injury using 

local acceleration responses at the crew location. If only global outputs are available, the 

method discussed in Appendix A to compute the local accelerations should be used. This 

method is used with the 7-parameter database. 

Brinkley contours are created from the local acceleration for the 4-parameter database and 

the 7-parameter database and compared as a function of the normal and tangential velocity on 

the impact plane. By using the normal and tangential velocities, as opposed to the initial 

velocities prescribed in the LS-DYNA key file, results for the two databases can now be 

compared. Figure 11 shows Brinkley contours at crew location 5 obtained with the 4-

parameter database shown in Fig. 11a (146 LS-DYNA runs) and the 7-parameter database 

(239 LS-DYNA runs) shown in Fig. 11b. Although there are some visible differences in the 

contour curvature and Brinkley values near the extreme velocities, as more LS-DYNA 

solutions are added to the 7-parameter database, results should move closer to those obtained 

12
with the 4-parameter data. This is because the 7-parameter RS model requires more LSDYNA solutions to properly capture the behavior of the system when compared to the 4-

parameter model. Unfortunately the question of how many LS-DYNA runs are needed does 

not have a simple answer. At best, users are encouraged to add solutions in those critical 

regions to see if there are significant changes in the results. 

Brinkley contours provide engineers with a portrait of those conditions that could cause 

injuries. To perform this assessment using a conventional approach, one would need to 

compute enough LS-DYNA solutions to populate the contour. However, this approach would 

only give the user a limited view of the functional relationship for the particular combination 

of variables. In contrast, by using the RS models, a more complete portrait of Brinkley values 

can be generated in which all the parameters are considered. For example, Figure 12 shows 

the Brinkley contours at crew location 5 for pitch angles of -28 degrees Fig. 12a (nominal), 

-18 degrees Fig. 12b, and -38 degree Fig. 12c. Aside from the fundamental differences in the 

curvatures, it should be apparent that at -18 degrees pitch there are many normal/tangential 

velocities that result in Brinkley values of 1. In contrast, results for the case when the wave 

slope is varied instead between -30 and 30 degrees (with the pitch angle held at -28 degrees) 

resulted in the Brinkley contours shown in Fig. 13. After examining Figs. 12 and 13, it clear 

that changes in the capsule pitch angles are not equivalent to changes in the wave slope angle. 

For completeness Figs. 12 and 13 show the computed impact angles, as defined in Appendix 

B, along with both pitch and the wave slope angles. Note that for a wave slope of 0 and a 

nominal capsule pitch angle of 28 degrees, the relative impact angle (angle between the water 

plane normal and capsule impact plane normal) is 28 degrees. On the other hand, a 30 

degrees wave slope corresponds to a 58 degrees impact angle. When studying Fig. 13, 

remember that the independent variables are the initial vertical and horizontal velocities and 

not the tangential and normal velocities. Hence, the normal and tangential velocities are 

derived quantities and therefore certain values are outside the range. 

7.7 Estimation of Critical Impact Conditions

A more important use of the RS model is to recover critical landing conditions. In Fig. 14 

is an illustration showing Brinkley contours with a fictitious probability distribution function 

superimposed for the normal velocity. If the probability distribution function is known, then it 

13
is possible to compute the most probable set of normal velocities that would cause Brinkley to 

exceed 1. This kind of information is precisely what the first and second order reliability 

methods provide. 

To demonstrate the use of RS models to estimate critical landing conditions, the 4-

parameter RS model is used with a genetic optimization algorithm to compute critical 

conditions. Because parameter distribution functions are not readily available, the optimizer 

is instructed to look for parameter sets that produce Brinkley values of 1. Table 4 shows one 

parameter set of the many sets computed by the optimizer. To verify that this is indeed a 

critical landing condition, LS-DYNA is executed with the parameters shown in the table and 

the time history from LS-DYNA is plotted in Fig. 15 in solid blue versus the prediction from 

the response surface model in dashed green. Although the RS prediction is not exact, it is very 

good, and more importantly it predicted a Brinkley index near 1. 

8. Concluding Remarks 

The paper discussed an approach that combines solutions from the commercial code LSDYNA with response surface techniques to estimate critical impact conditions of a capsule 

while landing on water. What makes this problem so challenging is the complexity of the 

water and capsule interaction, and the computational time it takes to generate enough 

solutions to understand the physics of the problem. Oftentimes, statistical analysis of 

problems this complex is not undertaken because the computational burden is significant. It is 

shown how to use a limited number of high fidelity LS-DYNA solutions with response 

surface techniques, like Extended Radial Basis Functions and Kriging, to conduct statistical 

analysis and to study conditions under which crew injuries are likely to occur. In addition, the 

approach also allows for sensitivity studies, using an approach developed by Sobol, to screen 

and rank parameters based on their contribution to the output. 

Initial results for the capsule landing on water show that the vertical velocity and the impact 

plane slope are the two most important parameters. Results also showed that although certain 

angle combinations, e.g. capsule pitch and wave slope, produce the same impact angles, the 

system response and Brinkley estimates can be significantly different. 

14
9. Acknowledgements 

The authors would like to thank James M. Corlis for his help during the problem definition 

phase, Dr. John T. Wang for providing the initial LS-DYNA water model that was the basis 

for the model used in the paper, and also special thanks to Gabrielle S. Snyder for her help 

setting up the computer cluster for the LS-DYNA runs. 

15
10. References 

1. Hartwig, H.H and Schmidt, J. L.: “Summary– Apollo Earth Landing System Impact 

Attenuation Trade-Off Study,” North American Aviation SID 62-897, July 1962. 

2. Stubbs, S.M.: “Investigation of the Skid-Rocker Landing Characteristics of Spacecraft 

Models," NASA TN D-1624, April 1963. 

3. Benson, H.E. and West, R.B.: “Results of Two One-Quarter Scale Apollo Model 

Impact Tests Utilizing Different Impact Attenuation Systems,” NASA TMX-65129, 

January 1965. 

4. “Water Impact of the Mercury Capsule Correlation of Analysis with NASA Tests,” 

Grumman Aircraft Engineering Corporation, NASA CR-55251, December 1962. 

5. Thompson, W.C.: “Dynamic Model Investigation of the Rough-Water Landing 

Characteristics of a Spacecraft,” NASA TN D-3774, April 1967 (Technical Film 

Supplement L-940). 

6. Reese, T. and Rose, J.D.: “Apollo Pad Abort Land Impact Tests at Kennedy Space 

Center,” NASA TM X-64365, April 1969. 

7. McCullough, J. E. and Lands, J. F. Jr.: “Apollo Command Module Land-Impact 

Tests,” NASA TN D-6979, October 1972. 

8. Baker, W. E., Westine, P. S., Garze, L. R., Hunter, P. A.: “Water Impact Studies of 

Model Apollo Command Module,” NASA CR-114949, August 1965. 

9. Baker, W.E and Westine, P.S.; “Model Tests for Structural Response of Apollo 

Command Module to Water Impact,” AIAA Journal of Spacecraft and Rockets, 0022-

4650, Vol. 4 No. 2, 1967. 

10. Stubbs, S.M.: “Dynamic Model Investigation of Water Pressures and Accelerations 

Encountered During Landings of the Apollo Spacecraft,” NASA TN D-3980, 

September 1967. (Technical Film Supplement L-960). 

11. Stubbs, S. M Hathaway, M. E.: “Effects of Bottom-Structure Flexibility on Water 

Landing Loads of Apollo Spacecraft Models,” NASA TN D-5108, March 1969. 

12. “Final Test Report Water Impact Test 92 Boilerplate 28 (Drop 3)”, North American 

Aviation SID 66-35, April 1967. 

13. “Final Test Report Water Impact Test 104 Spacecraft 2S-1 (Drop 2)”, North 

American Aviation SID 66-330A, April 1967. 

14. “Study of Apollo Water Impact Final Report: Volume 4 – Comparison with 

Experiments,” North American Aviation, NASA CR-92022, May 1967. 

15. “Apollo Applications Program Preliminary Evaluation Data – Apollo Landing and 

Recovery Subsystem Analyses,” North American Rockwell, SID-69-244, April 1969. 

16. Brinkley, J.W., and Moser, S.E.: "Development of Acceleration Exposure Limits for 

Advanced Escape Systems," AGARD CP-472, April, 1989. 

17. Anonymous: “LS-DYNA Keyword User’s Manual,” Version 971, Livermore 

Software Technology Company, Livermore, CA, August 2006. 

18. Sacks, J., Schiller, S.B, and Welch, W.J.: “Design of Computer Experiments,” 

Technometrics, February 1989, Vol. 31, No. 1 

19. Mullur, A. and Messac, A.: “Extended Radial Basis Functions: More Flexible and 

Effective Metamodeling.” AIAA Journal, Vol., 43, No. 6, June 2005. 

16
20. Mullur, A. and Messac, A.: “Metamodeling Using Extended Radial Basis Functions: 

A Comparative Approach.” Engineering with Computers (2006) 21: 203-217 DOI 

10.1007/s00366-0005-7 

21. Chipperfield, A., Fleming, P. J., Pohlheim, H. and Fonseca, C. M.: “Genetic 

Algorithm Toolbox for use with Matlab.” Technical Report No. 512, Department of 

Automatic Control and Systems Engineering, University of Sheffield, 1994 

22. Iman, R.J., and Shortencarier, M.J.: "A FORTRAN77 Program and User's Guide for 

Generation of Latin Hypercube and Random Samples for Use with Computer 

Models," NUREG/CR-3624, SAND83-2365, Sandia National Laboratories, 

Albuquerque, NM. 

23. Kalagnanam. J.R, and Diwekar, U.M.:" An efficient Sampling Technique for Off-line 

Quality Control." Technometrics, August 1997, Vol. 39, No. 3, pp 308-319. 

24. Montgomery, D.C: Design and Analysis of Experiments, 6th Edition, Copyright © 

2005 John Wiley and Sons, Inc. 

25. Halton, J.H., “On the Efficiency of Certain Quasi-Random Sequences of Points in 

Evaluating Multi-Dimensional Integrals,” Numerische Mathematik, Vol. 2, pp. 84-90, 

1960 

26. Robinson, D., and Atcitty, C., “Comparison of Quasi-and Pseudo-Monte Carlo 

Sampling for Reliability and Uncertainty Analysis,” American Institute of Aeronautics 

and Astronautics, A99-24882, and AIAA-99-1589. 

27. Sobol, I.M., Tarantola, S., Gatelli, D., Kucherenko, S.S, and Mauntz, W., "Estimating 

Approximation Error When Fixing Unessential Factors in Global Sensitivity 

Analysis," Reliability Engineering and System Safety 92 (2007) 957-960, © 2006 

Elsevier LTD 

28. Mullershon, H., and Liebsher, M, “Statistics and Non-Linear Sensitivity Analysis with 

LS-OPT and DSPEX.” 10th International LS-DYNA Users Conference, pp. 4-1,4-13 

29. Homma, T, and Saltelli, A, “Importance Measures in Global Sensitivity Analysis of 

Nonlinear Models.” Reliability Engineering and System Safety 52 (1996) 1-17, © 

1996 Elsevier LTD. 

30. Sudret, B, "Global Sensitivity Analysis Using Polynomial Chaos Expansion," 

Reliability Engineering and System Safety 93 (2008) 964-979, © 2007 Elsevier LTD 

31. Stander, N. and Goel, T.: “Metamodel Sensitivity to Sampling Strategies: A 

Crashworthiness Design Study.” 10th International LS-DYNA Users Conference, pp. 

4-154-29. 

32. Horta, L.G., Jackson, and K.E., Kellas, S: "A Computational Approach for Model 

Update of an LS-DYNA Energy Absorbing Cell." Proceedings of the 2008 AHS 64th

Annual Forum & Technology Display, April 29-May 1, Montreal Canada. 

33. Horta, L.G. and Reaves, M.C.:" On the Application of a Response Surface Technique 

to Analyze Roll-Over Stability of Capsules With Air Bags Using LS-DYNA." NASA 

TM 2008-215117, March 2008. 

34. Stander, N, Eggleston, T, Craig, K, and Roux, W: “LS-OPT User’s Manual A Design 

Optimization and Probabilistic Analysis Tool for the Engineering Analyst.” Version 2, 

Oct. 2003. 

17
11. Appendix A: Estimation of Response for Arbitrary points on Rigid Body 

With the complexity associated with LS-DYNA models, particularly when modeling 

fluids/structure interaction, the amount of data that is created can be overwhelming. For 

storage purposes, it is more efficient to store only a subset of the response data if other data 

can be recovered from it. The ability to recover information from a subset of data is also 

important in those cases where engineers, after completing a set of simulations, realized that 

they should have monitored other points not included in the stored output. Of course, one 

option is to re-run all simulations to store the data for the points left out, but a more efficient 

way, for rigid structures, is to store the least amount of information and then compute all other 

information from the stored data. This approach is, in fact, possible and more efficient 

provided that the initial information stored is judiciously selected to allow for general 3-

dimensional motion to be recovered. In the following, the reader is reminded of the equations 

to recover the responses at arbitrary points. 

Recall that from rigid body kinematics, the relative velocity of point c with respect to point 

m is given as 

c m cm/ / c m vv v R − = =× ω

G G G G G (11.1)

where ω

G is the angular velocity of the body and Rc m/

G is the relative position of c with respect 

to m. There are at least two ways that Eq. (11.1) can be used; 1) if the velocity of point c and 

m are known, the angular velocity of the body can be estimated knowing the body’s geometry, 

and/or 2) if the information for an arbitrary point d located at Rd m/

G is needed, given the body 

angular velocity, the relative velocity of d is also given by Eq. (11.1) . 

Similarly, if acceleration information is needed, the acceleration of point c with respect to 

m is given as 

Ac m cm m cm c = + = +× +× AA A v R / / ω α /m

GGG G G G G G (11.2)

where now the body angular acceleration vector α

G is introduced. As before, depending upon 

what information is stored with the rigid body simulation, accelerations for points anywhere 

in the body are easily computed using Eq. (11.2). 

18
For results in the present paper, the angular position, angular velocity, and angular 

accelerations are all available. Hence, Eqs. (11.1) and (11.2) are used with the location of 

new points for estimation of acceleration values. 

Although response extraction points on the body may seem futile, for cases where the 

location of the point at which the response is uncertain, this method provides a simple 

mechanism to examine the variation in the system response to establish uncertain bounds. 

Otherwise, new simulations would need to be executed when evaluating new points. 

19
12. Appendix B: Estimation of Impact Velocities and Angles 

In this study, many response quantities are computed directly within LS-DYNA but there 

are others like impact angles and velocities that must be derived from the LS-DYNA outputs. 

In the following, the equations used to compute the impact angles and velocities are provided. 

12.1 Estimation of Impact Angles

Consider, the capsule shown in Figure B.1. Assume that the capsule motion starts with an 

initial velocity and pitch orientation, falls under a gravity and impacts a water surface that is 

defined by a point Rp

G on the surface and its unit normal p n

G . Also consider a capsule with a 

moving reference frame attach to point m, which is located at the intersection of the capsule's 

center line and an imaginary plane containing the impact point IP. In the following, relative 

vector definitions will have a subscript “a/b” to indicate that it is a vector pointing from point 

b to point a, i.e. making the point b the reference. Vectors without an explicit reference are 

defined with respect to the inertial frame. Finally, if a vector quantity is used without the 

upper arrow, it is referring to the vector of components, for example / [ ] a b R = xyz ' . 

To estimate the impact angles of a capsule with the water, first the impact point IP on the 

water must be located at the intersection of the capsule and water planes. With the aid of Fig. 

B.1 let an impact plane on the capsule be defined by its normal c n

G and a center point m

located at Rm

G . As mentioned earlier, the water impact plane is defined by its normal p n

G and a 

point p located at Rp

G

c

. Now the vector difference between point c on the capsule impact 

plane (located at R

G ) and RmG must be orthogonal to c n

G , hence; 

( ) T

cc m nR R− = 0 (12.1) 

Similarly, for an arbitrary point on the water plane IP located at RIP

G , the vector difference 

between it and Rp

G must orthogonal to p n

G and thus satisfies the equation 

( ) T

p IP p nR R− = 0 (12.2)

The plane equations (12.1) and (12.2) are written in terms of vector differences of two 

arbitrary points on each of the respective planes. In addition to these two equations, the 

20
impact point, if it exists, is located along the direction B c p = n n × G G G , i.e., orthogonal to both 

normal vectors. For cases where both planes (capsule and water) are rectangular, B

G

points 

along the intersection line, but, if one of the impact planes is circular (which is our case), the 

intersection line reduces to a single point. 

Because the capsule impact plane normal c n

G is known or can be estimated from the capsule 

displacements, this normal can be used to express any point on the capsule impact plane as 

c m f =Uq R R = − (12.3)

where the matrix U is an 3 2x R constructed using the basis vectors from null space of , i.e., 

any orthogonal vector set perpendicular to ; whereas q is an 

c n

c n

2 1x R vector of arbitrary 

parameters to be computed later. Using the same logic, any point on the impact plane can be 

written as 

w Gz R R = = − IP p (12.4)

 where the matrix G is 3 2x R constructed from the null space of and z is an p n

2 1x R vector of 

arbitrary parameters. In Eqs. (12.3) and (12.4) each of the arbitrary vectors q and z has

dimensions reduced from 3 to 2. 

Since RIP

G is the impact point, the following vector summation must hold 

Rm IP + = f R (12.5)

and 

Rp I + = w R P (12.6)

Because the impact point has no projection along B

G

, this fact is represented as 

( ) T BR f m + = 0 (12.7)

or equivalently, 

( ) T BR w p + = 0 (12.8)

Combining Eqs. (12.3), (12.4), (12.5), (12.6), and (12.8), yields the following set of 4 

equations and four unknowns 

0

m p

T T

p

U G q R R

B G z B R

⎡ ⎤ ⎧ ⎫ ⎧⎪− + ⎫⎪ ⎨ ⎬⎨ ⎬ = ⎢ ⎥ ⎣ ⎦ ⎩ ⎭ ⎪⎩ ⎭ − ⎪ (12.9)

21
Solving for q and z, and substituting back into (12.5)or (12.6) provides the coordinates of the 

impact point RIP

G . 

Once the impact point location is known, the impact angles can be easily computed. 

Impact angleθ , which is independent of the impact point, is computed using the normal 

vectors as 

(12.10) 1 cos ( ) T

c p θ n n − =

Impact angle ψ is the angle between a traveling direction Rc m/

G and the impact point RIP m/G

(note that both vectors are now defined with respect to m). Using these definitions the angle 

is defined as 

1 / /

/ /

cos ( )

T

IP m c m

IP m c m

R R

R R ψ − = (12.11)

 (12.12) 

 As is always the case when defining angles, there is arbitrariness in the way the normals 

are defined and, as such, these definitions must be adjusted to conform to a particular sign 

convention. 

22
12.2 Estimation of Impact Velocities

Impact plane location and orientation is provided in terms of a point on the plane and its 

unit normal vector, but to compute impact velocities based on the capsule orientation and 

direction of travel, it is necessary to define the impact plane in terms of a local coordinate 

system. Using the impact plane normal unit vector p n

G a set of orthogonal axes can be defined 

in terms of the null space of . If the two bases vectors for the null space are represented as 

, one can use the vector

p n

G

1 2 Ggg = [ ] B c p = n n × G G G , which is normal to both planes, as one of 

the three components of the impact plane coordinate system. If the projection of the vector B

G

along and is given as 1 g 2 g 1 1 α = g D B and (where D means dot product), the 

traveling direction is defined as that for which 

2 2 α = g D B

αi is equal to zero. Since our bases vectors are 

not necessarily aligned with any particular axis, one can compute the component along the 

tangential direction (traveling direction) as one of these three vectors; 

m 

ψ

x 

z 

y 

Rp

G

Rm

G

Impact 

point IP 

c n

G

p n

Capsule Impact Plane G

Water Impact Plane 

p 

θ

C

Fig. B.1 Capsule and water plane parameter description 

23
1 1

1 1

1

1

1 1

1

2 2

2 1

2

( ) , 0, 1

, 0,

( ) , 0,

g Bg l g if o g B

g l if

g

g Bg l g if g B

α α

α α

α

= − ≠ ≠

= ≠

= − =

G D

G

G D

1 ,

1,

r

= or

w nl =− ×p

G G G

p ⎡ ⎤ −n lw ⎣ ⎦

 (12.13)

Using as the third component of the coordinate system, the coordinate system is 

now defined by the triad G G G . The negative sign on the impact plane normal is 

there to accommodate a particular sign convention. 

Finally, using these direction cosines, velocities along the normal and tangential 

components are computed using the dot product of the capsule velocities on impact and the 

impact plane triad. 

24
Table 1 Mass and inertial properties of capsule 

Mass, lbf s

2

 / in IXX, lbf s

2

 in IYY, lbf s

2

 in IZZ, lbf s

2

 in 

49.0536 216,252 136,104 177,684 

Table 2. Parameter definitions 

Parameter 

Number 

Parameter 

Name 

Upper 

Bound 

Nominal 

Lower 

Bound 

1 Vertical Speed 

Vx (in/sec) 

600 360 120 

2 Horizontal 

Speed Vz (in/sec) 

0 -480 -960 

3 Pitch Rotation 

θy (deg) 

-18 -28 -38 

4 Yaw Rotation 

θz (deg) 

10 0 -10 

5 Roll Rotation 

θx (deg) 

30 0 -30 

6 Wave Slope 

ψx (deg) 

30 0 -30 

7 Wave Roll 

ψy (deg) 

30 0 -30 

25
Table 3. Solid element dimensions for 4 mesh cases 

Mesh 

 Case 

Far-Field 

Air, in. 

Interface 

 Air, in. 

Interface 

 Water, in. 

Far-Field 

 Water, in. 

1 4 4 4 4 

2 16 16 4 4 

3 4 1 1 4 

4 4 4 1 4 

Table 4. Example of estimated critical landing parameters 

Parameter Nominal 

Critical 

Landing 

Condition 

Vertical Speed 

Vx (in/sec) 

360 501.9 

Horizontal Speed 

Vz (in/sec) 

-480 -47.4 

Pitch Rotation 

θy (deg) 

-28 -26 

Yaw Rotation 

θz (deg) 

0 0 

Roll Rotation 

θx (deg) 

0 0 

Wave Slope 

ψx (deg) 

0 -5 

Wave Roll 

ψy (deg) 

0 0 

26
 

Model 1 2 

Width, ft. (Y) 50.0 50.0

Length, ft. (Z) 50.0 75.0

Water Depth, ft. (X) 19.0 19.0

Air Depth, ft. (X) 14.0 14.3

Air Water/Air Dimensions 

Water 

z 

x 

Fig. 1 Model dimensions. 

Uphill

(Positive Slope Angle)

Heel-in

(Positive Pitch 

Angle)

Direction of Forward Velocity 

No 

Pitch

Downhill

(Negative Slope)

Toe-in

(Negative Pitch +Z 

+Y 

 3 4 

5 6 

 1 2

Direction of 

Forward Velocity

Fig. 2 Capsule and water orientation definitions with approximate crew location seating. 

27
0.00

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

0.50

0.00 0.02 0.04 0.06 0.08 0.10 0.12 0.14 0.16 0.18 0.20

Time, s

Case 1

Case 2

Case 3

Case 4

Brinkley 

Index 

Fig. 3 Mesh sensitivity results using Brinkley Index. 

Time 

To 

Peak 

Value 

(sec) 

4a) time to peak acceleration

Norm 

of 

Max 

Accel

(gs)

28

4b) norm of resultant acceleration

Crew Location 5 

Brinkley 

Index 

LS-DYNA Run Number 

4c) peak Brinkley value 

Fig. 4 Summary of LS-DYNA response for 7 parameter database. 
29

400

0 0.5 1 1.5 2 2.5 3

0

0.5

1

-300

-200

-100

0

100

200

300 LSDYNA

Solution

/s)

Time (sec)

bdaVel-x (in

Velocity 

Estrapolated

Solutions Interpolated

Solutions 

λ

Time (sec)

x 1 2 = +− λx x (1 ) λ

Parameter Variation 

Fig. 5 Example of response surface predictions of LS-DYNA time. 
30

0 0.05 0.1 0.15 0.2

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Time (sec)

Brinkley

VL 

Crew Location 4 Predictions

Estimate

Exact

LS-DYNA 

RS 

b) Parameters removed P=[540 in/s,-480 in/s,-28 deg,0 deg] 

0 0.05 0.1 0.15 0.2

0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

Time (sec)

B

Crew Location 4 Predictions

Estimate

Exact

rinkley

VL 

LS-DYNA 

RS 

Prediction Error=0.002 

Prediction Error=0.03 

Fig. 6 Verification of response surface and LS-DYNA time history predictions at 

crew location 4 after removing solution from the response surface. fitting 

a) Parameters removed P=[360 in/s,-480 in/s,-28 deg,0 deg] 

Brinkley 

Index 

Brinkley 

Index 
0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

-50 -40 -30 -20 -10 0 10 20

R

Vtan=40 ft/s, Vnorm=30 ft/s

Vtan=25 ft/s, Vnorm=43.3 ft/s

Vtan=80 ft/s, Vnorm=50 ft/s

31

210.8 210.9 211 211.1 211.2 211.3 211.4 211.5 211.6 211.7 211.8 -48.45

-48.4

-48.35

-48.3

-48.25

-48.2

-48.15

-48.1

-48.05

-48

-47.95

Max Accel. of Node 1

Total 

Variance

Single 

Parameter 

Variance 

Brinkley 

Index Crew 

Location 1 

Fig. 8 Example of analysis of variance. 

Brinkley Index Crew Location 2 

elative Pitch Angle (deg)

L BI ow

Brinkley 

Index 

Relative Pitch Angle (deg) 

Fig. 7 Single variable sensitivity analysis. 
1 2 3 4 5 6

0

0.2

0.4

0.6

0.8

1

O

Mean Values Brinkley Index Mean Values 

32

utput Number

Varian o ution ce C ntrib

Vx (in/s) 

Vz (in/s) 

Pitch (deg) 

0.470.380.400.380.470.40

W Slope (deg) 

Variance 

Contribution 

Crew Location 

Fig. 9 Variance analysis of Brinkley index using the 4 parameter RS model. 

Fig. 10 Variance of Brinkley index using the 7 parameter ERBF model. 

Brinkley Index Mean Values Mean Values

1 2 3 4 5 6

0

0.2

0.4

0.6

0.8

1

Output Number

Varian tion ce Contribu

Vx (in/s) 

Vz (in/s) 

Pitch (deg) 

Yaw (deg) 

Roll (deg) 

W Slope (deg) 

W Roll (deg) 

0.340.290.300.290.340.30

Variance 

Contribution 

Crew Location 
0.1

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

0.4

0.4

0.4

0.5

0.5

0.5

0.6

0.6

0.6

0.7

0.7

0.7

0.8

0.9

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

200 300 400 500 600

-900

-800

-700

-600

-500

-400

-300

-200

-100

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.2

0.2

0.2

0.3

0.3

0.3

0.3

0.4

0.4

0.4

0.5

0.5

0.5

0.6

0.6

0.7

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

33

200 300 400 500 600

-900

-800

-700

-600

-500

-400

-300

-200

-100

0

0.1

0.2

0.3

0.4

0.5

0.6

a) 4-Parameter database and 146 LS-DYNA runs 

0.7

b) 7-Parameter database and 239 LS-DYNA runs 

Fig. 11 Comparison of contour plots for the Brinkley index using the 4-parameter 

and a 7-parameter database. 
34

0.1

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

0.4

0.4

0.4 0.5

0.5

0.5

0.6

0.6

0.6

0.7

0.7

0.7

0.8

0.9

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

200 300 400 500 600

-900

-800

-700

-600

-500

-400

-300

-200

-100

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

0.4

0.4 0.4

0.5

0.5 0.5 0.6

0.6 0.6

0.7

0.7

0.7

0.8

0.8

0.8

0.9

0.9

0.9

1

1

1

1.1

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

200 300 400 500 600

-900

-800

-700

-600

-500

-400

-300

-200

-100

0

0.2

0.4

0.6

0.8

1

1.2

0.1

0.1

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

0.3

0.4

0.4

0.4

0.4

0.5

0.5

0.5

0.5

0.6

0.6

0.7

0.8

0.9

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

200 300 400 500 600

-900

-800

-700

-600

-500

-400

-300

-200

-100

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

a) -28 deg. Pitch, 

Impact angle θ=28 deg.

b) -18 deg. Pitch, 

Impact angle θ=18 deg.

c) -38 deg. Pitch, 

Impact angle θ=38 deg.

Fig. 12 Contour plots for Brinkley Index as a function of the impact tangential 

and normal velocity at the crew location 5. 
35

0.1

0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

0.4

0.4

0.4 0.5

0.5

0.5

0.6

0.6

0.6

0.7

0.7

0.7

0.8

0.9

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

200 300 400 500 600

-900

-800

-700

-600

-500

-400

-300

-200

-100

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

0.2

0.3

0.4

0.4

0.5

0.5

0.5

0.6

0.6

0.7

0.7

0.8

0.8

0.9 1

1.1

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

200 400 600 800

-700

-600

-500

-400

-300

-200

-100

0

100

200

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

0.2

0.2

0.2

0.2

0.2

0.2

0.4

0.4

0.6

0.6

0.8

1

1.2

1.4

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

-200 0 200 400

-1100

-1000

-900

-800

-700

-600

-500

-400

-300

-200

-100

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

a) 0 deg. Wave Slope, 

Impact angle θ=28 deg

b) 30 deg. Wave Slope, 

Impact angle θ=58 deg

c) -30 deg. Wave Slope, 

Impact angle θ=2 deg

Fig. 13 Contour plots for Brinkley index as a function of the impact tangential and 

normal velocity at the crew location 5. 
0.1

0.1

0.2

0.2

0.2

0.3

0.3

0.3

0.4

0.4

0.4

0.5

0.5

0.5

0.6

0.6

0.6

0.7

0.7

0.7

0.8

0.8

0.8

0.9

0.9

0.9

1

1

1

1.1

Normal Vel. (in/sec) 

Tangential Vel (in/sec)

Crew Location= 1 

0 0.05 0.1 0.15 0.2

0

0.2

0.4

0.6

0.8

1

1.2

1.4

Time (sec)

B

VL 

Crew Location 4 Predictions

rinkley

Estimate

Exact

LS-DYNA 

RS 

Fig. 15 Verification of predicted critical landing condition using the 

response surface model. 

Brinkley 

Index 

Fig. 14 Illustration of process to estimate most likely critical condition. 

200 300 400 500 600

-900

-800

-700

-600

-500

-400

-300

-200

-100

0 1.2

0.2

0.4

0.6

0.8

1

Nominal 

36
REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

2. REPORT TYPE 

Technical Memorandum

 4. TITLE AND SUBTITLE

A Computational Approach for Probabilistic Analysis of Water Impact 

Simulations 

5a. CONTRACT NUMBER

 6. AUTHOR(S)

Horta, Lucas G.; Mason, Brian H.; Lyle, Karen H.

 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-2199

 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

 8. PERFORMING ORGANIZATION

 REPORT NUMBER

L-19585

10. SPONSOR/MONITOR'S ACRONYM(S)

NASA

13. SUPPLEMENTARY NOTES

12. DISTRIBUTION/AVAILABILITY STATEMENT

Unclassified - Unlimited

Subject Category 18

Availability: NASA CASI (443) 757-5802

19a. NAME OF RESPONSIBLE PERSON

STI Help Desk (email: help@sti.nasa.gov)

14. ABSTRACT

NASA's development of new concepts for the Crew Exploration Vehicle Orion presents many similar challenges to those 

worked in the sixties during the Apollo program. However, with improved modeling capabilities, new challenges arise. For 

example, the use of the commercial code LS-DYNA, although widely used and accepted in the technical community, often 

involves high-dimensional, time consuming, and computationally intensive simulations. The challenge is to capture what is 

learned from a limited number of LS-DYNA simulations to develop models that allow users to conduct interpolation of 

solutions at a fraction of the computational time. This paper presents a description of the LS-DYNA model, a brief summary 

of the response surface techniques, the analysis of variance approach used in the sensitivity studies, equations used to estimate 

impact parameters, results showing conditions that might cause injuries, and concluding remarks. 

15. SUBJECT TERMS

Impact dynamics; Response surface; Uncertainty analysis; Water landing

18. NUMBER

 OF 

 PAGES

41

19b. TELEPHONE NUMBER (Include area code)

(443) 757-5802

a. REPORT

U

c. THIS PAGE

U

b. ABSTRACT

U

17. LIMITATION OF 

 ABSTRACT

UU

Prescribed by ANSI Std. Z39.18

Standard Form 298 (Rev. 8-98)

3. DATES COVERED (From - To)

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

5e. TASK NUMBER

5f. WORK UNIT NUMBER

644423.06.31.03.14.04.01

11. SPONSOR/MONITOR'S REPORT

 NUMBER(S)

NASA/TM-2009-215704

16. SECURITY CLASSIFICATION OF:

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing this burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and 

Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person 

shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

01- 04 - 2009