The actuation of a soft robot involves transforming its shape from an initial state to a desired operational state. To achieve task-specific design, it is necessary to map the shape between these two states to the robot's design parameters. This requires both a kinematic model of the soft robot and a shape-matching algorithm. However, existing kinematic models for soft robots are often limited in accuracy and generality due to the robot's flexibility and nonlinearity, and current shape-matching algorithms are not well-suited for 3D cases. To address this challenge, this paper presents a shape-matching design framework for bellow soft pneumatic actuators (SPAs) to expedite the actuator design process. First, a kinematic model of the bellow SPA is developed based on its novel modular design and a surrogate model, which is trained using an Artificial Neural Network and a dataset from Finite Element Method (FEM) simulations. Then, a 3D shape-matching algorithm, composed of a 3D piecewise-constant curvature segmentation and a bi-level Bayesian optimisation algorithm based on the surrogate model, is presented to find the optimal actuator design parameters that match the desired shape. An open-source design toolbox SPADA (Soft Pneumatic Actuator Design frAmework) is also developed to facilitate the use of the proposed design framework, including FEM simulation, shape-matching optimisation based on surrogate modelling, and automatic generation of the ready-to-print CAD file. Experimental results show an averaged root-mean-square error of 2.74 mm, validating the accuracy of the kinematics model. To demonstrate the proposed design framework, actuators are designed to match the predefined shapes in 2D and 3D space.