Laminate optimization¶
Simply put, optimization of a structure is the process of making the best or most effective use of a situation and/or a resource. Optimal solutions must generally satisfy a given set of design constraints and requirements that reflects both the limitations of the reality and the objectives for the optimization.
High-performance structural composites are most commonly found in light-weight design structures, suggesting that the most common objective for optimization of such designs is to minimize the weight.
Example 1: strength optimization¶
Constraints and requirements:
Ply material: Carbon/Epoxy(a)
Ply thickness: free variable
Possible orientations: 0, 90
Symmetric and balanced laminate
Load cases: Nx = 1000 N/mm and Ny = 2000 Nmm/mm
Stress exposure factor according to Maximum stress criterion shall be equal to 1.0 considering all plies and all positions
Objectives:
- Minimize the mass
Since the we are free to vary the thickness of individual plies, the task is rather straight-forward and an optimum solution can be found using the following steps:
- Create a laminate as simple as possible, which for the given constraints is either [0/90/0] or [90/0/90]
- Assume a given total thickness, for example unity $t=1$
- By brute force, solve for all combinations of thicknesses where $t_0=0...t$ while $t_{90} = t - t_0$.
- When the minimum exposure factor is found along with the optimum combination of thicknesses, scale the total thickness by the exposure factor (remember it is a linear system of relations).
Code:
import laminatelib, matlib, numpy
import matplotlib.pyplot as plt
m = matlib.get('Carbon/Epoxy(a)')
t=1 # total thickness
ts=numpy.linspace(0,t,100) # 100 values from 0 to t
fE=[] # to be filled with exposure factors for all values in ts
for t0 in ts:
t90=t-t0
layup = [ {'mat':m , 'ori': 0 , 'thi':t0/2},
{'mat':m , 'ori': 90 , 'thi':t90},
{'mat':m , 'ori': 0 , 'thi':t0/2}]
ABD=laminatelib.laminateStiffnessMatrix(layup)
load1,def1=laminatelib.solveLaminateLoadCase(ABD, Nx=1000, Ny=2000)
res=laminatelib.layerResults(layup,def1)
fE.append( max( res[0]['fail']['MS']['bot'], res[1]['fail']['MS']['bot']) ) # top layer results will be the same
plt.plot(ts,fE)
plt.show()
As shown in the figure, there is only one minimum for the exposure factor at approximatly $t_0=0.3$. Extracting the exact values:
minfE = numpy.amin(fE)
ind = numpy.where(fE == minfE)
t0=ts[ind[0]][0]
print(minfE,t0)
Scaling all thicknesses by the exposure factor:
t = t*minfE
t0=t0*minfE
t90=t-t0
print(t0,t90,t)
Solve with these values to verify that the exposure factor is equal to 1.0
layup = [ {'mat':m , 'ori': 0 , 'thi':t0/2},
{'mat':m , 'ori': 90 , 'thi':t90},
{'mat':m , 'ori': 0 , 'thi':t0/2}]
ABD=laminatelib.laminateStiffnessMatrix(layup)
load1,def1=laminatelib.solveLaminateLoadCase(ABD, Nx=1000, Ny=2000)
res=laminatelib.layerResults(layup,def1)
print( max( res[0]['fail']['MS']['bot'], res[1]['fail']['MS']['bot']) )
Example 2: stiffness optimization, fixed ply thickness¶
Constraints and requirements:
Ply material: Carbon/Epoxy(a)
Ply thickness: 0.2 mm
Possible orientations: 0 and/or 90
Symmetric and balanced laminate
Loads: Mx = 1000 Nmm/mm and My = 200 Nmm/mm
Allowed deformation: κx < 0.01 mm-1 and κy < 0.01 mm-1
Objectives:
- Minimize the mass
One major difference from the first example is the quantisized ply thickness.
For a laminate with a single 0-oriented layer:
layup = [ {'mat':m , 'ori': 0 , 'thi':0.2} ]
ABD=laminatelib.laminateStiffnessMatrix(layup)
load1,def1=laminatelib.solveLaminateLoadCase(ABD, Mx=1000, My=200)
print('[Kx, Ky, Kxy =]',def1[3:6])
The bending stiffness is clearly way to low (the resulting maximum curvature is almost 27 compared to the limit of 0.01).
A simple optimization technique is to add plies until the constraints are satisfied. For simplicity, let us just add more 0-oriented plies for this demonstration:
layup=[]
for i in range(0,100):
layup.append( {'mat':m , 'ori': 0 , 'thi':0.2} )
ABD=laminatelib.laminateStiffnessMatrix(layup)
load,deformation=laminatelib.solveLaminateLoadCase(ABD, Mx=1000, My=200)
if max(deformation[3:6])<0.01:
print('No. of plies=',i+1, 'where [Kx, Ky, Kxy =]',deformation[3:6])
break # satisfied, no reason to continue.
While the constraints and requirements are met with 14 plies, this is most likely not the optimum solution; some 90-oriented plies will probably reduce the required number of plies.
The brute force technique tests for all possible designs and will therefor find the absolut optimum solution.
The following code generates a list of all possible layup combinations for an even number of plies (6 plies for this example) where the number of combinations is 23 = 8
combs=[] # List of combinitions
angs=(0,90) # Possible orientations
for a1 in angs:
for a2 in angs:
for a3 in angs:
combs.append([a1,a2,a3,a3,a2,a1])
print(combs)
Alternative pythonic syntax of the same is:
angs=(0,90)
combs = [[a1,a2,a3,a3,a2,a1] for a1 in angs for a2 in angs for a3 in angs ]
print(combs)
Testing the performance of all eight layups and sort the result from best to worst:
maxdef=[]
for c in combs:
layup=[]
for a in c:
layup.append( {'mat':m , 'ori': a , 'thi':0.2} )
ABD=laminatelib.laminateStiffnessMatrix(layup)
load,deformation=laminatelib.solveLaminateLoadCase(ABD, Mx=1000, My=200)
maxdef.append(max(deformation[3:6]))
zipped = zip(maxdef, combs)
sorted_pairs = sorted(zipped)
for s in sorted_pairs:
print('Max curvature:',s[0],'Layup:',s[1])
Since none of the solutions satisfied the limit on deformation, let us try with 12 plies:
angs=(0,90)
combs = [[a1,a2,a3,a4,a5,a6,a6,a5,a4,a3,a2,a1]
for a1 in angs for a2 in angs for a3 in angs for a4 in angs for a5 in angs for a6 in angs]
print('No. of combinations:',len(combs))
Results for the laminates that satisfy the constraints and requirements:
maxdef=[]
for c in combs:
layup=[]
for a in c:
layup.append( {'mat':m , 'ori': a , 'thi':0.2} )
ABD=laminatelib.laminateStiffnessMatrix(layup)
load,deformation=laminatelib.solveLaminateLoadCase(ABD, Mx=1000, My=200)
maxdef.append(max(deformation[3:6]))
zipped = zip(maxdef, combs)
sorted_pairs = sorted(zipped)
for s in sorted_pairs:
if s[0]<0.01:
print('Max curvature:',s[0],'Layup:',s[1])
Out of 64 possible laminate designs, there are 17 acceptable solutions. However, an optimum could be found with fewer plies.
To go from 12 to 11 plies we only need to remove one of the instances of the paramter a6:
angs=(0,90)
combs = [[a1,a2,a3,a4,a5,a6,a5,a4,a3,a2,a1]
for a1 in angs for a2 in angs for a3 in angs for a4 in angs for a5 in angs for a6 in angs]
print('No. of combinations:',len(combs))
Results:
maxdef=[]
for c in combs:
layup=[]
for a in c:
layup.append( {'mat':m , 'ori': a , 'thi':0.2} )
ABD=laminatelib.laminateStiffnessMatrix(layup)
load,deformation=laminatelib.solveLaminateLoadCase(ABD, Mx=1000, My=200)
maxdef.append(max(deformation[3:6]))
zipped = zip(maxdef, combs)
sorted_pairs = sorted(zipped)
for s in sorted_pairs:
if s[0]<0.01:
print('Max curvature:',s[0],'Layup:',s[1])
Conclusion: 11 plies are sufficient and there are 4 equivalent layups with respect to weight. Out of the 4 solutions, we may argue that the first in the list above performs better (smallest deformation) and is therefor the optimum solution. This is however not explicitly stated for the objective (but could be included for some added value).
Example 3: strength optimization, fixed ply thickness¶
Constraints and requirements:
Ply material: Carbon/Epoxy(a)
Ply thickness: 0.2 mm
Possible orientations: 0, 90, 45 and -45
Symmetric and balanced laminate
Load cases, both should be considered:
- Case-1: Nxy = 500 N/mm
- Case-2: Mx = 1000 Nmm/mm
Stress exposure factor according to Maximum stress criterion shall be < 1.0 considering all plies and all positions
Objectives:
- Minimize the mass
For an odd number of plies (example: 8 plies):
angs=(0,90,45,-45)
combs = [[a1,a2,a3,a4,a4,a3,a2,a1]
for a1 in angs for a2 in angs for a3 in angs for a4 in angs ]
print('No. of combinations:',len(combs))
def isBalanced(c):
suma=0
for a in c:
if ((a !=0) and (a !=90)):
suma=suma+a
if suma == 0:
return True
else:
return False
balancedCombs=[]
for c in combs:
if isBalanced(c):
balancedCombs.append(c)
print('No. of combinations:',len(balancedCombs))
Thus, there are 256 combinations of symmetric layups but only 70 combinations of both symmetric and balanced laminates.
Another useful function: finding the greates exposure factor across all layers and positions based on the Maximum stress criterion:
def failExposure(LayerResults):
fe = []
for layer in LayerResults:
fe.append( layer['fail']['MS']['bot'] )
fe.append( layer['fail']['MS']['top'] )
return max(fe)
Solving and and reporting the five best solutions:
def solveStrength(combinations):
maxfe=[]
for c in combinations:
layup=[]
for a in c:
layup.append( {'mat':m , 'ori': a , 'thi':0.2} )
ABD=laminatelib.laminateStiffnessMatrix(layup)
loads,deformations=laminatelib.solveLaminateLoadCase(ABD, Nxy=500)
layerResults=laminatelib.layerResults(layup,deformations)
maxfe1=failExposure(layerResults)
loads,deformations=laminatelib.solveLaminateLoadCase(ABD, Mx=1000)
layerResults=laminatelib.layerResults(layup,deformations)
maxfe2=failExposure(layerResults)
maxfe.append(max(maxfe1,maxfe2))
zipped = zip(maxfe, combinations)
sorted_pairs = sorted(zipped)
return sorted_pairs
sp=solveStrength(balancedCombs)
for i in range(0,5):
print(sp[i][0],sp[i][1])
Number of plies must be increased to meet the requirements:
angs=(0,90,45,-45)
combs = [[a1,a2,a3,a4,a5,a6,a6,a5,a4,a3,a2,a1]
for a1 in angs for a2 in angs for a3 in angs for a4 in angs for a5 in angs for a6 in angs ]
print('No. of combinations:',len(combs))
balancedCombs=[]
for c in combs:
if isBalanced(c):
balancedCombs.append(c)
print('No. of combinations:',len(balancedCombs))
sp=solveStrength(balancedCombs)
for i in range(0,5):
print(sp[i][0],sp[i][1])
Try with one more ply:
angs=(0,90,45,-45)
combs = [[a1,a2,a3,a4,a5,a6,a7,a6,a5,a4,a3,a2,a1]
for a1 in angs for a2 in angs for a3 in angs for a4 in angs for a5 in angs for a6 in angs for a7 in angs ]
print('No. of combinations:',len(combs))
balancedCombs=[]
for c in combs:
if isBalanced(c):
balancedCombs.append(c)
print('No. of combinations:',len(balancedCombs))
sp=solveStrength(balancedCombs)
for i in range(0,5):
print(sp[i][0],sp[i][1])
Very close but one more ply is required:
angs=(0,90,45,-45)
combs = [[a1,a2,a3,a4,a5,a6,a7,a7,a6,a5,a4,a3,a2,a1]
for a1 in angs for a2 in angs for a3 in angs for a4 in angs for a5 in angs for a6 in angs for a7 in angs ]
print('No. of combinations:',len(combs))
balancedCombs=[]
for c in combs:
if isBalanced(c):
balancedCombs.append(c)
print('No. of combinations:',len(balancedCombs))
sp=solveStrength(balancedCombs)
for i in range(0,10):
print(sp[i][0],sp[i][1])
n=len(balancedCombs)
print(n)
Although the computational effort to solve all these laminates is manageable, some simple algorithms to identify and eliminate combinations that obviously seems to be less fitted can easily be made.
Let us randomly pick about 1% of the possible combinations:
from numpy.random import randint
samp=randint(low=0,high=n,size=int(n*0.01)) # about n/100 numbers betwwen 0 and total number of combinations (n)
print(samp)
sampCombs=[]
for i in samp:
sampCombs.append(balancedCombs[i])
Then, solve and sort from best to worst:
sp = solveStrength(sampCombs)
for c in sp:
print(c[0],c[1])
The resulting data can now be analyzed using various algorithms for pattern recognition. The basic question is: what characterizes good solutions and what characterizes bad solutions?
One low-hanging fruit is simply the difference between the composition of orientations. The simple code that follows computes the average number of plies for the different orientations for the 5 best and the 5 worst cases:
def averageNumberOf(sp,ori):
best=0
worst=0
for i in range(0,5):
best=best+sp[i][1].count(ori)
worst=worst+sp[len(sp)-1-i][1].count(ori)
return best/5, worst/5
print(' 0-plies, best/worst: ',averageNumberOf(sp,0))
print('90-plies, best/worst: ',averageNumberOf(sp,90))
print('45-plies, best/worst: ',averageNumberOf(sp,45))
Based on these values, we may do a relatively moderate filtering:
filteredCombs=[]
n=0
for c in balancedCombs:
if c.count(0)>4:
if c.count(90)<5:
if c.count(45)>2:
filteredCombs.append(c)
n=n+1
print('Remaining viable combinations:' ,n)
Additional random sampling from the remaining combinations could further reduce the options. In the following example, about 10% of the 210 combinations are randomly selected:
samp=randint(low=0,high=len(filteredCombs),size=20)
sampCombs=[]
for i in samp:
sampCombs.append(filteredCombs[i])
sp = solveStrength(sampCombs)
for c in sp:
print(c[0],c[1])
The obvious difference between best and worst is the orientation at outermost layers:
finalFilteredCombs=[]
n=0
for c in filteredCombs:
if (c[0]==0 and c[1]==0):
finalFilteredCombs.append(c)
n=n+1
print('Remaining viable combinations:' ,n)
sp = solveStrength(finalFilteredCombs)
for c in sp:
print(c[0],c[1])
