Files
김경종 4cc312954f
Tests / Hermetic test suite (push) Has been cancelled
Tests / Skill frontmatter validation (push) Has been cancelled
add wiki
2026-05-28 17:16:48 +09:00

399 lines
23 KiB
Markdown

<!-- source-page: 291 -->
EXAMPLE 4.30: Consider the two-node beam element shown in Fig. E4.30. Assume linear variations in the transverse displacement w and section rotation $\theta$ and a constant element transverse shear strain $\gamma$ . Establish the finite element equations.
$E =$ Young's modulus $G =$ shear modulus
![](images/page-291_6690cbc5c791f40c47c83e08f4fcd89dcfa4bd00b4531b527812babba30f4d72.jpg)
![](images/page-291_be6864772581dd3b47a75e07b9870430d79992fdfa822bdb16b76472682ada26.jpg)
<details>
<summary>text_image</summary>
w₁
θ₁
1
z, w
x, u
w₂
θ₂
2
L/2
L/2
</details>
Figure E4.30 Two-node beam element
We assume that the stresses are given by the strains, and so we can substitute $\pmb{\tau} = \mathbf{C}\pmb{\epsilon}$ into (4.114) and obtain
$$
\Pi_ {\mathrm{HR}} ^ {*} = \int_ {V} \left(- \frac {1}{2} \boldsymbol {\epsilon} ^ {T} \mathbf {C} \boldsymbol {\epsilon} + \boldsymbol {\epsilon} ^ {T} \mathbf {C} \boldsymbol {\partial} _ {\boldsymbol {\epsilon}} \mathbf {u} - \mathbf {u} ^ {T} \mathbf {f} ^ {B}\right) d V + \text { boundary terms } \tag {a}
$$
This variational indicator is also a Hellinger-Reissner functional, but comparing (a) with the functional in Exercise 4.36, we note that here strains and displacements are the independent variables (instead of the stresses and displacements in Exercise 4.36).
In our beam formulation the variables are u, w, and $\gamma_{xz}^{AS}$ (the superscript AS denotes the assumed constant value). Hence, the bending strain $\epsilon_{xx}$ is calculated from the displacement, and we can specialize (a) further:
$$
\tilde {\Pi} _ {H R} ^ {*} = \int_ {V} \left(\frac {1}{2} \epsilon_ {x x} E \epsilon_ {x x} - \frac {1}{2} \gamma_ {x z} ^ {A S} G \gamma_ {x z} ^ {A S} + \gamma_ {x z} ^ {A S} G \gamma_ {x z} - \mathbf {u} ^ {T} \mathbf {f} ^ {B}\right) d V + \text { boundary terms }
$$
where $\mathbf{u} = \begin{bmatrix} u \\ w \end{bmatrix}$ ; $\epsilon_{xx} = \frac{\partial u}{\partial x}$ ; $\gamma_{xz} = \frac{\partial w}{\partial x} + \frac{\partial u}{\partial z}$
Now invoking $\delta \tilde{\Pi}_{\mathrm{HR}}^{*} = 0$ , we obtain corresponding to $\delta \mathbf{u}$ , (not including boundary terms)
$$
\int_ {V} \left(\delta \epsilon_ {x x} E \epsilon_ {x x} + \delta \gamma_ {x z} G \gamma_ {x z} ^ {\mathrm{AS}}\right) d V = \int_ {V} \delta \mathbf {u} ^ {T} \mathbf {f} ^ {B} d V \tag {b}
$$
and corresponding to $\delta\gamma_{xz}^{AS}$ ,
$$
\int_ {V} \delta \gamma_ {x z} ^ {\mathrm{AS}} G (\gamma_ {x z} - \gamma_ {x z} ^ {\mathrm{AS}}) d V = 0 \tag {c}
$$
Let $\hat{\mathbf{u}} = \begin{bmatrix} w_1 \\ \theta_1 \\ w_2 \\ \theta_2 \end{bmatrix}$ ; $\hat{\boldsymbol{\epsilon}} = [\gamma^{\mathrm{AS}}]$
Then we can write
$$
\mathbf {u} = \mathbf {H} \hat {\mathbf {u}}; \quad \epsilon_ {x x} = \mathbf {B} _ {b} \hat {\mathbf {u}}
$$
$$
\gamma_ {x z} = \mathbf {B} _ {s} \hat {\mathbf {u}}; \quad \gamma_ {x z} ^ {\mathrm{AS}} = \mathbf {B} _ {s} ^ {\mathrm{AS}} \hat {\boldsymbol {\epsilon}}
$$
<!-- source-page: 292 -->
Substituting into (b) and (c), we obtain
$$
\left[ \begin{array}{l l} \mathbf {K} _ {u u} & \mathbf {K} _ {u \epsilon} \\ \mathbf {K} _ {u \epsilon} ^ {T} & \mathbf {K} _ {\epsilon \epsilon} \end{array} \right] \left[ \begin{array}{l} \hat {\mathbf {u}} \\ \hat {\boldsymbol {\epsilon}} \end{array} \right] = \left[ \begin{array}{c} \mathbf {R} _ {B} \\ \mathbf {0} \end{array} \right] \tag {d}
$$
where $\mathbf{K}_{uu} = \int_{V}\mathbf{B}_{b}^{T}E\mathbf{B}_{b}dV;$ $\mathbf{K}_{u\epsilon} = \int_{V}\mathbf{B}_{s}^{T}G\mathbf{B}_{s}^{\mathrm{AS}}dV$
$$
\mathbf {K} _ {\epsilon \epsilon} = - \int_ {V} \left(\mathbf {B} _ {s} ^ {A S}\right) ^ {T} G \mathbf {B} _ {s} ^ {A S} d V; \quad \mathbf {R} _ {B} = \int_ {V} \mathbf {H} ^ {T} \mathbf {f} ^ {B} d V
$$
We can now use static condensation on $\hat{e}$ to obtain the final element stiffness matrix:
$$
\mathbf {K} = \mathbf {K} _ {u u} - \mathbf {K} _ {u \epsilon} \mathbf {K} _ {\epsilon \epsilon} ^ {- 1} \mathbf {K} _ {u \epsilon} ^ {T}
$$
In our case, we have
$$
\mathbf {H} = \left[ \begin{array}{c c c c} 0 & - \frac {z}{L} \left(\frac {L}{2} - x\right) & 0 & - \frac {z}{L} \left(\frac {L}{2} + x\right) \\ \frac {1}{L} \left(\frac {L}{2} - x\right) & 0 & \frac {1}{L} \left(\frac {L}{2} + x\right) & 0 \end{array} \right]
$$
$$
\mathbf {B} _ {b} = \left[ \begin{array}{c c c c} 0 & \frac {z}{L} & 0 & - \frac {z}{L} \end{array} \right]
$$
$$
\mathbf {B} _ {s} = \left[ - \frac {1}{L} - \frac {1}{L} \left(\frac {L}{2} - x\right) \frac {1}{L} - \frac {1}{L} \left(\frac {L}{2} + x\right) \right]
$$
$$
\mathbf {B} _ {s} ^ {\mathrm{AS}} = [ 1 ]
$$
so that $\mathbf{K} = \begin{bmatrix} \frac{Gh}{L} & \frac{Gh}{2} & \frac{-Gh}{L} & \frac{Gh}{2} \\ \frac{Gh}{2} & \frac{GhL}{4} + \frac{Eh^3}{12L} & \frac{-Gh}{2} & \frac{GhL}{4} - \frac{Eh^3}{12L} \\ \frac{-Gh}{L} & \frac{-Gh}{2} & \frac{Gh}{L} & \frac{-Gh}{2} \\ \frac{Gh}{2} & \frac{GhL}{4} - \frac{Eh^3}{12L} & \frac{-Gh}{2} & \frac{GhL}{4} + \frac{Eh^3}{12L} \end{bmatrix}$ (e)
It is interesting to note that a pure displacement formulation would give a very similar stiffness matrix. The only difference is that the circled terms would be GhL/3 on the diagonal and GhL/6 in the off-diagonal locations. However, the element predictive capability of the pure displacement-based formulation is drastically different, displaying a behavior that is much too stiff when the element is thin (we discuss this phenomenon in Sections 4.5.7 and 5.4.1).
Note that if we assume a displacement vector corresponding to section rotations only,
$$
\hat {\mathbf {u}} = \left[ \begin{array}{c c c c} 0 & \alpha & 0 & - \alpha \end{array} \right]
$$
then using (e) the element displays bending stiffness only, whereas the pure displacement-based element shows an erroneous shear contribution.
Let us finally note that the stiffness matrix in (e) corresponds to the matrix obtained in the mixed interpolation approach discussed in detail in Section 5.4.1. Namely, if we use the last equation in (d), which corresponds to the equation (c), we obtain
$$
\gamma_ {X Z} ^ {\mathrm{AS}} = \frac {w _ {2} - w _ {1}}{L} - \frac {\theta_ {1} + \theta_ {2}}{2}
$$
<!-- source-page: 293 -->
which shows that the assumed shear strain value is equal to the shear strain value at the midpoint of the beam calculated from the nodal point displacements.
As pointed out above, the Hu-Washizu principle provides the basis for the derivation of various variational principles, and many different mixed finite element discretizations can be designed. However, whether a specific finite element discretization is effective for practical analysis depends on a number of factors, particularly on whether the method is general for a certain class of applications, whether the method is stable with a sufficiently high rate of convergence, how efficient the method is computationally, and how the method compares to alternative schemes. While mixed finite element discretizations can offer some advantages in certain analyses, compared to the standard displacement-based discretization, there are two large areas in which the use of mixed elements is much more efficient than the use of pure displacement-based elements. These two areas are the analysis of almost incompressible media and the analysis of plate and shell structures (see the following sections and Section 5.4).
# 4.4.3 Mixed Interpolation—Displacement/Pressure Formulations for Incompressible Analysis
The displacement-based finite element procedure described in Section 4.2 is very widely used because of its simplicity and general effectiveness. However, there are two problem areas in which the pure displacement-based finite elements are not sufficiently effective, namely, the analysis of incompressible (or almost incompressible) media and the analysis of plates and shells. In each of these cases, a mixed interpolation approach—which can be thought of as a special use of the Hu-Washizu variational principle (see Example 4.30)—is far more efficient.
We discuss the mixed interpolation for beam, plate, and shell analyses in Section 5.4, and we address here the analysis of incompressible media.
Although we are dealing with the solution of incompressible solid media, the same basic observations are also directly applicable to the analysis of incompressible fluids (see Section 7.4). For example, the elements summarized in Tables 4.6 and 4.7 (later in this section) are also used effectively in fluid flow solutions.
# The Basic Differential Equations for Incompressible Analysis
In the analysis of solids, it is frequently necessary to consider that the material is almost incompressible. For example, some rubberlike materials, and materials in inelastic conditions, may exhibit an almost incompressible response. Indeed, the compressibility effects may be so small that they could be neglected, in which case the material would be idealized as totally incompressible.
A basic observation in the analysis of almost incompressible media is that the pressure is difficult to predict accurately. Depending on how close the material is to being incompressible, the displacement-based finite element method may still provide accurate solutions, but the number of elements required to obtain a given solution accuracy is usually far greater than the number of elements required in a comparable analysis involving a compressible material.
<!-- source-page: 294 -->
To identify the basic difficulty in more detail, let us again consider the three-dimensional body in Fig. 4.1. The material of the body is isotropic and is described by Young's modulus $E$ and Poisson's ratio $\nu$ .
Using indicial notation, the governing differential equations for this body are (see Example 4.2)
$$
\tau_ {i j, j} + f _ {i} ^ {B} = 0 \quad \text { throughout the volume } V \text { of the body } \tag {4.122}
$$
$$
\tau_ {i j} n _ {j} = f _ {i} ^ {S _ {f}} \quad \text { on } S _ {f} \tag {4.123}
$$
$$
u _ {i} = u _ {i} ^ {S _ {u}} \quad \text { on } S _ {u} \tag {4.124}
$$
If the body is made of an almost incompressible material, we anticipate that the volumetric strains will be small in comparison to the deviatoric strains, and therefore we use the constitutive relations in the form (see Exercise 4.39)
$$
\tau_ {i j} = \kappa \epsilon_ {v} \delta_ {i j} + 2 G \epsilon_ {i j} ^ {\prime} \tag {4.125}
$$
where $\kappa$ is the bulk modulus,
$$
\kappa = \frac {E}{3 (1 - 2 \nu)} \tag {4.126}
$$
$\epsilon_{v}$ is the volumetric strain,
$$
\epsilon_ {V} = \epsilon_ {k k}
$$
$$
= \frac {\Delta V}{V} (= \epsilon_ {x x} + \epsilon_ {y y} + \epsilon_ {z z} \text { in Cartesian coordinates }) \tag {4.127}
$$
$\delta_{ij}$ is the Kronecker delta,
$$
\delta_ {i j} \left\{ \begin{array}{l l} = 1; & i = j \\ = 0; & i \neq j \end{array} \right. \tag {4.128}
$$
$\epsilon_{ij}^{\prime}$ are the deviatoric strain components,
$$
\epsilon_ {i j} ^ {\prime} = \epsilon_ {i j} - \frac {\epsilon_ {V}}{3} \delta_ {i j} \tag {4.129}
$$
and $G$ is the shear modulus,
$$
G = \frac {E}{2 (1 + \nu)} \tag {4.130}
$$
We also have for the pressure in the body,
$$
p = - \kappa \epsilon_ {V} \tag {4.131}
$$
where $p = -\frac{\tau_{kk}}{3}\left(= -\frac{\tau_{xx} + \tau_{yy} + \tau_{zz}}{3}\right)$ in Cartesian coordinates) (4.132)
<!-- source-page: 295 -->
Now let us gradually increase $\kappa$ (by increasing the Poisson ratio $\nu$ to approach 0.5). Then, as $\kappa$ increases, the volumetric strain $\epsilon_{\nu}$ decreases and becomes very small.
In fact, in total incompressibility $\nu$ is exactly equal to 0.5, the bulk modulus is infinite, the volumetric strain is zero, and the pressure is of course finite (and of the order of the applied boundary tractions). The stress components are then expressed as [see (4.125) and (4.131)]
$$
\tau_ {i j} = - p \delta_ {i j} + 2 G \epsilon_ {i j} ^ {\prime} \tag {4.133}
$$
and the solution of the governing differential equations (4.122) to (4.124) now involves using the displacements and the pressure as unknown variables.
In addition, special attention need also be given to the boundary conditions in (4.123) and (4.124) when material incompressibility is being considered and the displacements are prescribed on the complete surface of the body, i.e., when we have the special case $S_{u} = S$ , $S_{f} = 0$ . If the material is totally incompressible, a first condition is that the prescribed displacements $u_{i}^{S}$ must be compatible with the zero volumetric strain throughout the body. This physical observation is expressed as
$$
\epsilon_ {i i} = 0 \quad \text { throughout } V \tag {4.134}
$$
hence, $\int_{V}\epsilon_{ii}dV = \int_{S}\mathbf{u}^{s}\cdot \mathbf{n}dS = 0$ (4.135)
where we used the divergence theorem and n is the unit normal vector on the surface of the body. Hence, the displacements prescribed normal to the body surface must be such that the volume of the body is preserved. This condition will of course be automatically satisfied if the prescribed surface displacements are zero (the particles on the surface of the body are not displaced).
Assuming that the volumetric strain/boundary displacement compatibility is satisfied, for the case $S_{u} = S$ , the second condition is that the pressure must be prescribed at some point in the body. Otherwise, the pressure is not unique because an arbitrary constant pressure does not cause any deformations. Only when both these conditions are fulfilled is the problem well posed for solution.
Of course, the condition of prescribed displacements on the complete surface of the body is a somewhat special case in the analysis of solids, but we encounter an analogous situation frequently in fluid mechanics. Here the velocities may be prescribed on the complete boundary of the fluid domain (see Chapter 7).
Although we considered here a totally incompressible medium, it is clear that these considerations are also important when the material is only almost incompressible—a violation of the conditions discussed will lead to an ill-posed problem statement.
Of course, these observations also pertain to the use of the principle of virtual work. Let us consider the simple example shown in Fig. 4.19. Since only volumetric strain energy is present, the principle of virtual work gives for this case,
$$
\int_ {V} \overline {{\epsilon}} _ {V} \kappa \epsilon_ {V} d V = - \int_ {s _ {f}} \overline {{v}} ^ {s} p ^ {*} d S \tag {4.136}
$$
<!-- source-page: 296 -->
![](images/page-296_2abcdc2c16735bdbec8a702658e275e15c7ce443f3a249d6f0ea7a5b478e3cc4.jpg)
<details>
<summary>text_image</summary>
p*
v^s
S_f
Bulk modulus κ
Pressure p
y, v
L
x, u
</details>
Figure 4.19 Block of material in plane strain condition, subjected to uniform surface pressure $p^{*}$ ; G = 0
If the bulk modulus $\kappa$ is finite, we obtain directly from (4.136),
$$
v ^ {s} = - \frac {p ^ {*} L}{\kappa} \tag {4.137}
$$
and $p = p^{*}$ (4.138)
However, if $\kappa$ is infinite, we need to use instead of (4.136) the following form of the principle of virtual work, with the pressure $p$ unknown,
$$
\int_ {V} \overline {{\epsilon}} _ {V} (- p) d V = - \int_ {S _ {f}} \overline {{v}} ^ {S} p ^ {*} d S \tag {4.139}
$$
and we again obtain $p = p^{*}$ . Of course, the solution of (4.139) does not use the constitutive relation but only the equilibrium condition.
# The Finite Element Solution of Almost Incompressible Conditions
The preceding discussion indicates that when pursuing a pure displacement-based finite element analysis of an almost incompressible medium, significant difficulties must be expected. The very small volumetric strain, approaching zero in the limit of total incompressibility, is determined from derivatives of displacements, which are not as accurately predicted as the displacements themselves. Any error in the predicted volumetric strain will appear as a large error in the stresses, and this error will in turn also affect the displacement prediction since the external loads are balanced (using the principle of virtual work) by the stresses. In practice, therefore, a very fine finite element discretization may be required to obtain good solution accuracy.
<!-- source-page: 297 -->
![](images/page-297_f2e16a65759be099a7fa572a870952f4342ee234edd52626b7840a90fdd3f257.jpg)
<details>
<summary>text_image</summary>
p = 0.006 MPa
δ
10 mm
5 mm
10 mm
50 mm
10 mm
E = 55 MPa
v = 0.3, 0.499
Symmetry about
centerline of
bracket
</details>
![](images/page-297_4f64da14272a15017e998753a57db397fc9a4a008d92e30402462de3f0ac45b8.jpg)
<details>
<summary>natural_image</summary>
Pure geometric diagram of a T-shaped structure with dots and lines, no text or symbols present
</details>
(a) Geometry, material data, applied loading, and the coarse sixteen element mesh
Figure 4.20 Analysis of cantilever bracket in plane strain conditions. Nine-node displacement-based elements are used. The $16 \times 64 = 1024$ -element mesh is obtained by dividing each element of the 16-element mesh into 64 elements. Maximum principal stress results are shown using the band representation of Fig. 4.15. Also, $(\sigma_{1})_{\max}$ is the predicted maximum value of the maximum principal stress, and $\delta$ is defined in (a).
<!-- source-page: 298 -->
![](images/page-298_a90331681b5bff3fa2ef70bd4e9c9ddb1a7b455c1c5e36689b05aa24bcb0e653.jpg)
<details>
<summary>heatmap</summary>
| Value |
| ----------- |
| (σ₁)ₘₐₓ = 0.5050 |
| δ = 1.582 |
| (σ₁)ₘₐₓ = 0.6056 |
| δ = 1.669 |
</details>
(b) Displacement-based element solution results for the case Poisson's ratio $\nu = 0.30$ . Sixteen element and $16 \times 64$ element mesh results
Figure 4.20 (continued)
Figure 4.20 shows some results obtained in the analysis of a cantilever bracket subjected to pressure loading. We consider plane strain conditions and the cases of Poisson's ratio $\nu = 0.30$ and $\nu = 0.499$ . In all solutions, nine-node displacement-based elements have been used (with $3 \times 3$ Gauss integration; see Section 5.5.5). A coarse mesh and a very fine mesh are used, and Fig. 4.20(a) shows the coarse idealization using only 16 elements. The solution results for the maximum principal stress $\sigma_{1}$ are shown using the isoband representation discussed in Section 4.3.6. Here we have selected the bandwidth so as to be able to see the rather poor performance of the displacement-based element when the Poisson ratio is close to 0.5. Figure 4.20(b) shows that when $\nu = 0.30$ , the element stresses are reasonably smooth across boundaries for the coarse mesh and very smooth for the fine mesh. Indeed, the coarse idealization gives a quite reasonable stress prediction. However,
<!-- source-page: 299 -->
![](images/page-299_87728e35ba227ba6c28a49169163b079dd4e7a6838a16d21ab528d45cb4ea8a2.jpg)
<details>
<summary>area</summary>
| Stress Value | Max Stress (σ₁) |
| ------------ | --------------- |
| 1.955 | 1.044 |
</details>
![](images/page-299_4d574640efffd8fe770306299dd18cb92ecd33233060c6f280d43d2dc7c97522.jpg)
<details>
<summary>surface_3d</summary>
| Metric | Value |
|--------|-------|
| (σ₁)max | 1.343 |
| δ | 1.363 |
</details>
(c) Displacement-based element solution results for the case Poisson's ratio $\nu = 0.499$ . Sixteen element and $16 \times 64$ element mesh results
Figure 4.20 (continued)
when $\nu = 0.499$ , the same meshes of nine-node displacement-based elements result into poor stress predictions [see Fig. 4.20(c)]. Large stress fluctuations are seen in the individual elements of the coarse mesh and the fine mesh. $^{16}$ Hence, in summary, we see here that the displacement-based element used in the analysis is effective when $\nu = 0.3$ , but as $\nu$ approaches 0.5, the stress prediction becomes very inaccurate.
This discussion indicates what is very desirable, namely, a finite element formulation which gives essentially the same accuracy in results for a given mesh irrespective of what Poisson's ratio is used, even when $\nu$ is close to 0.5. Such behavior is observed if for the finite
$^{16}$ We discuss briefly in Section 5.5.6 the use of “reduced integration.” If in this analysis the reduced integration of $2 \times 2$ Gauss integration is attempted, the solution cannot be obtained because the resulting stiffness matrix is singular.
<!-- source-page: 300 -->
element formulation the predictive capability of displacements and stresses is independent of the bulk modulus used.
We refer to finite element formulations with this desirable behavior as nonlocking, whereas otherwise the finite elements are locking.
The term “locking” is based upon experiences in the analysis of beams, plates, and shells (see Section 5.4.1), where an inappropriate formulation—one that locks—results in displacements very much smaller than those intuitively expected for a given mesh (and calculated with an appropriate formulation; see, for example, Fig. 5.20). In the analysis of almost incompressible behavior, using a formulation that locks, the displacements are not necessarily that much in error but the stresses (the pressure) are very inaccurate. We note that the pure displacement formulation generally locks in almost incompressible analysis. These statements are discussed more precisely in Section 4.5.
Effective finite element formulations for the analysis of almost incompressible behavior that do not lock are obtained by interpolating displacements and pressure. Figure 4.21 shows the results obtained in the analysis of the cantilever bracket in Fig. 4.20 with a displacement/pressure formulation referred to as u/p formulation using the 9/3 element (see below for the explanation of the formulation and the element). We see that the isobands of the maximum principal stress have in all cases the desirable degree of smoothness and that the stress prediction does not deteriorate when Poisson's ratio v approaches 0.5.
To introduce the displacement/pressure formulations, we recall that in a pure displacement formulation, the evaluation of the pressure from the volumetric strain is difficult when $\kappa$ is large (in comparison to G) and that when a totally incompressible condition is considered, the pressure must be used as a solution variable [see (4.133)]. It therefore appears reasonable to work with the unknown displacements and pressure as solution variables when almost incompressible conditions are analyzed. Such analysis procedures, if properly formulated, should then also be directly applicable to the limit of incompressible conditions.
The basic approach of displacement/pressure finite element formulations is therefore to interpolate the displacements and the pressure. This requires that we express the principle of virtual work in terms of the independent variables u and p, which gives
$$
\boxed {\int_ {V} \overline {{{\epsilon}}} ^ {\prime T} \mathbf {S} d V - \int_ {V} \overline {{{\epsilon}}} _ {v} p d V = \mathcal {R}} \tag {4.140}
$$
where, as usual, the overbar indicates virtual quantities, R corresponds to the usual external virtual work [R is equal to the right-hand side of (4.7)], and S and $\epsilon'$ are the deviatoric stress and strain vectors,
$$
\mathbf {S} = \boldsymbol {\tau} + p \boldsymbol {\delta} \tag {4.141}
$$
$$
\epsilon^ {\prime} = \epsilon - \frac {1}{3} \epsilon_ {V} \delta \tag {4.142}
$$
where $\delta$ is a vector of the Kronecker delta symbol [see (4.128)].
Note that using the definition of p in (4.131), a uniform compressive stress gives a positive pressure and that in the simple example in Fig. 4.19, only the volumetric part of the internal virtual work contributed.
In (4.140) we have separated and then summed the deviatoric strain energy and the bulk strain energy. Since the displacements and pressure are considered independent vari-