NaN compliance matrix for more cambered airfoil

[tylercritchfield]

If I run the following example for composite material properties it results in a NaN-filled compliance matrix. I've taken this straight from the docs, only changing the airfoil coordinates. (Proposed solution below)

using GXBeam

# airfoil coordinates in the docs --> work great
# xaf = [1.00000000, 0.99619582, 0.98515158, 0.96764209, 0.94421447, 0.91510964, 0.88074158, 0.84177999, 0.79894110, 0.75297076, 0.70461763, 0.65461515, 0.60366461, 0.55242353, 0.50149950, 0.45144530, 0.40276150, 0.35589801, 0.31131449, 0.26917194, 0.22927064, 0.19167283, 0.15672257, 0.12469599, 0.09585870, 0.07046974, 0.04874337, 0.03081405, 0.01681379, 0.00687971, 0.00143518, 0.00053606, 0.00006572, 0.00001249, 0.00023032, 0.00079945, 0.00170287, 0.00354717, 0.00592084, 0.01810144, 0.03471169, 0.05589286, 0.08132751, 0.11073805, 0.14391397, 0.18067874, 0.22089879, 0.26433734, 0.31062190, 0.35933893, 0.40999990, 0.46204424, 0.51483073, 0.56767889, 0.61998250, 0.67114514, 0.72054815, 0.76758733, 0.81168064, 0.85227225, 0.88883823, 0.92088961, 0.94797259, 0.96977487, 0.98607009, 0.99640466, 1.00000000]
# yaf = [0.00000000, 0.00017047, 0.00100213, 0.00285474, 0.00556001, 0.00906779, 0.01357364, 0.01916802, 0.02580144, 0.03334313, 0.04158593, 0.05026338, 0.05906756, 0.06766426, 0.07571157, 0.08287416, 0.08882939, 0.09329359, 0.09592864, 0.09626763, 0.09424396, 0.09023579, 0.08451656, 0.07727756, 0.06875796, 0.05918984, 0.04880096, 0.03786904, 0.02676332, 0.01592385, 0.00647946, 0.00370956, 0.00112514, -0.00046881, -0.00191488, -0.00329201, -0.00470585, -0.00688469, -0.00912202, -0.01720842, -0.02488211, -0.03226730, -0.03908459, -0.04503763, -0.04986836, -0.05338180, -0.05551392, -0.05636585, -0.05605816, -0.05472399, -0.05254383, -0.04969990, -0.04637175, -0.04264894, -0.03859653, -0.03433153, -0.02996944, -0.02560890, -0.02134397, -0.01726049, -0.01343567, -0.00993849, -0.00679919, -0.00402321, -0.00180118, -0.00044469, 0.00000000]

# my airfoil coordinates
xaf = [1.0, 0.99669031, 0.98707731, 0.97177624, 0.95120013, 0.92554805, 0.89513234, 0.86042844, 0.82194932, 0.78024605, 0.73589923, 0.68948289, 0.64155561, 0.59265387, 0.5432852, 0.49390561, 0.44491485, 0.39667271, 0.34949257, 0.30376269, 0.25995386, 0.21850111, 0.17980252, 0.14421934, 0.11207729, 0.08364666, 0.05914156, 0.03873701, 0.02260069, 0.01080931, 0.00329342, 0.00180844, 0.00078264, 0.00016559, 2.4e-6, 9.545e-5, 0.0005152, 0.00124174, 0.00277989, 0.00480487, 0.0111091, 0.02423919, 0.04202841, 0.0645839, 0.09186189, 0.1236919, 0.15984299, 0.20001108, 0.24382591, 0.29085133, 0.34057636, 0.39242047, 0.4457723, 0.49999633, 0.55444002, 0.60844101, 0.66133467, 0.71246137, 0.7611783, 0.80686551, 0.84892722, 0.8867981, 0.91995456, 0.9479178, 0.97026156, 0.98661057, 0.99662097, 1.0]
yaf = [0.0, 0.00049003, 0.00207033, 0.00475426, 0.00825965, 0.01241627, 0.01723294, 0.0226708, 0.02862516, 0.034964, 0.0415088, 0.0480415, 0.05432933, 0.06013032, 0.06520439, 0.06931828, 0.07229087, 0.07399667, 0.0744182, 0.07365704, 0.07178323, 0.06883168, 0.064851, 0.05990381, 0.05406639, 0.04742811, 0.04011465, 0.03230144, 0.02418269, 0.01593784, 0.00791465, 0.00553261, 0.00337629, 0.0013778, 0.00015663, -0.00093355, -0.00196476, -0.0030302, -0.00467095, -0.00634666, -0.0102234, -0.01556325, -0.02006888, -0.02353976, -0.02598343, -0.02742948, -0.02793392, -0.02758763, -0.026503, -0.02481832, -0.02268597, -0.02024403, -0.01760693, -0.01488444, -0.01217351, -0.00956485, -0.00713548, -0.0049519, -0.0030634, -0.00151205, -0.00032506, 0.00048717, 0.00093779, 0.00105862, 0.00091038, 0.00057426, 0.00018251, 0.0]

chord = 1.9
twist = 0.0*pi/180
paxis = 0.4750 / chord

uni = Material(37.00e9, 9.00e9, 9.00e9, 4.00e9, 4.00e9, 4.00e9, 0.28, 0.28, 0.28, 1.86e3)
double = Material(10.30e9, 10.30e9, 10.30e9, 8.00e9, 8.00e9, 8.00e9, 0.30, 0.30, 0.30, 1.83e3)
gelcoat = Material(1e1, 1e1, 1e1, 1.0, 1.0, 1.0, 0.30, 0.30, 0.30, 1.83e3)
nexus = Material(10.30e9, 10.30e9, 10.30e9, 8.00e9, 8.00e9, 8.00e9, 0.30, 0.30, 0.30, 1.664e3)
balsa = Material(0.01e9, 0.01e9, 0.01e9, 2e5, 2e5, 2e5, 0.30, 0.30, 0.30, 0.128e3)
mat = [uni, double, gelcoat, nexus, balsa]

t = 0.001
theta = 20*pi/180
layer = Layer(balsa, t, theta)

idx = [3, 4, 2]  # material index
t = [0.000381, 0.00051, 18*0.00053]
theta = [0, 0, 20]*pi/180
layup1 = Layer.(mat[idx], t, theta)

idx = [3, 4, 2]
t = [0.000381, 0.00051, 33*0.00053]
theta = [0, 0, 20]*pi/180
layup2 = Layer.(mat[idx], t, theta)

idx = [3, 4, 2, 1, 5, 1, 2]
t = [0.000381, 0.00051, 17*0.00053, 38*0.00053, 1*0.003125, 37*0.00053, 16*0.00053]
theta = [0, 0, 20, 30, 0, 30, 20]*pi/180
layup3 = Layer.(mat[idx], t, theta)

idx = [3, 4, 2, 5, 2]
t = [0.000381, 0.00051, 17*0.00053, 0.003125, 16*0.00053]
theta = [0, 0, 20, 0, 0]*pi/180
layup4 = Layer.(mat[idx], t, theta)

segments = [layup1, layup2, layup3, layup4]

xbreak = [0.0, 0.0041, 0.1147, 0.5366, 1.0]

idx = [1, 5, 1]
t = [38*0.00053, 0.003125, 38*0.00053]
theta = [0, 0, 0]*pi/180
web = Layer.(mat[idx], t, theta)

webs = [web, web]

webloc = [0.15, 0.5]
nodes, elements = afmesh(xaf, yaf, chord, twist, paxis, xbreak, webloc, segments, webs)

S, sc, tc = compliance_matrix(nodes, elements)
M, mc = mass_matrix(nodes, elements)

I traced it to line 492 of afmesh.jl. There, yiu and yil are searched to find their intersection. Because the nodes don't match up exactly (I'm presuming this is the reason), the values are not compared to each other but rather to the value of zero. In the example, this works great. Below I've plotted yiu (blue) and yil (red), and sure enough they intersect when both surfaces cross zero.

But in my example, this isn't the case (see below). In fact, my second point in yil happens to be barely positive. An easy workaround for that is to limit the possible range to the back half or so of the surface. But even then, the intersection doesn't match up with where they cross zero...

I have an idea on how to get around this - I don't know if this would impact the rest of the code much, but I think it does what the current version is trying to do. The idea is to search backwards from the TE, as it appears that in both cases the nodes match up with each other (i.e. node 40 on top is approximately lined up with node 40 on the bottom, etc.), so instead of comparing both sets of nodes against zero you could compare to each other, such as

# index from the end for both top and bottom 
idx_intersection_from_TE = findlast(xt->xt<0, (reverse(yiu[end-20:end]) .- reverse(yil[end-20:end]))) 
# Note: need to specify a specific range so the number of elements is the same - instead of 20 it could 
# be half of the total length of one of them, I don't think it really matters as long as it's long enough 
# to include the intersection

# since they have different lengths, compute the index from the start separately for both surfaces
iu = length(yiu) - idx_intersection_from_TE + 1
il = length(yil) - idx_intersection_from_TE + 1

This might be overkill or perhaps not robust to every case, but it at least works for these two scenarios.

Thoughts @andrewning?

[andrewning]

I made an assumption that doesn't make any sense in retrospect. Will be pushing a fix shortly.