Friday, September 22, 2017

Scattering - RCS calculation with NEC2

Would you like to know the scattering Radar Cross Section (RCS) of your RC or UAV quadcopter toy?  If you are a Real, Card Carrying Geek, then of course you would...

Calculation of the RCS of an aircraft is difficult, due to the lack of good tools.  Even the commercial tools struggle with this problem, so before you waste money on a super expensive commercial tool, do try the Free ones first.

Here is some background information on the Numerical Electromagnetic Code (NEC), which is a method of moments type of successive approximation program needed to do the modelling:
http://www.aeronetworks.ca/2017/07/olde-skool-antenna-design-with-nec2-on.html

To calculate scattering of a plane wave, you need to make an electrical model of the aircraft.  Part3 of the NEC2 manual contains an example for calculation of the RCS of a stick model of an aircraft.  This looks absurdly simple, but there is method in the madness, since in many cases a stick model is good enough.  An aircraft with a cylindrical, conductive body (alluminium or carbon fibre composite) can be well enough approximated by simple thick wires.

However, if you want a better estimate, then you need to take a CAD model of the aircraft and convert that into a 1D mesh of short wire segments.  This can be done with the gmsh program, and then you need to convert the result into a format that your method of moments antenna modelling tool  can read.

Sticks and Stones

Here is a NEC2 stick model of a small aircraft for RCS calculation with xnec2c:

CM An approximation of a small aircraft Radar Cross Section
CM
CM AV dimensions:
CM The AV fuselage is about 3086 mm long and 680 mm high at the centre.
CM First approximation is a wire of 3086 mm, 340 mm radius.
CM Second approximation is a tapered wire 3086 long x 340 mm radius at hub,
CM tapering to 50 mm radius at the ends.
CM
CM The Gain is the scattering cross section in square wavelengths (sigma/lambda^2.
CM
CM Long range surveillance radar:
CM Fequency = 300 MHz:
CM Wave length = 1 m
CM 1/10th wave length = 100 mm
CM
CM Patch/Wire segment sizes:
CM Length: 3086 mm / 100 mm = 31 sections
CM Height: 680 mm / 100 mm = 7 sections
CM Tail: 2048 mm /100 mm = 21 sections
CM Nose: 1038 mm / 100 mm = 10 sections
CM Radius at hub = 340 mm
CM Taper radius at nose/tail: 50 mm
CM
CM Single wire:
CM GW 1,31,0,0,0,3.086,0,0,.340
CM
CM Two tapered wire segments:
CM GW 1,21,0,0,0,2.048,0,0,0
CM GC 0,0,1,.050,.340
CM GW 2,10,2.048,0,0,3.086,0,0,0
CM GC 0,0,1,.340,.050
CM
CM Excite:
CM EX 1,1,1,0,0
CM EX 1,1,1,0,90,30,-90
CM
CM Radiation Pattern examples for xnec2c:
CM RP 0,72,72,0,0,0,5,5,0,0,0
CM RP 0,90,90,0,0,0,4,4,0,0,0
CM
CE
GW 1,21,0,0,0,2.048,0,0,0
GC 0,0,1,.050,.340
GW 2,10,2.048,0,0,3.086,0,0,0
GC 0,0,1,.340,.050
GE
FR 0,1,0,0,300
EX 1,1,1,0,0
RP 0,90,90,0,0,0,4,4,0,0,0
EN


This model uses fat wire segments for the body, tapered towards the nose and tail, to make it a little more realistic.  The tail planes, wings and landing gear can be handled similarly, but there is diminishing returns - when you add the smaller features, the RCS doesn't change by much, since they are so much smaller than the main fuselage and wings.

CAD Solid Model to 1D Mesh

A CAD model STP file can be converted into a surface mesh with gmsh, but the default settings result in an enormous number (millions) of triangles and lines.  To be able to process the data on a general purpose computer, you need to reduce that to a few thousand elements.

Otherwise, you need to head over to a place like Amazon or Digital Ocean and rent a large virtual machine with many cores and oodles of RAM for a month, to do your calculation.

The gmsh program can read most CAD file types.  The output should be in a form that is easy to parse.  The following line will create an INRIA mesh file:

$ gmsh model.stp -1 clmin 20 -clmax 100 -format auto -o model-1d.mesh

The maximum edge length is set to 100 mm, which is one tenth of a typical long range 300 MHz radar with 1 m wave length.  The result for me was a file with about 30,000 edges.

Mesh to NEC2

Converting the mesh to something that NEC2 can handle required some scripting.  The following Bash script which I creatively named m2n, loops through the mesh file one line at a time to read each wire segment and then search the vertex file from the top down to find the start and end points x,y,z co-ordinates to generate a GW card for each wire.

Note that there may also be a z-axis problem - more on that down below.

#! /bin/bash
echo Convert a 1d INRIA mesh file to a NEC2 wire segment file
#
# Configuration:
# You have to set the file names down below
#
# Mesh generate example:
# ./Applications/gmsh model.stp -1 clmin 20 -clmax 100 -format auto -o model.mesh
#
# INRIA mesh format:
# MeshVersionFormatted 2
# Dimension
# 3
# Vertices
# 4
#     -303.41676307005          -79.958172002481           239.91717070429      1
#     -296.94332277259          -79.958171898826           240.50417216528      2
#     -251.72110263339          -87.385817660601           208.53092734841      3
#     -251.90752255856          -89.319412217105            210.5867628726      4
# Edges
# 4
# 1 2 1
# 2 3 1
# 3 4 2
# 4 1 2
#
# Vertex Format: x, y, z, vertex#
# Edge format: v1, v2, line#
#
# Therefore each edge line segment is defined by a record pointing to two x,y,z vertices
# The STP file is assumed to be in mm and is scaled to m with a GS card
#
# NEC Text File Format:
# The original Fortran NEC used punch cards with strict columnar data
# It is now done with lines in a text file
# A text line can use either spaces or commas as field delimeters
# Empty fields are padded with a zero
# The end of the line need not be padded with zeros
# Floating point numbers need to be in exponential format with decimal points
# Bash doesn't do exponential arithmetic (use bc), but the printf function can output it
# The locale settings in the printf statements ensure use of decimal points
#
#
# File Names:
#export mesh="short.mesh"
export vert="model.vert"
export nec="model.nec"
export mesh="model.mesh"
# Mesh file parse state:
# header, edges, edge, vertex1, vertex2, end
export state="header"
export edgestart="0"
export edgeend="0"
export edgenum="0"
export errmax=50000
export vertexnum=0
# NEC file data:
export wire="1"
export segments="1"
# Copy the vertices
$( cp "$mesh" "$vert" )
# Create the header of the NEC file
echo "CM NEC2 Card Stack Produced by m2n">"$nec"
echo "CE">>"$nec"
# Read lines from the mesh file
end_of_mesh_file=0
while [[ $end_of_mesh_file == 0 ]]
do
  read -r meshline
  # the last exit status is the
  # flag of the end of file
  end_of_mesh_file=$?
  #echo "meshline: $meshline"
  if [ "$meshline" == "End" ]
  then
    state="end"
    echo "state: $state"
    break
  fi
  if [ "$state" == "header" ]
  then
    #echo "state: $state"
    if [ "$meshline" == "Edges" ]
    then
      state="edges"
    fi
  elif [ "$state" == "edges" ]
  then
    #echo "state: $state"
    total="$meshline"
    state="edge"
  elif [ "$state" == "edge" ]
  then
    #echo "state: $state"
    edgeverts=($(echo $meshline | tr ' ' "\n"))
    # echo "edgestart: ${edgeverts[0]}"
    edgestart="${edgeverts[0]}"
    # echo "edgeend: ${edgeverts[1]}"
    edgeend="${edgeverts[1]}"
    # echo "edgenum: ${edgeverts[2]}"
    edgenum="${edgeverts[2]}"
    state="vertex1"
  fi
  if [ "$state" == "vertex1" ]
  then
    #echo "state: $state"
    # Read lines from the vertices file
    vertexnum=-5
    end_of_vertices_file=0
    while [[ $end_of_vertices_file == 0 ]]
    do
      read -r vertexline
      # the last exit status is the
      # flag of the end of file
      end_of_vertices_file=$?
      if [ "$vertexline" == "Edges" ]
      then
         echo "Error1: Vertex not found"
         exit 1
      fi
      if [ "$vertexnum" -ge "$errmax" ]
      then
         echo "Error11: Vertex not found"
         exit 1
      fi
      # Count the records manually
      (( vertexnum++ ))
      if [ "$vertexnum" == "$edgestart" ]
      then
        #echo "vertexline: $vertexline"
        vert1=($(echo $vertexline | tr ' ' "\n"))
        # echo "vertx1: ${vert1[0]}"
        vertx1="${vert1[0]}"
        # echo "verty1: ${vert1[1]}"
        verty1="${vert1[1]}"
         # echo "vertz1: ${vert1[2]}"
        vertz1="${vert1[2]}"
        (( wire++ ))
        break
      fi
    done < "$vert"
    state="vertex2"
    fi
    if [ "$state" == "vertex2" ]
    then
      #echo "state: $state"
      # Read lines from the vertices file
      vertexnum=-5
      end_of_vertices_file=0
      while [[ $end_of_vertices_file == 0 ]]
      do
        read -r vertexline
        # the last exit status is the
        # flag of the end of file
        end_of_vertices_file=$?
        if [ "$vertexline" == "Edges" ]
        then
           echo "Error2: Vertex not found"
           exit 1
        fi
        if [ "$vertexnum" -ge "$errmax" ]
        then
           echo "Error21: Vertex not found"
           exit 1
        fi
        # count the records manually
        (( vertexnum++ ))
        if [ "$vertexnum" == "$edgeend" ]
        then
          #echo "vertexline: $vertexline"
          vert2=($(echo $vertexline | tr ' ' "\n"))
          # echo "vertx2: ${vert2[0]}"
          vertx2="${vert2[0]}"
          # echo "verty2: ${vert2[1]}"
          verty2="${vert2[1]}"
          # echo "vertz2: ${vert2[2]}"
          vertz2="${vert2[2]}"
          LC_NUMERIC="en_US.UTF-8" printf "GW $wire $segments %.8E %.8E %.8E %.8E %.8E %.8E 1\n" \
            $vertx1 $verty1 $vertz1 $vertx2 $verty2 $vertz2
          LC_NUMERIC="en_US.UTF-8" printf "GW $wire $segments %.8E %.8E %.8E %.8E %.8E %.8E 1\n" \
            $vertx1 $verty1 $vertz1 $vertx2 $verty2 $vertz2 >>"$nec"
          break
        fi
      done < "$vert"
      state="edge"
    fi
  done < "$mesh"
# Create footer of NEC file
# Scale everything from mm to m
echo "GS 0 0 .001">>"$nec"
echo "GE">>"$nec"
echo "FR 0 1 0 0 300">>"$nec"
echo "EX 1 1 1 0 0">>"$nec"
echo "RP 0 90 90 0 0 0 4 4 0 0 0">>"$nec"
echo "EN">>"$nec"
exit 0


This script is sloooow...

Meshing a CAD file of a large model with gmsh may take 15 minutes, but converting the result to a NEC card file with this m2n script can take a whole day.  You have to test it with a much shortened file and you should make some errors in the file to see what the script does with them - having it exit due to an error after 24 hours of running on your personal two core 1 GHz super computer, would be disappointing.

Just commenting out most echo debug statements speeded the script up by about ten times, but the important thing is correctness, not speed, so I left one printf to the screen to see what it is doing.  However if someone would rewrite the program in C, please send me a copy!

Computing Troubles

All the meshing programs that I tried, have multiple bugs.  The INRIA mesh file created by gmsh has sequence number problems, but other meshers caused little bits and pieces of aircraft hanging in space around the main model, so gmsh is the least of all evils.

The vertex sequence number overflows and goes from 12076 back to 1 and then get stuck on 3 for a while and so on.   I therefore had to modify the above script to count the vertex record numbers by itself, since it cannot trust the values in the mesh file.  In the script I simply started the line counter at minus 5, to account for the 5 lines of overhead.

The edge record numbers are also doubled up: 1, 1, 2, 2, 3, 3..., but those are not used for anything here.

So with all meshing programs, you need to look at the data files with a text editor and see whether they make sense.  Avoid the binary file formats, since then you will never be able to figure out what went wrong.  Also note that if the resulting pattern plot looks like a hedgehog on a LSD trip, then your wire segments may simply be too long.  The wire segments should be about 1/10 the length of the incident wave.

Another problem is the speed of the calculations.  Much time is wasted accessing the disk drive and that can be alleviated with the use of a small RAM disk.  A RAM disk would be about 10 times faster than a SSD, which is also a good ten times faster than spinning rust.  On Linux systems, the tmpfs used for the /tmp directory is a RAM disk, so it is as simple as copying your files over there.

On a Mac, you need to manually make a RAM disk.  Here is a handy script that will make one and mount it inside your home directory.  Just remember that when you reboot, it goes to the great bit bucket in the sky:

#!/bin/bash 
echo Create a small RAM Disk in the user home directory
ramfs_size_mb=20 
mount_point=~/ramdisk 
ramfs_size_sectors=$((${ramfs_size_mb}*1024*1024/512)) 
ramdisk_dev=`hdid -nomount ram://${ramfs_size_sectors}` 
newfs_hfs -v 'Volatile' ${ramdisk_dev} 
mkdir -p ${mount_point} 
mount -o noatime -t hfs ${ramdisk_dev} ${mount_point} 


To process 30,000 edges, would require a very good super computer.  The complexity of the calculation increases exponentially with increased wire segment count.  On a small computer, you should be able to calculate scattering of a model with 3000 to 5000 elements.  It also looks like xnec2c currently has a maximum limit of 10,000 elements.  So if you only have a small computer, then you need to break your aircraft model up into multiple segments, to simplify your mesh models, but you will have to choose the segments wisely to get useful results.

Depending on the CAD model origin, one may also need to do some shifting on the z-axis, to ensure that the whole model is above the ground plane.  It is best to do this with a separate program to scan the whole file, find the minimum z and then shift everything up in a second pass.  This is hard to do since Bash doesn't handle floating point numbers, so one has to pipe the numbers into the calculator bc.  Here is my wonderfully named NEC Shift Up (nsu) script for this purpose:

#! /bin/bash
echo This script takes a NEC model file and shift it up in the z-axis,
echo to ensure that the whole model is 1 mm above zero.
#
# NEC Text File Format:
# The original Fortran NEC used punch cards with strict columnar data
# It is now done with lines in a text file
# A text line can use either spaces or commas as field delimeters
# Empty fields are padded with a zero
# The end of the line need not be padded with zeros
# Floating point numbers need to be in exponential format with decimal points
# Bash doesn't do exponential arithmetic (use bc), but the printf function can output it
# The locale settings in the printf statements ensure use of decimal points
#
# File Names:
export nec="model.nec"
export shift="shift.nec"
#
# Variables
export zmin="0"
export necline=""
export shiftline=""
#
# States: header, gw, head, shift, end
export state="header"
#
# Read lines from the NEC file to find the minimum z
end_of_nec_file=0
while [[ $end_of_nec_file == 0 ]]
do
  read -r necline
  # The last exit status is the flag of the end of file
  end_of_nec_file=$?
  #echo "necline: $necline"
  if [ "$necline" == "End" ]
  then
    echo "Error1: Premature EOF"
    exit 1
  fi
  if [ "$state" == "header" ]
  then
    echo "state: $state"
    if [ "$necline" == "CE" ]
    then
      state="gw"
      echo "state: $state"
    fi
  elif [ "$state" == "gw" ]
  then
    #echo "state: $state"
    if [ "$necline" == "GE" ]
    then
      state="head"
      break
    fi
    # Parse the GW x,y,z values into an array of 10 variables
    gw=($(echo $necline | tr ' ' "\n"))
    val="${gw[5]}"
    z=$( printf "%f" $val )
    #echo "zmin: $zmin, z: $z"
    if (( $(echo "$zmin > $z" |bc -l) ))
    then
      zmin="$z"
      echo "zmin: $zmin"
    fi
    val="${gw[8]}"
    z=$( printf "%.8f" $val )
    #echo "zmin: $zmin, z: $z"
    if (( $(echo "$zmin > $z" |bc -l) ))
    then
      zmin="$z"
      echo "zmin: $zmin"
    fi
  fi
done < "$nec"
echo "state: $state"
echo "zmin: $zmin"
# Read lines from the NEC file, write to a new file and scale up by minimum z
$( rm -f "shift" )
$( touch "shift" )
end_of_nec_file=0
while [[ $end_of_nec_file == 0 ]]
do
  read -r necline
  # The last exit status is the flag of the end of file
  end_of_nec_file=$?
  #echo "necline: $necline"
  if [ "$necline" == "End" ]
  then
    echo "Error1: Premature EOF"
    exit 1
  fi
  if [ "$state" == "head" ]
  then
    echo "state: $state"
    if [ "$necline" == "CE" ]
    then
      state="shift"
      echo "state: $state"
    fi
    echo "$necline">>$shift
  elif [ "$state" == "shift" ]
  then
    #echo "state: $state"
    # Parse the GW x,y,z values into an array of 10 variables
    gw=($(echo $necline | tr ' ' "\n"))
    gs="${gw[0]}"
    scale="${gw[3]}"
    if [ "$gs" == "GS" ]
    then
      echo "GS 0 0 $scale\n"
      echo "GS 0 0 $scale" >>$shift
      state="end"
      echo "state: $state"
    else
      wire="${gw[1]}"
      segments="${gw[2]}"
      x1="${gw[3]}"
      y1="${gw[4]}"
      val="${gw[5]}"
      z=$( printf "%f" $val )
      z1=$(echo "$z - $zmin + 1" |bc -l)
      x2="${gw[6]}"
      y2="${gw[7]}"
      val="${gw[8]}"
      z=$( printf "%.8f" $val )
      z2=$(echo "$z - $zmin + 1" |bc -l)
      printf "GW $wire $segments $x1 $y1 %.8E $x2 $y2 %.8E 1\n" $z1 $z2
      printf "GW $wire $segments $x1 $y1 %.8E $x2 $y2 %.8E 1\n" $z1 $z2 >>$shift
    fi
  elif [ "$state" == "end" ]
  then
    #echo "state: $state"
    echo "$necline"
    echo "$necline">>$shift
    if [ "$necline" == "EN" ]
    then
      state="done"
      break
    fi
  fi
done < "$nec"
echo "state: $state"
exit 0

The next problem is finding a computer with enough memory that can process all this without crashing. A Macbook Pro would certainly be a better candidate than a netbook.

Why do I do this with Bash?  I must be a masochist, but all workstations have Bash (even Windoze now has Bash) and a script allows one to debug text processing more rapidly than with C.

Thanks

Thank you to Neoklis Kyriazis http://www.qsl.net/5b4az/ for translating NEC2 to C.



Have fun!

Herman


No comments:

Post a Comment

On topic comments are welcome. Junk will be deleted.