Continuum Section
Lecture 5: Energy formulations and the Ritz Method

In this approach we don't start from the governing differential equation of the problem but from an energy minimization principle, for those problems in which such principle exists. As we saw in the case of "lumped" or discrete equilibrium problems, at equilibrium, the potential energy of the system is minimum (compare notes from first few lectures). This is known as the principle of minimum potential energy.

The definition of the potential energy of a continuum system has to be done for each system under consideration. It always leads to the concept of a functional:
	- the potential energy is a scalar,
	- as before, it is the sum of the internal energy and the potential of the external forces
	- it is a function of the primary continuous unknown of the problem
	- it involves and integral adding the local contributions of the continuous field defining the energy at each point of the continuum
	
Examples 
Potential energy of an elastic bar of possibly varying cross section:

The deformation of the bar is described by the axial displacement u(x).
The strain distribution in the bar is given by the derivative ϵ(x) = u'(x).
The stress distribution in the bar is σ(x) = E(x) ϵ(x), where E(x) is the elastic or Young modulus of the material, which can possibly vary along the axis "x".
The contribution to the internal (elastic strain) energy at each point is given by
		(1/2)σ(x) ϵ(x) = (1/2)E(x) ((ϵ(x))^2)= (1/2)E(x) u'(((x))^2)
The total internal energy is obtained by integrating this "local" contribution in the volume of the bar as:

		(∫_Ω)(1/2)E(x) u'(((x))^2)dΩ
Since the integrand does not vary in the cross section, only in "x", at a fixed "x" we can integrate in the area of the cross section at that point along the axis A(x), where the area of the cross section is allowed to vary along "x" without loss of generality. We obtain for the strain energy
		(∫_0^L)(1/2)E(x)A(x) u'(((x))^2)dx
The potential of the external forces depends on the type and point of application of the load. As a particular example, assume a single concentrated load P acting at the end of the bar (x=L). The contribution to the potential energy of this load is minus the work done by the force:
		- P u(L)
The total potential energy of the system is then:
		Π(u) = (∫_0^L)(1/2)E(x)A(x) u'(((x))^2)dx - P u(L)
We note that the bar will have different values of the potential energy for different values or expressions of the displacement distribution function u(x). That's why we make Π a function of u. In the lumped systems, the potential energy was a function of a discrete number of degrees of freedom. Now it is a "function of a function" or "functional".

According to the Principle of Minimum Potential Energy, the function u(x) which corresponds to equilibrium minimizes the potential energy of the system.

The mathematical tools needed to find the exact function u(x) which minimizes a functional, pertain to the theory of Calculus of Variations and is discussed in courses on energy and variational methods. 

Instead, we will adopt Ritz approach for finding an approximation to this solution:

Ritz method

1) Adopt and approximation to the exact solution
		u(x)~(Subscript[U, N])(x)= (∑_(j=1)^N)(Subscript[c, j])(Subscript[ϕ, j])(x)+(Subscript[ϕ, 0])(x)
2) Replace this approximate function in the continuous functional:
		Π((Subscript[U, N])) = (∫_0^L)(1/2)E(x)A(x) (Subscript[U, N])'(((x))^2)dx - P (Subscript[U, N])(L)
In effect, the functional Π is now only a function of the discrete number of algebraic coefficients (Subscript[c, j])
		Π((Subscript[c, 1]),(Subscript[c, 2]),..., (Subscript[c, N]))
thus replacing the functional with a "continuous" dependence on u(x) with a function of "N" discrete independent variables. 

3) The principle of minimum potential energy can now be applied just as it was done for a discrete system by the necessary conditions:
		((∂Π)/(∂(Subscript[c, 1])))((Subscript[c, 1]),(Subscript[c, 2]),...,(Subscript[c, N])) = 0
((∂Π)/(∂(Subscript[c, 2])))((Subscript[c, 1]),(Subscript[c, 2]),...,(Subscript[c, N])) = 0
...
((∂Π)/(∂(Subscript[c, N])))((Subscript[c, 1]),(Subscript[c, 2]),...,(Subscript[c, N])) = 0
which is a system of "N" "possibly coupled" algebraic equations with "N" unknowns.

Remark: note that in the case of the bar, the equations are going to be linear, as the strain energy is quadratic in the derivatives u'(x) and also in the coefficients (Subscript[c, j]), j=1,N. When differentiated with respect to (Subscript[c, j]), linear expressions will be obtained. The potential of the external forces in linear on u(L), and thus in (Subscript[c, j]), the derivative with respect to (Subscript[c, j]) will thus be a constant (right hand side of the system of linear equations.

Example of the bar:

Consider a bar with a constant Young's modulus "E",
a tapered cross sectional area: A(x) =(Subscript[A, 0])(1-(x/(2L))),
fixed at x=0: u(0) = 0, (essential boundary condition)

Use the approximation functions

The approximation function (Subscript[U, 2]) is:

(Subscript[ϕ, 0])= 0;
(Subscript[ϕ, 1])= x;
(Subscript[ϕ, 2])= (x^2);
(Subscript[U, 2])= (Subscript[ϕ, 0])+(Subscript[c, 1])(Subscript[ϕ, 1])+ (Subscript[c, 2])(Subscript[ϕ, 2])//Simplify

The Potential Energy is:

Π=(∫_0^L)(1/2)(Subscript[EA, 0])(1-(1/2)((((x/L)))^2)) ((((∂_x)(Subscript[U, 2])))^2)\[DifferentialD]x -P ((Subscript[U, 2])/. x->L)-(∫_0^L)(Subscript[q, 0])(Subscript[U, 2])\[DifferentialD]x //Expand

Minimize wrt (Subscript[c, 1]),(Subscript[c, 2]):

eq1 = ((∂_(Subscript[c, 1]))Π ==0)
eq2 = ((∂_(Subscript[c, 2]))Π ==0)
coeffs=Solve[{eq1,eq2},{(Subscript[c, 1]),(Subscript[c, 2])}]//Flatten

Reconstruct solution:

(Subscript[UU, 2])= (Subscript[U, 2])/. coeffs

Exact solution:

Clear[u]
exact=DSolve[{D[(Subscript[EA, 0])(1-(1/2)((((x/L)))^2)) u'[x],x] ==0,u[0]==0, ((Subscript[EA, 0])/2)u'[L]==P}, u[x],x]//Flatten;
u= u[x] /. exact

values = {(Subscript[EA, 0])->70 (10^9)(10^(-2)), L-> 3, P->1000,(Subscript[q, 0])-> 0};
u/. values
(Subscript[UU, 2])/. values

Plot[{u/. values, (Subscript[UU, 2])/. values}, {x,0,3L/. values},PlotStyle->{Thickness[0.007]}]

(Subscript[ϵ, 2]) = (∂_x)(Subscript[UU, 2])

(Subscript[σ, 2])= 70 (10^9) (Subscript[ϵ, 2])/. values

Plot[(Subscript[σ, 2]), {x,0,L/. values},PlotRange-> All]