Random Pressure Loads

16. Random Pressure Loads#

In a previous section we discussed random vibration input (see section 9). That section addresses a loading where the frequency content (or power spectral density) of the loading is known for a few points on the structure. In contrast, for hypersonic vehicles a random loading may occur at thousands of points on the surface. Many aspects of the loading are the same, but the specification is different, and for performance reasons, the solutions are performed differently.

The starting point for this analysis requires the following.

A surface sideset where the loading will be applied.
A temporal correlation function to apply on the surface. The temporal correlation function is the inverse Fourier transform of a power spectral density (PSD).
A spatial correlation relation. Currently, that relation may only be specified as a pair of exponential decay constants.

Details of the problem setup may be found in the User’s Manual. This section provides a simple example of the setup and an informal discussion of the sources of the data.

16.1. Example Problem Set-up#

For our example, we consider a cylinder in a flow field as shown in Figure 16.1 with the full input attached. The structure is a right circular cylinder of diameter 1 unit, and height 2 units. The flow is directed towards this cylinder in the \(X\) direction, and the PSD and corresponding temporal correlation function are shown in Figure 16.2.

Figure 16.1 Example Random Pressure Geometry#

Figure 16.2 Example Random Pressure PSD and Correlation Functions.#

We are interested in this example, in frequencies up to 500 Hz, so the cutoff frequency is 500 Hz. There is no point in adding energy above the desired cutoff frequency – it only complicates the procedure. [1] The PSD of the input thus controls much of the solution.

The spatial correlation is often more difficult to obtain. For our example, we require a decay constant of 2.0 units in the flow direction, and 5.0 perpendicular to the flow. One can think of corresponding decay distances of 0.5 and 0.2 respectively. Thus, down the flow, points more than about 1.5 units away will not be well correlated. [2] Perpendicular to the flow, correlations decay even faster.

One rarely has much experimental data about the spatial correlation. Some information is sometimes garnered from the temporal correlation. For example, if the correlation function has a characteristic time, \(\tau\), one would expect the spatial correlation length to be of the order of \(\delta = v\tau\), where \(v\) is the flow velocity. For a structure in a fluid, the dimensions of the turbulent layer also provide a bound on the spatial correlation.

16.2. Example: Input Specifications#

The physical quantities of the previous section can be interpreted and expressed as Sierra/SD input as follows.

The temporal correlation function of Figure 16.2 can be digitized as a Sierra/SD “time only” function. In input 16.1 we use a triangular pulse for simplicity. The correlation function should be symmetric about the origin, and it should have the value of 1.0 at the origin. The correlation_function is used in the load section, as shown in input 16.2.

In the loads section, we also define the following quantities.

cutoff_freq = 500
coordinate = 1 # to set flow direction
ntimes = 5 # varies from 3 to 20. Too small causes poor replication
           # of the temporal correlation function. Too large results
           # in ill conditioning and singularity.

Recall that the full correlation matrix is a tensor product of the spatial correlation with temporal components. The “NTimes” parameter controls the number of samples in the time domain.

All that remains is setting the spatial correlation decay constants in the loads section. The full text is shown in input 16.2.

correlation_length_z = 0.5
correlation_length_r = 0.2

Listing 16.1 Random Pressure Correlation Function#

FUNCTION 1
  type linear
    data   -0.001909859319285   0
    data   0   1
    data   0.0019098593192856   0
END

Listing 16.2 Random Pressure Load Section. Note that the “flow” direction is the \(Z\) coordinate direction of coordinate frame 1.#

LOADS
  sideset 1
    randompressure
    cutoff_freq = 500
    delta_t = 0.001
    correlation_function = 1
    ntimes   = 5
    correlation_length_z = 0.5
    correlation_length_r = 0.2
    coordinate = 1
END
BEGIN RECTANGULAR COORDINATE SYSTEM 1
  origin 0 0 0
  z point 1 0 0
  xz point 1 0 1
END

16.3. Example: Verifying the Load#

This is a fairly complex input, and it is advisable to verify the generated loads to ensure consistency. We examine four quantities.

average force on a node.
variance of the force on a node.
temporal force correlation on a single node.
cross correlation of forces between nodes.

All these quantities require output of the total input force, which is obtained by specifying force in the OUTPUTS section of the Sierra/SD input. We will use MATLAB tools to evaluate many of the results. Data can be read into MATLAB from the Exodus results using exo2mat or using other methods.

16.3.1. Average Nodal Force#

The average nodal force may be determined either by evaluating the MATLAB results directly, or using the statistics output from Sierra/SD. The built in statistical output is easiest. Supply the statistics keyword to the OUTPUTS section, and results will be written to an additional Exodus file. This has the added benefit that these results may be easily visualized using Paraview or Ensight. See Figure 16.3.

For long time integration, the average value of the nodal force should approach zero. Shorter time samples will have greater variation. The random variables depend on cutoff_freq. The number of random samples can be computed as,

\[ N_{samples} = Time_{analysis} \cdot \mbox{cutoff\_freq} \]

The fractional mean of the force should be within about \(3/\sqrt{N_{samples}}\). Or,

\[ Error_{mean} = \left|\frac{F_{mean}}{F_o}\right| < \frac{3}{\sqrt{N_{samples}}} \]

Here \(F_o\) is the force applied for a correlation function of 1. It involves the scale factors of the function, the sideset distribution factors and the effective area for each node. [3] See the comments section, 16.4 for, discussion on the effective area.

For the example in Figure 16.3, mean forces are of the order of 1/1000. In this example, we took 10,000 time steps, with each of 0.1 ms for a total time \(Time_{analysis}=1\,s\). With \(Delta\_T=1/\mbox{cutoff\_freq}=1\,ms\), the total number of random samples is \(N_{samples} = 1000\).

For nodes in the center of the loading area, the effective area is about 0.0098 square units. Because the sideset distribution factors are all one, we have \(F_o=0.0098\). Then,

\[ Error_{mean} = \frac{0.001}{0.0098} = 0.1 \]

which is greater than \(\frac{3}{\sqrt{1000}}= 0.095\). A distribution of the mean is shown in Figure 16.4.

Figure 16.3 Variation of Mean and Standard Deviation of Force Magnitude on the Surface.#

Figure 16.4 Distribution of Mean Forces on Surface.#

16.3.2. Variance of Nodal Force#

The standard deviation, which is the square root of the variance, is also available as an output from the analysis, and may be plotted on the structure using standard visualization tools. See Figure 16.3.

Again, the standard deviation is a statistical quantity, which is only meaningful for large numbers of samples. In the limit of large \(N\), the standard deviation should approach \(F_o\), as defined above, provided that the correlation function is 1 at time 0.

The plots show a value of \(F_{std} \approx 0.0085\) which is under the expected value of 0.0098. Because the averaging process tends to round out the correlation function, the measured values of the standard deviation are typically somewhat less than \(F_o\). The autocorrelation function analysis of the following section should make this more clear.

16.3.3. Temporal Nodal Force Autocorrelation#

The statistics of the loading on a single node should recover the initial input temporal correlation. Figure 16.5 shows the correlation function extracted from the raw time response data. The correlation function may be computed as,

\[ f_c(n) = \frac{1}{F_o^2} \sum_i w_i w_{i-n}. \]

Where \(w_i\) is the force on a node at time \(t_i\). This data can only be obtained using MATLAB or another external tool, i.e. it is not available as part of the statistical output. In MATLAB we get this with, C = xcorr(f1,f1), where f1 is the force time history on a node of the surface. We recover a correlation that is similar to the original triangle correlation in the input. Because of interpolation and finite sample length, we do not expect the same curves precisely.

The curves of Figure 16.5 should be considered “good enough” in a statistical sense. A temporal interpolation from multiples of Delta_T to the integration time step is being performed, which smooths the values of the correlation.

Figure 16.5 Nodal Force Autocorrelation.#

16.3.4. Spatial Cross Correlation#

The previous section discussed the autocorrelation function, i.e. the temporal correlation of signals on the same spatial location. We can also examine the cross correlation functions. We will only evaluate the functions at the peak.

This is more difficult. We use the MATLAB “find” method to get the indices of the nodes with \(x=-0.5\), and \(y=0\). We loop through these nodes, and compute the “xcorr” function between the node at the center and the other nodes. The peak value of this solution is then plotted versus the distance in Figure 16.6.

Figure 16.6 Nodal Force Spatial Cross Correlation.#

There are obvious differences between the measured loads and the target. The correlations for close distances are lower. This is understood to be generated by the temporal interpolation function. At large distances, the cross correlations never go to zero because of the finite length of the sample.

16.4. Random Pressure Comments#

16.4.1. Effective Area#

Random pressures are computed as force loads using a consistent pressure calculation. Pressures at the nodes are spread through the element shape functions to result in nodal forces. For a uniform mesh, this is similar to lumping the pressures from a fixed area onto the nodes with \(F=P\cdot Area\). In Figure 16.7 an element based mesh is shown along a corresponding effective area for the nodes. For a uniform quadrilateral mesh like the example above, the nodal effective area is the same as the area of an element face.

16.4.2. Temporal Interpolation#

To improve performance, the random pressure loading procedure computes random pressures at multiples of “Delta_T” and then interpolates to integration time steps. A piecewise linear interpolation introduces unacceptable errors; sinc interpolation is much better.

Interpolation can be avoided by choosing the integration and sampling times to be equal. In no case should the integration time be larger than the sampling time.

16.4.3. Singularities#

To compute the proper temporal and spatial correlations for the forces, we need to perform a Cholesky factorization of the correlation matrix. This factor will fail if the matrix is singular. Remember that the correlation matrix that we factor is a tensor product of temporal and spatial components, \(C=C_{spatial} \otimes C_{temporal}\). If either component is singular, the matrix \(C\) is singular. Several common issues can cause singularity of this matrix.

Taking NTimes too large or too small. For small Delta_T, NTimes must be large enough to sample the time correlation function. However, studies show that the condition number of the matrix grows exponentially with NTimes. The target value is 5. Values above 20 are not recommended; \(C_{temporal}\) is numerically singular.
Spatial degeneracy, leading to \(C_{spatial}=0\). We have only one means of entering the spatial correlation parameters, viz. the correlation_length variables pair. If either of these quantities are so large that \(\delta / \mbox{correlation\_length}\) is very close to zero (with \(\delta\) representing the distance from one node to another on the mesh), then the spatial portion of the matrix becomes singular. Effectively, these locations are no longer independent, but must apply the same load vector.
Using a Delta_T = 1/cutoff_freq and the default sinc function for a correlation function may generate a \(C_{temporal}\) singularity. [4] This is because we are now evaluating the correlation function at multiples of \(\pi\), where it is always zero.

16.4.4. Time Step#

The integration time step specified in the SOLUTION section must always be less than or equal to Delta_T.

16.4.5. Sinc Function#

The sinc function defined as \(\sin{(x)}/x\) is important in at least two places in the code. First, it is the only function available for the temporal interpolation function. Second, by default, we use the sinc function as the correlation function. In most cases, this use of the function should probably be replaced by another function. We use it as the default because it represents the Fourier transform of a flat PSD, which is the simplest loading.

16.5. Memory, Performance, Parallel and Anything Else of Interest#

The matrices generated for these operations are all square and dense. The matrix order is \(d=n_{spatial}\cdot n_{temporal}\). Here \(n_{spatial}\) is the number of points in the surface and \(n_{temporal}\) is the value of the NTimes parameter. Because memory requirements grow as the square of these variables, it is important to manage these carefully. Practically, models up to \(d=10^5\) are possible in parallel, but they take a lot of time.

The operation count for Cholesky factorization of a dense matrix is of order \(d^3\). Thus, the computational cost increases much faster than model size. The parallel solutions of the Cholesky system are not scalable. In a scalable problem, doubling the size of the problem, and also doubling the number of processors should not change the solution time. Although the sparse linear solvers for FE solution are scalable, the Cholesky factorization required to compute random pressure loads is not scalable.

The dense Cholesky factorization from the ScaLAPACK library is used. The parallel decomposition for this solve is completely different from the FEM decomposition, and is computed internally without user intervention. The user input for the parallel solution is exactly the same as the serial input. However, at this time, parallel solutions are limited to platforms built under the Intel compiler with MKL libraries. The solution will fail on other platforms.