# Nonlinear finite elements/Homework 4/Solutions

## Problem 1: Axially Loaded Bar with Nonlinear Modulus edit

** Given:**

Bar of uniform cross-section under axial body load and compressive axial end load ( ). Nonlinear material with constitutive relation

Finite element model uses quadratic elements with shape functions

where .

### Solution edit

#### Part 1 edit

Plot the shape functions for an element. Show that they add up to 1 at all points within the element.

A plot of the shape functions for an element is shown in Figure 1(a). The shape functions for a mesh are shown in Figure 1(b). If we add the shape functions we get

#### Part 2 edit

Plot the displacement, stress, strain, and the error norm when you use 10 elements to discretize the bar. Discuss and comment on your observations.

If we use an isoparametric form of the equations to simplify the integration of the stiffness matrix and load vector terms, the shape functions become

The map between and is

Therefore, the Jacobian of the deformation is

The terms of the stiffness matrix are given by

The element stiffness matrix is

The terms of the tangent stiffness matrix are given by

Therefore, the element tangent stiffness matrix is

The terms of the element load vector are given by

Therefore, the element load vector is

The stress-strain curve in compression is shown in Figure 2.

If we use a mesh containing equally sized elements, we get the solutions shown in Figure 3.

If we use a mesh containing equally sized elements, we get the solutions shown in Figure 4.

*Observations:*

- The two-node and three-node elements give the same solution.
- The solution using two elements shows that the stresses and strains are interpolated linearly within each element and are non longer constant.
- There is no jump in the stress/strain at the nodes.
- The solution converges quite rapidly - four iterations. This is due to the shape of the stress-strain curve in compression.

#### Part 2 edit

Do the same for 100 elements. What are your observations on the effect of mesh refinement on the solution?

If we use a mesh containing equally sized elements, we get the solutions shown in Figure 5.

*Observations:*

- The solution is approximated quite accurately with 10 elements. Using 100 elements does not lead to any significant improvement in the solution.
- The solution converges in four iterations unlike the case where a tensile load was applied.

#### Part 4 edit

Now, increase the applied load. At what load does the solution fail to converge? Explain why the solution fails to converge at this load.

The solution converges at all loads unlike in the tension problem. If the load is increased to , the elements invert and we get a negative Jacobian. Even though we get a solution, that solution is unphysical.

The Matlab code use to compute the results is shown below.

function HW4Prob1_2
E0 = 10;
R = -2;
a = 1;
W = 1;
H = 1;
A = 1;
L = 1;
tol = 1.0e-6;
nElem = 10;
h = L/nElem;
nNode = 2*nElem+1;
for i=1:nNode
xnode(i) = (i-1)*h/2;
end
for i=1:nElem
node1 = 2*(i-1)+1;
elem(i,:) = [node1 node1+1 node1+2];
end
u = zeros(nNode,1);
condition = 1;
condition_res = 1;
count = 0;
while (condition > tol)
count = count + 1
K = zeros(nNode);
T = zeros(nNode);
f = zeros(nNode,1);
for i=1:nElem
node1 = elem(i,1);
node2 = elem(i,2);
node3 = elem(i,3);
Ke = elementStiffness(A, E0, xnode, u, node1, node2, node3);
Te = elementTangentStiffness(A, E0, xnode, u, node1, node2, node3);
fe = elementLoad(a, xnode, node1, node2, node3);
K(node1:node3,node1:node3) = K(node1:node3,node1:node3) + Ke;
T(node1:node3,node1:node3) = T(node1:node3,node1:node3) + Te;
f(node1:node3) = f(node1:node3) + fe;
end
f(nNode) = f(nNode) + R;
u0 = 0;
Kred = K(2:nNode,2:nNode);
Tred = T(2:nNode,2:nNode);
fred = f(2:nNode);
ured = u(2:nNode);
r = Kred*ured - fred;
Tinv = inv(Tred);
u_new = ured - Tinv*r;
condition = norm(u_new - ured)/norm(u_new)
u = [u0; u_new];
iter(count) = count;
normres(count) = norm(r);
normu(count) = condition;
end
figure;
p01 = semilogy(iter, normres, 'ro-'); hold on;
p02 = semilogy(iter, normu, 'bo-');
set([p01 p02], 'LineWidth', 2);
xlabel('Iterations','FontName','palatino','FontSize',18);
ylabel('Error Norm','FontName','palatino','FontSize',18);
set(gca, 'LineWidth', 2, 'FontName','palatino','FontSize',18);
legend([p01 p02],'|r|','|u_{n+1}-u_n|/|u_{n+1}|');
axis square;
figure;
p0 = plotDisp(E0, A, L, R, a);
p = plot(xnode, u, 'ro-', 'LineWidth', 2); hold on;
xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18);
set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
legend([p0 p], 'Linear (Exact) ','Nonlinear (FEM)');
axis square;
for i=1:nElem
node1 = elem(i,1);
node2 = elem(i,2);
node3 = elem(i,3);
[xloc(i,:), eps(i,:), sig(i,:)] = ...
elementStrainStress(E0, xnode, u, node1, node2, node3);
end
figure;
p0 = plotStress(E0, A, L, R, a);
for i=1:nElem
p1 = plot(xloc(i,:), sig(i,:), 'r-o','LineWidth',2); hold on;
end
xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18);
set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
legend([p0 p1], 'Linear (Exact) ','Nonlinear (FEM)');
axis square;
figure;
p0 = plotStrain(E0, A, L, R, a);
for i=1:nElem
p1 = plot(xloc(i,:), eps(i,:), 'r-o','LineWidth',2); hold on;
end
xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
ylabel('\epsilon(x)', 'FontName', 'palatino', 'FontSize', 18);
set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
legend([p0 p1], 'Linear (Exact) ','Nonlinear (FEM)');
axis square;
function [Ke] = elementStiffness(A, E0, xnode, u, node1, node2, node3)
x1 = xnode(node1);
x2 = xnode(node2);
x3 = xnode(node3);
h = x3 - x1;
u1 = u(node1);
u2 = u(node2);
u3 = u(node3);
C1 = A*E0/(3*h^2);
Ke(1,1) = 7*h + u3 + 15*u1 - 16*u2;
Ke(1,2) = -8*(h + 2*u1 - 2*u2);
Ke(1,3) = (h - u3 + u1);
Ke(2,2) = 16*(h - u3 + u1);
Ke(2,3) = -8*(h - 2*u3 + 2*u2);
Ke(3,3) = (7*h - 15*u3 - u1 + 16*u2);
Ke(2,1) = Ke(1,2);
Ke(3,1) = Ke(1,3);
Ke(3,2) = Ke(2,3);
Ke = C1*Ke;
function [Te] = elementTangentStiffness(A, E0, xnode, u, ...
node1, node2, node3)
x1 = xnode(node1);
x2 = xnode(node2);
x3 = xnode(node3);
h = x3 - x1;
u1 = u(node1);
u2 = u(node2);
u3 = u(node3);
C1 = A*E0/(3*h^2);
Te(1,1) = 7*h + 2*u3 + 30*u1 - 32*u2;
Te(1,2) = -8*(h + 4*u1 - 4*u2);
Te(1,3) = (h - 2*u3 + 2*u1);
Te(2,2) = 16*(h - 2*u3 + 2*u1);
Te(2,3) = -8*(h - 4*u3 + 4*u2);
Te(3,3) = (7*h - 30*u3 - 2*u1 + 32*u2);
Te(2,1) = Te(1,2);
Te(3,1) = Te(1,3);
Te(3,2) = Te(2,3);
Te = C1*Te;
function [fe] = elementLoad(a, xnode, node1, node2, node3)
x1 = xnode(node1);
x2 = xnode(node2);
x3 = xnode(node3);
h = x3 - x1;
C1 = a*h/6;
fe(1,1) = x1;
fe(2,1) = 2*(x1+x3);
fe(3,1) = x3;
fe = C1*fe;
function [x, eps, sig] = elementStrainStress(E0, xnode, u, ...
node1, node2, node3)
x1 = xnode(node1);
x2 = xnode(node2);
x3 = xnode(node3);
h = x3 - x1;
u1 = u(node1);
u2 = u(node2);
u3 = u(node3);
ximin = -1;
ximax = 1;
nxi = 2;
dxi = (ximax - ximin)/nxi;
for i=1:nxi+1
xi = ximin + (i-1)*dxi;
J = h/2;
dN1_dxi = xi - 1/2;
dN2_dxi = -2*xi;
dN3_dxi = xi + 1/2;
B = [dN1_dxi/J dN2_dxi/J dN3_dxi/J];
u = [u1; u2; u3];
eps(i) = B*u;
sig(i) = E0*(1-eps(i))*eps(i);
F = 1 + eps(i);
if (F < 0)
display('Material has inverted!');
end
N1 = -0.5*(1-xi)*xi;
N2 = (1+xi)*(1-xi);
N3 = 0.5*(1+xi)*xi;
x(i) = x1*N1 + x2*N2 + x3*N3;
end
function [p] = plotDisp(E, A, L, R, a)
dx = 0.01;
nseg = L/dx;
for i=1:nseg+1
x(i) = (i-1)*dx;
u(i) = 1/(6*A*E)*(-a*x(i)^3 + (6*R + 3*a*L^2)*x(i));
end
p = plot(x, u, 'LineWidth', 3); hold on;
xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
ylabel('u(x)', 'FontName', 'palatino', 'FontSize', 18);
set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
function [p] = plotStrain(E, A, L, R, a)
dx = 0.01;
nseg = L/dx;
for i=1:nseg+1
x(i) = (i-1)*dx;
eps(i) = 1/(2*A*E)*(-a*x(i)^2 + (2*R + a*L^2));
end
p = plot(x, eps, 'LineWidth', 3); hold on;
xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
ylabel('\epsilon(x)', 'FontName', 'palatino', 'FontSize', 18);
set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);
function [p] = plotStress(E, A, L, R, a)
dx = 0.01;
nseg = L/dx;
for i=1:nseg+1
x(i) = (i-1)*dx;
sig(i) = 1/(2*A)*(-a*x(i)^2 + (2*R + a*L^2));
end
p = plot(x, sig, 'LineWidth', 3); hold on;
xlabel('x', 'FontName', 'palatino', 'FontSize', 18);
ylabel('\sigma(x)', 'FontName', 'palatino', 'FontSize', 18);
set(gca, 'LineWidth', 3, 'FontName', 'palatino', 'FontSize', 18);

You could alternatively write your own Maple code or use ANSYS for these calculations depending on your interests.

## Problem 2: Nonlinear Thermal Conduction edit

Consider the model chip system shown in Figure 6.

Assume that heat is being generated in the silicon carbide region of the chip at the rate of 100 W/cm . The cooling channels contain a fluid at a constant 250 K and (supposedly) can dissipate heat at the rate of 800 W/cm . The surroundings are at room temperature (300 K).

The thermal conductivity of the silicon carbide and the aluminum nitride are also given in Figure 6. The thermal conductivity of the epoxy resin is 1.5 W/m-K.

#### Part 1 edit

First run a linear simulation of the problem using room temperature thermal conductivities of the materials. This will act as a check of the model and setup and will help you determine any additional constraints that you will need. Do the cooling channels really dissipate heat at 800 W/cm ?

The following ANSYS input is used for this simulation. Due to symmetry, only half of the model will be meshed Figure 7. A constant temperature of 300K is applied around the edge Epoxy resin. Heat is being generated from areas made up by the Silicon Carbide at the rate of 100 W/cm . Each line forming the cooling channels has constant temperature of 250K.

/prep7
et,1,55
mp,kxx,1,28.634
mp,kxx,2,84.873
mp,kxx,3,1.5
radius = 0.00005
ewidth = 0.0001
eheight = 0.0001
i=0
k,i+1,0,0
k,i+2,radius,0
k,i+3,ewidth,0
k,i+4,ewidth,eheight
k,i+5,0,eheight
k,i+6,0,radius
local,11,1,0,0
k,i+7,radius,45
k,i+8,radius,-45
csys
k,i+9,0,-radius
k,i+10,0,-eheight
k,i+11,ewidth,-eheight
l,i+2,i+3,2
l,i+3,i+4,2
l,i+4,i+5,2
l,i+5,i+6,2
local,11,1,0,0
l,i+6,i+7,2
l,i+7,i+2,2
l,i+2,i+8,2
l,i+8,i+9,2
csys
l,i+9,i+10,2
l,i+10,i+11,2
l,i+11,i+3,2
l,i+7,i+4,2,2
l,i+8,i+11,2
a,i+2,i+3,i+4,i+7
a,i+7,i+4,i+5,i+6
a,i+2,i+3,i+11,i+8
a,i+8,i+11,i+10,i+9
AGEN,7,all,,,7*ewidth,0,0,,1
asel,s,area,,5,28
ARSYM,X,all, , , ,1,0
asel,s,area,,29,52
AGEN,,all,,,7*ewidth*7,,,,,1
asel,all
*do,i,1,37,7
blc4,ewidth*i,0,ewidth,eheight
blc4,ewidth*(i+1),0,ewidth,eheight
blc4,ewidth*(i+2),0,ewidth,eheight
blc4,ewidth*(i+3),0,ewidth,eheight
blc4,ewidth*(i+4),0,ewidth,eheight
blc4,ewidth*i,-eheight,ewidth,eheight
blc4,ewidth*(i+1),-eheight,ewidth,eheight
blc4,ewidth*(i+2),-eheight,ewidth,eheight
blc4,ewidth*(i+3),-eheight,ewidth,eheight
blc4,ewidth*(i+4),-eheight,ewidth,eheight
*enddo
blc4,ewidth*43,0,ewidth,eheight
blc4,ewidth*43,-eheight,ewidth,eheight
blc4,ewidth*44,0,ewidth,eheight
blc4,ewidth*44,-eheight,ewidth,eheight
*do,i,1,6
blc4,(i-1)*0.0001,0.0001,0.0001,0.0001
*enddo
*do,i,0,9
blc4,0.002+(0.0001*i),0.0001,0.0001,0.0001
*enddo
*do,i,0,14
blc4,0.0005+(i*0.0001),0.0001,0.0001,0.0001
*enddo
*do,i,0,14
blc4,0.003+(0.0001*i),0.0001,0.0001,0.0001
*enddo
blc4,0.0044,-0.0001,0.0001,0.0001
blc4,0.0044,0,0.0001,0.0001
blc4,0.0044,0.0001,0.0001,0.0001
*do,i,1,51
blc4,(i-1)*ewidth,2*eheight,ewidth,0.0003
blc4,(i-1)*ewidth,-eheight,ewidth,-0.0003
*enddo
*do,i,0,5
blc4,0.0045+(i*0.0001),eheight,ewidth,0.0001
blc4,0.0045+(i*0.0001),0,ewidth,0.0001
blc4,0.0045+(i*0.0001),-eheight,ewidth,0.0001
*enddo
aovlap,all
lesize,all,,,2
aatt,1
asel,s,loc,x,0,0.0045
asel,r,loc,y,-0.0001,0.0001
amesh,all
asel,s,loc,x,0,0.0005
asel,a,loc,x,0.002,0.003
asel,r,loc,y,0.0001,0.0002
amesh,all
asel,all
aatt,2
asel,s,loc,x,0.0005,0.002
asel,a,loc,x,0.003,0.0045
asel,r,loc,y,0.0001,0.0002
amesh,all
asel,all
aatt,3
asel,s,loc,x,0,0.005
asel,r,loc,y,-0.0001,-0.0004
amesh,all
asel,s,loc,x,0,0.005
asel,r,loc,y,0.0002,0.0005
amesh,all
asel,s,loc,x,0.0045,0.005
asel,r,loc,y,-0.0001,0.0002
amesh,all
asel,all
ESEL,S,MAT,,1
/color,elem,12,all
esel,all
ESEL,S,MAT,,2
/color,elem,8,all
esel,all
ESEL,S,MAT,,3
/color,elem,4,all
esel,all
nsel,s,loc,x,0,0
dsym,symm,x
nsel,all
ASEL,S,MAT,,2
bfa,all,hgen,1000000
asel,all
lsel,s,loc,x,0,0.005
lsel,r,loc,y,0.0005
DL,all, ,TEMP,300,1
lsel,s,loc,x,0,0.005
lsel,r,loc,y,-0.0004
DL,all, ,TEMP,300,1
lsel,s,loc,y,-0.0004,0.0005
lsel,r,loc,x,0.005
DL,all, ,TEMP,300,1
lsel,all
LSEL,S,LENGTH,,0.3927E-04
DL,all, ,TEMP,250,1
lsel,all
fini
/solu
solve
fini
/post1
plnsol,temp
plnsol,tfsum
fini

Shown in Figure 8 and 9 below are the contour plots of temperature and the vector sum of heat flux. From Figure 9 we see that the cooling channels do not dissipate heat at the rate suggested in the problem statement.

#### Part 2 edit

Fit curves to the thermal conductivity versus temperature data for Silicon Carbide and Aluminum Nitride.

Below are the curve-fits for each material.

- Silicon Carbide:
- Aluminum Nitride:

#### Part 3 edit

Perform a nonlinear simulation using these temperature dependent properties to find the equilibrium temperature inside the silicon carbide region. State all simplifying assumptions that you have made.

The same code that was used in the previous problem is used here. However, the thermal conductivities are tabulated according to their temperature.

/prep7
et,1,55
MPTEMP,,,,,,,,
MPTEMP,1,20
MPTEMP,2,60
MPTEMP,3,100
MPTEMP,4,140
MPTEMP,5,180
MPTEMP,6,220
MPTEMP,7,260
MPTEMP,8,300
MPDATA,KXX,1,,27
MPDATA,KXX,1,,24
MPDATA,KXX,1,,22
MPDATA,KXX,1,,21
MPDATA,KXX,1,,20
MPDATA,KXX,1,,20
MPDATA,KXX,1,,20
MPDATA,KXX,1,,20
MPDATA,KXX,2,,75
MPDATA,KXX,2,,59
MPDATA,KXX,2,,52
MPDATA,KXX,2,,49
MPDATA,KXX,2,,47
MPDATA,KXX,2,,45
MPDATA,KXX,2,,43
MPDATA,KXX,2,,41
mp,kxx,3,1.5
radius = 0.00005
ewidth = 0.0001
eheight = 0.0001
i=0
k,i+1,0,0
k,i+2,radius,0
k,i+3,ewidth,0
k,i+4,ewidth,eheight
k,i+5,0,eheight
k,i+6,0,radius
local,11,1,0,0
k,i+7,radius,45
k,i+8,radius,-45
csys
k,i+9,0,-radius
k,i+10,0,-eheight
k,i+11,ewidth,-eheight
l,i+2,i+3,2
l,i+3,i+4,2
l,i+4,i+5,2
l,i+5,i+6,2
local,11,1,0,0
l,i+6,i+7,2
l,i+7,i+2,2
l,i+2,i+8,2
l,i+8,i+9,2
csys
l,i+9,i+10,2
l,i+10,i+11,2
l,i+11,i+3,2
l,i+7,i+4,2,2
l,i+8,i+11,2
a,i+2,i+3,i+4,i+7
a,i+7,i+4,i+5,i+6
a,i+2,i+3,i+11,i+8
a,i+8,i+11,i+10,i+9
AGEN,7,all,,,7*ewidth,0,0,,1
asel,s,area,,5,28
ARSYM,X,all, , , ,1,0
asel,s,area,,29,52
AGEN,,all,,,7*ewidth*7,,,,,1
asel,all
*do,i,1,37,7
blc4,ewidth*i,0,ewidth,eheight
blc4,ewidth*(i+1),0,ewidth,eheight
blc4,ewidth*(i+2),0,ewidth,eheight
blc4,ewidth*(i+3),0,ewidth,eheight
blc4,ewidth*(i+4),0,ewidth,eheight
blc4,ewidth*i,-eheight,ewidth,eheight
blc4,ewidth*(i+1),-eheight,ewidth,eheight
blc4,ewidth*(i+2),-eheight,ewidth,eheight
blc4,ewidth*(i+3),-eheight,ewidth,eheight
blc4,ewidth*(i+4),-eheight,ewidth,eheight
*enddo
blc4,ewidth*43,0,ewidth,eheight
blc4,ewidth*43,-eheight,ewidth,eheight
blc4,ewidth*44,0,ewidth,eheight
blc4,ewidth*44,-eheight,ewidth,eheight
*do,i,1,6
blc4,(i-1)*0.0001,0.0001,0.0001,0.0001
*enddo
*do,i,0,9
blc4,0.002+(0.0001*i),0.0001,0.0001,0.0001
*enddo
*do,i,0,14
blc4,0.0005+(i*0.0001),0.0001,0.0001,0.0001
*enddo
*do,i,0,14
blc4,0.003+(0.0001*i),0.0001,0.0001,0.0001
*enddo
blc4,0.0044,-0.0001,0.0001,0.0001
blc4,0.0044,0,0.0001,0.0001
blc4,0.0044,0.0001,0.0001,0.0001
*do,i,1,51
blc4,(i-1)*ewidth,2*eheight,ewidth,0.0003
blc4,(i-1)*ewidth,-eheight,ewidth,-0.0003
*enddo
*do,i,0,5
blc4,0.0045+(i*0.0001),eheight,ewidth,0.0001
blc4,0.0045+(i*0.0001),0,ewidth,0.0001
blc4,0.0045+(i*0.0001),-eheight,ewidth,0.0001
*enddo
aovlap,all
lesize,all,,,2
aatt,1
asel,s,loc,x,0,0.0045
asel,r,loc,y,-0.0001,0.0001
amesh,all
asel,s,loc,x,0,0.0005
asel,a,loc,x,0.002,0.003
asel,r,loc,y,0.0001,0.0002
amesh,all
asel,all
aatt,2
asel,s,loc,x,0.0005,0.002
asel,a,loc,x,0.003,0.0045
asel,r,loc,y,0.0001,0.0002
amesh,all
asel,all
aatt,3
asel,s,loc,x,0,0.005
asel,r,loc,y,-0.0001,-0.0004
amesh,all
asel,s,loc,x,0,0.005
asel,r,loc,y,0.0002,0.0005
amesh,all
asel,s,loc,x,0.0045,0.005
asel,r,loc,y,-0.0001,0.0002
amesh,all
asel,all
ESEL,S,MAT,,1
/color,elem,12,all
esel,all
ESEL,S,MAT,,2
/color,elem,8,all
esel,all
ESEL,S,MAT,,3
/color,elem,4,all
esel,all
KBC,1
nsel,s,loc,x,0,0
dsym,symm,x
nsel,all
ASEL,S,MAT,,2
bfa,all,hgen,1000000
asel,all
lsel,s,loc,x,0,0.005
lsel,r,loc,y,0.0005
DL,all, ,TEMP,300,1
lsel,s,loc,x,0,0.005
lsel,r,loc,y,-0.0004
DL,all, ,TEMP,300,1
lsel,s,loc,y,-0.0004,0.0005
lsel,r,loc,x,0.005
DL,all, ,TEMP,300,1
lsel,all
LSEL,S,LENGTH,,0.3927E-04
DL,all, ,TEMP,250,1
lsel,all
fini
/solu
solve
fini
/post1
plnsol,temp
plnsol,tfsum
fini

Shown in Figure 10 and 11 below are the contour plots of temperature and the vector sum of heat flux.

There are many approaches to create the mesh for this problem. You will have to find the way the work best for you. These two codes listed are just two ways that the meshes can be manipulated.