We will solve a 1-d finite element problem for a tensioned cable using the method of Freytag. This derivation comes directly from the paper "Finite Element Analysis in situ", and is reproduced merely for learning purposes.

A cable with tension \(\lambda\) runs from \(x=a\) to \(x=b\) and have an applied load of \(q\). The governing equation is \(\lambda \frac{\partial^2 u}{\partial x^2} + q = 0\), where \(u\) gives the height of the cable. We have fixed boundary conditions \(u(a)=u_1\) and \(u(b)=u_2\).

In [392]:

a = .3; b = 3.5
lamb = 1.
u1 = 1.; u2 = 2.5

We make use of an approximate distance function \(\omega\) from the boundary. In this case, we can construct it using the distances from each end of the interval, \(\omega_1=x-a\) and \(\omega_2=b-x\). We take a suitable \(\omega = \omega_1^2 + \omega_2^2 - \sqrt{\omega_1^2 + \omega_2^2}\)

In [357]:

t=linspace(a,b,100); w1 = t-a; w2 = b-t; w = w1 + w2 - sqrt(w1*w1+w2*w2)
plt.plot(t,w1,c='r',label='$\omega_1$')
plt.plot(t,w2,c='r',label='$\omega_2$')
plt.plot(t,w,c='b',label='$\omega$')
plt.legend()

Out[357]:

<matplotlib.legend.Legend at 0x1137574d0>

We assume the solution has a form \(u=u_0+u^*=\sum_{i=1}^n C_i\eta_i + u^*\), where \(u_0\) satisfies homogeneous Dirichlet boundary conditions, and \(u^*\) is a function satisfying the nonhomogeneous boundary conditions. We use the distance functions to construct these fuctions.

\(u^* = \frac{\omega_1 u_2 + \omega_2 u_1}{\omega_1 + \omega_2}\).

\(\eta_i = \omega \chi_i\) where \(\chi_i\) is any shape function (that doesn't know about the boundary).

In [408]:

ustar = (w1*u2 + w2*u1)/(w1+w2)
plt.plot(t,w1,c='r',label='$\omega_1$');
plt.plot(t,w2,c='r',label='$\omega_2$');
plt.plot(t,ustar,c='b',label='$u^*$');

To define the shape functions, we use a fill the interval from \(x_i\) to \(x_f\), which are independent of \(a\) and \(b\), but should be chosen so \([x_i,x_f]\) contains \([a,b]\). We use \(n\) shape functions on this interval, resulting in a grid size of \(h=(x_f-x_i)/n\). Using linear hat shape functions

In [410]:

x0 = a-.1; x1 = b+.1; n=20; h=(x1-x0)/n
res_per_sf = 200
res = res_per_sf*n #numerical resolution
t=linspace(x0,x1,res); w1 = t-a; w2 = b-t; w = w1 + w2 - sqrt(w1*w1+w2*w2); ustar = (w1*u2 + w2*u1)/(w1+w2) #redefine distance functions in terms of res
x=linspace(x0,x1,n)
chi = zeros((n,res))
for i in range(n):
    if i!=0:
        interval = where( (x[i-1]<=t)&(t<=x[i]) )[0]
        chi[i,interval] = linspace(0,1,shape(interval)[0])
    if i!=n-1:
        interval = where( (x[i]<=t)&(t<=x[i+1]) )[0]
        chi[i,interval] = linspace(1,0,shape(interval)[0])
eta = w*chi;
for i in range(n):
    plt.plot(t,chi[i],c='r',label='$\chi_i$')
    plt.plot(t,eta[i],c='b',label='$\eta_i$')
plt.xlim(x0,x1);
#plt.savefig('eta.png',dpi=200)

Out[410]:

(0.19999999999999998, 3.6)

Using the assumed form of the solution, we plug this into the governing equation.

\(\int_a^b \left( \lambda\left(\frac{\partial^2 u_0}{\partial x^2} + \frac{\partial^2 u^*}{\partial x^2}\right) + q(x) \right)\eta_j(x) dx = 0, \hspace{10pt} j=1,...,n\)

\(-\sum_{i=1}^n C_i \int_a^b \frac{\partial \eta_j}{\partial x}\lambda \frac{\partial \eta_i}{\partial x}dx = - \int_a^b q(x) \eta_j(x) dx + \int_a^b \frac{\partial \eta_j}{\partial x}\lambda \frac{\partial u^*}{\partial x}dx \hspace{10pt} j=1,...,n\)

That is, we have a linear system \(Ax=b\) where

\(A_{ij} = -\int_a^b \frac{\partial \eta_j}{\partial x}\lambda \frac{\partial \eta_i}{\partial x}dx\)

\(b_i = \int_a^b \frac{\partial \eta_j}{\partial x}\lambda \frac{\partial u^*}{\partial x}dx - \int_a^b q(x) \eta_j(x) dx\)

\(x_i = C_i\).

Differentiating, we have \(\frac{\partial \eta_i}{\partial x} = \frac{\partial \omega}{\partial x}\chi_i + \omega \frac{\partial \chi_i}{\partial x}\).

Using \(\omega\) as above, we have \(\frac{\partial \omega}{\partial x} = \frac{\omega_1 - \omega_2}{\sqrt{\omega_1^2 +\omega_2^2}}\).

Similarly, \(\frac{\partial u^*}{\partial x} = \frac{u_2-u_1}{\omega_1+\omega_2} = \frac{u_2-u_1}{b-a}\)

In [387]:

dw = -(w1-w2)/sqrt(w1*w1+w2*w2)
dchi = zeros((n,res))
for i in range(n):
    if i!=0:
        dchi[i,where( (x[i-1]<t)&(t<x[i]) )[0]] = 1./h
    if i!=n-1:
        dchi[i,where( (x[i]<t)&(t<x[i+1]) )[0]] = -1./h
deta = dw*chi + w*dchi
for i in range(n):
    plt.plot(t,dw,c='m',label='dw')
    plt.plot(t,chi[i],c='c')
    plt.plot(t,dw*chi[i],c='g',label='$\dw*chi_i$')
    plt.plot(t,w*dchi[i],c='y',label='$\dw*chi_i$')
    plt.plot(t,dchi[i],c='r',label='$\dchi_i$')
    plt.plot(t,deta[i],c='b',label='$\eta_i$')
plt.xlim(x0,x1);
plt.ylim(-1/h*1.1,1/h*1.1);

In [402]:

#Ready to set up linear system.
#We need to discount the shape functions which lie completely outside the integration interval!
#For now, assume x0 and x1 are chosen close enough to a and b so all shape functions affect the integral.
q = -1.*ones_like(t) #load function
interval = where((a<=t)&(t<=b))[0] #integration interval
dx = (x1-x0)/res
A = -lamb*sum((deta.reshape(n,1,res)*deta)[...,interval],axis=-1)*dx
dustar = (u2-u1)/(b-a)
B = sum( (lamb*dustar*deta - q*eta)[...,interval] ,axis=-1)*dx